move sources to src directory

src/shapes/BoolShape.cc

@@ -0,0 +1,9 @@
+
+#include "BoolShape.h"
+
+void BoolShape::setMaterial(refptr<Material> material)
+{
+    m_material = material;
+    m_shape1->setMaterial(material);
+    m_shape2->setMaterial(material);
+}

src/shapes/BoolShape.h

@@ -0,0 +1,19 @@
+
+#ifndef BOOLSHAPE_H
+#define BOOLSHAPE_H BOOLSHAPE_H
+
+#include "Shape.h"
+
+class BoolShape : public Shape
+{
+    public:
+        virtual IntersectionList intersect(refptr<Shape> _this, const Ray & ray) = 0;
+        virtual void setMaterial(refptr<Material> material);
+        virtual refptr<Shape> clone() = 0;
+
+    protected:
+        refptr<Shape> m_shape1;
+        refptr<Shape> m_shape2;
+};
+
+#endif

src/shapes/Box.cc

@@ -0,0 +1,111 @@
+
+#include "Box.h"
+#include "util/Solver.h"
+#include <math.h>
+#include <iostream>
+using namespace std;
+
+#define FP_EQUAL(x,y) (fabs((x)-(y)) < 0.000001)
+
+Box::Box(refptr<Vector> size)
+{
+    m_size = *size;
+    m_size[0] = fabs(m_size[0]) / 2.0;
+    m_size[1] = fabs(m_size[1]) / 2.0;
+    m_size[2] = fabs(m_size[2]) / 2.0;
+}
+
+Shape::IntersectionList Box::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+
+    IntersectionList res;
+
+    /*** optimization ***/
+    if (ray_inv.getOrigin()[0] > m_size[0])
+    {
+        if (ray_inv.getDirection() % Vector::X > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[0] < -m_size[0])
+    {
+        if (ray_inv.getDirection() % Vector::X < 0)
+            return res;
+    }
+    if (ray_inv.getOrigin()[1] > m_size[1])
+    {
+        if (ray_inv.getDirection() % Vector::Y > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[1] < -m_size[1])
+    {
+        if (ray_inv.getDirection() % Vector::Y < 0)
+            return res;
+    }
+    if (ray_inv.getOrigin()[2] > m_size[2])
+    {
+        if (ray_inv.getDirection() % Vector::Z > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[2] < -m_size[2])
+    {
+        if (ray_inv.getDirection() % Vector::Z < 0)
+            return res;
+    }
+    /*** end optimization ***/
+
+    /*
+     * Ray equation:   R = R0 + tRd
+     *                 x = R0x + tRdx
+     *                 y = R0y + tRdy
+     *                 z = R0z + tRdz
+     * Side equation:  x - size_x = 0
+     * Combined:       R0x + (t)Rdx - size_x = 0
+     */
+    for (int dim = 0; dim < 3; dim++)
+    {
+        for (int side = -1; side <= 1; side += 2)
+        {
+            LinearSolver solver(ray_inv.getDirection()[dim],
+                                ray_inv.getOrigin()[dim]
+                                    + side * m_size[dim]);
+            Solver::Result solutions = solver.solve();
+            for (int i = 0; i < solutions.numResults; i++)
+            {
+                if (solutions.results[i] >= 0.0)
+                {
+                    Vector isect_point = ray_inv[solutions.results[i]];
+                    if ( (dim == 0 || fabs(isect_point[0]) <= m_size[0])
+                      && (dim == 1 || fabs(isect_point[1]) <= m_size[1])
+                      && (dim == 2 || fabs(isect_point[2]) <= m_size[2]) )
+                    {
+                        Vector normal(0, 0, -1);
+
+                        if      ( FP_EQUAL(isect_point[0],  m_size[0]) )
+                            normal = Vector(1, 0, 0);
+                        else if ( FP_EQUAL(isect_point[0], -m_size[0]) )
+                            normal = Vector(-1, 0, 0);
+                        else if ( FP_EQUAL(isect_point[1],  m_size[1]) )
+                            normal = Vector(0, 1, 0);
+                        else if ( FP_EQUAL(isect_point[1], -m_size[1]) )
+                            normal = Vector(0, -1, 0);
+                        else if ( FP_EQUAL(isect_point[2],  m_size[2]) )
+                            normal = Vector(0, 0, 1);
+
+                        res.add(Intersection(_this,
+                                    m_transform.transform_point(isect_point),
+                                    m_transform.transform_normal(normal))
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    return res;
+}
+
+refptr<Shape> Box::clone()
+{
+    return new Box(*this);
+}

src/shapes/Box.h

@@ -0,0 +1,19 @@
+
+#ifndef BOX_H
+#define BOX_H BOX_H
+
+#include "Shape.h"
+
+class Box : public Shape
+{
+    public:
+        Box(refptr<Vector> size);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        virtual refptr<Shape> clone();
+
+    protected:
+        Vector m_size;
+};
+
+#endif
+

src/shapes/Cyl.cc

@@ -0,0 +1,141 @@
+
+#include "Cyl.h"
+#include "util/Solver.h"
+#include <math.h>
+#include <iostream>
+using namespace std;
+
+#define FP_EQUAL(x,y) (fabs((x)-(y)) < 0.000001)
+
+Cyl::Cyl(double bottom_radius, double top_radius, double height)
+{
+    m_bottom_radius = fabs(bottom_radius);
+    m_bottom_radius_2 = bottom_radius * bottom_radius;
+    m_top_radius = fabs(top_radius);
+    m_top_radius_2 = top_radius * top_radius;
+    m_height = fabs(height);
+    if (m_height == 0.0)
+        m_height = 1.0;
+    m_slope = (m_top_radius - m_bottom_radius) / m_height; /* rise over run */
+    if ( ! FP_EQUAL(m_slope, 0.0) )
+        m_inv_slope = -1.0 / m_slope;
+}
+
+Shape::IntersectionList Cyl::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+    IntersectionList res;
+
+    /* First intersect with the bottom plane, if it has positive area */
+    if (m_bottom_radius > 0.0)
+    {
+        LinearSolver solver(-ray_inv.getDirection()[2],
+                            -ray_inv.getOrigin()[2]);
+        Solver::Result solutions = solver.solve();
+        if (solutions.numResults > 0 && solutions.results[0] > 0.0)
+        {
+            Vector isect_point = ray_inv[solutions.results[0]];
+            if (isect_point[0]*isect_point[0] + isect_point[1]*isect_point[1]
+                    <= m_bottom_radius_2)
+            {
+                res.add(Intersection(_this,
+                             m_transform.transform_point(isect_point),
+                             m_transform.transform_normal(Vector(0, 0, -1))));
+            }
+        }
+    }
+
+    /* Same for the top plane */
+    if (m_top_radius > 0.0)
+    {
+        LinearSolver solver(ray_inv.getDirection()[2],
+                            ray_inv.getOrigin()[2] - m_height);
+        Solver::Result solutions = solver.solve();
+        if (solutions.numResults > 0 && solutions.results[0] > 0.0)
+        {
+            Vector isect_point = ray_inv[solutions.results[0]];
+            if (isect_point[0]*isect_point[0] + isect_point[1]*isect_point[1]
+                    <= m_top_radius_2)
+            {
+                res.add(Intersection(_this,
+                             m_transform.transform_point(isect_point),
+                             m_transform.transform_normal(Vector(0, 0, 1))));
+            }
+        }
+    }
+
+    /*
+     * Now see if the ray hit the side of the cylinder/cone thingy
+     * Ray equation:   R = R0 + tRd
+     *                 x = R0x + tRdx
+     *                 y = R0y + tRdy
+     *                 z = R0z + tRdz
+     * Side equation:  x^2 + y^2 = (m*z + b)^2
+     * Combined:       (R0x+t*Rdx)^2 + (R0y+t*Rdy)^2 = (m*(R0z+t*Rdz) + b)^2
+     * Expanded:       (R0x*R0x + (t*t)*Rdx*Rdx + 2*R0x*t*Rdx)
+     *                 + (R0y*R0y + (t*t)*Rdy*Rdy + 2*R0y*t*Rdy)
+     *                 = ((m*R0z+m*t*Rdz)^2 + b*b + 2*m*(R0z+t*Rdz)*b)
+     *                 = ((m*R0z*m*R0z + m*m*t*t*Rdz*Rdz + 2*m*R0z*m*t*Rdz)
+     *                  + b*b + 2*m*R0z*b + 2*m*t*Rdz*b)
+     * Quadratic form: (t*t)(Rdx*Rdx + Rdy*Rdy - m*m*Rdz*Rdz)
+     *                 + (t)(2*R0x*Rdx + 2*R0y*Rdy - 2*m*R0z*m*Rdz - 2*m*Rdz*b)
+     *                 + (R0x*R0x + R0y*R0y - m*m*R0z*R0z - b*b - 2*m*R0z*b)
+     *                 = 0
+     */
+    double Rdx = ray_inv.getDirection()[0];
+    double Rdy = ray_inv.getDirection()[1];
+    double Rdz = ray_inv.getDirection()[2];
+    double R0x = ray_inv.getOrigin()[0];
+    double R0y = ray_inv.getOrigin()[1];
+    double R0z = ray_inv.getOrigin()[2];
+    double m = m_slope;
+    double m2 = m*m;
+    double b = m_bottom_radius;
+    QuadraticSolver solver(Rdx * Rdx + Rdy * Rdy - m2 * Rdz * Rdz,
+                           2.0 * (R0x*Rdx + R0y*Rdy - m2*R0z*Rdz - m*Rdz*b),
+                           R0x*R0x + R0y*R0y - m2*R0z*R0z - b*b - 2*m*R0z*b);
+
+    Solver::Result solutions = solver.solve();
+    for (int i = 0; i < solutions.numResults; i++)
+    {
+        if (solutions.results[i] >= 0.0)
+        {
+            Vector isect_point = ray_inv[solutions.results[i]];
+            if (isect_point[2] > 0.0 && isect_point[2] < m_height)
+            {
+                Vector normal;
+
+                if ( FP_EQUAL(m_slope, 0.0) )
+                {
+                    normal = Vector(isect_point[0], isect_point[1], 0.0);
+                }
+                else
+                {
+                    double x = isect_point[0];
+                    double y = isect_point[1];
+                    double dist = sqrt(  isect_point[0] * isect_point[0]
+                                       + isect_point[1] * isect_point[1] );
+                    if (dist > 0.0)
+                    {
+                        double scale = 1.0 / dist;
+                        x *= scale;
+                        y *= scale;
+                    }
+                    normal = Vector(x, y, m_inv_slope);
+                }
+                normal.normalize();
+
+                res.add(Intersection(_this,
+                                     m_transform.transform_point(isect_point),
+                                     m_transform.transform_normal(normal)));
+            }
+        }
+    }
+
+    return res;
+}
+
+refptr<Shape> Cyl::clone()
+{
+    return new Cyl(*this);
+}

src/shapes/Cyl.h

@@ -0,0 +1,25 @@
+
+#ifndef CYL_H
+#define CYL_H CYL_H
+
+#include "Shape.h"
+
+class Cyl : public Shape
+{
+    public:
+        Cyl(double bottom_radius, double top_radius, double height);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        virtual refptr<Shape> clone();
+
+    protected:
+        double m_bottom_radius;
+        double m_bottom_radius_2;
+        double m_top_radius;
+        double m_top_radius_2;
+        double m_height;
+        double m_slope;
+        double m_inv_slope;
+};
+
+#endif
+

src/shapes/Extrude.cc

@@ -0,0 +1,267 @@
+
+#include <math.h>
+
+#include <iostream>
+#include <algorithm>                /* sort() */
+
+#include "Extrude.h"
+#include "util/Polygon.h"
+#include "util/Solver.h"
+
+using namespace std;
+
+#define FP_EQUAL(x,y) (fabs((x)-(y)) < 0.000001)
+
+Extrude::Extrude()
+{
+}
+
+class IntersectListComparator
+{
+    public:
+        IntersectListComparator(Vector start) : m_start(start) {}
+        bool operator()(const Shape::Intersection & i1,
+                const Shape::Intersection & i2) const
+        {
+            double d1 = (i1.position - m_start).mag2();
+            double d2 = (i2.position - m_start).mag2();
+            return d1 < d2;
+        }
+    protected:
+        Vector m_start;
+};
+
+Shape::IntersectionList Extrude::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+    IntersectionList res;
+
+    /*** optimization ***/
+    if (ray_inv.getOrigin()[0] > m_aabb.max[0])
+    {
+        if (ray_inv.getDirection() % Vector::X > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[0] < m_aabb.min[0])
+    {
+        if (ray_inv.getDirection() % Vector::X < 0)
+            return res;
+    }
+    if (ray_inv.getOrigin()[1] > m_aabb.max[1])
+    {
+        if (ray_inv.getDirection() % Vector::Y > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[1] < m_aabb.min[1])
+    {
+        if (ray_inv.getDirection() % Vector::Y < 0)
+            return res;
+    }
+    if (ray_inv.getOrigin()[2] > m_aabb.max[2])
+    {
+        if (ray_inv.getDirection() % Vector::Z > 0)
+            return res;
+    }
+    else if (ray_inv.getOrigin()[2] < m_aabb.min[2])
+    {
+        if (ray_inv.getDirection() % Vector::Z < 0)
+            return res;
+    }
+    /*** end optimization ***/
+
+    int n_polygons = m_polygons.size();
+    int n_offsets = m_offsets.size();
+
+    double distance = 0.0;
+    Vector scale(1.0, 1.0, 1.0);
+    Vector shift(0.0, 0.0, 0.0);
+    for (int p = 0; p < n_polygons; p++)
+    {
+        refptr<Polygon> polygon = m_polygons[p];
+        for (int o = 0; o < n_offsets; o++)
+        {
+            Offset & offset = m_offsets[o];
+            for (int pt = 0, n_pts = polygon->size(); pt < n_pts; pt++)
+            {
+                Vector p1 = scale.mult(*(*polygon)[pt]) + shift;
+                Vector p2 = scale.mult(*(*polygon)[(pt+1) % n_pts]) + shift;
+                Vector p3 = scale.mult(offset.scale)
+                    .mult(*(*polygon)[(pt+1) % n_pts])
+                    + shift + offset.shift;
+                Vector p4 = scale.mult(offset.scale).mult(*(*polygon)[pt])
+                    + shift + offset.shift;
+                p1[2] += distance;
+                p2[2] += distance;
+                p3[2] += distance + offset.distance;
+                p4[2] += distance + offset.distance;
+                Vector e1 = p2 - p1;
+                Vector e2 = p3 - p2;
+                Vector normal = (e1 * e2).normalize();
+                double a = normal[0], b = normal[1], c = normal[2];
+                double d = -(a * p1[0] + b * p1[1] + c * p1[2]);
+                LinearSolver solver(  a * ray_inv.getDirection()[0]
+                                    + b * ray_inv.getDirection()[1]
+                                    + c * ray_inv.getDirection()[2],
+                                      a * ray_inv.getOrigin()[0]
+                                    + b * ray_inv.getOrigin()[1]
+                                    + c * ray_inv.getOrigin()[2]
+                                    + d);
+                Solver::Result solutions = solver.solve();
+                if (solutions.numResults > 0 && solutions.results[0] > 0.0)
+                {
+                    Vector ipoint = ray_inv[solutions.results[0]];
+                    Polygon quad;
+                    quad.add(p1).add(p2).add(p3).add(p4);
+                    if (quad.containsPointConvex(ipoint))
+                    {
+                        res.add(Intersection(_this,
+                                    m_transform.transform_point(ipoint),
+                                    m_transform.transform_normal(normal)));
+                    }
+                }
+            }
+            distance += offset.distance;
+            scale = scale.mult(offset.scale);
+            shift += offset.shift;
+        }
+    }
+
+    double a = 0, b = 0, c = -1.0, d = 0.0;
+    LinearSolver solver(  a * ray_inv.getDirection()[0]
+                        + b * ray_inv.getDirection()[1]
+                        + c * ray_inv.getDirection()[2],
+                          a * ray_inv.getOrigin()[0]
+                        + b * ray_inv.getOrigin()[1]
+                        + c * ray_inv.getOrigin()[2]
+                        + d);
+    Solver::Result solutions = solver.solve();
+    if (solutions.numResults > 0 && solutions.results[0] > 0.0)
+    {
+        Vector ipoint = ray_inv[solutions.results[0]];
+        for (int p = 0; p < n_polygons; p++)
+        {
+            refptr<Polygon> polygon = m_polygons[p];
+            if (polygon->containsPoint2D(ipoint))
+            {
+                res.add(Intersection(_this,
+                            m_transform.transform_point(ipoint),
+                            m_transform.transform_normal(Vector(a, b, c))));
+            }
+        }
+    }
+    if (scale[0] > 0.0 && scale[1] > 0.0)
+    {
+        a = 0, b = 0, c = 1.0, d = -distance;
+        LinearSolver solver(  a * ray_inv.getDirection()[0]
+                            + b * ray_inv.getDirection()[1]
+                            + c * ray_inv.getDirection()[2],
+                              a * ray_inv.getOrigin()[0]
+                            + b * ray_inv.getOrigin()[1]
+                            + c * ray_inv.getOrigin()[2]
+                            + d);
+        Solver::Result solutions = solver.solve();
+        if (solutions.numResults > 0 && solutions.results[0] > 0.0)
+        {
+            Vector ipoint = ray_inv[solutions.results[0]];
+            for (int p = 0; p < n_polygons; p++)
+            {
+                Polygon tp = *m_polygons[p];
+                for (int i = 0, sz = tp.size(); i < sz; i++)
+                {
+                    tp[i] = new Vector(tp[i]->mult(scale) + shift);
+                    (*tp[i])[2] += distance;
+                }
+                if (tp.containsPoint2D(ipoint))
+                {
+                    res.add(Intersection(_this,
+                                m_transform.transform_point(ipoint),
+                                m_transform.transform_normal(Vector(a, b, c))));
+                    break;
+                }
+            }
+        }
+    }
+
+    sort(res.begin(), res.end(), IntersectListComparator(ray.getOrigin()));
+
+    return res;
+}
+
+void Extrude::addPolygon(refptr<Polygon> polygon)
+{
+    m_polygons.push_back(polygon);
+    updateAABB();
+}
+
+void Extrude::addOffset(double distance,
+        const Vector & scale, const Vector & shift)
+{
+    m_offsets.push_back(Offset(distance, scale, shift));
+    updateAABB();
+}
+
+void Extrude::updateAABB()
+{
+    int n_polygons = m_polygons.size();
+    int n_offsets = m_offsets.size();
+    bool first = true;
+    double distance = 0.0;
+    Vector scale(1.0, 1.0, 1.0);
+    Vector shift(0.0, 0.0, 0.0);
+    for (int p = 0; p < n_polygons; p++)
+    {
+        refptr<Polygon> polygon = m_polygons[p];
+        for (int o = 0; o < n_offsets; o++)
+        {
+            Offset & offset = m_offsets[o];
+            for (int pt = 0, n_pts = polygon->size(); pt < n_pts; pt++)
+            {
+                Vector p1 = scale.mult(*(*polygon)[pt]) + shift;
+                Vector p2 = scale.mult(*(*polygon)[(pt+1) % n_pts]) + shift;
+                Vector p3 = scale.mult(offset.scale)
+                    .mult(*(*polygon)[(pt+1) % n_pts])
+                    + shift + offset.shift;
+                Vector p4 = scale.mult(offset.scale).mult(*(*polygon)[pt])
+                    + shift + offset.shift;
+                p1[2] += distance;
+                p2[2] += distance;
+                p3[2] += distance + offset.distance;
+                p4[2] += distance + offset.distance;
+                if (first)
+                {
+                    m_aabb.min = p1;
+                    m_aabb.max = p1;
+                    first = false;
+                }
+                expandAABB(p1);
+                expandAABB(p2);
+                expandAABB(p3);
+                expandAABB(p4);
+            }
+            distance += offset.distance;
+            scale = scale.mult(offset.scale);
+            shift += offset.shift;
+        }
+    }
+}
+
+void Extrude::expandAABB(const Vector & v)
+{
+    if (v[0] < m_aabb.min[0])
+        m_aabb.min[0] = v[0];
+    if (v[0] > m_aabb.max[0])
+        m_aabb.max[0] = v[0];
+    if (v[1] < m_aabb.min[1])
+        m_aabb.min[1] = v[1];
+    if (v[1] > m_aabb.max[1])
+        m_aabb.max[1] = v[1];
+    if (v[2] < m_aabb.min[2])
+        m_aabb.min[2] = v[2];
+    if (v[2] > m_aabb.max[2])
+        m_aabb.max[2] = v[2];
+}
+
+refptr<Shape> Extrude::clone()
+{
+    return new Extrude(*this);
+}

src/shapes/Extrude.h

@@ -0,0 +1,43 @@
+
+#ifndef EXTRUDE_H
+#define EXTRUDE_H EXTRUDE_H
+
+#include <vector>
+
+#include "util/Polygon.h"
+#include "util/AABB.h"
+
+#include "Shape.h"
+
+class Extrude : public Shape
+{
+    public:
+        Extrude();
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        void addPolygon(refptr<Polygon> polygon);
+        void addOffset(double distance,
+                const Vector & scale, const Vector & shift);
+        class Offset
+        {
+            public:
+                Offset(double d, const Vector & s, const Vector & p)
+                    : distance(d), scale(s), shift(p)
+                {
+                }
+                double distance;
+                Vector scale;
+                Vector shift;
+        };
+        virtual refptr<Shape> clone();
+
+    protected:
+        std::vector< refptr<Polygon> > m_polygons;
+        std::vector<Offset> m_offsets;
+        AABB m_aabb;
+
+        void updateAABB();
+        void expandAABB(const Vector & v);
+};
+
+#endif
+

src/shapes/Intersect.cc

@@ -0,0 +1,115 @@
+
+#include "Intersect.h"
+#include <iostream>
+using namespace std;
+
+Intersect::Intersect(const vector< refptr<Shape> > & shapes)
+{
+    int num_shapes = shapes.size();
+    if (num_shapes > 2)
+    {
+        m_shape2 = shapes[num_shapes - 1];
+        vector< refptr<Shape> > rest = shapes;
+        rest.pop_back();
+        m_shape1 = new Intersect(rest);
+    }
+    else if (num_shapes == 2)
+    {
+        m_shape1 = shapes[0];
+        m_shape2 = shapes[1];
+    }
+    else
+    {
+        cerr << __FILE__ << ": " << __LINE__
+             << ": error: Intersect::Intersect() called with only "
+             << num_shapes
+             << " sub-shapes!"
+             << endl;
+        exit(4);
+    }
+}
+
+Shape::IntersectionList Intersect::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+    IntersectionList res1 = m_shape1->intersect(m_shape1, ray_inv);
+
+    /*** optimization ***/
+    if (res1.size() == 0)
+        return res1;
+    /*** end optimization ***/
+
+    IntersectionList res2 = m_shape2->intersect(m_shape2, ray_inv);
+    BoolIntersectionList merged(res1, res2, ray_inv.getOrigin());
+
+    IntersectionList res;
+    bool in1 = false, in2 = false;
+
+    /* initially go through the merged intersections to see whether
+     * the ray started inside one of the sub-objects */
+    for (int i = 0, sz = merged.size(), saw1 = 0, saw2 = 0;
+         i < sz && (!saw1 || !saw2);
+         i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool back = dot < 0.0;
+        bool left = merged[i].left;
+        if (back)
+        {
+            if (left && !saw1)
+                in1 = true;
+            else if (!left && !saw2)
+                in2 = true;
+        }
+        if (left)
+            saw1 = 1;
+        else
+            saw2 = 1;
+    }
+
+    bool in_bool = in1 && in2;
+    for (int i = 0, sz = merged.size(); i < sz; i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool front = dot > 0.0;
+        bool left = merged[i].left;
+        if (front)
+        {
+            if (left)
+                in1 = true;
+            else
+                in2 = true;
+            if (!in_bool && in1 && in2)
+            {
+                /* we found an intersection point with the boolean object */
+                in_bool = true;
+                res.add( merged[i].intersection.transform(m_transform) );
+            }
+        }
+        else
+        {
+            if (in_bool && in1 && in2)
+            {
+                /* we found an intersection point with the boolean object */
+                res.add( merged[i].intersection.transform(m_transform) );
+            }
+            if (left)
+                in1 = false;
+            else
+                in2 = false;
+            in_bool = false;
+        }
+    }
+
+    return res;
+}
+
+refptr<Shape> Intersect::clone()
+{
+    Intersect * s = new Intersect(*this);
+    s->m_shape1 = m_shape1->clone();
+    s->m_shape2 = m_shape2->clone();
+    return refptr<Shape>(s);
+}

src/shapes/Intersect.h

@@ -0,0 +1,17 @@
+
+#ifndef INTERSECT_H
+#define INTERSECT_H INTERSECT_H
+
+#include "BoolShape.h"
+#include <vector>
+
+class Intersect : public BoolShape
+{
+    public:
+        Intersect(const std::vector< refptr<Shape> > & shapes);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        refptr<Shape> clone();
+};
+
+#endif
+

src/shapes/Plane.cc

@@ -0,0 +1,57 @@
+
+#include "Plane.h"
+#include "util/Solver.h"
+#include <math.h>
+#include <iostream>
+using namespace std;
+
+Plane::Plane(double a, double b, double c, double d)
+{
+    m_a = a;
+    m_b = b;
+    m_c = c;
+    m_d = d;
+    m_normal = Vector(m_a, m_b, m_c);
+    m_normal.normalize();
+}
+
+Shape::IntersectionList Plane::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+
+    IntersectionList res;
+    /*
+     * Plane equation: ax + by + cz + d = 0
+     * Ray equation:   R = R0 + tRd
+     *                 x = R0x + tRdx
+     *                 y = R0y + tRdy
+     *                 z = R0z + tRdz
+     * Combined: a(R0x + tRdx) + b(R0y + tRdy) + c(R0z + tRdz) + d = 0
+     *           aR0x + (t)aRdx + bR0y + (t)bRdy + cR0z + (t)cRdz + d = 0
+     *           (t)(aRdx + bRdy + cRdz) + aR0x + bR0y + cR0z + d = 0
+     */
+    LinearSolver solver(  m_a * ray_inv.getDirection()[0]
+                        + m_b * ray_inv.getDirection()[1]
+                        + m_c * ray_inv.getDirection()[2],
+                          m_a * ray_inv.getOrigin()[0]
+                        + m_b * ray_inv.getOrigin()[1]
+                        + m_c * ray_inv.getOrigin()[2]
+                        + m_d);
+    Solver::Result solutions = solver.solve();
+    if (solutions.numResults > 0)
+    {
+        if (solutions.results[0] > 0.0)
+        {
+            Vector isect_point = ray_inv[solutions.results[0]];
+            res.add(Intersection(_this,
+                                 m_transform.transform_point(isect_point),
+                                 m_transform.transform_normal(m_normal)));
+        }
+    }
+    return res;
+}
+
+refptr<Shape> Plane::clone()
+{
+    return new Plane(*this);
+}

src/shapes/Plane.h

@@ -0,0 +1,20 @@
+
+#ifndef PLANE_H
+#define PLANE_H PLANE_H
+
+#include "Shape.h"
+
+class Plane : public Shape
+{
+    public:
+        Plane(double a, double b, double c, double d);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        virtual refptr<Shape> clone();
+
+    protected:
+        double m_a, m_b, m_c, m_d;
+        Vector m_normal;
+};
+
+#endif
+

src/shapes/Shape.cc

@@ -0,0 +1,70 @@
+
+#include "Shape.h"
+#include <algorithm>                /* sort() */
+#include <functional>               /* binary_function */
+#include <utility>
+using namespace std;
+
+static refptr<Material> default_material;
+static bool default_material_initialized = false;
+
+Shape::Shape()
+{
+    if (default_material_initialized == false)
+    {
+        default_material = new Material();
+        default_material_initialized = true;
+    }
+    m_material = default_material;
+}
+
+Shape::~Shape()
+{
+}
+
+void Shape::setMaterial(refptr<Material> material)
+{
+    m_material = material;
+}
+
+class BoolIntersectionComparator
+    : public std::binary_function<Shape::BoolIntersection,
+                                  Shape::BoolIntersection,
+                                  bool>
+{
+    public:
+        BoolIntersectionComparator(const Vector & refPoint)
+        {
+            m_refPoint = refPoint;
+        }
+        bool operator()(const Shape::BoolIntersection & i1,
+                        const Shape::BoolIntersection & i2) const
+        {
+            return (  m_refPoint.dist_to(i1.intersection.position)
+                    < m_refPoint.dist_to(i2.intersection.position) );
+        }
+    protected:
+        Vector m_refPoint;
+};
+
+Shape::BoolIntersectionList::BoolIntersectionList(const IntersectionList & l1,
+                                                  const IntersectionList & l2,
+                                                  const Vector & startPoint)
+{
+    for (size_t i = 0, sz = l1.size();
+         i < sz;
+         i++)
+    {
+        m_intersections.push_back( BoolIntersection(l1[i], true) );
+    }
+    for (size_t i = 0, sz = l2.size();
+         i < sz;
+         i++)
+    {
+        m_intersections.push_back( BoolIntersection(l2[i], false) );
+    }
+
+    sort(m_intersections.begin(),
+         m_intersections.end(),
+         BoolIntersectionComparator(startPoint));
+}

src/shapes/Shape.h

@@ -0,0 +1,131 @@
+
+#ifndef SHAPE_H
+#define SHAPE_H SHAPE_H
+
+#include "util/Solver.h"
+#include "util/Ray.h"
+#include "util/Vector.h"
+#include "util/Transform.h"
+#include "util/Material.h"
+#include "util/refptr.h"
+#include <vector>
+#include <utility>
+
+class Shape
+{
+    public:
+        class Intersection
+        {
+            public:
+                Intersection() {};
+                Intersection(refptr<Shape> shape,
+                             const Vector & position,
+                             const Vector & normal)
+                {
+                    this->shape = shape;
+                    this->position = position;
+                    this->normal = normal;
+                }
+                refptr<Shape> shape;
+                Vector position;
+                Vector normal;
+                Intersection transform(Transform & t)
+                {
+                    return Intersection(shape,
+                            t.transform_point(position),
+                            t.transform_normal(normal));
+                }
+        };
+        class BoolIntersection
+        {
+            public:
+                BoolIntersection() {};
+                BoolIntersection(const Intersection & intersection, bool left)
+                {
+                    this->intersection = intersection;
+                    this->left = left;
+                }
+                Intersection intersection;
+                bool left;
+        };
+        class IntersectionList
+        {
+            public:
+                void add(const Intersection & i)
+                {
+                    m_intersections.push_back(i);
+                }
+                Intersection & operator[](int i)
+                {
+                    return m_intersections[i];
+                }
+                const Intersection & operator[](int i) const
+                {
+                    return m_intersections[i];
+                }
+                size_t size() const { return m_intersections.size(); }
+                std::vector<Intersection>::iterator begin()
+                {
+                    return m_intersections.begin();
+                }
+                std::vector<Intersection>::iterator end()
+                {
+                    return m_intersections.end();
+                }
+            protected:
+                std::vector< Intersection > m_intersections;
+        };
+        class BoolIntersectionList
+        {
+            public:
+                BoolIntersectionList(const IntersectionList & l1,
+                                     const IntersectionList & l2,
+                                     const Vector & startPoint);
+                BoolIntersection & operator[](int i)
+                {
+                    return m_intersections[i];
+                }
+                const BoolIntersection & operator[](int i) const
+                {
+                    return m_intersections[i];
+                }
+                size_t size() const { return m_intersections.size(); }
+                std::vector<BoolIntersection>::iterator begin()
+                {
+                    return m_intersections.begin();
+                }
+                std::vector<BoolIntersection>::iterator end()
+                {
+                    return m_intersections.end();
+                }
+
+            protected:
+                std::vector< BoolIntersection > m_intersections;
+        };
+
+        Shape();
+        virtual ~Shape();
+        virtual IntersectionList intersect(refptr<Shape> _this,
+                                           const Ray & ray) = 0;
+        virtual refptr<Shape> clone() = 0;
+
+        void setTransform(const Transform & t)
+        {
+            m_transform = t;
+            m_inverse = m_transform.getInverse();
+        }
+        Transform & getTransform() { return m_transform; }
+
+        virtual void setMaterial(refptr<Material> material);
+        refptr<Material> getMaterial() const { return m_material; }
+
+    protected:
+        Transform m_transform;
+        Transform m_inverse;
+        refptr<Material> m_material;
+};
+
+typedef refptr<Shape> ShapeRef;
+
+#endif
+

src/shapes/Sphere.cc

@@ -0,0 +1,49 @@
+
+#include "Sphere.h"
+#include "util/Solver.h"
+#include <math.h>
+#include <iostream>
+using namespace std;
+
+Sphere::Sphere(double radius)
+{
+    m_radius = radius;
+    m_radius2 = radius * radius;
+}
+
+Shape::IntersectionList Sphere::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+
+    IntersectionList res;
+
+    QuadraticSolver solver(1.0,
+                           2 * ( ray_inv.getOrigin()[0] * ray_inv.getDirection()[0]
+                               + ray_inv.getOrigin()[1] * ray_inv.getDirection()[1]
+                               + ray_inv.getOrigin()[2] * ray_inv.getDirection()[2] ),
+                             ray_inv.getOrigin()[0] * ray_inv.getOrigin()[0]
+                           + ray_inv.getOrigin()[1] * ray_inv.getOrigin()[1]
+                           + ray_inv.getOrigin()[2] * ray_inv.getOrigin()[2]
+                           - m_radius2);
+    Solver::Result quadSolutions = solver.solve();
+    for (int i = 0; i < quadSolutions.numResults; i++)
+    {
+        if (quadSolutions.results[i] >= 0.0)
+        {
+            Vector isect_point = ray_inv[quadSolutions.results[i]];
+            Vector normal = isect_point;
+            normal.normalize();
+            res.add(
+                Intersection(_this,
+                             m_transform.transform_point(isect_point),
+                             m_transform.transform_normal(normal))
+            );
+        }
+    }
+    return res;
+}
+
+refptr<Shape> Sphere::clone()
+{
+    return new Sphere(*this);
+}

src/shapes/Sphere.h

@@ -0,0 +1,20 @@
+
+#ifndef SPHERE_H
+#define SPHERE_H SPHERE_H
+
+#include "Shape.h"
+
+class Sphere : public Shape
+{
+    public:
+        Sphere(double radius);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        virtual refptr<Shape> clone();
+
+    protected:
+        double m_radius;
+        double m_radius2;
+};
+
+#endif
+

src/shapes/Subtract.cc

@@ -0,0 +1,115 @@
+
+#include "Subtract.h"
+#include <iostream>
+using namespace std;
+
+Subtract::Subtract(const vector< refptr<Shape> > & shapes)
+{
+    int num_shapes = shapes.size();
+    if (num_shapes > 2)
+    {
+        m_shape2 = shapes[num_shapes - 1];
+        vector< refptr<Shape> > rest = shapes;
+        rest.pop_back();
+        m_shape1 = new Subtract(rest);
+    }
+    else if (num_shapes == 2)
+    {
+        m_shape1 = shapes[0];
+        m_shape2 = shapes[1];
+    }
+    else
+    {
+        cerr << __FILE__ << ": " << __LINE__
+             << ": error: Subtract::Subtract() called with only "
+             << num_shapes
+             << " sub-shapes!"
+             << endl;
+        exit(4);
+    }
+}
+
+Shape::IntersectionList Subtract::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+    IntersectionList res1 = m_shape1->intersect(m_shape1, ray_inv);
+
+    /*** optimization ***/
+    if (res1.size() == 0)
+        return res1;
+    /*** end optimization ***/
+
+    IntersectionList res2 = m_shape2->intersect(m_shape2, ray_inv);
+    BoolIntersectionList merged(res1, res2, ray_inv.getOrigin());
+
+    IntersectionList res;
+    bool in1 = false, in2 = false;
+
+    /* initially go through the merged intersections to see whether
+     * the ray started inside one of the sub-objects */
+    for (int i = 0, sz = merged.size(), saw1 = 0, saw2 = 0;
+         i < sz && (!saw1 || !saw2);
+         i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool back = dot < 0.0;
+        bool left = merged[i].left;
+        if (back)
+        {
+            if (left && !saw1)
+                in1 = true;
+            else if (!left && !saw2)
+                in2 = true;
+        }
+        if (left)
+            saw1 = 1;
+        else
+            saw2 = 1;
+    }
+
+    bool in_bool = in1 && !in2;
+    for (int i = 0, sz = merged.size(); i < sz; i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool front = dot > 0.0;
+        bool left = merged[i].left;
+        if (left)
+            in1 = front;
+        else
+            in2 = front;
+        if (!in_bool && in1 && !in2)
+        {
+            /* we found an intersection point
+             * to get into the boolean object */
+            in_bool = true;
+            BoolIntersection bi = merged[i];
+            Intersection i = bi.intersection;
+            if ( ! left )   /* if this point came from object B (A - B) */
+                i.normal = - i.normal;
+            res.add(i.transform(m_transform));
+        }
+        else if (in_bool && !(in1 && !in2))
+        {
+            /* we found an intersection point
+             * to get out of the boolean object */
+            BoolIntersection bi = merged[i];
+            Intersection i = bi.intersection;
+            if ( ! left )   /* if this point came from object B (A - B) */
+                i.normal = - i.normal;
+            res.add(i.transform(m_transform));
+            in_bool = false;
+        }
+    }
+
+    return res;
+}
+
+refptr<Shape> Subtract::clone()
+{
+    Subtract * s = new Subtract(*this);
+    s->m_shape1 = m_shape1->clone();
+    s->m_shape2 = m_shape2->clone();
+    return refptr<Shape>(s);
+}

src/shapes/Subtract.h

@@ -0,0 +1,17 @@
+
+#ifndef SUBTRACT_H
+#define SUBTRACT_H SUBTRACT_H
+
+#include "BoolShape.h"
+#include <vector>
+
+class Subtract : public BoolShape
+{
+    public:
+        Subtract(const std::vector< refptr<Shape> > & shapes);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        refptr<Shape> clone();
+};
+
+#endif
+

src/shapes/Union.cc

@@ -0,0 +1,110 @@
+
+#include "Union.h"
+#include <iostream>
+#include <vector>
+using namespace std;
+
+Union::Union(const vector< refptr<Shape> > & shapes)
+{
+    int num_shapes = shapes.size();
+    if (num_shapes > 2)
+    {
+        m_shape2 = shapes[num_shapes - 1];
+        vector< refptr<Shape> > rest = shapes;
+        rest.pop_back();
+        m_shape1 = new Union(rest);
+    }
+    else if (num_shapes == 2)
+    {
+        m_shape1 = shapes[0];
+        m_shape2 = shapes[1];
+    }
+    else
+    {
+        cerr << __FILE__ << ": " << __LINE__
+             << ": error: Union::Union() called with only "
+             << num_shapes
+             << " sub-shapes!"
+             << endl;
+        exit(4);
+    }
+}
+
+Shape::IntersectionList Union::intersect(refptr<Shape> _this, const Ray & ray)
+{
+    Ray ray_inv = m_inverse.transform_ray(ray);
+    IntersectionList res1 = m_shape1->intersect(m_shape1, ray_inv);
+    IntersectionList res2 = m_shape2->intersect(m_shape2, ray_inv);
+    BoolIntersectionList merged(res1, res2, ray_inv.getOrigin());
+
+    IntersectionList res;
+    bool in1 = false, in2 = false;
+
+    /* initially go through the merged intersections to see whether
+     * the ray started inside one of the sub-objects */
+    for (int i = 0, sz = merged.size(), saw1 = 0, saw2 = 0;
+         i < sz && (!saw1 || !saw2);
+         i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool back = dot < 0.0;
+        bool left = merged[i].left;
+        if (back)
+        {
+            if (left && !saw1)
+                in1 = true;
+            else if (!left && !saw2)
+                in2 = true;
+        }
+        if (left)
+            saw1 = 1;
+        else
+            saw2 = 1;
+    }
+
+    bool in_bool = in1 || in2;
+    for (int i = 0, sz = merged.size(); i < sz; i++)
+    {
+        Vector normal = merged[i].intersection.normal;
+        double dot = - (ray_inv.getDirection() % normal);
+        bool front = dot > 0.0;
+        bool left = merged[i].left;
+        if (front)
+        {
+            if (left)
+                in1 = true;
+            else
+                in2 = true;
+            if (!in_bool && (in1 || in2))
+            {
+                /* we found an intersection point with the boolean object */
+                in_bool = true;
+                res.add( merged[i].intersection.transform(m_transform) );
+            }
+        }
+        else
+        {
+            if (left)
+                in1 = false;
+            else
+                in2 = false;
+            if (in_bool && !(in1 || in2))
+            {
+                /* we found an intersection point with the boolean object */
+                res.add( merged[i].intersection.transform(m_transform) );
+                in_bool = false;
+            }
+        }
+    }
+
+    return res;
+}
+
+refptr<Shape> Union::clone()
+{
+    Union * s = new Union(*this);
+    s->m_shape1 = m_shape1->clone();
+    s->m_shape2 = m_shape2->clone();
+    return refptr<Shape>(s);
+}

src/shapes/Union.h

@@ -0,0 +1,17 @@
+
+#ifndef UNION_H
+#define UNION_H UNION_H
+
+#include "BoolShape.h"
+#include <vector>
+
+class Union : public BoolShape
+{
+    public:
+        Union(const std::vector< refptr<Shape> > & shapes);
+        IntersectionList intersect(refptr<Shape> _this, const Ray & ray);
+        refptr<Shape> clone();
+};
+
+#endif
+

src/shapes/shapes.h

@@ -0,0 +1,10 @@
+
+#include "BoolShape.h"
+#include "Box.h"
+#include "Cyl.h"
+#include "Extrude.h"
+#include "Intersect.h"
+#include "Plane.h"
+#include "Sphere.h"
+#include "Subtract.h"
+#include "Union.h"