Index: bspline.h
===================================================================
--- bspline.h	(revision 0)
+++ bspline.h	(revision 0)
@@ -0,0 +1,52 @@
+/*                           B S P L I N E . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file bspline.h
+ *
+ *
+ */
+
+#ifndef __BSPLINE_H__
+#define __BSPLINE_H__
+
+#include <vector>
+
+#include "nurbs.h"
+
+class BSpline {
+public:
+    int p; // degree
+    int m; // num_knots-1
+    int n; // num_samples-1 (aka number of control points)
+    std::vector<double> params;
+    std::vector<double> knots;
+    ON_2dPointArray controls;
+};
+
+#endif /* __BSPLINE_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: PullbackCurve.cpp
===================================================================
--- PullbackCurve.cpp	(revision 53826)
+++ PullbackCurve.cpp	(working copy)
@@ -42,24 +42,14 @@
 /* interface header */
 #include "PullbackCurve.h"
 
+#include "bspline.h"
 
+
 #define RANGE_HI 0.55
 #define RANGE_LO 0.45
 #define UNIVERSAL_SAMPLE_COUNT 1001
 
 
-/* FIXME: duplicated with opennurbs_ext.cpp */
-class BSpline {
-public:
-    int p; // degree
-    int m; // num_knots-1
-    int n; // num_samples-1 (aka number of control points)
-    std::vector<double> params;
-    std::vector<double> knots;
-    ON_2dPointArray controls;
-};
-
-
 bool
 isFlat(const ON_2dPoint& p1, const ON_2dPoint& m, const ON_2dPoint& p2, double flatness)
 {
Index: gcpdata.cpp
===================================================================
--- gcpdata.cpp	(revision 0)
+++ gcpdata.cpp	(revision 0)
@@ -0,0 +1,181 @@
+/*                           G C P D A T A . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file gcpdata.cpp
+ *
+ *
+ */
+#include "gcpdata.h"
+
+#include "common.h"
+
+#include "bio.h"
+#include <assert.h>
+
+#include "tnt.h"
+#include "jama_lu.h"
+
+#include "vmath.h"
+#include "bu.h"
+
+#include "brep.h"
+#include "dvec.h"
+
+#include "nurbs.h"
+
+#include "bspline.h"
+#include "pbcdata.h"
+
+namespace brlcad {
+    bool
+    gcp_gradient(pt2d_t out_grad, GCPData& data, pt2d_t uv)
+    {
+        bool evaluated = data.surf->Ev2Der(uv[0],
+                                           uv[1],
+                                           data.S,
+                                           data.du,
+                                           data.dv,
+                                           data.duu,
+                                           data.duv,
+                                           data.dvv); // calc S(u, v) dS/du dS/dv d2S/du2 d2S/dv2 d2S/dudv
+        if (!evaluated) return false;
+        out_grad[0] = 2 * (data.du * (data.S - data.pt));
+        out_grad[1] = 2 * (data.dv * (data.S - data.pt));
+        return true;
+    }
+
+
+    bool
+    gcp_newton_iteration(pt2d_t out_uv, GCPData& data, pt2d_t grad, pt2d_t in_uv)
+    {
+        ON_3dVector delta = data.S - data.pt;
+        double g1du = 2 * (data.duu * delta) + 2 * (data.du * data.du);
+        double g2du = 2 * (data.duv * delta) + 2 * (data.dv * data.du);
+        double g1dv = g2du;
+        double g2dv = 2 * (data.dvv * delta) + 2 * (data.dv * data.dv);
+        mat2d_t jacob = { g1du, g1dv,
+                          g2du, g2dv };
+        mat2d_t inv_jacob;
+        if (mat2d_inverse(inv_jacob, jacob)) {
+            pt2d_t tmp;
+            mat2d_pt2d_mul(tmp, inv_jacob, grad);
+            pt2dsub(out_uv, in_uv, tmp);
+            return true;
+        } else {
+            std::cerr << "inverse failed!" << std::endl; // XXX fix the error handling
+            return false;
+        }
+    }
+
+
+    bool
+    get_closest_point(ON_2dPoint& outpt,
+                      ON_BrepFace* face,
+                      const ON_3dPoint& point,
+                      brlcad::SurfaceTree* tree,
+                      double tolerance)
+    {
+        int try_count = 0;
+        bool delete_tree = false;
+        bool found = false;
+        double d_last = DBL_MAX;
+        pt2d_t curr_grad = {0.0, 0.0};
+        pt2d_t new_uv = {0.0, 0.0};
+        GCPData data;
+        data.surf = face->SurfaceOf();
+        data.pt = point;
+
+        TRACE("get_closest_point: " << PT(point));
+
+        // get initial estimate
+        brlcad::SurfaceTree* a_tree = tree;
+        if (a_tree == NULL) {
+            a_tree = new brlcad::SurfaceTree(face);
+            delete_tree = true;
+        }
+        ON_Interval u, v;
+        ON_2dPoint est = a_tree->getClosestPointEstimate(point, u, v);
+        pt2d_t uv = {est[0], est[1]};
+
+        // do the newton iterations
+        // terminates on 1 of 3 conditions:
+        // 1. if the gradient falls below an epsilon (preferred :-)
+        // 2. if the gradient diverges
+        // 3. iterated MAX_FCP_ITERATIONS
+    try_again:
+        int diverge_count = 0;
+        for (int i = 0; i < BREP_MAX_FCP_ITERATIONS; i++) {
+            assert(gcp_gradient(curr_grad, data, uv));
+
+            ON_3dPoint p = data.surf->PointAt(uv[0], uv[1]);
+            double d = p.DistanceTo(point);
+            TRACE("dist: " << d);
+
+            if (NEAR_EQUAL(d, d_last, tolerance)) {
+                found = true; break;
+            } else if (d > d_last) {
+                TRACE("diverged!");
+                diverge_count++;
+            }
+            if (gcp_newton_iteration(new_uv, data, curr_grad, uv)) {
+                move(uv, new_uv);
+            } else {
+                // iteration failed, nudge diagonally
+                uv[0] += VUNITIZE_TOL;
+                uv[1] += VUNITIZE_TOL;
+            }
+            d_last = d;
+        }
+        if (found) {
+            // check to see if we've left the surface domain
+            double l, h;
+            data.surf->GetDomain(0, &l, &h);
+            if (uv[0] < l) uv[0] = l; // clamp if out of range!
+            if (uv[0] > h) uv[0] = h;
+            data.surf->GetDomain(1, &l, &h);
+            if (uv[1] < l) uv[1] = l;
+            if (uv[1] > h) uv[1] = h;
+
+            outpt[0] = uv[0];
+            outpt[1] = uv[1];
+        } else {
+            TRACE("FAILED TO FIND CLOSEST POINT!");
+            // XXX: try the mid point of the domain -- HACK! but it seems to work!?
+            if (try_count == 0) {
+                uv[0] = u.Mid();
+                uv[1] = v.Mid();
+                ++try_count;
+                goto try_again;
+            }
+        }
+
+        if (delete_tree) delete a_tree;
+        return found;
+    }
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: pbcdata.h
===================================================================
--- pbcdata.h	(revision 0)
+++ pbcdata.h	(revision 0)
@@ -0,0 +1,51 @@
+/*                           P B C D A T A . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file pbcdata.h
+ *
+ *
+ */
+
+#ifndef __PBCDATA_H__
+#define __PBCDATA_H__
+
+namespace brlcad {
+    class PBCData {
+    public:
+        double tolerance;
+        double flatness;
+        const ON_Curve* curve;
+        const ON_Surface* surf;
+        ON_BrepFace* face;
+        SurfaceTree* tree;
+        ON_2dPointArray samples;
+    };
+}
+
+#endif /* __PBCDATA_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: triangle.cpp
===================================================================
--- triangle.cpp	(revision 0)
+++ triangle.cpp	(revision 0)
@@ -0,0 +1,132 @@
+/*                           T R I A N G L E . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file triangle.h
+ *
+ *
+ */
+
+#include "triangle.h"
+
+bool
+triangle_intersection(const struct Triangle &TriA, const struct Triangle &TriB, ON_3dPoint &center)
+{
+    ON_Plane plane_a(TriA.A, TriA.B, TriA.C);
+    ON_Plane plane_b(TriB.A, TriB.B, TriB.C);
+    ON_Line intersect;
+    if (!ON_Intersect(plane_a, plane_b, intersect))
+	return false;
+    ON_3dVector line_normal = ON_CrossProduct(plane_a.Normal(), intersect.Direction());
+
+    // dpi - >0: one side of the line, <0: another side, ==0: on the line
+    double dp1 = ON_DotProduct(TriA.A - intersect.from, line_normal);
+    double dp2 = ON_DotProduct(TriA.B - intersect.from, line_normal);
+    double dp3 = ON_DotProduct(TriA.C - intersect.from, line_normal);
+
+    int points_on_one_side = (dp1 >= 0) + (dp2 >= 0) + (dp3 >= 0);
+    if (points_on_one_side == 0 || points_on_one_side == 3)
+	return false;
+
+    double dp4 = ON_DotProduct(TriB.A - intersect.from, line_normal);
+    double dp5 = ON_DotProduct(TriB.B - intersect.from, line_normal);
+    double dp6 = ON_DotProduct(TriB.C - intersect.from, line_normal);
+    points_on_one_side = (dp4 >= 0) + (dp5 >= 0) + (dp6 >= 0);
+    if (points_on_one_side == 0 || points_on_one_side == 3)
+	return false;
+
+    double t[4];
+    int count = 0;
+    double t1, t2;
+
+    // dpi*dpj < 0 - the corresponding points are on different sides
+    // - the line segment between them are intersected by the plane-plane
+    // intersection line
+    if (dp1*dp2 < 0) {
+	intersect.ClosestPointTo(TriA.A, &t1);
+	intersect.ClosestPointTo(TriA.B, &t2);
+	double d1 = TriA.A.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriA.B.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (dp1*dp3 < 0) {
+	intersect.ClosestPointTo(TriA.A, &t1);
+	intersect.ClosestPointTo(TriA.C, &t2);
+	double d1 = TriA.A.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriA.C.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (dp2*dp3 < 0 && count != 2) {
+	intersect.ClosestPointTo(TriA.B, &t1);
+	intersect.ClosestPointTo(TriA.C, &t2);
+	double d1 = TriA.B.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriA.C.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (count != 2)
+	return false;
+    if (dp4*dp5 < 0) {
+	intersect.ClosestPointTo(TriB.A, &t1);
+	intersect.ClosestPointTo(TriB.B, &t2);
+	double d1 = TriB.A.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriB.B.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (dp4*dp6 < 0) {
+	intersect.ClosestPointTo(TriB.A, &t1);
+	intersect.ClosestPointTo(TriB.C, &t2);
+	double d1 = TriB.A.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriB.C.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (dp5*dp6 < 0 && count != 4) {
+	intersect.ClosestPointTo(TriB.B, &t1);
+	intersect.ClosestPointTo(TriB.C, &t2);
+	double d1 = TriB.B.DistanceTo(intersect.PointAt(t1));
+	double d2 = TriB.C.DistanceTo(intersect.PointAt(t2));
+	if (!ZERO(d1 + d2))
+	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
+    }
+    if (count != 4)
+	return false;
+    if (t[0] > t[1])
+	std::swap(t[1], t[0]);
+    if (t[2] > t[3])
+	std::swap(t[3], t[2]);
+    double left = std::max(t[0], t[2]);
+    double right = std::min(t[1], t[3]);
+    if (left > right)
+	return false;
+    center = intersect.PointAt((left+right)/2);
+    return true;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: subsurface.h
===================================================================
--- subsurface.h	(revision 0)
+++ subsurface.h	(revision 0)
@@ -0,0 +1,129 @@
+/*                           S U B S U R F A C E . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file subsurface.h
+ *
+ *
+ */
+
+#ifndef __SUBSURFACE_H__
+#define __SUBSURFACE_H__
+
+#include "nurbs.h"
+
+struct Subsurface {
+private:
+    ON_BoundingBox m_node;
+public:
+    ON_Surface *m_surf;
+    ON_Interval m_u, m_v;
+    Subsurface *m_children[4];
+
+    Subsurface() : m_surf(NULL)
+    {
+	m_children[0] = m_children[1] = m_children[2] = m_children[3] = NULL;
+    }
+    Subsurface(const Subsurface &_ssurf)
+    {
+	m_surf = _ssurf.m_surf->Duplicate();
+	m_u = _ssurf.m_u;
+	m_v = _ssurf.m_v;
+	m_children[0] = m_children[1] = m_children[2] = m_children[3] = NULL;
+	SetBBox(_ssurf.m_node);
+    }
+    ~Subsurface()
+    {
+	for (int i = 0; i < 4; i++) {
+	    if (m_children[i])
+		delete m_children[i];
+	}
+	delete m_surf;
+    }
+    int Split()
+    {
+	for (int i = 0; i < 4; i++)
+	    m_children[i] = new Subsurface();
+	ON_Surface *temp_surf1 = NULL, *temp_surf2 = NULL;
+	ON_BOOL32 ret = true;
+	ret = m_surf->Split(0, m_u.Mid(), temp_surf1, temp_surf2);
+	if (!ret) {
+	    delete temp_surf1;
+	    delete temp_surf2;
+	    return -1;
+	}
+
+	ret = temp_surf1->Split(1, m_v.Mid(), m_children[0]->m_surf, m_children[1]->m_surf);
+	delete temp_surf1;
+	if (!ret) {
+	    delete temp_surf2;
+	    return -1;
+	}
+	m_children[0]->m_u = ON_Interval(m_u.Min(), m_u.Mid());
+	m_children[0]->m_v = ON_Interval(m_v.Min(), m_v.Mid());
+	m_children[0]->SetBBox(m_children[0]->m_surf->BoundingBox());
+	m_children[1]->m_u = ON_Interval(m_u.Min(), m_u.Mid());
+	m_children[1]->m_v = ON_Interval(m_v.Mid(), m_v.Max());
+	m_children[1]->SetBBox(m_children[1]->m_surf->BoundingBox());
+
+	ret = temp_surf2->Split(1, m_v.Mid(), m_children[2]->m_surf, m_children[3]->m_surf);
+	delete temp_surf2;
+	if (!ret)
+	    return -1;
+	m_children[2]->m_u = ON_Interval(m_u.Mid(), m_u.Max());
+	m_children[2]->m_v = ON_Interval(m_v.Min(), m_v.Mid());
+	m_children[2]->SetBBox(m_children[2]->m_surf->BoundingBox());
+	m_children[3]->m_u = ON_Interval(m_u.Mid(), m_u.Max());
+	m_children[3]->m_v = ON_Interval(m_v.Mid(), m_v.Max());
+	m_children[3]->SetBBox(m_children[3]->m_surf->BoundingBox());
+
+	return 0;
+    }
+    void GetBBox(ON_3dPoint &min, ON_3dPoint &max)
+    {
+	min = m_node.m_min;
+	max = m_node.m_max;
+    }
+    void SetBBox(const ON_BoundingBox &bbox)
+    {
+	m_node = bbox;
+	/* Make sure that each dimension of the bounding box is greater than
+	 * ON_ZERO_TOLERANCE.
+	 * It does the same work as building the surface tree when ray tracing
+	 */
+	for (int i = 0; i < 3; i++) {
+	    double d = m_node.m_max[i] - m_node.m_min[i];
+	    if (ON_NearZero(d, ON_ZERO_TOLERANCE)) {
+		m_node.m_min[i] -= 0.001;
+		m_node.m_max[i] += 0.001;
+	    }
+	}
+    }
+};
+
+#endif /* __BSPLINE_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: gcpdata.h
===================================================================
--- gcpdata.h	(revision 0)
+++ gcpdata.h	(revision 0)
@@ -0,0 +1,88 @@
+/*                           G C P D A T A . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file gcpdata.h
+ *
+ *
+ */
+
+#ifndef __GCPDATA_H__
+#define __GCPDATA_H__
+
+#include "common.h"
+
+#include "bio.h"
+#include <assert.h>
+#include <vector>
+
+#include "tnt.h"
+#include "jama_lu.h"
+
+#include "vmath.h"
+#include "bu.h"
+
+#include "brep.h"
+#include "dvec.h"
+
+#include "nurbs.h"
+
+#include "bspline.h"
+#include "pbcdata.h"
+
+namespace brlcad {
+    typedef struct _gcp_data {
+        const ON_Surface* surf;
+        ON_3dPoint pt;
+
+        // for newton iteration
+        ON_3dPoint S;
+        ON_3dVector du;
+        ON_3dVector dv;
+        ON_3dVector duu;
+        ON_3dVector dvv;
+        ON_3dVector duv;
+    } GCPData;
+
+
+    bool
+    gcp_gradient(pt2d_t out_grad, GCPData& data, pt2d_t uv);
+
+
+    bool
+    gcp_newton_iteration(pt2d_t out_uv, GCPData& data, pt2d_t grad, pt2d_t in_uv);
+
+
+    bool
+    get_closest_point(ON_2dPoint& outpt,
+                      ON_BrepFace* face,
+                      const ON_3dPoint& point,
+                      brlcad::SurfaceTree* tree,
+                      double tolerance);
+}
+#endif /* __GCPDATA_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: pointpair.h
===================================================================
--- pointpair.h	(revision 0)
+++ pointpair.h	(revision 0)
@@ -0,0 +1,46 @@
+/*                           P O I N T P A I R . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file pointpair.h
+ *
+ *
+ */
+
+#ifndef __POINTPAIR_H__
+#define __POINTPAIR_H__
+
+struct PointPair {
+    int indexA, indexB;
+    double distance;
+    bool operator < (const PointPair &_pp) const {
+	return distance < _pp.distance;
+    }
+};
+
+#endif /* __POINTPAIR_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: triangle.h
===================================================================
--- triangle.h	(revision 0)
+++ triangle.h	(revision 0)
@@ -0,0 +1,65 @@
+/*                           T R I A N G L E . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2012 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file triangle.h
+ *
+ *
+ */
+
+#ifndef __TRIANGLE_H__
+#define __TRIANGLE_H__
+
+#include "nurbs.h"
+
+struct Triangle {
+    ON_3dPoint A, B, C;
+    inline void CreateFromPoints(ON_3dPoint &p_A, ON_3dPoint &p_B, ON_3dPoint &p_C) {
+	A = p_A;
+	B = p_B;
+	C = p_C;
+    }
+    ON_3dPoint BarycentricCoordinate(ON_3dPoint &pt)
+    {
+	ON_Matrix M(3, 3);
+	M[0][0] = A[0], M[0][1] = B[0], M[0][2] = C[0];
+	M[1][0] = A[1], M[1][1] = B[1], M[1][2] = C[1];
+	M[2][0] = A[2], M[2][1] = B[2], M[2][2] = C[2];
+	M.Invert(0.0);
+	ON_Matrix MX(3, 1);
+	MX[0][0] = pt[0], MX[1][0] = pt[1], MX[2][0] = pt[2];
+	M.Multiply(M, MX);
+	return ON_3dPoint(M[0][0], M[1][0], M[2][0]);
+    }
+};
+
+
+bool
+triangle_intersection(const struct Triangle &TriA, const struct Triangle &TriB, ON_3dPoint &center);
+
+#endif /* __TRIANGLE_H__ */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: opennurbs_ext.cpp
===================================================================
--- opennurbs_ext.cpp	(revision 53826)
+++ opennurbs_ext.cpp	(working copy)
@@ -40,7 +40,14 @@
 
 #include "nurbs.h"
 
+#include "bspline.h"
+#include "pbcdata.h"
+#include "gcpdata.h"
+#include "subsurface.h"
+#include "pointpair.h"
+#include "triangle.h"
 
+
 #define RANGE_HI 0.55
 #define RANGE_LO 0.45
 #define UNIVERSAL_SAMPLE_COUNT 1001
@@ -2499,174 +2506,9 @@
 
 
 //--------------------------------------------------------------------------------
-// get_closest_point implementation
-
-typedef struct _gcp_data {
-    const ON_Surface* surf;
-    ON_3dPoint pt;
-
-    // for newton iteration
-    ON_3dPoint S;
-    ON_3dVector du;
-    ON_3dVector dv;
-    ON_3dVector duu;
-    ON_3dVector dvv;
-    ON_3dVector duv;
-} GCPData;
-
-
-bool
-gcp_gradient(pt2d_t out_grad, GCPData& data, pt2d_t uv)
-{
-    bool evaluated = data.surf->Ev2Der(uv[0],
-				       uv[1],
-				       data.S,
-				       data.du,
-				       data.dv,
-				       data.duu,
-				       data.duv,
-				       data.dvv); // calc S(u, v) dS/du dS/dv d2S/du2 d2S/dv2 d2S/dudv
-    if (!evaluated) return false;
-    out_grad[0] = 2 * (data.du * (data.S - data.pt));
-    out_grad[1] = 2 * (data.dv * (data.S - data.pt));
-    return true;
-}
-
-
-bool
-gcp_newton_iteration(pt2d_t out_uv, GCPData& data, pt2d_t grad, pt2d_t in_uv)
-{
-    ON_3dVector delta = data.S - data.pt;
-    double g1du = 2 * (data.duu * delta) + 2 * (data.du * data.du);
-    double g2du = 2 * (data.duv * delta) + 2 * (data.dv * data.du);
-    double g1dv = g2du;
-    double g2dv = 2 * (data.dvv * delta) + 2 * (data.dv * data.dv);
-    mat2d_t jacob = { g1du, g1dv,
-		      g2du, g2dv };
-    mat2d_t inv_jacob;
-    if (mat2d_inverse(inv_jacob, jacob)) {
-	pt2d_t tmp;
-	mat2d_pt2d_mul(tmp, inv_jacob, grad);
-	pt2dsub(out_uv, in_uv, tmp);
-	return true;
-    } else {
-	std::cerr << "inverse failed!" << std::endl; // XXX fix the error handling
-	return false;
-    }
-}
-
-
-bool
-get_closest_point(ON_2dPoint& outpt,
-		  ON_BrepFace* face,
-		  const ON_3dPoint& point,
-		  SurfaceTree* tree,
-		  double tolerance)
-{
-    int try_count = 0;
-    bool delete_tree = false;
-    bool found = false;
-    double d_last = DBL_MAX;
-    pt2d_t curr_grad = {0.0, 0.0};
-    pt2d_t new_uv = {0.0, 0.0};
-    GCPData data;
-    data.surf = face->SurfaceOf();
-    data.pt = point;
-
-    TRACE("get_closest_point: " << PT(point));
-
-    // get initial estimate
-    SurfaceTree* a_tree = tree;
-    if (a_tree == NULL) {
-	a_tree = new SurfaceTree(face);
-	delete_tree = true;
-    }
-    ON_Interval u, v;
-    ON_2dPoint est = a_tree->getClosestPointEstimate(point, u, v);
-    pt2d_t uv = {est[0], est[1]};
-
-    // do the newton iterations
-    // terminates on 1 of 3 conditions:
-    // 1. if the gradient falls below an epsilon (preferred :-)
-    // 2. if the gradient diverges
-    // 3. iterated MAX_FCP_ITERATIONS
-try_again:
-    int diverge_count = 0;
-    for (int i = 0; i < BREP_MAX_FCP_ITERATIONS; i++) {
-	assert(gcp_gradient(curr_grad, data, uv));
-
-	ON_3dPoint p = data.surf->PointAt(uv[0], uv[1]);
-	double d = p.DistanceTo(point);
-	TRACE("dist: " << d);
-
-	if (NEAR_EQUAL(d, d_last, tolerance)) {
-	    found = true; break;
-	} else if (d > d_last) {
-	    TRACE("diverged!");
-	    diverge_count++;
-	}
-	if (gcp_newton_iteration(new_uv, data, curr_grad, uv)) {
-	    move(uv, new_uv);
-	} else {
-	    // iteration failed, nudge diagonally
-	    uv[0] += VUNITIZE_TOL;
-	    uv[1] += VUNITIZE_TOL;
-	}
-	d_last = d;
-    }
-    if (found) {
-	// check to see if we've left the surface domain
-	double l, h;
-	data.surf->GetDomain(0, &l, &h);
-	if (uv[0] < l) uv[0] = l; // clamp if out of range!
-	if (uv[0] > h) uv[0] = h;
-	data.surf->GetDomain(1, &l, &h);
-	if (uv[1] < l) uv[1] = l;
-	if (uv[1] > h) uv[1] = h;
-
-	outpt[0] = uv[0];
-	outpt[1] = uv[1];
-    } else {
-	TRACE("FAILED TO FIND CLOSEST POINT!");
-	// XXX: try the mid point of the domain -- HACK! but it seems to work!?
-	if (try_count == 0) {
-	    uv[0] = u.Mid();
-	    uv[1] = v.Mid();
-	    ++try_count;
-	    goto try_again;
-	}
-    }
-
-    if (delete_tree) delete a_tree;
-    return found;
-}
-
-
-//--------------------------------------------------------------------------------
 // pullback_curve implementation
-class PBCData {
-public:
-    double tolerance;
-    double flatness;
-    const ON_Curve* curve;
-    const ON_Surface* surf;
-    ON_BrepFace* face;
-    SurfaceTree* tree;
-    ON_2dPointArray samples;
-};
 
 
-class BSpline {
-public:
-    int p; // degree
-    int m; // num_knots-1
-    int n; // num_samples-1 (aka number of control points)
-    std::vector<double> params;
-    std::vector<double> knots;
-    ON_2dPointArray controls;
-};
-
-
 bool
 isFlat(const ON_2dPoint& p1, const ON_2dPoint& m, const ON_2dPoint& p2, double flatness)
 {
@@ -3000,225 +2842,6 @@
 }
 
 
-struct Triangle {
-    ON_3dPoint A, B, C;
-    inline void CreateFromPoints(ON_3dPoint &p_A, ON_3dPoint &p_B, ON_3dPoint &p_C) {
-	A = p_A;
-	B = p_B;
-	C = p_C;
-    }
-    ON_3dPoint BarycentricCoordinate(ON_3dPoint &pt)
-    {
-	ON_Matrix M(3, 3);
-	M[0][0] = A[0], M[0][1] = B[0], M[0][2] = C[0];
-	M[1][0] = A[1], M[1][1] = B[1], M[1][2] = C[1];
-	M[2][0] = A[2], M[2][1] = B[2], M[2][2] = C[2];
-	M.Invert(0.0);
-	ON_Matrix MX(3, 1);
-	MX[0][0] = pt[0], MX[1][0] = pt[1], MX[2][0] = pt[2];
-	M.Multiply(M, MX);
-	return ON_3dPoint(M[0][0], M[1][0], M[2][0]);
-    }
-};
-
-
-bool
-triangle_intersection(const struct Triangle &TriA, const struct Triangle &TriB, ON_3dPoint &center)
-{
-    ON_Plane plane_a(TriA.A, TriA.B, TriA.C);
-    ON_Plane plane_b(TriB.A, TriB.B, TriB.C);
-    ON_Line intersect;
-    if (!ON_Intersect(plane_a, plane_b, intersect))
-	return false;
-    ON_3dVector line_normal = ON_CrossProduct(plane_a.Normal(), intersect.Direction());
-
-    // dpi - >0: one side of the line, <0: another side, ==0: on the line
-    double dp1 = ON_DotProduct(TriA.A - intersect.from, line_normal);
-    double dp2 = ON_DotProduct(TriA.B - intersect.from, line_normal);
-    double dp3 = ON_DotProduct(TriA.C - intersect.from, line_normal);
-
-    int points_on_one_side = (dp1 >= 0) + (dp2 >= 0) + (dp3 >= 0);
-    if (points_on_one_side == 0 || points_on_one_side == 3)
-	return false;
-
-    double dp4 = ON_DotProduct(TriB.A - intersect.from, line_normal);
-    double dp5 = ON_DotProduct(TriB.B - intersect.from, line_normal);
-    double dp6 = ON_DotProduct(TriB.C - intersect.from, line_normal);
-    points_on_one_side = (dp4 >= 0) + (dp5 >= 0) + (dp6 >= 0);
-    if (points_on_one_side == 0 || points_on_one_side == 3)
-	return false;
-
-    double t[4];
-    int count = 0;
-    double t1, t2;
-
-    // dpi*dpj < 0 - the corresponding points are on different sides
-    // - the line segment between them are intersected by the plane-plane
-    // intersection line
-    if (dp1*dp2 < 0) {
-	intersect.ClosestPointTo(TriA.A, &t1);
-	intersect.ClosestPointTo(TriA.B, &t2);
-	double d1 = TriA.A.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriA.B.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (dp1*dp3 < 0) {
-	intersect.ClosestPointTo(TriA.A, &t1);
-	intersect.ClosestPointTo(TriA.C, &t2);
-	double d1 = TriA.A.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriA.C.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (dp2*dp3 < 0 && count != 2) {
-	intersect.ClosestPointTo(TriA.B, &t1);
-	intersect.ClosestPointTo(TriA.C, &t2);
-	double d1 = TriA.B.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriA.C.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (count != 2)
-	return false;
-    if (dp4*dp5 < 0) {
-	intersect.ClosestPointTo(TriB.A, &t1);
-	intersect.ClosestPointTo(TriB.B, &t2);
-	double d1 = TriB.A.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriB.B.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (dp4*dp6 < 0) {
-	intersect.ClosestPointTo(TriB.A, &t1);
-	intersect.ClosestPointTo(TriB.C, &t2);
-	double d1 = TriB.A.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriB.C.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (dp5*dp6 < 0 && count != 4) {
-	intersect.ClosestPointTo(TriB.B, &t1);
-	intersect.ClosestPointTo(TriB.C, &t2);
-	double d1 = TriB.B.DistanceTo(intersect.PointAt(t1));
-	double d2 = TriB.C.DistanceTo(intersect.PointAt(t2));
-	if (!ZERO(d1 + d2))
-	    t[count++] = t1 + d1/(d1+d2)*(t2-t1);
-    }
-    if (count != 4)
-	return false;
-    if (t[0] > t[1])
-	std::swap(t[1], t[0]);
-    if (t[2] > t[3])
-	std::swap(t[3], t[2]);
-    double left = std::max(t[0], t[2]);
-    double right = std::min(t[1], t[3]);
-    if (left > right)
-	return false;
-    center = intersect.PointAt((left+right)/2);
-    return true;
-}
-
-
-struct PointPair {
-    int indexA, indexB;
-    double distance;
-    bool operator < (const PointPair &_pp) const {
-	return distance < _pp.distance;
-    }
-};
-
-
-struct Subsurface {
-private:
-    ON_BoundingBox m_node;
-public:
-    ON_Surface *m_surf;
-    ON_Interval m_u, m_v;
-    Subsurface *m_children[4];
-
-    Subsurface() : m_surf(NULL)
-    {
-	m_children[0] = m_children[1] = m_children[2] = m_children[3] = NULL;
-    }
-    Subsurface(const Subsurface &_ssurf)
-    {
-	m_surf = _ssurf.m_surf->Duplicate();
-	m_u = _ssurf.m_u;
-	m_v = _ssurf.m_v;
-	m_children[0] = m_children[1] = m_children[2] = m_children[3] = NULL;
-	SetBBox(_ssurf.m_node);
-    }
-    ~Subsurface()
-    {
-	for (int i = 0; i < 4; i++) {
-	    if (m_children[i])
-		delete m_children[i];
-	}
-	delete m_surf;
-    }
-    int Split()
-    {
-	for (int i = 0; i < 4; i++)
-	    m_children[i] = new Subsurface();
-	ON_Surface *temp_surf1 = NULL, *temp_surf2 = NULL;
-	ON_BOOL32 ret = true;
-	ret = m_surf->Split(0, m_u.Mid(), temp_surf1, temp_surf2);
-	if (!ret) {
-	    delete temp_surf1;
-	    delete temp_surf2;
-	    return -1;
-	}
-
-	ret = temp_surf1->Split(1, m_v.Mid(), m_children[0]->m_surf, m_children[1]->m_surf);
-	delete temp_surf1;
-	if (!ret) {
-	    delete temp_surf2;
-	    return -1;
-	}
-	m_children[0]->m_u = ON_Interval(m_u.Min(), m_u.Mid());
-	m_children[0]->m_v = ON_Interval(m_v.Min(), m_v.Mid());
-	m_children[0]->SetBBox(m_children[0]->m_surf->BoundingBox());
-	m_children[1]->m_u = ON_Interval(m_u.Min(), m_u.Mid());
-	m_children[1]->m_v = ON_Interval(m_v.Mid(), m_v.Max());
-	m_children[1]->SetBBox(m_children[1]->m_surf->BoundingBox());
-
-	ret = temp_surf2->Split(1, m_v.Mid(), m_children[2]->m_surf, m_children[3]->m_surf);
-	delete temp_surf2;
-	if (!ret)
-	    return -1;
-	m_children[2]->m_u = ON_Interval(m_u.Mid(), m_u.Max());
-	m_children[2]->m_v = ON_Interval(m_v.Min(), m_v.Mid());
-	m_children[2]->SetBBox(m_children[2]->m_surf->BoundingBox());
-	m_children[3]->m_u = ON_Interval(m_u.Mid(), m_u.Max());
-	m_children[3]->m_v = ON_Interval(m_v.Mid(), m_v.Max());
-	m_children[3]->SetBBox(m_children[3]->m_surf->BoundingBox());
-
-	return 0;
-    }
-    void GetBBox(ON_3dPoint &min, ON_3dPoint &max)
-    {
-	min = m_node.m_min;
-	max = m_node.m_max;
-    }
-    void SetBBox(const ON_BoundingBox &bbox)
-    {
-	m_node = bbox;
-	/* Make sure that each dimension of the bounding box is greater than
-	 * ON_ZERO_TOLERANCE.
-	 * It does the same work as building the surface tree when ray tracing
-	 */
-	for (int i = 0; i < 3; i++) {
-	    double d = m_node.m_max[i] - m_node.m_min[i];
-	    if (ON_NearZero(d, ON_ZERO_TOLERANCE)) {
-		m_node.m_min[i] -= 0.001;
-		m_node.m_max[i] += 0.001;
-	    }
-	}
-    }
-};
-
-
 #define INTERSECT_MAX_DEPTH 8
 int
 surface_surface_intersection(const ON_Surface* surfA,
Index: CMakeLists.txt
===================================================================
--- CMakeLists.txt	(revision 53826)
+++ CMakeLists.txt	(working copy)
@@ -22,6 +22,8 @@
   opennurbs_ext.cpp
   #opennurbs_fit.cpp
   PullbackCurve.cpp
+  gcpdata.cpp
+  triangle.cpp
   )
 
 set(libnurbs_ignored_files
