brlcad-commits

[brlcad-commits] SF.net SVN: brlcad:[55433] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55433
          http://sourceforge.net/p/brlcad/code/55433
Author:   phoenixyjll
Date:     2013-05-15 17:55:05 +0000 (Wed, 15 May 2013)
Log Message:
-----------
Use max_dis_2dA and max_dis_2dB for merging polylines.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-05-15 17:38:44 UTC (rev 55432)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-05-15 17:55:05 UTC (rev 55433)
@@ -192,9 +192,9 @@
 
 struct PointPair {
     int indexA, indexB;
-    double distance;
+    double distance3d, distanceA, distanceB;
     bool operator < (const PointPair &_pp) const {
-	return distance < _pp.distance;
+	return distance3d < _pp.distance3d;
     }
 };
 
@@ -504,8 +504,8 @@
 	max_dis = pow(surfA->BoundingBox().Volume()*surfB->BoundingBox().Volume(), 1.0/6.0) * 0.2;
     }
 
-    double max_dis_2dA = ON_2dVector(surfA->Domain(0).Length(), surfA->Domain(1).Length()).Length() * 0.01;
-    double max_dis_2dB = ON_2dVector(surfB->Domain(0).Length(), surfB->Domain(1).Length()).Length() * 0.01;
+    double max_dis_2dA = ON_2dVector(surfA->Domain(0).Length(), surfA->Domain(1).Length()).Length() * 0.1;
+    double max_dis_2dB = ON_2dVector(surfB->Domain(0).Length(), surfB->Domain(1).Length()).Length() * 0.1;
     bu_log("max_dis: %lf\n", max_dis);
     bu_log("max_dis_2dA: %lf\n", max_dis_2dA);
     bu_log("max_dis_2dB: %lf\n", max_dis_2dB);
@@ -515,8 +515,10 @@
     for (int i = 0; i < curvept.Count(); i++) {
 	for (int j = i + 1; j < curvept.Count(); j++) {
 	    PointPair ppair;
-	    ppair.distance = curvept[i].DistanceTo(curvept[j]);
-	    if (ppair.distance < max_dis) {
+	    ppair.distance3d = curvept[i].DistanceTo(curvept[j]);
+	    ppair.distanceA = curveuv[i].DistanceTo(curveuv[j]);
+	    ppair.distanceB = curvest[i].DistanceTo(curvest[j]);
+	    if (ppair.distanceA < max_dis_2dA && ppair.distanceB < max_dis_2dB) {
 		ppair.indexA = i;
 		ppair.indexB = j;
 		ptpairs.push_back(ppair);
@@ -682,3 +684,12 @@
 }
 
 #endif
+
+// Local Variables:
+// tab-width: 8
+// mode: C++
+// c-basic-offset: 4
+// indent-tabs-mode: t
+// c-file-style: "stroustrup"
+// End:
+// ex: shiftwidth=4 tabstop=8

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55637] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55637
          http://sourceforge.net/p/brlcad/code/55637
Author:   phoenixyjll
Date:     2013-06-03 13:48:45 +0000 (Mon, 03 Jun 2013)
Log Message:
-----------
remove the unused variable names to quell compiler warnings.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-03 13:43:43 UTC (rev 55636)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-03 13:48:45 UTC (rev 55637)
@@ -64,11 +64,11 @@
  * Point-curve intersections (PCI)
  */
 bool
-ON_Intersect(const ON_3dPoint& pointA,
-	     const ON_Curve& curveB,
-	     ON_ClassArray<ON_PX_EVENT>& x,
-	     double tolerance,
-	     const ON_Interval* curveB_domain)
+ON_Intersect(const ON_3dPoint&,
+	     const ON_Curve&,
+	     ON_ClassArray<ON_PX_EVENT>&,
+	     double,
+	     const ON_Interval*)
 {
     // Implement later.
     return false;
@@ -78,12 +78,12 @@
  * Point-surface intersections (PSI)
  */
 bool
-ON_Intersect(const ON_3dPoint& pointA,
-	     const ON_Surface& surfaceB,
-	     ON_ClassArray<ON_PX_EVENT>& x,
-	     double tolerance,
-	     const ON_Interval* surfaceB_udomain,
-	     const ON_Interval* surfaceB_vdomain)
+ON_Intersect(const ON_3dPoint&,
+	     const ON_Surface&,
+	     ON_ClassArray<ON_PX_EVENT>&,
+	     double,
+	     const ON_Interval*,
+	     const ON_Interval*)
 {
     // Implement later.
     return false;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55682] brlcad/trunk/src/libbrep/intersect.cpp

[<br...@us...>]

Revision: 55682
          http://sourceforge.net/p/brlcad/code/55682
Author:   brlcad
Date:     2013-06-06 15:59:17 +0000 (Thu, 06 Jun 2013)
Log Message:
-----------
bio.h is a system header wrapper but should still probably come after the <> headers

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-05 22:32:42 UTC (rev 55681)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-06 15:59:17 UTC (rev 55682)
@@ -25,16 +25,17 @@
 
 #include "common.h"
 
-#include "bio.h"
 #include <assert.h>
 #include <vector>
 #include <algorithm>
+#include "bio.h"
 
 #include "vmath.h"
 #include "bu.h"
 
 #include "brep.h"
 
+
 /**
  * Point-point intersections (PPI)
  */

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55692] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55692
          http://sourceforge.net/p/brlcad/code/55692
Author:   phoenixyjll
Date:     2013-06-07 16:59:29 +0000 (Fri, 07 Jun 2013)
Log Message:
-----------
Try to add PS support using get_closest_point().

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-07 15:29:06 UTC (rev 55691)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-07 16:59:29 UTC (rev 55692)
@@ -38,6 +38,14 @@
 
 /**
  * Point-point intersections (PPI)
+ *
+ * approach:
+ *
+ * - Calculate the distance of the two points.
+ *
+ * - If the distance is within the tolerance, the two points
+ *   intersect. Thd mid-point of them is the intersection point,
+ *   and half of the distance is the uncertainty radius.
  */
 bool
 ON_Intersect(const ON_3dPoint& pointA,
@@ -79,14 +87,49 @@
  * Point-surface intersections (PSI)
  */
 bool
-ON_Intersect(const ON_3dPoint&,
-	     const ON_Surface&,
-	     ON_ClassArray<ON_PX_EVENT>&,
-	     double,
-	     const ON_Interval*,
-	     const ON_Interval*)
+ON_Intersect(const ON_3dPoint& pointA,
+	     const ON_Surface& surfaceB,
+	     ON_ClassArray<ON_PX_EVENT>& x,
+	     double tolerance,
+	     const ON_Interval* surfaceB_udomain,
+	     const ON_Interval* surfaceB_vdomain)
 {
-    // Implement later.
+    if (tolerance <= 0.0)
+	tolerance = 0.01;
+
+    ON_Interval u_domain, v_domain;
+    if (surfaceB_udomain == 0)
+	u_domain = surfaceB.Domain(0);
+    else
+	u_domain = *surfaceB_udomain;
+    if (surfaceB_vdomain == 0)
+	v_domain = surfaceB.Domain(1);
+    else
+	v_domain = *surfaceB_vdomain;
+
+    // We need ON_BrepFace for get_closest_point().
+    // TODO: Use routines like SubSurface in SSI to reduce computation.
+    ON_Brep *brep = ON_Brep::New();
+    brep->AddSurface(surfaceB.Duplicate());
+    brep->NewFace(0);
+    ON_2dPoint closest_point_uv;
+    brlcad::get_closest_point(closest_point_uv, brep->Face(0), pointA);
+    delete brep;
+
+    ON_3dPoint closest_point_3d = surfaceB.PointAt(closest_point_uv.x, closest_point_uv.y);
+    bu_log("%lf %lf %lf\n", pointA.x, pointA.y, pointA.z);
+    bu_log("%lf %lf %lf\n", closest_point_3d.x, closest_point_3d.y, closest_point_3d.z);
+    if (pointA.DistanceTo(closest_point_3d) <= tolerance) {
+	ON_PX_EVENT Event;
+	Event.m_type = ON_PX_EVENT::psx_point;
+	Event.m_A = pointA;
+	Event.m_b = closest_point_uv;
+	Event.m_B = closest_point_3d;
+	Event.m_Mid = (pointA + Event.m_B) * 0.5;
+	Event.m_radius = pointA.DistanceTo(Event.m_B) * 0.5;
+	x.Append(Event);
+	return true;
+    }
     return false;
 }
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55715] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55715
          http://sourceforge.net/p/brlcad/code/55715
Author:   phoenixyjll
Date:     2013-06-12 07:12:47 +0000 (Wed, 12 Jun 2013)
Log Message:
-----------
Calculate point-curve intersection using sub-division.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-11 21:55:33 UTC (rev 55714)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-12 07:12:47 UTC (rev 55715)
@@ -36,6 +36,168 @@
 #include "brep.h"
 
 
+/* Sub-division support for a curve.
+ * It's similar to generating the bounding box tree, when the Split()
+ * method is called, the curve is splitted into two parts, whose bounding
+ * boxes become the children of the original curve's bbox.
+ */
+struct Subcurve {
+private:
+    ON_BoundingBox m_node;
+public:
+    ON_Curve *m_curve;
+    ON_Interval m_t;
+    Subcurve *m_children[2];
+    ON_BOOL32 m_islinear;
+
+    Subcurve() : m_curve(NULL), m_islinear(false)
+    {
+	m_children[0] = m_children[1] = NULL;
+    }
+    Subcurve(const Subcurve &_scurve)
+    {
+	m_islinear = _scurve.m_islinear;
+	m_curve = _scurve.m_curve->Duplicate();
+	m_t = _scurve.m_t;
+	m_children[0] = m_children[1] = NULL;
+	SetBBox(_scurve.m_node);
+    }
+    ~Subcurve()
+    {
+	for (int i = 0; i < 2; i++) {
+	    if (m_children[i])
+		delete m_children[i];
+	}
+	delete m_curve;
+    }
+    int Split()
+    {
+	for (int i = 0; i < 2; i++)
+	    m_children[i] = new Subcurve();
+
+	if (m_curve->Split(m_t.Mid(), m_children[0]->m_curve, m_children[1]->m_curve) == false)
+	    return -1;
+
+	m_children[0]->m_t = ON_Interval(m_t.Min(), m_t.Mid());
+	m_children[0]->SetBBox(m_children[0]->m_curve->BoundingBox());
+	m_children[0]->m_islinear = m_children[0]->m_curve->IsLinear();
+	m_children[1]->m_t = ON_Interval(m_t.Mid(), m_t.Max());
+	m_children[1]->SetBBox(m_children[1]->m_curve->BoundingBox());
+	m_children[1]->m_islinear = m_children[1]->m_curve->IsLinear();
+
+	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;
+    }
+    bool IsPointIn(const ON_3dPoint &pt)
+    {
+	return m_node.IsPointIn(pt);
+    }
+};
+
+
+/* Sub-division support for a surface.
+ * It's similar to generating the bounding box tree, when the Split()
+ * method is called, the surface is splitted into two parts, whose bounding
+ * boxes become the children of the original surface's bbox.
+ */
+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;
+	    }
+	}
+    }
+};
+
+
 /**
  * Point-point intersections (PPI)
  *
@@ -69,20 +231,100 @@
     return false;
 }
 
+
+#define MAX_PCI_DEPTH 8
 /**
  * Point-curve intersections (PCI)
  */
 bool
-ON_Intersect(const ON_3dPoint&,
-	     const ON_Curve&,
-	     ON_ClassArray<ON_PX_EVENT>&,
-	     double,
-	     const ON_Interval*)
+ON_Intersect(const ON_3dPoint& pointA,
+	     const ON_Curve& curveB,
+	     ON_ClassArray<ON_PX_EVENT>& x,
+	     double tolerance,
+	     const ON_Interval* curveB_domain)
 {
-    // Implement later.
+    bu_log("PCI called.\n");
+    Subcurve root;
+    if (curveB_domain == NULL) {
+	root.m_curve = curveB.Duplicate();
+	root.m_t = curveB.Domain();
+    }
+    else {
+	// Use ON_Curve::Split() to get the sub-curve inside curveB_domain
+	ON_Curve *temp_curve1 = NULL, *temp_curve2 = NULL, *temp_curve3 = NULL;
+	if (!curveB.Split(curveB_domain->Min(), temp_curve1, temp_curve2))
+	    return false;
+	delete temp_curve1;
+	temp_curve1 = NULL;
+	if (!temp_curve2->Split(curveB_domain->Max(), temp_curve1, temp_curve3))
+	    return false;
+	delete temp_curve1;
+	delete temp_curve2;
+	root.m_curve = temp_curve3;
+	root.m_t = *curveB_domain;
+    }
+
+    root.SetBBox(root.m_curve->BoundingBox());
+    root.m_islinear = root.m_curve->IsLinear();
+
+    if (!root.IsPointIn(pointA))
+	return false;
+
+    std::vector<Subcurve*> candidates, next_candidates;
+    candidates.push_back(&root);
+
+    // Find the sub-curves that are possibly intersected with the point.
+    for (int i = 0; i < MAX_PCI_DEPTH; i++) {
+	next_candidates.clear();
+	for (unsigned int j = 0; j < candidates.size(); j++) {
+	    if (candidates[j]->m_islinear) {
+		next_candidates.push_back(candidates[j]);
+	    } else {
+		if (candidates[j]->Split() == 0) {
+		    if (candidates[j]->m_children[0]->IsPointIn(pointA))
+			next_candidates.push_back(candidates[j]->m_children[0]);
+		    if (candidates[j]->m_children[1]->IsPointIn(pointA))
+			next_candidates.push_back(candidates[j]->m_children[1]);
+		} else
+		    next_candidates.push_back(candidates[j]);
+	    }
+	}
+	candidates = next_candidates;
+    }
+
+    for (unsigned int i = 0; i < candidates.size(); i++) {
+	// Use linear approximation to get the intersection point
+	ON_Line line(candidates[i]->m_curve->PointAtStart(), candidates[i]->m_curve->PointAtEnd());
+	double t, dis;
+	line.ClosestPointTo(pointA, &t);
+	// Get the closest point on the line *segment* to the point
+	double closest_point_t;
+	if (t < 0)
+	    closest_point_t = 0;
+	else if (t > 1)
+	    closest_point_t = 1;
+	else
+	    closest_point_t = t;
+	ON_3dPoint closest_point = line.PointAt(closest_point_t);
+	dis = pointA.DistanceTo(closest_point);
+	
+	if (dis <= tolerance) {
+	    ON_PX_EVENT *Event = new ON_PX_EVENT;
+	    Event->m_type = ON_PX_EVENT::pcx_point;
+	    Event->m_A = pointA;
+	    Event->m_B = closest_point;
+	    Event->m_b[0] = candidates[i]->m_t.Min()+candidates[i]->m_t.Length()*closest_point_t;
+	    Event->m_Mid = (pointA + Event->m_B) * 0.5;
+	    Event->m_radius = dis * 0.5;
+	    x.Append(*Event);
+	    return true;
+	}
+    }
+    // All candidates have no intersection
     return false;
 }
 
+
 /**
  * Point-surface intersections (PSI)
  */
@@ -296,97 +538,7 @@
 };
 
 
-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
+#define MAX_SSI_DEPTH 8
 int
 ON_Intersect(const ON_Surface* surfA,
 	     const ON_Surface* surfB,
@@ -429,7 +581,7 @@
      * method of splitting a NURBS surface.
      * So finally only a small subset of the surface tree is created.
      */
-    for (int h = 0; h <= INTERSECT_MAX_DEPTH; h++) {
+    for (int h = 0; h <= MAX_SSI_DEPTH; h++) {
 	if (nodepairs.empty())
 	    break;
 	NodePairs tmp_pairs;
@@ -444,7 +596,7 @@
 		if (ret1) {
 		    if (ret2) { /* both splits failed */
 			tmp_pairs.push_back(*i);
-			h = INTERSECT_MAX_DEPTH;
+			h = MAX_SSI_DEPTH;
 		    } else { /* the first failed */
 			for (int j = 0; j < 4; j++)
 			    tmp_pairs.push_back(std::make_pair((*i).first, (*i).second->m_children[j]));
@@ -469,7 +621,7 @@
 	    (*i).first->GetBBox(box_a.m_min, box_a.m_max);
 	    (*i).second->GetBBox(box_b.m_min, box_b.m_max);
 	    if (box_intersect.Intersection(box_a, box_b)) {
-		if (h == INTERSECT_MAX_DEPTH) {
+		if (h == MAX_SSI_DEPTH) {
 		    // We have arrived at the bottom of the trees.
 		    // Get an estimate of the intersection point lying on the intersection curve
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55716] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55716
          http://sourceforge.net/p/brlcad/code/55716
Author:   phoenixyjll
Date:     2013-06-12 07:16:14 +0000 (Wed, 12 Jun 2013)
Log Message:
-----------
ws.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-12 07:12:47 UTC (rev 55715)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-12 07:16:14 UTC (rev 55716)
@@ -307,7 +307,7 @@
 	    closest_point_t = t;
 	ON_3dPoint closest_point = line.PointAt(closest_point_t);
 	dis = pointA.DistanceTo(closest_point);
-	
+
 	if (dis <= tolerance) {
 	    ON_PX_EVENT *Event = new ON_PX_EVENT;
 	    Event->m_type = ON_PX_EVENT::pcx_point;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55793] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55793
          http://sourceforge.net/p/brlcad/code/55793
Author:   phoenixyjll
Date:     2013-06-17 14:38:53 +0000 (Mon, 17 Jun 2013)
Log Message:
-----------
The tolerance should be considered in the IsPointIn() test.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-15 17:28:43 UTC (rev 55792)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-17 14:38:53 UTC (rev 55793)
@@ -96,9 +96,11 @@
     {
 	m_node = bbox;
     }
-    bool IsPointIn(const ON_3dPoint &pt)
+    bool IsPointIn(const ON_3dPoint &pt, double tolerance = 0.0)
     {
-	return m_node.IsPointIn(pt);
+	ON_3dVector vtol(tolerance,tolerance,tolerance);
+	ON_BoundingBox new_bbox(m_node.m_min-vtol, m_node.m_max+vtol);
+	return new_bbox.IsPointIn(pt);
     }
 };
 
@@ -243,7 +245,9 @@
 	     double tolerance,
 	     const ON_Interval* curveB_domain)
 {
-    bu_log("PCI called.\n");
+    if (tolerance <= 0.0)
+	tolerance = 0.01;
+
     Subcurve root;
     if (curveB_domain == NULL) {
 	root.m_curve = curveB.Duplicate();
@@ -267,7 +271,7 @@
     root.SetBBox(root.m_curve->BoundingBox());
     root.m_islinear = root.m_curve->IsLinear();
 
-    if (!root.IsPointIn(pointA))
+    if (!root.IsPointIn(pointA, tolerance))
 	return false;
 
     std::vector<Subcurve*> candidates, next_candidates;
@@ -281,9 +285,9 @@
 		next_candidates.push_back(candidates[j]);
 	    } else {
 		if (candidates[j]->Split() == 0) {
-		    if (candidates[j]->m_children[0]->IsPointIn(pointA))
+		    if (candidates[j]->m_children[0]->IsPointIn(pointA, tolerance))
 			next_candidates.push_back(candidates[j]->m_children[0]);
-		    if (candidates[j]->m_children[1]->IsPointIn(pointA))
+		    if (candidates[j]->m_children[1]->IsPointIn(pointA, tolerance))
 			next_candidates.push_back(candidates[j]->m_children[1]);
 		} else
 		    next_candidates.push_back(candidates[j]);

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55795] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55795
          http://sourceforge.net/p/brlcad/code/55795
Author:   phoenixyjll
Date:     2013-06-18 08:40:42 +0000 (Tue, 18 Jun 2013)
Log Message:
-----------
Use macros to represent default tolerance and change it to 0.001 (the same as the default tolerance of curve/curve, curve/surface, surface/surface defined by openNURBS).

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-17 14:43:25 UTC (rev 55794)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-18 08:40:42 UTC (rev 55795)
@@ -218,7 +218,7 @@
 	     double tolerance)
 {
     if (tolerance <= 0.0)
-	tolerance = 0.01;
+	tolerance = ON_ZERO_TOLERANCE;
 
     if (pointA.DistanceTo(pointB) <= tolerance) {
 	ON_PX_EVENT Event;
@@ -234,10 +234,18 @@
 }
 
 
-#define MAX_PCI_DEPTH 8
 /**
  * Point-curve intersections (PCI)
  */
+
+// The maximal depth for subdivision - trade-off between accurancy and
+// performance
+#define MAX_PCI_DEPTH 8
+
+// We make the default tolerance for PSI the same as that of curve-curve
+// intersections defined by openNURBS (see other/openNURBS/opennurbs_curve.h)
+#define PCI_DEFAULT_TOLERANCE 0.001
+
 bool
 ON_Intersect(const ON_3dPoint& pointA,
 	     const ON_Curve& curveB,
@@ -246,7 +254,7 @@
 	     const ON_Interval* curveB_domain)
 {
     if (tolerance <= 0.0)
-	tolerance = 0.01;
+	tolerance = PCI_DEFAULT_TOLERANCE;
 
     Subcurve root;
     if (curveB_domain == NULL) {
@@ -332,6 +340,11 @@
 /**
  * Point-surface intersections (PSI)
  */
+
+// We make the default tolerance for PSI the same as that of curve-surface
+// intersections defined by openNURBS (see other/openNURBS/opennurbs_curve.h)
+#define PSI_DEFAULT_TOLERANCE 0.001
+
 bool
 ON_Intersect(const ON_3dPoint& pointA,
 	     const ON_Surface& surfaceB,
@@ -341,7 +354,7 @@
 	     const ON_Interval* surfaceB_vdomain)
 {
     if (tolerance <= 0.0)
-	tolerance = 0.01;
+	tolerance = PSI_DEFAULT_TOLERANCE;
 
     ON_Interval u_domain, v_domain;
     if (surfaceB_udomain == 0)
@@ -542,7 +555,10 @@
 };
 
 
+// The maximal depth for subdivision - trade-off between accurancy and
+// performance
 #define MAX_SSI_DEPTH 8
+
 int
 ON_Intersect(const ON_Surface* surfA,
 	     const ON_Surface* surfB,

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55796] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55796
          http://sourceforge.net/p/brlcad/code/55796
Author:   phoenixyjll
Date:     2013-06-18 08:48:54 +0000 (Tue, 18 Jun 2013)
Log Message:
-----------
The input u_domain and v_domain should be considered.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-18 08:40:42 UTC (rev 55795)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-18 08:48:54 UTC (rev 55796)
@@ -257,7 +257,7 @@
 	tolerance = PCI_DEFAULT_TOLERANCE;
 
     Subcurve root;
-    if (curveB_domain == NULL) {
+    if (curveB_domain == NULL || *curveB_domain == curveB.Domain()) {
 	root.m_curve = curveB.Duplicate();
 	root.m_t = curveB.Domain();
     }
@@ -378,7 +378,9 @@
     ON_3dPoint closest_point_3d = surfaceB.PointAt(closest_point_uv.x, closest_point_uv.y);
     bu_log("%lf %lf %lf\n", pointA.x, pointA.y, pointA.z);
     bu_log("%lf %lf %lf\n", closest_point_3d.x, closest_point_3d.y, closest_point_3d.z);
-    if (pointA.DistanceTo(closest_point_3d) <= tolerance) {
+    if (pointA.DistanceTo(closest_point_3d) <= tolerance
+	&& u_domain.Includes(closest_point_uv.x)
+	&& v_domain.Includes(closest_point_uv.y)) {
 	ON_PX_EVENT Event;
 	Event.m_type = ON_PX_EVENT::psx_point;
 	Event.m_A = pointA;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55806] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55806
          http://sourceforge.net/p/brlcad/code/55806
Author:   phoenixyjll
Date:     2013-06-20 12:57:54 +0000 (Thu, 20 Jun 2013)
Log Message:
-----------
Use a smaller depth for PSI to improve performance.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-20 12:37:29 UTC (rev 55805)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-20 12:57:54 UTC (rev 55806)
@@ -345,6 +345,12 @@
 // intersections defined by openNURBS (see other/openNURBS/opennurbs_curve.h)
 #define PSI_DEFAULT_TOLERANCE 0.001
 
+// The default maximal depth for creating a surfacee tree (8) is way too
+// much - killing performance. Since it's only used for getting an
+// estimation, and we use Newton iterations afterwards, it's reasonable
+// to use a smaller depth in this step.
+#define MAX_PSI_DEPTH 4
+
 bool
 ON_Intersect(const ON_3dPoint& pointA,
 	     const ON_Surface& surfaceB,
@@ -372,16 +378,17 @@
     brep->AddSurface(surfaceB.Duplicate());
     brep->NewFace(0);
     ON_2dPoint closest_point_uv;
-    if (brlcad::get_closest_point(closest_point_uv, brep->Face(0), pointA) == false) {
+    brlcad::SurfaceTree *tree = new brlcad::SurfaceTree(brep->Face(0), true, MAX_PSI_DEPTH);
+    if (brlcad::get_closest_point(closest_point_uv, brep->Face(0), pointA, tree) == false) {
 	delete brep;
+	delete tree;
 	return false;
     }
 
     delete brep;
+    delete tree;
 
     ON_3dPoint closest_point_3d = surfaceB.PointAt(closest_point_uv.x, closest_point_uv.y);
-    bu_log("%lf %lf %lf\n", pointA.x, pointA.y, pointA.z);
-    bu_log("%lf %lf %lf\n", closest_point_3d.x, closest_point_3d.y, closest_point_3d.z);
     if (pointA.DistanceTo(closest_point_3d) <= tolerance
 	&& u_domain.Includes(closest_point_uv.x)
 	&& v_domain.Includes(closest_point_uv.y)) {

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55810] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55810
          http://sourceforge.net/p/brlcad/code/55810
Author:   phoenixyjll
Date:     2013-06-21 10:37:25 +0000 (Fri, 21 Jun 2013)
Log Message:
-----------
Add Newton-Raphson iteration in PCI to improve accuracy (use the one from subdivision and linear approximation as a starting point)

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-21 00:09:55 UTC (rev 55809)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-21 10:37:25 UTC (rev 55810)
@@ -246,6 +246,11 @@
 // intersections defined by openNURBS (see other/openNURBS/opennurbs_curve.h)
 #define PCI_DEFAULT_TOLERANCE 0.001
 
+// Newton-Raphson iteration needs a good start. Although we can almost get
+// a good starting point every time, but in case of some unfortunate cases,
+// we define a maximum iteration, preventing from infinite loop.
+#define PCI_MAX_ITERATIONS 100
+
 bool
 ON_Intersect(const ON_3dPoint& pointA,
 	     const ON_Curve& curveB,
@@ -305,29 +310,60 @@
     }
 
     for (unsigned int i = 0; i < candidates.size(); i++) {
-	// Use linear approximation to get the intersection point
+	// Use linear approximation to get an estimated intersection point
 	ON_Line line(candidates[i]->m_curve->PointAtStart(), candidates[i]->m_curve->PointAtEnd());
 	double t, dis;
 	line.ClosestPointTo(pointA, &t);
-	// Get the closest point on the line *segment* to the point
-	double closest_point_t;
+
+	// Make sure line_t belongs to [0, 1]
+	double line_t;
 	if (t < 0)
-	    closest_point_t = 0;
+	    line_t = 0;
 	else if (t > 1)
-	    closest_point_t = 1;
+	    line_t = 1;
 	else
-	    closest_point_t = t;
-	ON_3dPoint closest_point = line.PointAt(closest_point_t);
-	dis = pointA.DistanceTo(closest_point);
+	    line_t = t;
 
+	double closest_point_t = candidates[i]->m_t.Min() + candidates[i]->m_t.Length()*line_t;
+	ON_3dPoint closest_point = curveB.PointAt(closest_point_t);
+	dis = closest_point.DistanceTo(pointA);
+
+	// Use Newton-Raphson method with the starting point from linear
+	// approximation.
+	// It's similar to point projection (or point inversion, or get
+	// closest point on curve). We calculate the root of:
+	//       (curve(t) - pointA) \dot tangent(t) = 0 (*)
+	// Then curve(t) may be the closest point on the curve to pointA.
+	//
+	// In some cases, the closest point (intersection point) is the
+	// endpoints of the curve, it does not satisfy (*), but that's
+	// OK because it already comes out with the linear approximation
+	// above.
+
+	double last_t = DBL_MAX;
+	int iterations = 0;
+	while (!ON_NearZero(last_t-closest_point_t)
+	    && dis > tolerance
+	    && iterations++ < PCI_MAX_ITERATIONS) {
+	    ON_3dVector tangent, s_deriv;
+	    last_t = closest_point_t;
+	    curveB.Ev2Der(closest_point_t, closest_point, tangent, s_deriv);
+	    double F = ON_DotProduct(closest_point-pointA, tangent);
+	    double derivF = ON_DotProduct(closest_point-pointA, s_deriv) + ON_DotProduct(tangent, tangent);
+	    if (!ON_NearZero(derivF))
+		closest_point_t -= F/derivF;
+	    dis = closest_point.DistanceTo(pointA);
+	}
+	closest_point = curveB.PointAt(closest_point_t);
+
 	if (dis <= tolerance) {
 	    ON_PX_EVENT *Event = new ON_PX_EVENT;
 	    Event->m_type = ON_PX_EVENT::pcx_point;
 	    Event->m_A = pointA;
 	    Event->m_B = closest_point;
-	    Event->m_b[0] = candidates[i]->m_t.Min()+candidates[i]->m_t.Length()*closest_point_t;
+	    Event->m_b[0] = closest_point_t;
 	    Event->m_Mid = (pointA + Event->m_B) * 0.5;
-	    Event->m_radius = dis * 0.5;
+	    Event->m_radius = closest_point.DistanceTo(pointA) * 0.5;
 	    x.Append(*Event);
 	    return true;
 	}

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55822] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55822
          http://sourceforge.net/p/brlcad/code/55822
Author:   phoenixyjll
Date:     2013-06-24 15:38:39 +0000 (Mon, 24 Jun 2013)
Log Message:
-----------
We have to set removeTrimmed to false, otherwise the surface tree cannot be built correctly.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-24 13:13:46 UTC (rev 55821)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-24 15:38:39 UTC (rev 55822)
@@ -414,7 +414,7 @@
     brep->AddSurface(surfaceB.Duplicate());
     brep->NewFace(0);
     ON_2dPoint closest_point_uv;
-    brlcad::SurfaceTree *tree = new brlcad::SurfaceTree(brep->Face(0), true, MAX_PSI_DEPTH);
+    brlcad::SurfaceTree *tree = new brlcad::SurfaceTree(brep->Face(0), false, MAX_PSI_DEPTH);
     if (brlcad::get_closest_point(closest_point_uv, brep->Face(0), pointA, tree) == false) {
 	delete brep;
 	delete tree;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55834] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55834
          http://sourceforge.net/p/brlcad/code/55834
Author:   phoenixyjll
Date:     2013-06-25 07:32:26 +0000 (Tue, 25 Jun 2013)
Log Message:
-----------
overlap_tolerance is used now.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 06:48:22 UTC (rev 55833)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 07:32:26 UTC (rev 55834)
@@ -542,7 +542,7 @@
 	     const ON_Curve* curveB,
 	     ON_SimpleArray<ON_X_EVENT>& x,
 	     double intersection_tolerance,
-	     double UNUSED(overlap_tolerance),
+	     double overlap_tolerance,
 	     const ON_Interval* curveA_domain,
 	     const ON_Interval* curveB_domain)
 {

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55836] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55836
          http://sourceforge.net/p/brlcad/code/55836
Author:   phoenixyjll
Date:     2013-06-25 08:54:23 +0000 (Tue, 25 Jun 2013)
Log Message:
-----------
Fix a fatal typo.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 07:43:03 UTC (rev 55835)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 08:54:23 UTC (rev 55836)
@@ -576,10 +576,10 @@
 		splittedA.push_back((*i).first->m_children[1]);
 	    }
 	    if ((*i).second->m_islinear || (*i).second->Split() != 0) {
-		splittedA.push_back((*i).second);
+		splittedB.push_back((*i).second);
 	    } else {
-		splittedA.push_back((*i).second->m_children[0]);
-		splittedA.push_back((*i).second->m_children[1]);
+		splittedB.push_back((*i).second->m_children[0]);
+		splittedB.push_back((*i).second->m_children[1]);
 	    }
 	    for (unsigned int j = 0; j < splittedA.size(); j++)
 		for (unsigned int k = 0; k < splittedB.size(); k++)
@@ -643,6 +643,7 @@
 	    ON_3dPoint pointB2 = curveB->PointAt(t_b2);
 	    double distance1 = pointA1.DistanceTo(pointB1);
 	    double distance2 = pointA2.DistanceTo(pointB2);
+
 	    // Check the validity of the solution
 	    if (distance1 < intersection_tolerance && distance2 < intersection_tolerance) {
 		ON_X_EVENT Event;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55837] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55837
          http://sourceforge.net/p/brlcad/code/55837
Author:   phoenixyjll
Date:     2013-06-25 09:00:25 +0000 (Tue, 25 Jun 2013)
Log Message:
-----------
Add logic in case that one side may fail.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 08:54:23 UTC (rev 55836)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 09:00:25 UTC (rev 55837)
@@ -669,6 +669,24 @@
 		else
 		    Event.m_type = ON_X_EVENT::ccx_overlap;
 		x.Append(Event);
+	    } else if (distance1 < intersection_tolerance) {
+		// in case that the second one was not correct
+		ON_X_EVENT Event;
+		Event.m_A[0] = pointA1;
+		Event.m_B[0] = pointB1;
+		Event.m_a[0] = t_a1;
+		Event.m_b[0] = t_b1;
+		Event.m_type = ON_X_EVENT::ccx_point;
+		x.Append(Event);
+	    } else if (distance2 < intersection_tolerance) {
+		// in case that the first one was not correct
+		ON_X_EVENT Event;
+		Event.m_A[0] = pointA2;
+		Event.m_B[0] = pointB2;
+		Event.m_a[0] = t_a2;
+		Event.m_b[0] = t_b2;
+		Event.m_type = ON_X_EVENT::ccx_point;
+		x.Append(Event);
 	    }
 	}
     }

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55838] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55838
          http://sourceforge.net/p/brlcad/code/55838
Author:   phoenixyjll
Date:     2013-06-25 09:40:14 +0000 (Tue, 25 Jun 2013)
Log Message:
-----------
Fix an opposite logic, and use pointers for ON_X_EVENTs.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 09:00:25 UTC (rev 55837)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 09:40:14 UTC (rev 55838)
@@ -524,7 +524,7 @@
 	J[1][1] = -derivB.y;
 	F[0][0] = pointA.x - pointB.x;
 	F[1][0] = pointA.y - pointB.y;
-	if (J.Invert(0.0)) {
+	if (!J.Invert(0.0)) {
 	    // FIXME: More elegant error handling.
 	    bu_log("Inverse failed.\n");
 	    continue;
@@ -626,13 +626,13 @@
 	    double distance = pointA.DistanceTo(pointB);
 	    // Check the validity of the solution
 	    if (distance < intersection_tolerance) {
-		ON_X_EVENT Event;
-		Event.m_A[0] = pointA;
-		Event.m_B[0] = pointB;
-		Event.m_a[0] = t_a1;
-		Event.m_b[0] = t_b1;
-		Event.m_type = ON_X_EVENT::ccx_point;
-		x.Append(Event);
+		ON_X_EVENT *Event = new ON_X_EVENT;
+		Event->m_A[0] = pointA;
+		Event->m_B[0] = pointB;
+		Event->m_a[0] = t_a1;
+		Event->m_b[0] = t_b1;
+		Event->m_type = ON_X_EVENT::ccx_point;
+		x.Append(*Event);
 	    }
 	} else {
 	    // Check overlap
@@ -646,15 +646,15 @@
 
 	    // Check the validity of the solution
 	    if (distance1 < intersection_tolerance && distance2 < intersection_tolerance) {
-		ON_X_EVENT Event;
-		Event.m_A[0] = pointA1;
-		Event.m_A[1] = pointA2;
-		Event.m_B[0] = pointB1;
-		Event.m_B[1] = pointB2;
-		Event.m_a[0] = t_a1;
-		Event.m_a[1] = t_a2;
-		Event.m_b[0] = t_b1;
-		Event.m_b[1] = t_b2;
+		ON_X_EVENT *Event = new ON_X_EVENT;
+		Event->m_A[0] = pointA1;
+		Event->m_A[1] = pointA2;
+		Event->m_B[0] = pointB1;
+		Event->m_B[1] = pointB2;
+		Event->m_a[0] = t_a1;
+		Event->m_a[1] = t_a2;
+		Event->m_b[0] = t_b1;
+		Event->m_b[1] = t_b2;
 		int j;
 		for (j = 1; j < CCI_OVERLAP_TEST_POINTS; j++) {
 		    double strike = 1.0/CCI_OVERLAP_TEST_POINTS;
@@ -665,28 +665,28 @@
 			break;
 		}
 		if (j != CCI_OVERLAP_TEST_POINTS)
-		    Event.m_type = ON_X_EVENT::ccx_point;
+		    Event->m_type = ON_X_EVENT::ccx_point;
 		else
-		    Event.m_type = ON_X_EVENT::ccx_overlap;
-		x.Append(Event);
+		    Event->m_type = ON_X_EVENT::ccx_overlap;
+		x.Append(*Event);
 	    } else if (distance1 < intersection_tolerance) {
 		// in case that the second one was not correct
-		ON_X_EVENT Event;
-		Event.m_A[0] = pointA1;
-		Event.m_B[0] = pointB1;
-		Event.m_a[0] = t_a1;
-		Event.m_b[0] = t_b1;
-		Event.m_type = ON_X_EVENT::ccx_point;
-		x.Append(Event);
+		ON_X_EVENT *Event = new ON_X_EVENT;
+		Event->m_A[0] = pointA1;
+		Event->m_B[0] = pointB1;
+		Event->m_a[0] = t_a1;
+		Event->m_b[0] = t_b1;
+		Event->m_type = ON_X_EVENT::ccx_point;
+		x.Append(*Event);
 	    } else if (distance2 < intersection_tolerance) {
 		// in case that the first one was not correct
-		ON_X_EVENT Event;
-		Event.m_A[0] = pointA2;
-		Event.m_B[0] = pointB2;
-		Event.m_a[0] = t_a2;
-		Event.m_b[0] = t_b2;
-		Event.m_type = ON_X_EVENT::ccx_point;
-		x.Append(Event);
+		ON_X_EVENT *Event = new ON_X_EVENT;
+		Event->m_A[0] = pointA2;
+		Event->m_B[0] = pointB2;
+		Event->m_a[0] = t_a2;
+		Event->m_b[0] = t_b2;
+		Event->m_type = ON_X_EVENT::ccx_point;
+		x.Append(*Event);
 	    }
 	}
     }

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55847] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55847
          http://sourceforge.net/p/brlcad/code/55847
Author:   phoenixyjll
Date:     2013-06-26 02:59:09 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
Add tolerance in the bounding box intersections.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-25 21:52:09 UTC (rev 55846)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 02:59:09 UTC (rev 55847)
@@ -102,10 +102,12 @@
 	ON_BoundingBox new_bbox(m_node.m_min-vtol, m_node.m_max+vtol);
 	return new_bbox.IsPointIn(pt);
     }
-    bool Intersect(const Subcurve& other) const
+    bool Intersect(const Subcurve& other, double tolerance = 0.0) const
     {
+	ON_3dVector vtol(tolerance,tolerance,tolerance);
+	ON_BoundingBox new_bbox(m_node.m_min-vtol, m_node.m_max+vtol);
 	ON_BoundingBox intersection;
-	return intersection.Intersection(m_node, other.m_node);
+	return intersection.Intersection(new_bbox, other.m_node);
     }
 };
 
@@ -559,7 +561,7 @@
 
     typedef std::vector<std::pair<Subcurve*, Subcurve*> > NodePairs;
     NodePairs candidates, next_candidates;
-    if (rootA.Intersect(rootB))
+    if (rootA.Intersect(rootB, intersection_tolerance))
 	candidates.push_back(std::make_pair(&rootA, &rootB));
 
     // Use sub-division and bounding box intersections first.
@@ -583,7 +585,7 @@
 	    }
 	    for (unsigned int j = 0; j < splittedA.size(); j++)
 		for (unsigned int k = 0; k < splittedB.size(); k++)
-		    if (splittedA[j]->Intersect(*splittedB[k]))
+		    if (splittedA[j]->Intersect(*splittedB[k], intersection_tolerance))
 			next_candidates.push_back(std::make_pair(splittedA[j], splittedB[k]));
 	}
 	candidates = next_candidates;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55848] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55848
          http://sourceforge.net/p/brlcad/code/55848
Author:   phoenixyjll
Date:     2013-06-26 03:10:22 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
Check duplication before appending to the array x.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 02:59:09 UTC (rev 55847)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 03:10:22 UTC (rev 55848)
@@ -591,6 +591,7 @@
 	candidates = next_candidates;
     }
 
+    ON_SimpleArray<ON_X_EVENT> tmp_x;
     // For intersected bounding boxes, we calculate an accurate intersection
     // point.
     for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
@@ -634,7 +635,7 @@
 		Event->m_a[0] = t_a1;
 		Event->m_b[0] = t_b1;
 		Event->m_type = ON_X_EVENT::ccx_point;
-		x.Append(*Event);
+		tmp_x.Append(*Event);
 	    }
 	} else {
 	    // Check overlap
@@ -670,7 +671,7 @@
 		    Event->m_type = ON_X_EVENT::ccx_point;
 		else
 		    Event->m_type = ON_X_EVENT::ccx_overlap;
-		x.Append(*Event);
+		tmp_x.Append(*Event);
 	    } else if (distance1 < intersection_tolerance) {
 		// in case that the second one was not correct
 		ON_X_EVENT *Event = new ON_X_EVENT;
@@ -679,7 +680,7 @@
 		Event->m_a[0] = t_a1;
 		Event->m_b[0] = t_b1;
 		Event->m_type = ON_X_EVENT::ccx_point;
-		x.Append(*Event);
+		tmp_x.Append(*Event);
 	    } else if (distance2 < intersection_tolerance) {
 		// in case that the first one was not correct
 		ON_X_EVENT *Event = new ON_X_EVENT;
@@ -688,11 +689,21 @@
 		Event->m_a[0] = t_a2;
 		Event->m_b[0] = t_b2;
 		Event->m_type = ON_X_EVENT::ccx_point;
-		x.Append(*Event);
+		tmp_x.Append(*Event);
 	    }
 	}
     }
 
+    for (int i = 0; i < tmp_x.Count(); i++) {
+	int j;
+	for (j = 0; j < x.Count(); j++) {
+	    if (ON_NearZero(tmp_x[i].m_a[0] - x[j].m_a[0]) && ON_NearZero(tmp_x[i].m_a[1] - x[j].m_a[1])
+		&& ON_NearZero(tmp_x[i].m_b[0] - x[j].m_b[0]) && ON_NearZero(tmp_x[i].m_b[1] - x[j].m_b[1]))
+		break;
+	}
+	if (j == x.Count())
+	    x.Append(tmp_x[i]);
+    }
     return x.Count();
 }
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55849] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55849
          http://sourceforge.net/p/brlcad/code/55849
Author:   phoenixyjll
Date:     2013-06-26 04:22:06 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
If the inverse fails, we try another two directions.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 03:10:22 UTC (rev 55848)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 04:22:06 UTC (rev 55849)
@@ -505,7 +505,7 @@
     //   y_a(t_a) - y_b(t_b) = 0
     //   z_a(t_a) - z_b(t_b) = 0
     // It's an over-determined system.
-    // We use Newton-Raphson iterations to solve the first two equations,
+    // We use Newton-Raphson iterations to solve two equations of the three,
     // and use the third for checking.
     double last_t_a = DBL_MAX*.5, last_t_b = DBL_MAX*.5;
     ON_3dPoint pointA = curveA->PointAt(t_a);
@@ -527,9 +527,26 @@
 	F[0][0] = pointA.x - pointB.x;
 	F[1][0] = pointA.y - pointB.y;
 	if (!J.Invert(0.0)) {
-	    // FIXME: More elegant error handling.
-	    bu_log("Inverse failed.\n");
-	    continue;
+	    // bu_log("Inverse failed.\n");
+	    J[0][0] = derivA.x;
+	    J[0][1] = -derivB.x;
+	    J[1][0] = derivA.z;
+	    J[1][1] = -derivB.z;
+	    F[0][0] = pointA.x - pointB.x;
+	    F[1][0] = pointA.z - pointB.z;
+	    if (!J.Invert(0.0)) {
+		// bu_log("Inverse failed again.\n");
+		J[0][0] = derivA.y;
+		J[0][1] = -derivB.y;
+		J[1][0] = derivA.z;
+		J[1][1] = -derivB.z;
+		F[0][0] = pointA.y - pointB.y;
+		F[1][0] = pointA.z - pointB.z;
+		if (!J.Invert(0.0)) {
+		    // bu_log("Inverse failed and never try again.\n");
+		    continue;
+		}
+	    }
 	}
 	ON_Matrix Delta;
 	Delta.Multiply(J, F);

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55851] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55851
          http://sourceforge.net/p/brlcad/code/55851
Author:   phoenixyjll
Date:     2013-06-26 09:17:26 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
More work on the tolerance value to get a more accurate and correct result.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 09:09:46 UTC (rev 55850)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 09:17:26 UTC (rev 55851)
@@ -498,7 +498,7 @@
 }
 
 void
-newton_cci(double& t_a, double& t_b, const ON_Curve* curveA, const ON_Curve* curveB, double intersection_tolerance)
+newton_cci(double& t_a, double& t_b, const ON_Curve* curveA, const ON_Curve* curveB)
 {
     // Equations:
     //   x_a(t_a) - x_b(t_b) = 0
@@ -513,7 +513,6 @@
     double distance = pointA.DistanceTo(pointB);
     int iteration = 0;
     while (fabs(last_t_a - t_a) + fabs(last_t_b - t_b) > ON_ZERO_TOLERANCE
-	&& distance > intersection_tolerance
 	&& iteration++ < CCI_MAX_ITERATIONS) {
 	last_t_a = t_a, last_t_b = t_b;
 	ON_3dVector derivA, derivB;
@@ -634,12 +633,16 @@
 	// If they converge to one point, it's considered an intersection
 	// point, otherwise it's considered an overlap event.
 	double t_a1 = i->first->m_t.Min(), t_b1 = i->second->m_t.Min();
-	newton_cci(t_a1, t_b1, curveA, curveB, intersection_tolerance);
+	newton_cci(t_a1, t_b1, curveA, curveB);
 	double t_a2 = i->first->m_t.Max(), t_b2 = i->second->m_t.Max();
-	newton_cci(t_a2, t_b2, curveA, curveB, intersection_tolerance);
+	newton_cci(t_a2, t_b2, curveA, curveB);
 
-	if (fabs(t_a1 - t_a2) + fabs(t_b1 - t_b2) < ON_ZERO_TOLERANCE*2) {
-	    // the errors may double, so we use ON_ZERO_TOLERANCE*2
+	ON_3dPoint pointA1 = curveA->PointAt(t_a1);
+	ON_3dPoint pointB1 = curveB->PointAt(t_b1);
+	ON_3dPoint pointA2 = curveA->PointAt(t_a2);
+	ON_3dPoint pointB2 = curveB->PointAt(t_b2);
+	if (pointA1.DistanceTo(pointA2) < intersection_tolerance
+	    && pointB1.DistanceTo(pointB2) < intersection_tolerance) {
 	    // it's considered the same point
 	    ON_3dPoint pointA = curveA->PointAt(t_a1);
 	    ON_3dPoint pointB = curveB->PointAt(t_b1);
@@ -656,11 +659,7 @@
 	    }
 	} else {
 	    // Check overlap
-	    bu_log("Maybe overlap.\n");
-	    ON_3dPoint pointA1 = curveA->PointAt(t_a1);
-	    ON_3dPoint pointB1 = curveB->PointAt(t_b1);
-	    ON_3dPoint pointA2 = curveA->PointAt(t_a2);
-	    ON_3dPoint pointB2 = curveB->PointAt(t_b2);
+	    // bu_log("Maybe overlap.\n");
 	    double distance1 = pointA1.DistanceTo(pointB1);
 	    double distance2 = pointA2.DistanceTo(pointB2);
 
@@ -714,8 +713,10 @@
     for (int i = 0; i < tmp_x.Count(); i++) {
 	int j;
 	for (j = 0; j < x.Count(); j++) {
-	    if (ON_NearZero(tmp_x[i].m_a[0] - x[j].m_a[0]) && ON_NearZero(tmp_x[i].m_a[1] - x[j].m_a[1])
-		&& ON_NearZero(tmp_x[i].m_b[0] - x[j].m_b[0]) && ON_NearZero(tmp_x[i].m_b[1] - x[j].m_b[1]))
+	    if (tmp_x[i].m_A[0].DistanceTo(x[j].m_A[0]) < intersection_tolerance
+		&& tmp_x[i].m_A[1].DistanceTo(x[j].m_A[1]) < intersection_tolerance
+		&& tmp_x[i].m_B[0].DistanceTo(x[j].m_B[0]) < intersection_tolerance
+		&& tmp_x[i].m_B[1].DistanceTo(x[j].m_B[1]) < intersection_tolerance)
 		break;
 	}
 	if (j == x.Count())

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55852] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55852
          http://sourceforge.net/p/brlcad/code/55852
Author:   phoenixyjll
Date:     2013-06-26 09:21:02 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
Remove the set-but-not-used variable distance.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 09:17:26 UTC (rev 55851)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 09:21:02 UTC (rev 55852)
@@ -510,7 +510,6 @@
     double last_t_a = DBL_MAX*.5, last_t_b = DBL_MAX*.5;
     ON_3dPoint pointA = curveA->PointAt(t_a);
     ON_3dPoint pointB = curveB->PointAt(t_b);
-    double distance = pointA.DistanceTo(pointB);
     int iteration = 0;
     while (fabs(last_t_a - t_a) + fabs(last_t_b - t_b) > ON_ZERO_TOLERANCE
 	&& iteration++ < CCI_MAX_ITERATIONS) {
@@ -551,7 +550,6 @@
 	Delta.Multiply(J, F);
 	t_a -= Delta[0][0];
 	t_b -= Delta[1][0];
-	distance = pointA.DistanceTo(pointB);
     }
 }
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55857] brlcad/trunk/src/libbrep/intersect.cpp

[<br...@us...>]

Revision: 55857
          http://sourceforge.net/p/brlcad/code/55857
Author:   brlcad
Date:     2013-06-26 16:01:37 +0000 (Wed, 26 Jun 2013)
Log Message:
-----------
style conformance

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 15:58:30 UTC (rev 55856)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 16:01:37 UTC (rev 55857)
@@ -98,13 +98,13 @@
     }
     bool IsPointIn(const ON_3dPoint &pt, double tolerance = 0.0)
     {
-	ON_3dVector vtol(tolerance,tolerance,tolerance);
+	ON_3dVector vtol(tolerance, tolerance, tolerance);
 	ON_BoundingBox new_bbox(m_node.m_min-vtol, m_node.m_max+vtol);
 	return new_bbox.IsPointIn(pt);
     }
     bool Intersect(const Subcurve& other, double tolerance = 0.0) const
     {
-	ON_3dVector vtol(tolerance,tolerance,tolerance);
+	ON_3dVector vtol(tolerance, tolerance, tolerance);
 	ON_BoundingBox new_bbox(m_node.m_min-vtol, m_node.m_max+vtol);
 	ON_BoundingBox intersection;
 	return intersection.Intersection(new_bbox, other.m_node);
@@ -497,6 +497,7 @@
     return true;
 }
 
+
 void
 newton_cci(double& t_a, double& t_b, const ON_Curve* curveA, const ON_Curve* curveB)
 {
@@ -553,6 +554,7 @@
     }
 }
 
+
 int
 ON_Intersect(const ON_Curve* curveA,
 	     const ON_Curve* curveB,
@@ -963,18 +965,22 @@
 		if ((*i).second->m_children[0] == NULL)
 		    ret2 = (*i).second->Split();
 		if (ret1) {
-		    if (ret2) { /* both splits failed */
+		    if (ret2) {
+			/* both splits failed */
 			tmp_pairs.push_back(*i);
 			h = MAX_SSI_DEPTH;
-		    } else { /* the first failed */
+		    } else {
+			/* the first failed */
 			for (int j = 0; j < 4; j++)
 			    tmp_pairs.push_back(std::make_pair((*i).first, (*i).second->m_children[j]));
 		    }
 		} else {
-		    if (ret2) { /* the second failed */
+		    if (ret2) {
+			/* the second failed */
 			for (int j = 0; j < 4; j++)
 			    tmp_pairs.push_back(std::make_pair((*i).first->m_children[j], (*i).second));
-		    } else { /* both success */
+		    } else {
+			/* both success */
 			for (int j = 0; j < 4; j++)
 			    for (int k = 0; k < 4; k++)
 				tmp_pairs.push_back(std::make_pair((*i).first->m_children[j], (*i).second->m_children[k]));
@@ -1294,6 +1300,7 @@
     return x.Count();
 }
 
+
 // Local Variables:
 // tab-width: 8
 // mode: C++

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55866] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55866
          http://sourceforge.net/p/brlcad/code/55866
Author:   phoenixyjll
Date:     2013-06-27 05:32:46 +0000 (Thu, 27 Jun 2013)
Log Message:
-----------
Merge the overlap events that are continuous, and eliminate the intersection points that are inside the overlap events.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-26 19:07:29 UTC (rev 55865)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-27 05:32:46 UTC (rev 55866)
@@ -552,10 +552,23 @@
 	t_a -= Delta[0][0];
 	t_b -= Delta[1][0];
     }
+
+    // Make sure the solution is inside the domains
+    t_a = std::min(t_a, curveA->Domain().Max());
+    t_a = std::max(t_a, curveA->Domain().Min());
+    t_b = std::min(t_b, curveB->Domain().Max());
+    t_b = std::max(t_b, curveB->Domain().Min());
 }
 
 
 int
+compare_by_m_a0(const ON_X_EVENT* a, const ON_X_EVENT* b)
+{
+    return a->m_a[0] - b->m_a[0];
+}
+
+
+int
 ON_Intersect(const ON_Curve* curveA,
 	     const ON_Curve* curveB,
 	     ON_SimpleArray<ON_X_EVENT>& x,
@@ -666,14 +679,26 @@
 	    // Check the validity of the solution
 	    if (distance1 < intersection_tolerance && distance2 < intersection_tolerance) {
 		ON_X_EVENT *Event = new ON_X_EVENT;
-		Event->m_A[0] = pointA1;
-		Event->m_A[1] = pointA2;
-		Event->m_B[0] = pointB1;
-		Event->m_B[1] = pointB2;
-		Event->m_a[0] = t_a1;
-		Event->m_a[1] = t_a2;
-		Event->m_b[0] = t_b1;
-		Event->m_b[1] = t_b2;
+		// We make sure that m_a[0] <= m_a[1]
+		if (t_a1 <= t_a2) {
+		    Event->m_A[0] = pointA1;
+		    Event->m_A[1] = pointA2;
+		    Event->m_B[0] = pointB1;
+		    Event->m_B[1] = pointB2;
+		    Event->m_a[0] = t_a1;
+		    Event->m_a[1] = t_a2;
+		    Event->m_b[0] = t_b1;
+		    Event->m_b[1] = t_b2;
+		} else {
+		    Event->m_A[0] = pointA2;
+		    Event->m_A[1] = pointA1;
+		    Event->m_B[0] = pointB2;
+		    Event->m_B[1] = pointB1;
+		    Event->m_a[0] = t_a2;
+		    Event->m_a[1] = t_a1;
+		    Event->m_b[0] = t_b2;
+		    Event->m_b[1] = t_b1;
+		}
 		int j;
 		for (j = 1; j < CCI_OVERLAP_TEST_POINTS; j++) {
 		    double strike = 1.0/CCI_OVERLAP_TEST_POINTS;
@@ -710,18 +735,66 @@
 	}
     }
 
+    ON_SimpleArray<ON_X_EVENT> points, overlap, pending;
+    // Use an O(n^2) naive approach to eliminate duplications
     for (int i = 0; i < tmp_x.Count(); i++) {
 	int j;
-	for (j = 0; j < x.Count(); j++) {
-	    if (tmp_x[i].m_A[0].DistanceTo(x[j].m_A[0]) < intersection_tolerance
-		&& tmp_x[i].m_A[1].DistanceTo(x[j].m_A[1]) < intersection_tolerance
-		&& tmp_x[i].m_B[0].DistanceTo(x[j].m_B[0]) < intersection_tolerance
-		&& tmp_x[i].m_B[1].DistanceTo(x[j].m_B[1]) < intersection_tolerance)
+	if (tmp_x[i].m_type == ON_X_EVENT::ccx_overlap) {
+	    overlap.Append(tmp_x[i]);
+	    continue;
+	}
+	// ccx_point
+	for (j = 0; j < points.Count(); j++) {
+	    if (tmp_x[i].m_A[0].DistanceTo(points[j].m_A[0]) < intersection_tolerance
+		&& tmp_x[i].m_A[1].DistanceTo(points[j].m_A[1]) < intersection_tolerance
+		&& tmp_x[i].m_B[0].DistanceTo(points[j].m_B[0]) < intersection_tolerance
+		&& tmp_x[i].m_B[1].DistanceTo(points[j].m_B[1]) < intersection_tolerance)
 		break;
 	}
-	if (j == x.Count())
-	    x.Append(tmp_x[i]);
+	if (j == points.Count()) {
+	    points.Append(tmp_x[i]);
+	}	
     }
+
+    // Merge the overlap events that are continuous
+    overlap.QuickSort(compare_by_m_a0);
+    for (int i = 0; i < overlap.Count(); i++) {
+	bool merged = false;
+	for (int j = 0; j < pending.Count(); j++) {
+	    if (pending[j].m_a[1] < overlap[i].m_a[0] - intersection_tolerance) {
+		x.Append(pending[j]);
+		pending.Remove(j);
+		j--;
+		continue;
+	    }
+	    if (overlap[i].m_a[0] < pending[j].m_a[1] + intersection_tolerance
+		&& overlap[i].m_b[0] < pending[j].m_b[1] + intersection_tolerance) {
+		// Need to merge
+		merged = true;
+		pending[j].m_a[1] = std::max(overlap[i].m_a[1], pending[j].m_a[1]);
+		pending[j].m_b[1] = std::max(overlap[i].m_b[1], pending[j].m_b[1]);
+		break;
+	    }
+	}
+	if (merged == false)
+	    pending.Append(overlap[i]);
+    }
+    for (int i = 0; i < pending.Count(); i++)
+	x.Append(pending[i]);
+
+    // The intersection points shouldn't be inside the overlapped parts.
+    int overlap_events = x.Count();
+    for (int i = 0; i < points.Count(); i++) {
+	int j;
+	for (j = 0; j < overlap_events; j++) {
+	    if (points[i].m_a[0] > x[j].m_a[0] - intersection_tolerance
+		&& points[i].m_a[0] < x[j].m_a[1] + intersection_tolerance)
+		break;
+	}
+	if (j == overlap_events)
+	    x.Append(points[i]);
+    }
+
     return x.Count();
 }
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55869] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55869
          http://sourceforge.net/p/brlcad/code/55869
Author:   phoenixyjll
Date:     2013-06-27 09:30:55 +0000 (Thu, 27 Jun 2013)
Log Message:
-----------
Some special handling for linear curves.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-27 08:35:08 UTC (rev 55868)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-27 09:30:55 UTC (rev 55869)
@@ -624,113 +624,160 @@
     // For intersected bounding boxes, we calculate an accurate intersection
     // point.
     for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
-	/*
-	// First, use linear approximation to get a starting point
-	ON_Line lineA(curveA->PointAt(i->first->m_t.Min()), curveA->PointAt(i->first->m_t.Max()));
-	ON_Line lineB(curveB->PointAt(i->second->m_t.Min()), curveB->PointAt(i->second->m_t.Max()));
-	double t_lineA, t_lineB;
-	double t_a, t_b;
-	if (ON_IntersectLineLine(lineA, lineB, &t_lineA, &t_lineB, ON_ZERO_TOLERANCE, true)) {
-	// The line segments intersect
-	t_a = i->first->m_t.ParameterAt(t_lineA);
-	t_b = i->second->m_t.ParameterAt(t_lineB);
-	} else {
-	// Sometimes the approximated line segments do not intersect,
-	// but the curves DO intersect. We use the mid-points of the
-	// sub-curves as the starting point.
-	t_a = i->first->m_t.Mid();
-	t_b = i->second->m_t.Mid();
-	}*/
+	if (i->first->m_islinear && i->second->m_islinear) {
+	    // Both of them are linear, we use ON_IntersectLineLine
+	    ON_Line lineA(curveA->PointAt(i->first->m_t.Min()), curveA->PointAt(i->first->m_t.Max()));
+	    ON_Line lineB(curveB->PointAt(i->second->m_t.Min()), curveB->PointAt(i->second->m_t.Max()));
+	    if (lineA.Direction().IsParallelTo(lineB.Direction())) {
+		// They are parallel
+		if (lineA.MinimumDistanceTo(lineB) < intersection_tolerance) {
+		    // we report a ccx_overlap event
+		    double startB_on_A, endB_on_A, startA_on_B, endA_on_B;
+		    lineA.ClosestPointTo(lineB.from, &startB_on_A);
+		    lineA.ClosestPointTo(lineB.to, &endB_on_A);
+		    lineB.ClosestPointTo(lineA.from, &startA_on_B);
+		    lineB.ClosestPointTo(lineA.to, &endA_on_B);
 
-	// Use two different start points - the two end-points of the interval
-	// If they converge to one point, it's considered an intersection
-	// point, otherwise it's considered an overlap event.
-	double t_a1 = i->first->m_t.Min(), t_b1 = i->second->m_t.Min();
-	newton_cci(t_a1, t_b1, curveA, curveB);
-	double t_a2 = i->first->m_t.Max(), t_b2 = i->second->m_t.Max();
-	newton_cci(t_a2, t_b2, curveA, curveB);
+		    if (startB_on_A > 1+intersection_tolerance || endB_on_A < -intersection_tolerance
+			|| startA_on_B > 1+intersection_tolerance || endA_on_B < -intersection_tolerance)
+			continue;
 
-	ON_3dPoint pointA1 = curveA->PointAt(t_a1);
-	ON_3dPoint pointB1 = curveB->PointAt(t_b1);
-	ON_3dPoint pointA2 = curveA->PointAt(t_a2);
-	ON_3dPoint pointB2 = curveB->PointAt(t_b2);
-	if (pointA1.DistanceTo(pointA2) < intersection_tolerance
-	    && pointB1.DistanceTo(pointB2) < intersection_tolerance) {
-	    // it's considered the same point
-	    ON_3dPoint pointA = curveA->PointAt(t_a1);
-	    ON_3dPoint pointB = curveB->PointAt(t_b1);
-	    double distance = pointA.DistanceTo(pointB);
-	    // Check the validity of the solution
-	    if (distance < intersection_tolerance) {
-		ON_X_EVENT *Event = new ON_X_EVENT;
-		Event->m_A[0] = pointA;
-		Event->m_B[0] = pointB;
-		Event->m_a[0] = t_a1;
-		Event->m_b[0] = t_b1;
-		Event->m_type = ON_X_EVENT::ccx_point;
-		tmp_x.Append(*Event);
+		    double t_a1, t_a2, t_b1, t_b2;
+		    if (startB_on_A > 0.0)
+			t_a1 = startB_on_A, t_b1 = 0.0;
+		    else
+			t_a1 = 0.0, t_b1 = startA_on_B;
+		    if (endB_on_A < 1.0)
+			t_a2 = endB_on_A, t_b2 = 1.0;
+		    else
+			t_a2 = 1.0, t_b2 = endA_on_B;
+
+		    ON_X_EVENT* Event = new ON_X_EVENT;
+		    Event->m_A[0] = lineA.PointAt(t_a1);
+		    Event->m_A[1] = lineA.PointAt(t_a2);
+		    Event->m_B[0] = lineB.PointAt(t_b1);
+		    Event->m_B[1] = lineB.PointAt(t_b2);
+		    Event->m_a[0] = i->first->m_t.ParameterAt(t_a1);
+		    Event->m_a[1] = i->first->m_t.ParameterAt(t_a2);
+		    Event->m_b[0] = i->second->m_t.ParameterAt(t_b1);
+		    Event->m_b[1] = i->second->m_t.ParameterAt(t_b2);
+		    Event->m_type = ON_X_EVENT::TYPE::ccx_overlap;
+		    tmp_x.Append(*Event);
+		}
+	    } else {
+		// They are not parallel, check intersection point
+		double t_lineA, t_lineB;
+		double t_a, t_b;
+		if (ON_IntersectLineLine(lineA, lineB, &t_lineA, &t_lineB, ON_ZERO_TOLERANCE, true)) {
+		    // The line segments intersect
+		    t_a = i->first->m_t.ParameterAt(t_lineA);
+		    t_b = i->second->m_t.ParameterAt(t_lineB);
+
+		    ON_X_EVENT* Event = new ON_X_EVENT;
+		    Event->m_A[0] = lineA.PointAt(t_lineA);
+		    Event->m_B[0] = lineB.PointAt(t_lineB);
+		    Event->m_a[0] = t_a;
+		    Event->m_b[0] = t_b;
+		    Event->m_type = ON_X_EVENT::TYPE::ccx_point;
+		    tmp_x.Append(*Event);
+		}
 	    }
 	} else {
-	    // Check overlap
-	    // bu_log("Maybe overlap.\n");
-	    double distance1 = pointA1.DistanceTo(pointB1);
-	    double distance2 = pointA2.DistanceTo(pointB2);
+	    // They are not both linear.
 
-	    // Check the validity of the solution
-	    if (distance1 < intersection_tolerance && distance2 < intersection_tolerance) {
-		ON_X_EVENT *Event = new ON_X_EVENT;
-		// We make sure that m_a[0] <= m_a[1]
-		if (t_a1 <= t_a2) {
+	    // Use two different start points - the two end-points of the interval
+	    // If they converge to one point, it's considered an intersection
+	    // point, otherwise it's considered an overlap event.
+	    // FIXME: Find a better machanism to check overlapping, because this method
+	    // may miss some overlap cases. (Overlap events can also converge to one
+	    // point)
+	    double t_a1 = i->first->m_t.Min(), t_b1 = i->second->m_t.Min();
+	    newton_cci(t_a1, t_b1, curveA, curveB);
+	    double t_a2 = i->first->m_t.Max(), t_b2 = i->second->m_t.Max();
+	    newton_cci(t_a2, t_b2, curveA, curveB);
+
+	    ON_3dPoint pointA1 = curveA->PointAt(t_a1);
+	    ON_3dPoint pointB1 = curveB->PointAt(t_b1);
+	    ON_3dPoint pointA2 = curveA->PointAt(t_a2);
+	    ON_3dPoint pointB2 = curveB->PointAt(t_b2);
+	    if (pointA1.DistanceTo(pointA2) < intersection_tolerance
+		&& pointB1.DistanceTo(pointB2) < intersection_tolerance) {
+		// it's considered the same point
+		ON_3dPoint pointA = curveA->PointAt(t_a1);
+		ON_3dPoint pointB = curveB->PointAt(t_b1);
+		double distance = pointA.DistanceTo(pointB);
+		// Check the validity of the solution
+		if (distance < intersection_tolerance) {
+		    ON_X_EVENT *Event = new ON_X_EVENT;
+		    Event->m_A[0] = pointA;
+		    Event->m_B[0] = pointB;
+		    Event->m_a[0] = t_a1;
+		    Event->m_b[0] = t_b1;
+		    Event->m_type = ON_X_EVENT::ccx_point;
+		    tmp_x.Append(*Event);
+		}
+	    } else {
+		// Check overlap
+		// bu_log("Maybe overlap.\n");
+		double distance1 = pointA1.DistanceTo(pointB1);
+		double distance2 = pointA2.DistanceTo(pointB2);
+
+		// Check the validity of the solution
+		if (distance1 < intersection_tolerance && distance2 < intersection_tolerance) {
+		    ON_X_EVENT *Event = new ON_X_EVENT;
+		    // We make sure that m_a[0] <= m_a[1]
+		    if (t_a1 <= t_a2) {
+			Event->m_A[0] = pointA1;
+			Event->m_A[1] = pointA2;
+			Event->m_B[0] = pointB1;
+			Event->m_B[1] = pointB2;
+			Event->m_a[0] = t_a1;
+			Event->m_a[1] = t_a2;
+			Event->m_b[0] = t_b1;
+			Event->m_b[1] = t_b2;
+		    } else {
+			Event->m_A[0] = pointA2;
+			Event->m_A[1] = pointA1;
+			Event->m_B[0] = pointB2;
+			Event->m_B[1] = pointB1;
+			Event->m_a[0] = t_a2;
+			Event->m_a[1] = t_a1;
+			Event->m_b[0] = t_b2;
+			Event->m_b[1] = t_b1;
+		    }
+		    int j;
+		    for (j = 1; j < CCI_OVERLAP_TEST_POINTS; j++) {
+			double strike = 1.0/CCI_OVERLAP_TEST_POINTS;
+			ON_3dPoint test_point = curveA->PointAt(t_a1 + (t_a2-t_a1)*j*strike);
+			ON_ClassArray<ON_PX_EVENT> pci_x;
+			// Use point-curve intersection
+			if (!ON_Intersect(test_point, *curveA, pci_x, overlap_tolerance))
+			    break;
+		    }
+		    if (j != CCI_OVERLAP_TEST_POINTS)
+			Event->m_type = ON_X_EVENT::ccx_point;
+		    else
+			Event->m_type = ON_X_EVENT::ccx_overlap;
+		    tmp_x.Append(*Event);
+		} else if (distance1 < intersection_tolerance) {
+		    // in case that the second one was not correct
+		    ON_X_EVENT *Event = new ON_X_EVENT;
 		    Event->m_A[0] = pointA1;
-		    Event->m_A[1] = pointA2;
 		    Event->m_B[0] = pointB1;
-		    Event->m_B[1] = pointB2;
 		    Event->m_a[0] = t_a1;
-		    Event->m_a[1] = t_a2;
 		    Event->m_b[0] = t_b1;
-		    Event->m_b[1] = t_b2;
-		} else {
+		    Event->m_type = ON_X_EVENT::ccx_point;
+		    tmp_x.Append(*Event);
+		} else if (distance2 < intersection_tolerance) {
+		    // in case that the first one was not correct
+		    ON_X_EVENT *Event = new ON_X_EVENT;
 		    Event->m_A[0] = pointA2;
-		    Event->m_A[1] = pointA1;
 		    Event->m_B[0] = pointB2;
-		    Event->m_B[1] = pointB1;
 		    Event->m_a[0] = t_a2;
-		    Event->m_a[1] = t_a1;
 		    Event->m_b[0] = t_b2;
-		    Event->m_b[1] = t_b1;
+		    Event->m_type = ON_X_EVENT::ccx_point;
+		    tmp_x.Append(*Event);
 		}
-		int j;
-		for (j = 1; j < CCI_OVERLAP_TEST_POINTS; j++) {
-		    double strike = 1.0/CCI_OVERLAP_TEST_POINTS;
-		    ON_3dPoint test_point = curveA->PointAt(t_a1 + (t_a2-t_a1)*j*strike);
-		    ON_ClassArray<ON_PX_EVENT> pci_x;
-		    // Use point-curve intersection
-		    if (!ON_Intersect(test_point, *curveA, pci_x, overlap_tolerance))
-			break;
-		}
-		if (j != CCI_OVERLAP_TEST_POINTS)
-		    Event->m_type = ON_X_EVENT::ccx_point;
-		else
-		    Event->m_type = ON_X_EVENT::ccx_overlap;
-		tmp_x.Append(*Event);
-	    } else if (distance1 < intersection_tolerance) {
-		// in case that the second one was not correct
-		ON_X_EVENT *Event = new ON_X_EVENT;
-		Event->m_A[0] = pointA1;
-		Event->m_B[0] = pointB1;
-		Event->m_a[0] = t_a1;
-		Event->m_b[0] = t_b1;
-		Event->m_type = ON_X_EVENT::ccx_point;
-		tmp_x.Append(*Event);
-	    } else if (distance2 < intersection_tolerance) {
-		// in case that the first one was not correct
-		ON_X_EVENT *Event = new ON_X_EVENT;
-		Event->m_A[0] = pointA2;
-		Event->m_B[0] = pointB2;
-		Event->m_a[0] = t_a2;
-		Event->m_b[0] = t_b2;
-		Event->m_type = ON_X_EVENT::ccx_point;
-		tmp_x.Append(*Event);
 	    }
 	}
     }

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[brlcad-commits] SF.net SVN: brlcad:[55870] brlcad/trunk/src/libbrep/intersect.cpp

[<pho...@us...>]

Revision: 55870
          http://sourceforge.net/p/brlcad/code/55870
Author:   phoenixyjll
Date:     2013-06-27 09:39:22 +0000 (Thu, 27 Jun 2013)
Log Message:
-----------
remove trailing ws.

Modified Paths:
--------------
    brlcad/trunk/src/libbrep/intersect.cpp

Modified: brlcad/trunk/src/libbrep/intersect.cpp
===================================================================
--- brlcad/trunk/src/libbrep/intersect.cpp	2013-06-27 09:30:55 UTC (rev 55869)
+++ brlcad/trunk/src/libbrep/intersect.cpp	2013-06-27 09:39:22 UTC (rev 55870)
@@ -800,7 +800,7 @@
 	}
 	if (j == points.Count()) {
 	    points.Append(tmp_x[i]);
-	}	
+	}
     }
 
     // Merge the overlap events that are continuous

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.