Add constructive solid geometry ops. These work to some extent, but

not very well; I'm doing a b-rep, where the boundaries are complex
polygons, and there's too many special cases. I should probably
replace this with a triangle mesh solution.

[git-p4: depot-paths = "//depot/solvespace/": change = 1731]

dsc.h

@@ -40,6 +40,8 @@ public:
     double x, y, z;
     
     static Vector MakeFrom(double x, double y, double z);
+    static Vector AtIntersectionOfPlanes(Vector n1, double d1,
+                                         Vector n2, double d2);
 
     bool Equals(Vector v);
     Vector Plus(Vector b);

entity.cpp

@@ -528,6 +528,7 @@ void Entity::LineDrawOrGetDistanceOrEdge(Vector a, Vector b) {
     LineDrawOrGetDistance(a, b);
     if(dogd.edges && !construction) {
         SEdge edge;
+        edge.tag = 0;
         edge.a = a; edge.b = b;
         dogd.edges->l.Add(&edge);
     }
@@ -762,7 +763,7 @@ void Entity::DrawOrGetDistance(int order) {
             double thetaa, thetab, dtheta;
             ArcGetAngles(&thetaa, &thetab, &dtheta);
 
-            int i, n = (int)((40*dtheta)/(2*PI));
+            int i, n = (int)((20*dtheta)/(2*PI));
             Vector prev = pa;
             for(i = 1; i <= n; i++) {
                 double theta = thetaa + (dtheta*i)/n;
@@ -783,7 +784,7 @@ void Entity::DrawOrGetDistance(int order) {
             Vector center = SS.GetEntity(point[0])->PointGetNum();
             Vector u = q.RotationU(), v = q.RotationV();
 
-            int i, c = 40;
+            int i, c = 20;
             Vector prev = u.ScaledBy(r).Plus(center);
             for(i = 1; i <= c; i++) {
                 double phi = (2*PI*i)/c;

glhelper.cpp

@@ -158,7 +158,60 @@ void glxFillPolygon(SPolygon *p)
     gluDeleteTess(gt);
 }
 
-void glxMarkPolygonNormal(SPolygon *p) {
+void glxDebugPolygon(SPolygon *p)
+{
+    int i, j;
+    glLineWidth(2);
+    glPointSize(7);
+    glDisable(GL_DEPTH_TEST);
+    for(i = 0; i < p->l.n; i++) {
+        SContour *sc = &(p->l.elem[i]);
+        for(j = 0; j < (sc->l.n-1); j++) {
+            Vector a = (sc->l.elem[j]).p;
+            Vector b = (sc->l.elem[j+1]).p;
+
+            glxLockColorTo(0, 0, 1);
+            Vector d = (a.Minus(b)).WithMagnitude(-0);
+            glBegin(GL_LINES);
+                glxVertex3v(a.Plus(d));
+                glxVertex3v(b.Minus(d));
+            glEnd();
+            glxLockColorTo(1, 0, 0);
+            glBegin(GL_POINTS);
+                glxVertex3v(a.Plus(d));
+                glxVertex3v(b.Minus(d));
+            glEnd();
+        }
+    }
+}
+
+void glxDebugEdgeList(SEdgeList *el)
+{
+    int i;
+    glLineWidth(2);
+    glPointSize(7);
+    glDisable(GL_DEPTH_TEST);
+    for(i = 0; i < el->l.n; i++) {
+        SEdge *se = &(el->l.elem[i]);
+        if(se->tag) continue;
+        Vector a = se->a, b = se->b;
+
+        glxLockColorTo(0, 1, 0);
+        Vector d = (a.Minus(b)).WithMagnitude(-0);
+        glBegin(GL_LINES);
+            glxVertex3v(a.Plus(d));
+            glxVertex3v(b.Minus(d));
+        glEnd();
+        glxLockColorTo(0, 0, 1);
+        glBegin(GL_POINTS);
+            glxVertex3v(a.Plus(d));
+            glxVertex3v(b.Minus(d));
+        glEnd();
+    }
+}
+
+void glxMarkPolygonNormal(SPolygon *p)
+{
     Vector tail = Vector::MakeFrom(0, 0, 0);
     int i, j, cnt = 0;
     // Choose some reasonable center point.

graphicswin.cpp

@@ -112,9 +112,8 @@ void GraphicsWindow::Init(void) {
     showNormals = true;
     showPoints = true;
     showConstraints = true;
-    showSolids = true;
     showHdnLines = false;
-    showSolids = true;
+    showSolids = false;
 
     solving = SOLVE_ALWAYS;
 
@@ -1181,5 +1180,36 @@ void GraphicsWindow::Paint(int w, int h) {
     for(i = 0; i < MAX_SELECTED; i++) {
         selection[i].Draw();
     }
+
+    if(SS.group.n >= 5) {
+        if(1) {
+            SPolyhedron p; ZERO(&p);
+
+            (SS.group.elem[2].polyh).Boolean(&p, 0, 
+                &(SS.group.elem[4].polyh));
+
+            glEnable(GL_LIGHTING);
+            GLfloat vec[] = { 0.3f, 0.3f, 0.3f, 1.0 };
+            glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, vec);
+            GLfloat vec2[] = { 1.0f, 0.3f, 0.3f, 1.0 };
+            glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, vec2);
+            glEnable(GL_LIGHTING);
+            glEnable(GL_DEPTH_TEST);
+            for(i = 0; i < p.l.n; i++) {
+                glxFillPolygon(&(p.l.elem[i]));
+            }
+            glDisable(GL_LIGHTING);
+            glDisable(GL_DEPTH_TEST);
+
+            p.Clear();
+        } else {
+            SPolygon p; ZERO(&p);
+
+            (SS.group.elem[1].polyg).Boolean(&p, SEdgeList::INTERSECT, 
+                &(SS.group.elem[2].polyg));
+            glxDebugPolygon(&p);
+            p.Clear();
+        }
+    }
 }
 

polygon.cpp

@@ -1,8 +1,19 @@
 #include "solvespace.h"
 
+static int ByDouble(const void *av, const void *bv) {
+    const double *a = (const double *)av;
+    const double *b = (const double *)bv;
+    if(*a == *b) {
+        return 0;
+    } else if(*a > *b) {
+        return 1;
+    } else {
+        return -1;
+    }
+}
+
 bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt) {
     dest->Clear();
-    l.ClearTags();
 
     for(;;) {
         Vector first, last;
@@ -42,8 +53,10 @@ bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt) {
             }
             if(i >= l.n) {
                 // Couldn't assemble a closed contour; mark where.
+                if(errorAt) {
                     errorAt->a = first;
                     errorAt->b = last;
+                }
                 return false;
             }
 
@@ -51,6 +64,109 @@ bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt) {
     }
 }
 
+void SEdgeList::CopyBreaking(SEdgeList *dest) {
+    int i, j, k;
+    
+    for(i = 0; i < l.n; i++) {
+        SEdge *ei = &(l.elem[i]);
+        Vector p0i = ei->a;
+        Vector dpi = (ei->b).Minus(ei->a);
+
+        double inter[100];
+        int inters = 0;
+        for(j = 0; j < l.n; j++) {
+            if(i == j) continue;
+            SEdge *ej = &(l.elem[j]);
+            Vector p0j = ej->a;
+            Vector dpj = (ej->b).Minus(ej->a);
+
+            // Find the intersection, if any
+            Vector dn = dpi.Cross(dpj);
+            if(dn.Magnitude() < 0.001) continue; // parallel, non-intersecting
+            Vector dni = dn.Cross(dpi);
+            Vector dnj = dn.Cross(dpj);
+            double tj =  ((p0i.Minus(p0j)).Dot(dni))/(dpj.Dot(dni));
+            double ti = -((p0i.Minus(p0j)).Dot(dnj))/(dpi.Dot(dnj));
+            // could also test for skew, but assume it's all in plane so not
+
+            if(ti <= 0  || ti >= 1) continue;
+            if(tj < -0.001 || tj > 1.001) continue;
+
+            inter[inters++] = ti;
+        }
+        inter[inters++] = 0;
+        inter[inters++] = 1;
+        qsort(inter, inters, sizeof(inter[0]), ByDouble);
+
+        for(k = 1; k < inters; k++) {
+            SEdge ne;
+            ne.tag = 0;
+            ne.a = p0i.Plus(dpi.ScaledBy(inter[k-1]));
+            ne.b = p0i.Plus(dpi.ScaledBy(inter[k]));
+            dest->l.Add(&ne);
+        }
+    }
+}
+
+void SEdgeList::CullDuplicates(void) {
+    int i, j;
+    for(i = 0; i < l.n; i++) {
+        SEdge *se = &(l.elem[i]);
+        if(se->tag) continue;
+
+        if(((se->a).Minus(se->b)).Magnitude() < 0.01) {
+            se->tag = 1;
+            continue;
+        }
+
+        for(j = i+1; j < l.n; j++) {
+            SEdge *st = &(l.elem[j]);
+            if(st->tag) continue;
+
+            if(((se->a).Equals(st->a) && (se->b).Equals(st->b)) ||
+               ((se->a).Equals(st->b) && (se->b).Equals(st->a)))
+            {
+                // This is an exact duplicate, so mark it as unused now.
+                st->tag = 1;
+                break;
+            }
+        }
+    }
+}
+
+bool SEdgeList::BooleanOp(int op, bool inA, bool inB) {
+    if(op == UNION) {
+        return inA || inB;
+    } else if(op == DIFF) {
+        return inA && (!inB);
+    } else if(op == INTERSECT) {
+        return inA && inB;
+    } else oops();
+}
+
+void SEdgeList::CullForBoolean(int op, SPolygon *a, SPolygon *b) {
+    int i;
+    for(i = 0; i < l.n; i++) {
+        SEdge *se = &(l.elem[i]);
+        if(se->tag) continue;
+
+        Vector tp = ((se->a).Plus(se->b)).ScaledBy(0.5);
+        Vector nudge = ((se->a).Minus(se->b)).Cross(a->normal);
+        nudge = nudge.WithMagnitude(.01);
+        Vector tp1 = tp.Plus(nudge);
+        Vector tp2 = tp.Minus(nudge);
+       
+        bool inf1 = BooleanOp(op, a->ContainsPoint(tp1), b->ContainsPoint(tp1));
+        bool inf2 = BooleanOp(op, a->ContainsPoint(tp2), b->ContainsPoint(tp2));
+
+        if((inf1 && inf2) || (!inf1 && !inf2)) {
+            // The "in polygon" state doesn't change as you cross the edge;
+            // so it doesn't lie on the output polygon.
+            se->tag = 1;
+        }
+    }
+}
+
 void SPolygon::Clear(void) {
     int i;
     for(i = 0; i < l.n; i++) {
@@ -88,6 +204,92 @@ Vector SPolygon::ComputeNormal(void) {
     return (l.elem[0]).ComputeNormal();
 }
 
+bool SPolygon::ContainsPoint(Vector p) {
+    bool inside = false;
+    int i;
+    for(i = 0; i < l.n; i++) {
+        SContour *sc = &(l.elem[i]);
+        if(sc->ContainsPointProjdToNormal(normal, p)) {
+            inside = !inside;
+        }
+    }
+    return inside;
+}
+
+void SPolygon::IntersectAgainstPlane(SEdgeList *dest, Vector p0, Vector n) {
+    if(l.n == 0 || (l.elem[0].l.n == 0)) return;
+    double od = normal.Dot(l.elem[0].l.elem[0].p);
+    double d = n.Dot(p0);
+
+    Vector u = (normal.Cross(n));
+    if(u.Magnitude() < 0.001) {
+        if(n.Dot(normal) < 0) od = -od;
+        if(fabs(od - d) < 0.001) {
+            // The planes are coincident; so the intersection is a copy of
+            // this polygon.
+            MakeEdgesInto(dest);
+        }
+        return;
+    }
+
+    u = u.WithMagnitude(1);
+    Vector v = normal.Cross(u);
+
+    Vector lp = Vector::AtIntersectionOfPlanes(n, d, normal, od);
+    double vp = v.Dot(lp);
+
+    double inter[100];
+    int inters = 0;
+    int i;
+    for(i = 0; i < l.n; i++) {
+        SContour *sc = &(l.elem[i]);
+        // The 0.01 is because I mishandle the case where the intersection
+        // plane goes through a vertex
+        sc->IntersectAgainstPlane(inter, &inters, u, v, vp);
+    }
+    qsort(inter, inters, sizeof(inter[0]), ByDouble);
+    for(i = 0; i < inters; i += 2) {
+        SEdge se;
+        se.tag = 0;
+        se.a = lp.Plus(u.ScaledBy(inter[i]));
+        se.b = lp.Plus(u.ScaledBy(inter[i+1]));
+        dest->l.Add(&se);
+    }
+}
+
+void SPolygon::CopyBreaking(SPolyhedron *dest, SPolyhedron *against, int how) {
+    if(l.n == 0 || (l.elem[0].l.n == 0)) return;
+    Vector p0 = l.elem[0].l.elem[0].p;
+
+    SEdgeList el; ZERO(&el);
+    int i;
+    for(i = 0; i < against->l.n; i++) {
+        SPolygon *pb = &(against->l.elem[i]);
+        pb->IntersectAgainstPlane(&el, p0, normal);
+    }
+    el.CullDuplicates();
+
+    SPolygon inter; ZERO(&inter);
+    bool worked = el.AssemblePolygon(&inter, NULL);
+    inter.normal = normal;
+
+    SPolygon res; ZERO(&res);
+    if(how == 0) {
+        this->Boolean(&res, SEdgeList::DIFF, &inter);
+        res.normal = normal;
+    } else if(how == 1) {
+        this->Boolean(&res, SEdgeList::INTERSECT, &inter);
+        res.normal = normal.ScaledBy(-1);
+    } else oops();
+
+    if(res.l.n > 0) {
+        dest->l.Add(&res);
+    }
+
+    el.l.Clear();
+    inter.Clear();
+}
+
 void SPolygon::FixContourDirections(void) {
     // Outside curve looks counterclockwise, projected against our normal.
     int i, j;
@@ -112,10 +314,37 @@ void SPolygon::FixContourDirections(void) {
     }
 }
 
+bool SPolygon::Boolean(SPolygon *dest, int op, SPolygon *b) {
+    SEdgeList el;
+    ZERO(&el);
+    this->MakeEdgesInto(&el);
+    b->MakeEdgesInto(&el);
+    
+    SEdgeList br;
+    ZERO(&br);
+    el.CopyBreaking(&br);
+
+    br.CullDuplicates();
+    br.CullForBoolean(op, this, b);
+
+    SEdge e;
+    bool ret = br.AssemblePolygon(dest, &e);
+    if(!ret) {
+        br.l.ClearTags();
+        br.CullDuplicates();
+        br.CullForBoolean(op, this, b);
+    }
+
+    br.l.Clear();
+    el.l.Clear();
+    return ret;
+}
+
 void SContour::MakeEdgesInto(SEdgeList *el) {
     int i;
     for(i = 0; i < (l.n-1); i++) {
         SEdge e;
+        e.tag = 0;
         e.a = l.elem[i].p;
         e.b = l.elem[i+1].p;
         el->l.Add(&e);
@@ -133,7 +362,7 @@ Vector SContour::ComputeNormal(void) {
             n = nt;
         }
     }
-    return n;
+    return n.WithMagnitude(1);
 }
 
 bool SContour::IsClockwiseProjdToNormal(Vector n) {
@@ -169,11 +398,31 @@ bool SContour::ContainsPointProjdToNormal(Vector n, Vector p) {
     for(int i = 0; i < (l.n - 1); i++) {
         double ua = (l.elem[i  ].p).Dot(u);
         double va = (l.elem[i  ].p).Dot(v);
-        double ub = (l.elem[i+1].p).Dot(u);
+        // The curve needs to be exactly closed; approximation is death.
-        double vb = (l.elem[i+1].p).Dot(v);
+        double ub = (l.elem[(i+1)%(l.n-1)].p).Dot(u);
+        double vb = (l.elem[(i+1)%(l.n-1)].p).Dot(v);
+
+        if ((((va <= vp) && (vp < vb)) ||
+             ((vb <= vp) && (vp < va))) &&
+            (up < (ub - ua) * (vp - va) / (vb - va) + ua))
+        {
+          inside = !inside;
+        }
+    }
+
+    return inside;
+}
+
+
+void SContour::IntersectAgainstPlane(double *inter, int *inters,
+                                     Vector u, Vector v, double vp)
+{
+    for(int i = 0; i < (l.n - 1); i++) {
+        double ua = (l.elem[i  ].p).Dot(u);
+        double va = (l.elem[i  ].p).Dot(v);
+        double ub = (l.elem[(i+1)%(l.n-1)].p).Dot(u);
+        double vb = (l.elem[(i+1)%(l.n-1)].p).Dot(v);
 
-        // Write the parametric equation of the line, standardized so that
-        // t = 0 has smaller v than t = 1
         double u0, v0, du, dv;
 
         if(va < vb) {
@@ -184,14 +433,15 @@ bool SContour::ContainsPointProjdToNormal(Vector n, Vector p) {
             du = (ua - ub); dv = (va - vb);
         }
 
-        if(dv == 0) continue; // intersects our horiz ray either 0 or 2 times
+        if(dv == 0) continue; 
 
         double t = (vp - v0)/dv;
+        if(t >= 0 && t < 1) {
             double ui = u0 + t*du;
-        if(ui > up && t >= 0 && t < 1) inside = !inside;
+            // Our line v = vp intersects the edge; record the u value
+            inter[(*inters)++] = ui;
+        }
     }
-
-    return inside;
 }
 
 void SContour::Reverse(void) {
@@ -203,3 +453,41 @@ void SContour::Reverse(void) {
         l.elem[i] = t;
     }
 }
+
+void SPolyhedron::AddFace(SPolygon *p) {
+    l.Add(p);
+}
+
+void SPolyhedron::Clear(void) {
+    int i;
+    for(i = 0; i < l.n; i++) {
+        (l.elem[i]).Clear();
+    }
+    l.Clear();
+}
+
+void SPolyhedron::IntersectAgainstPlane(SEdgeList *d, Vector p0, Vector n) {
+    int i;
+    for(i = 0; i < l.n; i++) {
+        SPolygon *sp = &(l.elem[i]);
+        sp->IntersectAgainstPlane(d, p0, n);
+    }
+}
+
+bool SPolyhedron::Boolean(SPolyhedron *dest, int op, SPolyhedron *b) {
+    int i;
+    dbp(">>>");
+
+    for(i = 0; i < l.n; i++) {
+        SPolygon *sp = &(l.elem[i]);
+        sp->CopyBreaking(dest, b, 0);
+    }
+
+    for(i = 0; i < b->l.n; i++) {
+        SPolygon *sp = &(b->l.elem[i]);
+        sp->CopyBreaking(dest, this, 1);
+    }
+
+    return true;
+}
+

polygon.h

@@ -3,6 +3,7 @@
 #define __POLYGON_H
 
 class SPolygon;
+class SPolyhedron;
 
 template <class T>
 class SList {
@@ -44,7 +45,11 @@ public:
     SList<SEdge>    l;
 
     bool AssemblePolygon(SPolygon *dest, SEdge *errorAt);
-    void BreakEdgesInto(SEdgeList *el);
+    void CopyBreaking(SEdgeList *dest);
+    static const int UNION = 0, DIFF = 1, INTERSECT = 2;
+    bool BooleanOp(int op, bool inA, bool inB);
+    void CullForBoolean(int op, SPolygon *a, SPolygon *b);
+    void CullDuplicates(void);
 };
 
 class SPoint {
@@ -62,6 +67,8 @@ public:
     Vector ComputeNormal(void);
     bool IsClockwiseProjdToNormal(Vector n);
     bool ContainsPointProjdToNormal(Vector n, Vector p);
+    void IntersectAgainstPlane(double *inter, int *inters,
+                               Vector u, Vector v, double vp);
 };
 
 class SPolygon {
@@ -72,14 +79,26 @@ public:
     Vector ComputeNormal(void);
     void AddEmptyContour(void);
     void AddPoint(Vector p);
+    bool ContainsPoint(Vector p);
     void MakeEdgesInto(SEdgeList *el);
     void FixContourDirections(void);
     void Clear(void);
+
+    bool Boolean(SPolygon *dest, int op, SPolygon *b);
+
+    void CopyBreaking(SPolyhedron *dest, SPolyhedron *against, int how);
+    void IntersectAgainstPlane(SEdgeList *dest, Vector p0, Vector n);
 };
 
 class SPolyhedron {
-    SList<SPolygon> l;
 public:
+    SList<SPolygon> l;
+
+    void AddFace(SPolygon *p);
+    void Clear(void);
+
+    void IntersectAgainstPlane(SEdgeList *dest, Vector p0, Vector n);
+    bool Boolean(SPolyhedron *dest, int op, SPolyhedron *b);
 };
 
 #endif

sketch.cpp

@@ -383,12 +383,13 @@ void Group::CopyEntity(hEntity in, int a, hParam dx, hParam dy, hParam dz,
 
 void Group::MakePolygons(void) {
     int i;
-    for(i = 0; i < faces.n; i++) {
+    polyh.Clear();
-        (faces.elem[i]).Clear();
+    polyg.Clear();
-    }
+
-    faces.Clear();
+    SEdgeList edges;
+    ZERO(&edges);
+
     if(type == DRAWING_3D || type == DRAWING_WORKPLANE) {
-        edges.l.Clear();
         int i;
         for(i = 0; i < SS.entity.n; i++) {
             Entity *e = &(SS.entity.elem[i]);
@@ -396,17 +397,14 @@ void Group::MakePolygons(void) {
 
             e->GenerateEdges(&edges);
         }
-        SPolygon poly;
-        memset(&poly, 0, sizeof(poly));
         SEdge error;
-        if(edges.AssemblePolygon(&poly, &error)) {
+        if(edges.AssemblePolygon(&polyg, &error)) {
             polyError.yes = false;
-            poly.normal = poly.ComputeNormal();
+            polyg.normal = polyg.ComputeNormal();
-            faces.Add(&poly);
         } else {
             polyError.yes = true;
             polyError.notClosedAt = error;
-            poly.Clear();
+            polyg.Clear();
         }
     } else if(type == EXTRUDE) {
         Vector translate;
@@ -423,57 +421,55 @@ void Group::MakePolygons(void) {
         // Get the source polygon to extrude, and break it down to edges
         edges.l.Clear();
         Group *src = SS.GetGroup(opA);
-        if(src->faces.n != 1) return;
 
-        (src->faces.elem[0]).MakeEdgesInto(&edges);
+        (src->polyg).MakeEdgesInto(&edges);
         for(i = 0; i < edges.l.n; i++) {
             SEdge *edge = &(edges.l.elem[i]);
             edge->a = (edge->a).Plus(t0);
             edge->b = (edge->b).Plus(t0);
         }
 
-        SPolygon poly;
+        SPolygon np;
-        SEdge error;
+        memset(&np, 0, sizeof(np));
-
         // The bottom
-        memset(&poly, 0, sizeof(poly));
+        if(!edges.AssemblePolygon(&np, NULL)) oops();
-        if(!edges.AssemblePolygon(&poly, &error)) oops();
+        Vector n = np.ComputeNormal();
-        Vector n = poly.ComputeNormal();
         if(translate.Dot(n) > 0) {
             n = n.ScaledBy(-1);
         }
-        poly.normal = n;
+        np.normal = n;
-        poly.FixContourDirections();
+        np.FixContourDirections();
-        faces.Add(&poly);
+        polyh.AddFace(&np);
 
         // Regenerate the edges, with the contour directions fixed up.
         edges.l.Clear();
-        poly.MakeEdgesInto(&edges);
+        np.MakeEdgesInto(&edges);
 
         // The sides
         int i;
         for(i = 0; i < edges.l.n; i++) {
             SEdge *edge = &(edges.l.elem[i]);
-            memset(&poly, 0, sizeof(poly));
+            memset(&np, 0, sizeof(np));
-            poly.AddEmptyContour();
+            np.AddEmptyContour();
-            poly.AddPoint(edge->a);
+            np.AddPoint(edge->a);
-            poly.AddPoint(edge->b);
+            np.AddPoint(edge->b);
-            poly.AddPoint((edge->b).Plus(dt));
+            np.AddPoint((edge->b).Plus(dt));
-            poly.AddPoint((edge->a).Plus(dt));
+            np.AddPoint((edge->a).Plus(dt));
-            poly.AddPoint(edge->a);
+            np.AddPoint(edge->a);
-            poly.normal = ((edge->a).Minus(edge->b).Cross(n)).WithMagnitude(1);
+            np.normal = ((edge->a).Minus(edge->b).Cross(n)).WithMagnitude(1);
-            faces.Add(&poly);
+            polyh.AddFace(&np);
 
             edge->a = (edge->a).Plus(dt);
             edge->b = (edge->b).Plus(dt);
         }
 
         // The top
-        memset(&poly, 0, sizeof(poly));
+        memset(&np, 0, sizeof(np));
-        if(!edges.AssemblePolygon(&poly, &error)) oops();
+        if(!edges.AssemblePolygon(&np, NULL)) oops();
-        poly.normal = n.ScaledBy(-1);
+        np.normal = n.ScaledBy(-1);
-        faces.Add(&poly);
+        polyh.AddFace(&np);
     }
+    edges.l.Clear();
 }
 
 void Group::Draw(void) {
@@ -500,18 +496,24 @@ void Group::Draw(void) {
         glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, vec);
         GLfloat vec2[] = { 1.0f, 0.3f, 0.3f, 1.0 };
         glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, vec2);
-        for(i = 0; i < faces.n; i++) {
+        for(i = 0; i < polyh.l.n; i++) {
-            glxFillPolygon(&(faces.elem[i]));
+            glxFillPolygon(&(polyh.l.elem[i]));
 #if 0
             // Debug stuff to show normals to the faces on-screen
             glDisable(GL_LIGHTING);
             glDisable(GL_DEPTH_TEST);
-            glxMarkPolygonNormal(&(faces.elem[i]));
+            glxMarkPolygonNormal(&(polyh.l.elem[i]));
             glEnable(GL_LIGHTING);
             glEnable(GL_DEPTH_TEST);
 #endif
-
         }
+
+        GLfloat vec3[] = { 0.3f, 1.0f, 0.3f, 0.5 };
+        glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, vec3);
+        GLfloat vec4[] = { 1.0f, 0.3f, 0.3f, 0.5 };
+        glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, vec4);
+        glxFillPolygon(&polyg);
+
         glDisable(GL_LIGHTING);
     }
 }

sketch.h

@@ -106,8 +106,8 @@ public:
         bool        negateV;
     } wrkpl;
 
-    SEdgeList       edges;
+    SPolygon        polyg;
-    SList<SPolygon> faces;
+    SPolyhedron     polyh;
     struct {
         SEdge           notClosedAt;
         bool            yes;

solvespace.h

@@ -13,6 +13,7 @@
 #endif
 
 #define SWAP(T, a, b) do { T temp = (a); (a) = (b); (b) = temp; } while(0)
+#define ZERO(v) memset((v), 0, sizeof(*(v)))
 
 #define isforname(c) (isalnum(c) || (c) == '_' || (c) == '-' || (c) == '#')
 
@@ -70,6 +71,8 @@ void vl(void); // debug function to validate
 // Utility functions that are provided in the platform-independent code.
 void glxVertex3v(Vector u);
 void glxFillPolygon(SPolygon *p);
+void glxDebugPolygon(SPolygon *p);
+void glxDebugEdgeList(SEdgeList *l);
 void glxMarkPolygonNormal(SPolygon *p);
 void glxWriteText(char *str);
 void glxWriteTextRefCenter(char *str);

util.cpp

@@ -178,7 +178,7 @@ Vector Vector::MakeFrom(double x, double y, double z) {
 }
 
 bool Vector::Equals(Vector v) {
-    double tol = 0.1;
+    double tol = 0.00001;
     if(fabs(x - v.x) > tol) return false;
     if(fabs(y - v.y) > tol) return false;
     if(fabs(z - v.z) > tol) return false;
@@ -358,6 +358,17 @@ Vector Vector::ClosestOrtho(void) {
     } else oops();
 }
 
+Vector Vector::AtIntersectionOfPlanes(Vector n1, double d1,
+                                      Vector n2, double d2)
+{
+    double det = (n1.Dot(n1))*(n2.Dot(n2)) - 
+                 (n1.Dot(n2))*(n1.Dot(n2));
+    double c1 = (d1*n2.Dot(n2) - d2*n1.Dot(n2))/det;
+    double c2 = (d2*n1.Dot(n1) - d1*n1.Dot(n2))/det;
+    
+    return (n1.ScaledBy(c1)).Plus(n2.ScaledBy(c2));
+}
+
 Point2d Point2d::Plus(Point2d b) {
     Point2d r;
     r.x = x + b.x;