Revision: 73092
http://sourceforge.net/p/brlcad/code/73092
Author: brlcad
Date: 2019-05-20 06:25:40 +0000 (Mon, 20 May 2019)
Log Message:
-----------
ws indent cleanup
Modified Paths:
--------------
brlcad/trunk/src/libbg/QuickHull.cxx
brlcad/trunk/src/libbg/bg_private.h
brlcad/trunk/src/libbg/chull.c
brlcad/trunk/src/libbg/chull3d.cxx
brlcad/trunk/src/libbg/clipper.cxx
brlcad/trunk/src/libbg/polygon.c
brlcad/trunk/src/libbg/polygon_point_in.c
brlcad/trunk/src/libbg/polygon_triangulate.cxx
brlcad/trunk/src/libbg/spsr.c
brlcad/trunk/src/libbg/tests/chull.c
brlcad/trunk/src/libbg/tests/obr.c
brlcad/trunk/src/libbg/tests/polygon_triangulate.c
brlcad/trunk/src/libbg/tri_ray.c
brlcad/trunk/src/libbg/trimesh.c
brlcad/trunk/src/libbg/trimesh_isect.cxx
brlcad/trunk/src/libbg/util.c
Modified: brlcad/trunk/src/libbg/QuickHull.cxx
===================================================================
--- brlcad/trunk/src/libbg/QuickHull.cxx 2019-05-20 06:20:41 UTC (rev
73091)
+++ brlcad/trunk/src/libbg/QuickHull.cxx 2019-05-20 06:25:40 UTC (rev
73092)
@@ -24,494 +24,494 @@
namespace quickhull {
- template<>
- const float QuickHull<float>::Epsilon = 0.0001f;
+template<>
+const float QuickHull<float>::Epsilon = 0.0001f;
- template<>
- const double QuickHull<double>::Epsilon = 0.0000001;
+template<>
+const double QuickHull<double>::Epsilon = 0.0000001;
- /*
- * Implementation of the algorithm
- */
+/*
+ * Implementation of the algorithm
+ */
- template<typename T>
- ConvexHull<T> QuickHull<T>::getConvexHull(const
std::vector<Vector3<T>>& pointCloud, bool CCW, bool useOriginalIndices, T
epsilon) {
- VertexDataSource<T> vertexDataSource(pointCloud);
- return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
- }
+template<typename T>
+ConvexHull<T> QuickHull<T>::getConvexHull(const std::vector<Vector3<T>>&
pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
+ VertexDataSource<T> vertexDataSource(pointCloud);
+ return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+}
- template<typename T>
- ConvexHull<T> QuickHull<T>::getConvexHull(const Vector3<T>* vertexData,
size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
- VertexDataSource<T> vertexDataSource(vertexData,vertexCount);
- return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
- }
+template<typename T>
+ConvexHull<T> QuickHull<T>::getConvexHull(const Vector3<T>* vertexData, size_t
vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
+ VertexDataSource<T> vertexDataSource(vertexData,vertexCount);
+ return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+}
- template<typename T>
- ConvexHull<T> QuickHull<T>::getConvexHull(const T* vertexData, size_t
vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
- VertexDataSource<T> vertexDataSource((const
vec3*)vertexData,vertexCount);
- return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
- }
+template<typename T>
+ConvexHull<T> QuickHull<T>::getConvexHull(const T* vertexData, size_t
vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
+ VertexDataSource<T> vertexDataSource((const vec3*)vertexData,vertexCount);
+ return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+}
- template<typename FloatType>
- HalfEdgeMesh<FloatType, IndexType>
QuickHull<FloatType>::getConvexHullAsMesh(const FloatType* vertexData, size_t
vertexCount, bool CCW, FloatType epsilon) {
- VertexDataSource<FloatType> vertexDataSource((const
vec3*)vertexData,vertexCount);
- buildMesh(vertexDataSource, CCW, false, epsilon);
- return HalfEdgeMesh<FloatType, IndexType>(m_mesh, m_vertexData);
- }
+template<typename FloatType>
+HalfEdgeMesh<FloatType, IndexType>
QuickHull<FloatType>::getConvexHullAsMesh(const FloatType* vertexData, size_t
vertexCount, bool CCW, FloatType epsilon) {
+ VertexDataSource<FloatType> vertexDataSource((const
vec3*)vertexData,vertexCount);
+ buildMesh(vertexDataSource, CCW, false, epsilon);
+ return HalfEdgeMesh<FloatType, IndexType>(m_mesh, m_vertexData);
+}
- template<typename T>
- void QuickHull<T>::buildMesh(const VertexDataSource<T>& pointCloud,
bool UNUSED(CCW), bool UNUSED(useOriginalIndices), T epsilon) {
- if (pointCloud.size()==0) {
- m_mesh = MeshBuilder<T>();
- return;
- }
- m_vertexData = pointCloud;
+template<typename T>
+void QuickHull<T>::buildMesh(const VertexDataSource<T>& pointCloud, bool
UNUSED(CCW), bool UNUSED(useOriginalIndices), T epsilon) {
+ if (pointCloud.size()==0) {
+ m_mesh = MeshBuilder<T>();
+ return;
+ }
+ m_vertexData = pointCloud;
- // Very first: find extreme values and use them to compute the
scale of the point cloud.
- m_extremeValues = getExtremeValues();
- m_scale = getScale(m_extremeValues);
+ // Very first: find extreme values and use them to compute the scale of
the point cloud.
+ m_extremeValues = getExtremeValues();
+ m_scale = getScale(m_extremeValues);
- // Epsilon we use depends on the scale
- m_epsilon = epsilon*m_scale;
- m_epsilonSquared = m_epsilon*m_epsilon;
+ // Epsilon we use depends on the scale
+ m_epsilon = epsilon*m_scale;
+ m_epsilonSquared = m_epsilon*m_epsilon;
- m_planar = false; // The planar case happens when all the points
appear to lie on a two dimensional subspace of R^3.
- createConvexHalfEdgeMesh();
- if (m_planar) {
- const size_t extraPointIndex = m_planarPointCloudTemp.size()-1;
- for (auto& he : m_mesh.m_halfEdges) {
- if (he.m_endVertex == extraPointIndex) {
- he.m_endVertex = 0;
- }
- }
- m_vertexData = pointCloud;
- m_planarPointCloudTemp.clear();
+ m_planar = false; // The planar case happens when all the points appear to
lie on a two dimensional subspace of R^3.
+ createConvexHalfEdgeMesh();
+ if (m_planar) {
+ const size_t extraPointIndex = m_planarPointCloudTemp.size()-1;
+ for (auto& he : m_mesh.m_halfEdges) {
+ if (he.m_endVertex == extraPointIndex) {
+ he.m_endVertex = 0;
}
}
+ m_vertexData = pointCloud;
+ m_planarPointCloudTemp.clear();
+ }
+}
- template<typename T>
- ConvexHull<T> QuickHull<T>::getConvexHull(const VertexDataSource<T>&
pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
- buildMesh(pointCloud,CCW,useOriginalIndices,epsilon);
- return ConvexHull<T>(m_mesh,m_vertexData, CCW, useOriginalIndices);
+template<typename T>
+ConvexHull<T> QuickHull<T>::getConvexHull(const VertexDataSource<T>&
pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
+ buildMesh(pointCloud,CCW,useOriginalIndices,epsilon);
+ return ConvexHull<T>(m_mesh,m_vertexData, CCW, useOriginalIndices);
+}
+
+template<typename T>
+void QuickHull<T>::createConvexHalfEdgeMesh() {
+ // Temporary variables used during iteration
+ std::vector<IndexType> visibleFaces;
+ std::vector<IndexType> horizonEdges;
+ struct FaceData {
+ IndexType m_faceIndex;
+ IndexType m_enteredFromHalfEdge; // If the face turns out not to be
visible, this half edge will be marked as horizon edge
+ FaceData(IndexType fi, IndexType he) :
m_faceIndex(fi),m_enteredFromHalfEdge(he) {}
+ };
+ std::vector<FaceData> possiblyVisibleFaces;
+
+ // Compute base tetrahedron
+ m_mesh = getInitialTetrahedron();
+ assert(m_mesh.m_faces.size()==4);
+
+ // Init face stack with those faces that have points assigned to them
+ std::deque<IndexType> faceList;
+ for (size_t i=0;i < 4;i++) {
+ auto& f = m_mesh.m_faces[i];
+ if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size()>0) {
+ faceList.push_back(i);
+ f.m_inFaceStack = 1;
}
+ }
- template<typename T>
- void QuickHull<T>::createConvexHalfEdgeMesh() {
- // Temporary variables used during iteration
- std::vector<IndexType> visibleFaces;
- std::vector<IndexType> horizonEdges;
- struct FaceData {
- IndexType m_faceIndex;
- IndexType m_enteredFromHalfEdge; // If the face turns out not
to be visible, this half edge will be marked as horizon edge
- FaceData(IndexType fi, IndexType he) :
m_faceIndex(fi),m_enteredFromHalfEdge(he) {}
- };
- std::vector<FaceData> possiblyVisibleFaces;
+ // Process faces until the face list is empty.
+ size_t iter = 0;
+ while (!faceList.empty()) {
+ iter++;
+ if (iter == std::numeric_limits<size_t>::max()) {
+ // Visible face traversal marks visited faces with iteration
counter (to mark that the face has been visited on this iteration) and the max
value represents unvisited faces. At this point we have to reset iteration
counter. This shouldn't be an
+ // issue on 64 bit machines.
+ iter = 0;
+ }
- // Compute base tetrahedron
- m_mesh = getInitialTetrahedron();
- assert(m_mesh.m_faces.size()==4);
+ const IndexType topFaceIndex = faceList.front();
+ faceList.pop_front();
- // Init face stack with those faces that have points assigned to
them
- std::deque<IndexType> faceList;
- for (size_t i=0;i < 4;i++) {
- auto& f = m_mesh.m_faces[i];
- if (f.m_pointsOnPositiveSide &&
f.m_pointsOnPositiveSide->size()>0) {
- faceList.push_back(i);
- f.m_inFaceStack = 1;
- }
- }
+ auto& tf = m_mesh.m_faces[topFaceIndex];
+ tf.m_inFaceStack = 0;
- // Process faces until the face list is empty.
- size_t iter = 0;
- while (!faceList.empty()) {
- iter++;
- if (iter == std::numeric_limits<size_t>::max()) {
- // Visible face traversal marks visited faces with
iteration counter (to mark that the face has been visited on this iteration)
and the max value represents unvisited faces. At this point we have to reset
iteration counter. This shouldn't be an
- // issue on 64 bit machines.
- iter = 0;
- }
+ assert(!tf.m_pointsOnPositiveSide ||
tf.m_pointsOnPositiveSide->size()>0);
+ if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) {
+ continue;
+ }
- const IndexType topFaceIndex = faceList.front();
- faceList.pop_front();
+ // Pick the most distant point to this triangle plane as the point to
which we extrude
+ const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint];
+ const size_t activePointIndex = tf.m_mostDistantPoint;
- auto& tf = m_mesh.m_faces[topFaceIndex];
- tf.m_inFaceStack = 0;
+ // Find out the faces that have our active point on their positive side
(these are the "visible faces"). The face on top of the stack of course is one
of them. At the same time, we create a list of horizon edges.
+ horizonEdges.clear();
+ possiblyVisibleFaces.clear();
+ visibleFaces.clear();
+
possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits<size_t>::max());
+ while (possiblyVisibleFaces.size()) {
+ const auto faceData = possiblyVisibleFaces.back();
+ possiblyVisibleFaces.pop_back();
+ auto& pvf = m_mesh.m_faces[faceData.m_faceIndex];
+ assert(!pvf.isDisabled());
- assert(!tf.m_pointsOnPositiveSide ||
tf.m_pointsOnPositiveSide->size()>0);
- if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) {
+ if (pvf.m_visibilityCheckedOnIteration == iter) {
+ if (pvf.m_isVisibleFaceOnCurrentIteration) {
continue;
}
-
- // Pick the most distant point to this triangle plane as the
point to which we extrude
- const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint];
- const size_t activePointIndex = tf.m_mostDistantPoint;
-
- // Find out the faces that have our active point on their
positive side (these are the "visible faces"). The face on top of the stack of
course is one of them. At the same time, we create a list of horizon edges.
- horizonEdges.clear();
- possiblyVisibleFaces.clear();
- visibleFaces.clear();
-
possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits<size_t>::max());
- while (possiblyVisibleFaces.size()) {
- const auto faceData = possiblyVisibleFaces.back();
- possiblyVisibleFaces.pop_back();
- auto& pvf = m_mesh.m_faces[faceData.m_faceIndex];
- assert(!pvf.isDisabled());
-
- if (pvf.m_visibilityCheckedOnIteration == iter) {
- if (pvf.m_isVisibleFaceOnCurrentIteration) {
- continue;
+ }
+ else {
+ const Plane<T>& P = pvf.m_P;
+ pvf.m_visibilityCheckedOnIteration = iter;
+ const T d = P.m_N.dotProduct(activePoint)+P.m_D;
+ if (d>0) {
+ pvf.m_isVisibleFaceOnCurrentIteration = 1;
+ pvf.m_horizonEdgesOnCurrentIteration = 0;
+ visibleFaces.push_back(faceData.m_faceIndex);
+ for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) {
+ if (m_mesh.m_halfEdges[heIndex].m_opp !=
faceData.m_enteredFromHalfEdge) {
+ possiblyVisibleFaces.emplace_back(
m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex );
}
}
- else {
- const Plane<T>& P = pvf.m_P;
- pvf.m_visibilityCheckedOnIteration = iter;
- const T d = P.m_N.dotProduct(activePoint)+P.m_D;
- if (d>0) {
- pvf.m_isVisibleFaceOnCurrentIteration = 1;
- pvf.m_horizonEdgesOnCurrentIteration = 0;
- visibleFaces.push_back(faceData.m_faceIndex);
- for (auto heIndex :
m_mesh.getHalfEdgeIndicesOfFace(pvf)) {
- if (m_mesh.m_halfEdges[heIndex].m_opp !=
faceData.m_enteredFromHalfEdge) {
- possiblyVisibleFaces.emplace_back(
m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex );
- }
- }
- continue;
- }
- assert(faceData.m_faceIndex != topFaceIndex);
- }
-
- // The face is not visible. Therefore, the halfedge we came
from is part of the horizon edge.
- pvf.m_isVisibleFaceOnCurrentIteration = 0;
- horizonEdges.push_back(faceData.m_enteredFromHalfEdge);
- // Store which half edge is the horizon edge. The other
half edges of the face will not be part of the final mesh so their data slots
can by recycled.
- const auto halfEdges =
m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]);
- const std::int8_t ind =
(halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 :
(halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2);
-
m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration
|= (1<<ind);
- }
- const size_t horizonEdgeCount = horizonEdges.size();
-
- // Order horizon edges so that they form a loop. This may fail
due to numerical instability in which case we give up trying to solve horizon
edge for this point and accept a minor degeneration in the convex hull.
- if (!reorderHorizonEdges(horizonEdges)) {
- std::cerr << "Failed to solve horizon edge." << std::endl;
- auto it =
std::find(tf.m_pointsOnPositiveSide->begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex);
- tf.m_pointsOnPositiveSide->erase(it);
- if (tf.m_pointsOnPositiveSide->size()==0) {
- reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide);
- }
continue;
}
+ assert(faceData.m_faceIndex != topFaceIndex);
+ }
- // Except for the horizon edges, all half edges of the visible
faces can be marked as disabled. Their data slots will be reused.
- // The faces will be disabled as well, but we need to remember
the points that were on the positive side of them - therefore
- // we save pointers to them.
- m_newFaceIndices.clear();
- m_newHalfEdgeIndices.clear();
- m_disabledFacePointVectors.clear();
- size_t disableCounter = 0;
- for (auto faceIndex : visibleFaces) {
- auto& disabledFace = m_mesh.m_faces[faceIndex];
- auto halfEdges =
m_mesh.getHalfEdgeIndicesOfFace(disabledFace);
- for (size_t j=0;j<3;j++) {
- if ((disabledFace.m_horizonEdgesOnCurrentIteration &
(1<<j)) == 0) {
- if (disableCounter < horizonEdgeCount*2) {
- // Use on this iteration
- m_newHalfEdgeIndices.push_back(halfEdges[j]);
- disableCounter++;
- }
- else {
- // Mark for reusal on later iteration step
- m_mesh.disableHalfEdge(halfEdges[j]);
- }
- }
+ // The face is not visible. Therefore, the halfedge we came from is
part of the horizon edge.
+ pvf.m_isVisibleFaceOnCurrentIteration = 0;
+ horizonEdges.push_back(faceData.m_enteredFromHalfEdge);
+ // Store which half edge is the horizon edge. The other half edges
of the face will not be part of the final mesh so their data slots can by
recycled.
+ const auto halfEdges =
m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]);
+ const std::int8_t ind =
(halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 :
(halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2);
+
m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration
|= (1<<ind);
+ }
+ const size_t horizonEdgeCount = horizonEdges.size();
+
+ // Order horizon edges so that they form a loop. This may fail due to
numerical instability in which case we give up trying to solve horizon edge for
this point and accept a minor degeneration in the convex hull.
+ if (!reorderHorizonEdges(horizonEdges)) {
+ std::cerr << "Failed to solve horizon edge." << std::endl;
+ auto it =
std::find(tf.m_pointsOnPositiveSide->begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex);
+ tf.m_pointsOnPositiveSide->erase(it);
+ if (tf.m_pointsOnPositiveSide->size()==0) {
+ reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide);
+ }
+ continue;
+ }
+
+ // Except for the horizon edges, all half edges of the visible faces
can be marked as disabled. Their data slots will be reused.
+ // The faces will be disabled as well, but we need to remember the
points that were on the positive side of them - therefore
+ // we save pointers to them.
+ m_newFaceIndices.clear();
+ m_newHalfEdgeIndices.clear();
+ m_disabledFacePointVectors.clear();
+ size_t disableCounter = 0;
+ for (auto faceIndex : visibleFaces) {
+ auto& disabledFace = m_mesh.m_faces[faceIndex];
+ auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace);
+ for (size_t j=0;j<3;j++) {
+ if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1<<j)) ==
0) {
+ if (disableCounter < horizonEdgeCount*2) {
+ // Use on this iteration
+ m_newHalfEdgeIndices.push_back(halfEdges[j]);
+ disableCounter++;
}
- // Disable the face, but retain pointer to the points that
were on the positive side of it. We need to assign those points
- // to the new faces we create shortly.
- auto t = std::move(m_mesh.disableFace(faceIndex));
- if (t) {
- assert(t->size()); // Because we should not assign
point vectors to faces unless needed...
- m_disabledFacePointVectors.push_back(std::move(t));
+ else {
+ // Mark for reusal on later iteration step
+ m_mesh.disableHalfEdge(halfEdges[j]);
}
}
- if (disableCounter < horizonEdgeCount*2) {
- const size_t newHalfEdgesNeeded =
horizonEdgeCount*2-disableCounter;
- for (size_t i=0;i<newHalfEdgesNeeded;i++) {
- m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge());
- }
- }
+ }
+ // Disable the face, but retain pointer to the points that were on
the positive side of it. We need to assign those points
+ // to the new faces we create shortly.
+ auto t = std::move(m_mesh.disableFace(faceIndex));
+ if (t) {
+ assert(t->size()); // Because we should not assign point
vectors to faces unless needed...
+ m_disabledFacePointVectors.push_back(std::move(t));
+ }
+ }
+ if (disableCounter < horizonEdgeCount*2) {
+ const size_t newHalfEdgesNeeded = horizonEdgeCount*2-disableCounter;
+ for (size_t i=0;i<newHalfEdgesNeeded;i++) {
+ m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge());
+ }
+ }
- // Create new faces using the edgeloop
- for (size_t i = 0; i < horizonEdgeCount; i++) {
- const IndexType AB = horizonEdges[i];
+ // Create new faces using the edgeloop
+ for (size_t i = 0; i < horizonEdgeCount; i++) {
+ const IndexType AB = horizonEdges[i];
- auto horizonEdgeVertexIndices =
m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]);
- IndexType A,B,C;
- A = horizonEdgeVertexIndices[0];
- B = horizonEdgeVertexIndices[1];
- C = activePointIndex;
+ auto horizonEdgeVertexIndices =
m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]);
+ IndexType A,B,C;
+ A = horizonEdgeVertexIndices[0];
+ B = horizonEdgeVertexIndices[1];
+ C = activePointIndex;
- const IndexType newFaceIndex = m_mesh.addFace();
- m_newFaceIndices.push_back(newFaceIndex);
+ const IndexType newFaceIndex = m_mesh.addFace();
+ m_newFaceIndices.push_back(newFaceIndex);
- const IndexType CA = m_newHalfEdgeIndices[2*i+0];
- const IndexType BC = m_newHalfEdgeIndices[2*i+1];
+ const IndexType CA = m_newHalfEdgeIndices[2*i+0];
+ const IndexType BC = m_newHalfEdgeIndices[2*i+1];
- m_mesh.m_halfEdges[AB].m_next = BC;
- m_mesh.m_halfEdges[BC].m_next = CA;
- m_mesh.m_halfEdges[CA].m_next = AB;
+ m_mesh.m_halfEdges[AB].m_next = BC;
+ m_mesh.m_halfEdges[BC].m_next = CA;
+ m_mesh.m_halfEdges[CA].m_next = AB;
- m_mesh.m_halfEdges[BC].m_face = newFaceIndex;
- m_mesh.m_halfEdges[CA].m_face = newFaceIndex;
- m_mesh.m_halfEdges[AB].m_face = newFaceIndex;
+ m_mesh.m_halfEdges[BC].m_face = newFaceIndex;
+ m_mesh.m_halfEdges[CA].m_face = newFaceIndex;
+ m_mesh.m_halfEdges[AB].m_face = newFaceIndex;
- m_mesh.m_halfEdges[CA].m_endVertex = A;
- m_mesh.m_halfEdges[BC].m_endVertex = C;
+ m_mesh.m_halfEdges[CA].m_endVertex = A;
+ m_mesh.m_halfEdges[BC].m_endVertex = C;
- auto& newFace = m_mesh.m_faces[newFaceIndex];
+ auto& newFace = m_mesh.m_faces[newFaceIndex];
- const Vector3<T> planeNormal =
mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint);
- newFace.m_P = Plane<T>(planeNormal,activePoint);
- newFace.m_he = AB;
+ const Vector3<T> planeNormal =
mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint);
+ newFace.m_P = Plane<T>(planeNormal,activePoint);
+ newFace.m_he = AB;
- m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ?
i*2-1 : 2*horizonEdgeCount-1];
- m_mesh.m_halfEdges[BC].m_opp =
m_newHalfEdgeIndices[((i+1)*2) % (horizonEdgeCount*2)];
- }
+ m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ? i*2-1 :
2*horizonEdgeCount-1];
+ m_mesh.m_halfEdges[BC].m_opp = m_newHalfEdgeIndices[((i+1)*2) %
(horizonEdgeCount*2)];
+ }
- // Assign points that were on the positive side of the disabled
faces to the new faces.
- for (auto& disabledPoints : m_disabledFacePointVectors) {
- assert(disabledPoints);
- for (const auto& point : *(disabledPoints)) {
- if (point == activePointIndex) {
- continue;
- }
- for (size_t j=0;j<horizonEdgeCount;j++) {
- if
(addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]], point)) {
- break;
- }
- }
- }
- // The points are no longer needed: we can move them to the
vector pool for reuse.
- reclaimToIndexVectorPool(disabledPoints);
+ // Assign points that were on the positive side of the disabled faces
to the new faces.
+ for (auto& disabledPoints : m_disabledFacePointVectors) {
+ assert(disabledPoints);
+ for (const auto& point : *(disabledPoints)) {
+ if (point == activePointIndex) {
+ continue;
}
-
- // Increase face stack size if needed
- for (const auto newFaceIndex : m_newFaceIndices) {
- auto& newFace = m_mesh.m_faces[newFaceIndex];
- if (newFace.m_pointsOnPositiveSide) {
- assert(newFace.m_pointsOnPositiveSide->size()>0);
- if (!newFace.m_inFaceStack) {
- faceList.push_back(newFaceIndex);
- newFace.m_inFaceStack = 1;
- }
+ for (size_t j=0;j<horizonEdgeCount;j++) {
+ if (addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]],
point)) {
+ break;
}
}
}
-
- // Cleanup
- m_indexVectorPool.clear();
+ // The points are no longer needed: we can move them to the vector
pool for reuse.
+ reclaimToIndexVectorPool(disabledPoints);
}
- /*
- * Private helper functions
- */
-
- template <typename T>
- std::array<IndexType,6> QuickHull<T>::getExtremeValues() {
- std::array<IndexType,6> outIndices{{0,0,0,0,0,0}};
- T extremeVals[6] =
{m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z};
- const size_t vCount = m_vertexData.size();
- for (size_t i=1;i<vCount;i++) {
- const Vector3<T>& pos = m_vertexData[i];
- if (pos.x>extremeVals[0]) {
- extremeVals[0]=pos.x;
- outIndices[0]=(IndexType)i;
+ // Increase face stack size if needed
+ for (const auto newFaceIndex : m_newFaceIndices) {
+ auto& newFace = m_mesh.m_faces[newFaceIndex];
+ if (newFace.m_pointsOnPositiveSide) {
+ assert(newFace.m_pointsOnPositiveSide->size()>0);
+ if (!newFace.m_inFaceStack) {
+ faceList.push_back(newFaceIndex);
+ newFace.m_inFaceStack = 1;
}
- else if (pos.x<extremeVals[1]) {
- extremeVals[1]=pos.x;
- outIndices[1]=(IndexType)i;
- }
- if (pos.y>extremeVals[2]) {
- extremeVals[2]=pos.y;
- outIndices[2]=(IndexType)i;
- }
- else if (pos.y<extremeVals[3]) {
- extremeVals[3]=pos.y;
- outIndices[3]=(IndexType)i;
- }
- if (pos.z>extremeVals[4]) {
- extremeVals[4]=pos.z;
- outIndices[4]=(IndexType)i;
- }
- else if (pos.z<extremeVals[5]) {
- extremeVals[5]=pos.z;
- outIndices[5]=(IndexType)i;
- }
}
- return outIndices;
}
+ }
- template<typename T>
- bool QuickHull<T>::reorderHorizonEdges(std::vector<IndexType>&
horizonEdges) {
- const size_t horizonEdgeCount = horizonEdges.size();
- for (size_t i=0;i<horizonEdgeCount-1;i++) {
- const IndexType endVertex = m_mesh.m_halfEdges[ horizonEdges[i]
].m_endVertex;
- bool foundNext = false;
- for (size_t j=i+1;j<horizonEdgeCount;j++) {
- const IndexType beginVertex = m_mesh.m_halfEdges[
m_mesh.m_halfEdges[horizonEdges[j]].m_opp ].m_endVertex;
- if (beginVertex == endVertex) {
- std::swap(horizonEdges[i+1],horizonEdges[j]);
- foundNext = true;
- break;
- }
- }
- if (!foundNext) {
- return false;
- }
- }
- assert(m_mesh.m_halfEdges[ horizonEdges[horizonEdges.size()-1]
].m_endVertex == m_mesh.m_halfEdges[ m_mesh.m_halfEdges[horizonEdges[0]].m_opp
].m_endVertex);
- return true;
+ // Cleanup
+ m_indexVectorPool.clear();
+}
+
+/*
+ * Private helper functions
+ */
+
+template <typename T>
+std::array<IndexType,6> QuickHull<T>::getExtremeValues() {
+ std::array<IndexType,6> outIndices{{0,0,0,0,0,0}};
+ T extremeVals[6] =
{m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z};
+ const size_t vCount = m_vertexData.size();
+ for (size_t i=1;i<vCount;i++) {
+ const Vector3<T>& pos = m_vertexData[i];
+ if (pos.x>extremeVals[0]) {
+ extremeVals[0]=pos.x;
+ outIndices[0]=(IndexType)i;
}
+ else if (pos.x<extremeVals[1]) {
+ extremeVals[1]=pos.x;
+ outIndices[1]=(IndexType)i;
+ }
+ if (pos.y>extremeVals[2]) {
+ extremeVals[2]=pos.y;
+ outIndices[2]=(IndexType)i;
+ }
+ else if (pos.y<extremeVals[3]) {
+ extremeVals[3]=pos.y;
+ outIndices[3]=(IndexType)i;
+ }
+ if (pos.z>extremeVals[4]) {
+ extremeVals[4]=pos.z;
+ outIndices[4]=(IndexType)i;
+ }
+ else if (pos.z<extremeVals[5]) {
+ extremeVals[5]=pos.z;
+ outIndices[5]=(IndexType)i;
+ }
+ }
+ return outIndices;
+}
- template <typename T>
- T QuickHull<T>::getScale(const std::array<IndexType,6>& extremeValues) {
- T s = 0;
- for (size_t i=0;i<6;i++) {
- const T* v = (const T*)(&m_vertexData[extremeValues[i]]);
- v += i/2;
- auto a = std::abs(*v);
- if (a>s) {
- s = a;
- }
+template<typename T>
+bool QuickHull<T>::reorderHorizonEdges(std::vector<IndexType>& horizonEdges) {
+ const size_t horizonEdgeCount = horizonEdges.size();
+ for (size_t i=0;i<horizonEdgeCount-1;i++) {
+ const IndexType endVertex = m_mesh.m_halfEdges[ horizonEdges[i]
].m_endVertex;
+ bool foundNext = false;
+ for (size_t j=i+1;j<horizonEdgeCount;j++) {
+ const IndexType beginVertex = m_mesh.m_halfEdges[
m_mesh.m_halfEdges[horizonEdges[j]].m_opp ].m_endVertex;
+ if (beginVertex == endVertex) {
+ std::swap(horizonEdges[i+1],horizonEdges[j]);
+ foundNext = true;
+ break;
}
- return s;
}
+ if (!foundNext) {
+ return false;
+ }
+ }
+ assert(m_mesh.m_halfEdges[ horizonEdges[horizonEdges.size()-1]
].m_endVertex == m_mesh.m_halfEdges[ m_mesh.m_halfEdges[horizonEdges[0]].m_opp
].m_endVertex);
+ return true;
+}
- template <typename T>
- MeshBuilder<T> QuickHull<T>::getInitialTetrahedron() {
- const size_t vertexCount = m_vertexData.size();
+template <typename T>
+T QuickHull<T>::getScale(const std::array<IndexType,6>& extremeValues) {
+ T s = 0;
+ for (size_t i=0;i<6;i++) {
+ const T* v = (const T*)(&m_vertexData[extremeValues[i]]);
+ v += i/2;
+ auto a = std::abs(*v);
+ if (a>s) {
+ s = a;
+ }
+ }
+ return s;
+}
- // If we have at most 4 points, just return a degenerate
tetrahedron:
- if (vertexCount <= 4) {
- IndexType v[4] =
{0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)};
- const Vector3<T> N =
mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]);
- const Plane<T> trianglePlane(N,m_vertexData[v[0]]);
- if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) {
- std::swap(v[0],v[1]);
- }
- return MeshBuilder<T>(v[0],v[1],v[2],v[3]);
- }
+template <typename T>
+MeshBuilder<T> QuickHull<T>::getInitialTetrahedron() {
+ const size_t vertexCount = m_vertexData.size();
- // Find two most distant extreme points.
- T maxD = m_epsilonSquared;
- std::pair<IndexType,IndexType> selectedPoints;
- for (size_t i=0;i<6;i++) {
- for (size_t j=i+1;j<6;j++) {
- const T d = m_vertexData[ m_extremeValues[i]
].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] );
- if (d > maxD) {
- maxD=d;
- selectedPoints={m_extremeValues[i],m_extremeValues[j]};
- }
- }
- }
- if (EQUAL(maxD, m_epsilonSquared)) {
- // A degenerate case: the point cloud seems to consists of a
single point
- return
MeshBuilder<T>(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1));
- }
- assert(selectedPoints.first != selectedPoints.second);
+ // If we have at most 4 points, just return a degenerate tetrahedron:
+ if (vertexCount <= 4) {
+ IndexType v[4] =
{0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)};
+ const Vector3<T> N =
mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]);
+ const Plane<T> trianglePlane(N,m_vertexData[v[0]]);
+ if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) {
+ std::swap(v[0],v[1]);
+ }
+ return MeshBuilder<T>(v[0],v[1],v[2],v[3]);
+ }
- // Find the most distant point to the line between the two chosen
extreme points.
- const Ray<T> r(m_vertexData[selectedPoints.first],
(m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first]));
- maxD = m_epsilonSquared;
- size_t maxI=std::numeric_limits<size_t>::max();
- const size_t vCount = m_vertexData.size();
- for (size_t i=0;i<vCount;i++) {
- const T distToRay =
mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i],r);
- if (distToRay > maxD) {
- maxD=distToRay;
- maxI=i;
- }
+ // Find two most distant extreme points.
+ T maxD = m_epsilonSquared;
+ std::pair<IndexType,IndexType> selectedPoints;
+ for (size_t i=0;i<6;i++) {
+ for (size_t j=i+1;j<6;j++) {
+ const T d = m_vertexData[ m_extremeValues[i]
].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] );
+ if (d > maxD) {
+ maxD=d;
+ selectedPoints={m_extremeValues[i],m_extremeValues[j]};
}
- if (EQUAL(maxD, m_epsilonSquared)) {
- // It appears that the point cloud belongs to a 1 dimensional
- // subspace of R^3: convex hull has no volume => return a thin
- // triangle Pick any point other than selectedPoints.first and
- // selectedPoints.second as the third point of the triangle
- auto it =
std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
- return ve != m_vertexData[selectedPoints.first] && ve
!= m_vertexData[selectedPoints.second];
- });
- const IndexType thirdPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
- it =
std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
- return ve != m_vertexData[selectedPoints.first] && ve
!= m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint];
- });
- const IndexType fourthPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
- return
MeshBuilder<T>(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint);
- }
+ }
+ }
+ if (EQUAL(maxD, m_epsilonSquared)) {
+ // A degenerate case: the point cloud seems to consists of a single
point
+ return
MeshBuilder<T>(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1));
+ }
+ assert(selectedPoints.first != selectedPoints.second);
- // These three points form the base triangle for our tetrahedron.
- assert(selectedPoints.first != maxI && selectedPoints.second !=
maxI);
- std::array<size_t,3> baseTriangle{{selectedPoints.first,
selectedPoints.second, maxI}};
- const Vector3<T> baseTriangleVertices[]={
m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]],
m_vertexData[baseTriangle[2]] };
+ // Find the most distant point to the line between the two chosen extreme
points.
+ const Ray<T> r(m_vertexData[selectedPoints.first],
(m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first]));
+ maxD = m_epsilonSquared;
+ size_t maxI=std::numeric_limits<size_t>::max();
+ const size_t vCount = m_vertexData.size();
+ for (size_t i=0;i<vCount;i++) {
+ const T distToRay =
mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i],r);
+ if (distToRay > maxD) {
+ maxD=distToRay;
+ maxI=i;
+ }
+ }
+ if (EQUAL(maxD, m_epsilonSquared)) {
+ // It appears that the point cloud belongs to a 1 dimensional
+ // subspace of R^3: convex hull has no volume => return a thin
+ // triangle Pick any point other than selectedPoints.first and
+ // selectedPoints.second as the third point of the triangle
+ auto it =
std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
+
return ve != m_vertexData[selectedPoints.first] && ve !=
m_vertexData[selectedPoints.second];
+ });
+ const IndexType thirdPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
+ it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const
vec3& ve) {
+ return ve
!= m_vertexData[selectedPoints.first] && ve !=
m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint];
+ });
+ const IndexType fourthPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
+ return
MeshBuilder<T>(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint);
+ }
- // Next step is to find the 4th vertex of the tetrahedron. We
- // naturally choose the point farthest away from the triangle
- // plane.
- maxD=m_epsilon;
- maxI=0;
- const Vector3<T> N =
mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]);
- Plane<T> trianglePlane(N,baseTriangleVertices[0]);
- for (size_t i=0;i<vCount;i++) {
- const T d =
std::abs(mathutils::getSignedDistanceToPlane(m_vertexData[i],trianglePlane));
- if (d>maxD) {
- maxD=d;
- maxI=i;
- }
- }
- if (EQUAL(maxD, m_epsilon)) {
- // All the points seem to lie on a 2D subspace of R^3. How to
handle this? Well, let's add one extra point to the point cloud so that the
convex hull will have volume.
- m_planar = true;
- const vec3 N1 =
mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]);
- m_planarPointCloudTemp.clear();
-
m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end());
- const vec3 extraPoint = N1 + m_vertexData[0];
- m_planarPointCloudTemp.push_back(extraPoint);
- maxI = m_planarPointCloudTemp.size()-1;
- m_vertexData = VertexDataSource<T>(m_planarPointCloudTemp);
- }
+ // These three points form the base triangle for our tetrahedron.
+ assert(selectedPoints.first != maxI && selectedPoints.second != maxI);
+ std::array<size_t,3> baseTriangle{{selectedPoints.first,
selectedPoints.second, maxI}};
+ const Vector3<T> baseTriangleVertices[]={ m_vertexData[baseTriangle[0]],
m_vertexData[baseTriangle[1]], m_vertexData[baseTriangle[2]] };
- // Enforce CCW orientation (if user prefers clockwise orientation,
swap two vertices in each triangle when final mesh is created)
- const Plane<T> triPlane(N,baseTriangleVertices[0]);
- if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) {
- std::swap(baseTriangle[0],baseTriangle[1]);
- }
+ // Next step is to find the 4th vertex of the tetrahedron. We
+ // naturally choose the point farthest away from the triangle
+ // plane.
+ maxD=m_epsilon;
+ maxI=0;
+ const Vector3<T> N =
mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]);
+ Plane<T> trianglePlane(N,baseTriangleVertices[0]);
+ for (size_t i=0;i<vCount;i++) {
+ const T d =
std::abs(mathutils::getSignedDistanceToPlane(m_vertexData[i],trianglePlane));
+ if (d>maxD) {
+ maxD=d;
+ maxI=i;
+ }
+ }
+ if (EQUAL(maxD, m_epsilon)) {
+ // All the points seem to lie on a 2D subspace of R^3. How to handle
this? Well, let's add one extra point to the point cloud so that the convex
hull will have volume.
+ m_planar = true;
+ const vec3 N1 =
mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]);
+ m_planarPointCloudTemp.clear();
+
m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end());
+ const vec3 extraPoint = N1 + m_vertexData[0];
+ m_planarPointCloudTemp.push_back(extraPoint);
+ maxI = m_planarPointCloudTemp.size()-1;
+ m_vertexData = VertexDataSource<T>(m_planarPointCloudTemp);
+ }
- // Create a tetrahedron half edge mesh and compute planes defined
by each triangle
- MeshBuilder<T>
mesh(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
- for (auto& f : mesh.m_faces) {
- auto v = mesh.getVertexIndicesOfFace(f);
- const Vector3<T>& va = m_vertexData[v[0]];
- const Vector3<T>& vb = m_vertexData[v[1]];
- const Vector3<T>& vc = m_vertexData[v[2]];
- const Vector3<T> N1 = mathutils::getTriangleNormal(va, vb, vc);
- const Plane<T> trianglePlane1(N1,va);
- f.m_P = trianglePlane1;
- }
+ // Enforce CCW orientation (if user prefers clockwise orientation, swap
two vertices in each triangle when final mesh is created)
+ const Plane<T> triPlane(N,baseTriangleVertices[0]);
+ if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) {
+ std::swap(baseTriangle[0],baseTriangle[1]);
+ }
- // Finally we assign a face for each vertex outside the tetrahedron
(vertices inside the tetrahedron have no role anymore)
- for (size_t i=0;i<vCount;i++) {
- for (auto& face : mesh.m_faces) {
- if (addPointToFace(face, i)) {
- break;
- }
- }
+ // Create a tetrahedron half edge mesh and compute planes defined by each
triangle
+ MeshBuilder<T> mesh(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
+ for (auto& f : mesh.m_faces) {
+ auto v = mesh.getVertexIndicesOfFace(f);
+ const Vector3<T>& va = m_vertexData[v[0]];
+ const Vector3<T>& vb = m_vertexData[v[1]];
+ const Vector3<T>& vc = m_vertexData[v[2]];
+ const Vector3<T> N1 = mathutils::getTriangleNormal(va, vb, vc);
+ const Plane<T> trianglePlane1(N1,va);
+ f.m_P = trianglePlane1;
+ }
+
+ // Finally we assign a face for each vertex outside the tetrahedron
(vertices inside the tetrahedron have no role anymore)
+ for (size_t i=0;i<vCount;i++) {
+ for (auto& face : mesh.m_faces) {
+ if (addPointToFace(face, i)) {
+ break;
}
- return mesh;
}
+ }
+ return mesh;
+}
- /*
- * Explicit template specifications for float and double
- */
+/*
+ * Explicit template specifications for float and double
+ */
- template class QuickHull<float>;
- template class QuickHull<double>;
+template class QuickHull<float>;
+template class QuickHull<double>;
}
// Local Variables:
Modified: brlcad/trunk/src/libbg/bg_private.h
===================================================================
--- brlcad/trunk/src/libbg/bg_private.h 2019-05-20 06:20:41 UTC (rev 73091)
+++ brlcad/trunk/src/libbg/bg_private.h 2019-05-20 06:25:40 UTC (rev 73092)
@@ -66,8 +66,8 @@
*/
extern int
coplanar_3d_to_2d(point2d_t **points_2d, const point_t *origin_pnt,
- const vect_t *u_axis, const vect_t *v_axis,
- const point_t *points_3d, int n);
+ const vect_t *u_axis, const vect_t *v_axis,
+ const point_t *points_3d, int n);
/**
* @brief
@@ -83,8 +83,8 @@
*/
extern int
coplanar_2d_to_3d(point_t **points_3d, const point_t *origin_pnt,
- const vect_t *u_axis, const vect_t *v_axis,
- const point2d_t *points_2d, int n);
+ const vect_t *u_axis, const vect_t *v_axis,
+ const point2d_t *points_2d, int n);
#endif /* LIBGM_LIBGM_PRIVATE_H */
Modified: brlcad/trunk/src/libbg/chull.c
===================================================================
--- brlcad/trunk/src/libbg/chull.c 2019-05-20 06:20:41 UTC (rev 73091)
+++ brlcad/trunk/src/libbg/chull.c 2019-05-20 06:25:40 UTC (rev 73092)
@@ -74,7 +74,7 @@
/* test if next vertex is inside the deque hull */
if ((isLeft(D[bot], D[bot+1], polyline[i]) > 0) &&
(isLeft(D[top-1], D[top], polyline[i]) > 0) )
- continue; /* skip an interior vertex */
+ continue; /* skip an interior vertex */
/* incrementally add an exterior vertex to the deque hull
get the rightmost tangent at the deque bot */
Modified: brlcad/trunk/src/libbg/chull3d.cxx
===================================================================
--- brlcad/trunk/src/libbg/chull3d.cxx 2019-05-20 06:20:41 UTC (rev 73091)
+++ brlcad/trunk/src/libbg/chull3d.cxx 2019-05-20 06:25:40 UTC (rev 73092)
@@ -45,7 +45,7 @@
extern "C" int
bg_3d_chull(int **faces, int *num_faces, point_t **vertices, int *num_vertices,
- const point_t *input_points_3d, int num_input_pnts)
+ const point_t *input_points_3d, int num_input_pnts)
{
int pnt_cnt = 0;
int f_ind = 0;
Modified: brlcad/trunk/src/libbg/clipper.cxx
===================================================================
--- brlcad/trunk/src/libbg/clipper.cxx 2019-05-20 06:20:41 UTC (rev 73091)
+++ brlcad/trunk/src/libbg/clipper.cxx 2019-05-20 06:25:40 UTC (rev 73092)
@@ -50,11 +50,11 @@
namespace ClipperLib {
- static long64 const loRange = 1518500249; //sqrt(2^63 -1)/2
- static long64 const hiRange = 6521908912666391106LL; //sqrt(2^127 -1)/2
- static double const pi = 3.141592653589793238;
- enum Direction { dRightToLeft, dLeftToRight };
- enum RangeTest { rtLo, rtHi, rtError };
+static long64 const loRange = 1518500249; //sqrt(2^63 -1)/2
+static long64 const hiRange = 6521908912666391106LL; //sqrt(2^127 -1)/2
+static double const pi = 3.141592653589793238;
+enum Direction { dRightToLeft, dLeftToRight };
+enum RangeTest { rtLo, rtHi, rtError };
#define HORIZONTAL (-1.0E+40)
#define TOLERANCE (1.0e-20)
@@ -61,3233 +61,3233 @@
#define NEAR_ZERO(val) (((val) > -TOLERANCE) && ((val) < TOLERANCE))
#define NEAR_EQUAL(a, b) NEAR_ZERO((a) - (b))
- inline long64 Abs(long64 val)
+inline long64 Abs(long64 val)
+{
+ if (val < 0) return -val; else return val;
+}
+//------------------------------------------------------------------------------
+
+//------------------------------------------------------------------------------
+// Int128 class (enables safe math on signed 64bit integers)
+// eg Int128 val1((long64)9223372036854775807); //ie 2^63 -1
+// Int128 val2((long64)9223372036854775807);
+// Int128 val3 = val1 * val2;
+// val3.AsString => "85070591730234615847396907784232501249" (8.5e+37)
+//------------------------------------------------------------------------------
+
+class Int128
+{
+public:
+
+ Int128(long64 _lo = 0)
{
- if (val < 0) return -val; else return val;
+ hi = 0;
+ if (_lo < 0) {
+ lo = -_lo;
+ Negate(*this);
+ } else
+ lo = _lo;
}
-
//------------------------------------------------------------------------------
-
//------------------------------------------------------------------------------
- // Int128 class (enables safe math on signed 64bit integers)
- // eg Int128 val1((long64)9223372036854775807); //ie 2^63 -1
- // Int128 val2((long64)9223372036854775807);
- // Int128 val3 = val1 * val2;
- // val3.AsString => "85070591730234615847396907784232501249" (8.5e+37)
-
//------------------------------------------------------------------------------
+ Int128(const Int128 &val): hi(val.hi), lo(val.lo){}
- class Int128
+ long64 operator = (const long64 &val)
{
- public:
+ hi = 0;
+ lo = Abs(val);
+ if (val < 0) Negate(*this);
+ return val;
+ }
- Int128(long64 _lo = 0)
- {
- hi = 0;
- if (_lo < 0) {
- lo = -_lo;
- Negate(*this);
- } else
- lo = _lo;
- }
+ bool operator == (const Int128 &val) const
+ {return (hi == val.hi && lo == val.lo);}
- Int128(const Int128 &val): hi(val.hi), lo(val.lo){}
+ bool operator != (const Int128 &val) const { return !(*this == val);}
- long64 operator = (const long64 &val)
- {
- hi = 0;
- lo = Abs(val);
- if (val < 0) Negate(*this);
- return val;
- }
+ bool operator > (const Int128 &val) const
+ {
+ if (hi > val.hi) return true;
+ else if (hi < val.hi) return false;
+ else return ulong64(lo) > ulong64(val.lo);
+ }
- bool operator == (const Int128 &val) const
- {return (hi == val.hi && lo == val.lo);}
+ bool operator < (const Int128 &val) const
+ {
+ if (hi < val.hi) return true;
+ else if (hi > val.hi) return false;
+ else return ulong64(lo) < ulong64(val.lo);
+ }
- bool operator != (const Int128 &val) const { return !(*this ==
val);}
+ Int128& operator += (const Int128 &rhs)
+ {
+ hi += rhs.hi;
+ lo += rhs.lo;
+ if (ulong64(lo) < ulong64(rhs.lo)) hi++;
+ return *this;
+ }
- bool operator > (const Int128 &val) const
- {
- if (hi > val.hi) return true;
- else if (hi < val.hi) return false;
- else return ulong64(lo) > ulong64(val.lo);
- }
+ Int128 operator + (const Int128 &rhs) const
+ {
+ Int128 result(*this);
+ result+= rhs;
+ return result;
+ }
- bool operator < (const Int128 &val) const
- {
- if (hi < val.hi) return true;
- else if (hi > val.hi) return false;
- else return ulong64(lo) < ulong64(val.lo);
- }
+ Int128& operator -= (const Int128 &rhs)
+ {
+ Int128 tmp(rhs);
+ Negate(tmp);
+ *this += tmp;
+ return *this;
+ }
- Int128& operator += (const Int128 &rhs)
- {
- hi += rhs.hi;
- lo += rhs.lo;
- if (ulong64(lo) < ulong64(rhs.lo)) hi++;
- return *this;
- }
+ Int128 operator - (const Int128 &rhs) const
+ {
+ Int128 result(*this);
+ result-= rhs;
+ return result;
+ }
- Int128 operator + (const Int128 &rhs) const
- {
- Int128 result(*this);
- result+= rhs;
- return result;
- }
+ Int128 operator * (const Int128 &rhs) const {
+ if ( !(hi == 0 || hi == -1) || !(rhs.hi == 0 || rhs.hi == -1))
+ throw "Int128 operator*: overflow error";
+ bool negate = (hi < 0) != (rhs.hi < 0);
- Int128& operator -= (const Int128 &rhs)
- {
- Int128 tmp(rhs);
- Negate(tmp);
- *this += tmp;
- return *this;
- }
+ Int128 tmp(*this);
+ if (tmp.hi < 0) Negate(tmp);
+ ulong64 int1Hi = ulong64(tmp.lo) >> 32;
+ ulong64 int1Lo = ulong64(tmp.lo & 0xFFFFFFFF);
- Int128 operator - (const Int128 &rhs) const
- {
- Int128 result(*this);
- result-= rhs;
- return result;
- }
+ tmp = rhs;
+ if (tmp.hi < 0) Negate(tmp);
+ ulong64 int2Hi = ulong64(tmp.lo) >> 32;
+ ulong64 int2Lo = ulong64(tmp.lo & 0xFFFFFFFF);
- Int128 operator * (const Int128 &rhs) const {
- if ( !(hi == 0 || hi == -1) || !(rhs.hi == 0 || rhs.hi == -1))
- throw "Int128 operator*: overflow error";
- bool negate = (hi < 0) != (rhs.hi < 0);
+ //nb: see comments in clipper.pas
+ ulong64 a = int1Hi * int2Hi;
+ ulong64 b = int1Lo * int2Lo;
+ ulong64 c = int1Hi * int2Lo + int1Lo * int2Hi;
- Int128 tmp(*this);
- if (tmp.hi < 0) Negate(tmp);
- ulong64 int1Hi = ulong64(tmp.lo) >> 32;
- ulong64 int1Lo = ulong64(tmp.lo & 0xFFFFFFFF);
+ tmp.hi = long64(a + (c >> 32));
+ tmp.lo = long64(c << 32);
+ tmp.lo += long64(b);
+ if (ulong64(tmp.lo) < b) tmp.hi++;
+ if (negate) Negate(tmp);
+ return tmp;
+ }
- tmp = rhs;
- if (tmp.hi < 0) Negate(tmp);
- ulong64 int2Hi = ulong64(tmp.lo) >> 32;
- ulong64 int2Lo = ulong64(tmp.lo & 0xFFFFFFFF);
+ Int128 operator/ (const Int128 &rhs) const
+ {
+ if (rhs.lo == 0 && rhs.hi == 0)
+ throw "Int128 operator/: divide by zero";
+ bool negate = (rhs.hi < 0) != (hi < 0);
+ Int128 result(*this), denom(rhs);
+ if (result.hi < 0) Negate(result);
+ if (denom.hi < 0) Negate(denom);
+ if (denom > result) return Int128(0); //result is only a fraction of 1
+ Negate(denom);
- //nb: see comments in clipper.pas
- ulong64 a = int1Hi * int2Hi;
- ulong64 b = int1Lo * int2Lo;
- ulong64 c = int1Hi * int2Lo + int1Lo * int2Hi;
-
- tmp.hi = long64(a + (c >> 32));
- tmp.lo = long64(c << 32);
- tmp.lo += long64(b);
- if (ulong64(tmp.lo) < b) tmp.hi++;
- if (negate) Negate(tmp);
- return tmp;
- }
-
- Int128 operator/ (const Int128 &rhs) const
- {
- if (rhs.lo == 0 && rhs.hi == 0)
- throw "Int128 operator/: divide by zero";
- bool negate = (rhs.hi < 0) != (hi < 0);
- Int128 result(*this), denom(rhs);
- if (result.hi < 0) Negate(result);
- if (denom.hi < 0) Negate(denom);
- if (denom > result) return Int128(0); //result is only a
fraction of 1
- Negate(denom);
-
- Int128 p(0);
- for (int i = 0; i < 128; ++i)
- {
- p.hi = p.hi << 1;
- if (p.lo < 0) p.hi++;
- p.lo = long64(p.lo) << 1;
- if (result.hi < 0) p.lo++;
- result.hi = result.hi << 1;
- if (result.lo < 0) result.hi++;
- result.lo = long64(result.lo) << 1;
- Int128 p2(p);
- p += denom;
- if (p.hi < 0) p = p2;
- else result.lo++;
- }
- if (negate) Negate(result);
- return result;
- }
-
- double AsDouble() const
- {
- const double shift64 = 18446744073709551616.0; //2^64
- const double bit64 = 9223372036854775808.0;
- if (hi < 0)
- {
- Int128 tmp(*this);
- Negate(tmp);
- if (tmp.lo < 0)
- return (double)tmp.lo - bit64 - tmp.hi * shift64;
- else
- return -(double)tmp.lo - tmp.hi * shift64;
- }
- else if (lo < 0)
- return -(double)lo + bit64 + hi * shift64;
- else
- return (double)lo + (double)hi * shift64;
- }
-
- //for bug testing ...
- std::string AsString() const
- {
- std::string result;
- unsigned char r = 0;
- Int128 tmp(0), val(*this);
- if (hi < 0) Negate(val);
- result.resize(50);
- std::string::size_type i = result.size() -1;
- while (val.hi != 0 || val.lo != 0)
- {
- Div10(val, tmp, r);
- result[i--] = char('0' + r);
- val = tmp;
- }
- if (hi < 0) result[i--] = '-';
- result.erase(0,i+1);
- if (result.size() == 0) result = "0";
- return result;
- }
-
- private:
- long64 hi;
- long64 lo;
-
- static void Negate(Int128 &val)
- {
- if (val.lo == 0)
- {
- if( val.hi == 0) return;
- val.lo = ~val.lo;
- val.hi = ~val.hi +1;
- }
- else
- {
- val.lo = ~val.lo +1;
- val.hi = ~val.hi;
- }
- }
-
- //debugging only ...
- void Div10(const Int128 val, Int128& result, unsigned char &
remainder) const
- {
- remainder = 0;
- result = 0;
- for (int i = 63; i >= 0; --i)
- {
- if ((val.hi & ((long64)1 << i)) != 0)
- remainder = char((remainder * 2) + 1); else
- remainder *= char(2);
- if (remainder >= 10)
- {
- result.hi += ((long64)1 << i);
- remainder -= char(10);
- }
- }
- for (int i = 63; i >= 0; --i)
- {
- if ((val.lo & ((long64)1 << i)) != 0)
- remainder = char((remainder * 2) + 1); else
- remainder *= char(2);
- if (remainder >= 10)
- {
- result.lo += ((long64)1 << i);
- remainder -= char(10);
- }
- }
- }
- };
-
-
//------------------------------------------------------------------------------
-
//------------------------------------------------------------------------------
-
- RangeTest TestRange(const Polygon &pts)
- {
- RangeTest result = rtLo;
- for (Polygon::size_type i = 0; i < pts.size(); ++i)
+ Int128 p(0);
+ for (int i = 0; i < 128; ++i)
{
- if (Abs(pts[i].X) > hiRange || Abs(pts[i].Y) > hiRange)
- return rtError;
- else if (Abs(pts[i].X) > loRange || Abs(pts[i].Y) > loRange)
- result = rtHi;
+ p.hi = p.hi << 1;
+ if (p.lo < 0) p.hi++;
+ p.lo = long64(p.lo) << 1;
+ if (result.hi < 0) p.lo++;
+ result.hi = result.hi << 1;
+ if (result.lo < 0) result.hi++;
+ result.lo = long64(result.lo) << 1;
+ Int128 p2(p);
+ p += denom;
+ if (p.hi < 0) p = p2;
+ else result.lo++;
}
+ if (negate) Negate(result);
return result;
}
-
//------------------------------------------------------------------------------
- bool Orientation(const Polygon &poly)
+ double AsDouble() const
{
- int highI = (int)poly.size() -1;
- if (highI < 2) return false;
- bool UseFullInt64Range = false;
-
- int j = 0, jplus, jminus;
- for (int i = 0; i <= highI; ++i)
+ const double shift64 = 18446744073709551616.0; //2^64
+ const double bit64 = 9223372036854775808.0;
+ if (hi < 0)
{
- if (Abs(poly[i].X) > hiRange || Abs(poly[i].Y) > hiRange)
- throw "Coordinate exceeds range bounds.";
- if (Abs(poly[i].X) > loRange || Abs(poly[i].Y) > loRange)
- UseFullInt64Range = true;
- if (poly[i].Y < poly[j].Y) continue;
- if ((poly[i].Y > poly[j].Y || poly[i].X < poly[j].X)) j = i;
- };
- if (j == highI) jplus = 0;
- else jplus = j +1;
- if (j == 0) jminus = highI;
- else jminus = j -1;
-
- IntPoint vec1, vec2;
- //get cross product of vectors of the edges adjacent to highest point
...
- vec1.X = poly[j].X - poly[jminus].X;
- vec1.Y = poly[j].Y - poly[jminus].Y;
- vec2.X = poly[jplus].X - poly[j].X;
- vec2.Y = poly[jplus].Y - poly[j].Y;
-
- if (UseFullInt64Range)
- {
- Int128 cross = Int128(vec1.X) * Int128(vec2.Y) -
- Int128(vec2.X) * Int128(vec1.Y);
- return cross > 0;
+ Int128 tmp(*this);
+ Negate(tmp);
+ if (tmp.lo < 0)
+ return (double)tmp.lo - bit64 - tmp.hi * shift64;
+ else
+ return -(double)tmp.lo - tmp.hi * shift64;
}
+ else if (lo < 0)
+ return -(double)lo + bit64 + hi * shift64;
else
- {
- return (vec1.X * vec2.Y - vec2.X * vec1.Y) > 0;
- }
+ return (double)lo + (double)hi * shift64;
}
-
//------------------------------------------------------------------------------
- bool Orientation(OutRec *outRec, bool UseFullInt64Range)
+ //for bug testing ...
+ std::string AsString() const
{
- OutPt *opBottom = outRec->pts, *op = outRec->pts->next;
- while (op != outRec->pts)
+ std::string result;
+ unsigned char r = 0;
+ Int128 tmp(0), val(*this);
+ if (hi < 0) Negate(val);
+ result.resize(50);
+ std::string::size_type i = result.size() -1;
+ while (val.hi != 0 || val.lo != 0)
{
- if (op->pt.Y >= opBottom->pt.Y)
- {
- if (op->pt.Y > opBottom->pt.Y || op->pt.X < opBottom->pt.X)
- opBottom = op;
- }
- op = op->next;
+ Div10(val, tmp, r);
+ result[i--] = char('0' + r);
+ val = tmp;
}
-
- IntPoint vec1, vec2;
- vec1.X = op->pt.X - op->prev->pt.X;
- vec1.Y = op->pt.Y - op->prev->pt.Y;
- vec2.X = op->next->pt.X - op->pt.X;
- vec2.Y = op->next->pt.Y - op->pt.Y;
-
- if (UseFullInt64Range)
- {
- Int128 cross = Int128(vec1.X) * Int128(vec2.Y) - Int128(vec2.X) *
Int128(vec1.Y);
- return cross > 0;
- }
- else
- {
- return (vec1.X * vec2.Y - vec2.X * vec1.Y) > 0;
- }
+ if (hi < 0) result[i--] = '-';
+ result.erase(0,i+1);
+ if (result.size() == 0) result = "0";
+ return result;
}
-
//------------------------------------------------------------------------------
- inline bool PointsEqual( const IntPoint &pt1, const IntPoint &pt2)
- {
- return ( pt1.X == pt2.X && pt1.Y == pt2.Y );
- }
-
//------------------------------------------------------------------------------
+private:
+ long64 hi;
+ long64 lo;
- double Area(const Polygon &poly)
+ static void Negate(Int128 &val)
{
- int highI = (int)poly.size() -1;
- if (highI < 2) return 0;
- bool UseFullInt64Range;
- RangeTest rt = TestRange(poly);
- switch (rt) {
- case rtLo:
- UseFullInt64Range = false;
- break;
- case rtHi:
- UseFullInt64Range = true;
- break;
- default:
- throw "Coordinate exceeds range bounds.";
+ if (val.lo == 0)
+ {
+ if( val.hi == 0) return;
+ val.lo = ~val.lo;
+ val.hi = ~val.hi +1;
}
-
- if (UseFullInt64Range) {
- Int128 a(0);
- a = (Int128(poly[highI].X) * Int128(poly[0].Y)) -
- Int128(poly[0].X) * Int128(poly[highI].Y);
- for (int i = 0; i < highI; ++i)
- a += Int128(poly[i].X) * Int128(poly[i+1].Y) -
- Int128(poly[i+1].X) * Int128(poly[i].Y);
- return a.AsDouble() / 2;
- }
else
{
- double a;
- a = (double)poly[highI].X * poly[0].Y - (double)poly[0].X *
poly[highI].Y;
- for (int i = 0; i < highI; ++i)
- a += (double)poly[i].X * poly[i+1].Y - (double)poly[i+1].X *
poly[i].Y;
- return a/2;
+ val.lo = ~val.lo +1;
+ val.hi = ~val.hi;
}
}
-
//------------------------------------------------------------------------------
- bool PointIsVertex(const IntPoint &pt, OutPt *pp)
+ //debugging only ...
+ void Div10(const Int128 val, Int128& result, unsigned char & remainder)
const
{
- OutPt *pp2 = pp;
- do
+ remainder = 0;
+ result = 0;
+ for (int i = 63; i >= 0; --i)
{
- if (PointsEqual(pp2->pt, pt)) return true;
- pp2 = pp2->next;
- }
- while (pp2 != pp);
- return false;
- }
-
//------------------------------------------------------------------------------
-
- bool PointInPolygon(const IntPoint &pt, OutPt *pp, bool UseFullInt64Range)
- {
- OutPt *pp2 = pp;
- bool result = false;
- if (UseFullInt64Range) {
- do
+ if ((val.hi & ((long64)1 << i)) != 0)
+ remainder = char((remainder * 2) + 1); else
+ remainder *= char(2);
+ if (remainder >= 10)
{
- if ((((pp2->pt.Y <= pt.Y) && (pt.Y < pp2->prev->pt.Y)) ||
- ((pp2->prev->pt.Y <= pt.Y) && (pt.Y < pp2->pt.Y)))
&&
- Int128(pt.X - pp2->pt.X) < (Int128(pp2->prev->pt.X -
pp2->pt.X) *
- Int128(pt.Y - pp2->pt.Y)) / Int128(pp2->prev->pt.Y
- pp2->pt.Y))
- result = !result;
- pp2 = pp2->next;
+ result.hi += ((long64)1 << i);
+ remainder -= char(10);
}
- while (pp2 != pp);
}
- else
+ for (int i = 63; i >= 0; --i)
{
- do
+ if ((val.lo & ((long64)1 << i)) != 0)
+ remainder = char((remainder * 2) + 1); else
+ remainder *= char(2);
+ if (remainder >= 10)
{
- if ((((pp2->pt.Y <= pt.Y) && (pt.Y < pp2->prev->pt.Y)) ||
- ((pp2->prev->pt.Y <= pt.Y) && (pt.Y < pp2->pt.Y)))
&&
- (pt.X < (pp2->prev->pt.X - pp2->pt.X) * (pt.Y -
pp2->pt.Y) /
- (pp2->prev->pt.Y - pp2->pt.Y) + pp2->pt.X )) result =
!result;
- pp2 = pp2->next;
+ result.lo += ((long64)1 << i);
+ remainder -= char(10);
}
- while (pp2 != pp);
}
- return result;
}
-
//------------------------------------------------------------------------------
+};
- bool SlopesEqual(TEdge &e1, TEdge &e2, bool UseFullInt64Range)
+//------------------------------------------------------------------------------
+//------------------------------------------------------------------------------
+
+RangeTest TestRange(const Polygon &pts)
+{
+ RangeTest result = rtLo;
+ for (Polygon::size_type i = 0; i < pts.size(); ++i)
{
- if (e1.ybot == e1.ytop) return (e2.ybot == e2.ytop);
- else if (e1.xbot == e1.xtop) return (e2.xbot == e2.xtop);
- else if (UseFullInt64Range)
- return Int128(e1.ytop - e1.ybot) * Int128(e2.xtop - e2.xbot) ==
- Int128(e1.xtop - e1.xbot) * Int128(e2.ytop - e2.ybot);
- else return (e1.ytop - e1.ybot)*(e2.xtop - e2.xbot) ==
- (e1.xtop - e1.xbot)*(e2.ytop - e2.ybot);
+ if (Abs(pts[i].X) > hiRange || Abs(pts[i].Y) > hiRange)
+ return rtError;
+ else if (Abs(pts[i].X) > loRange || Abs(pts[i].Y) > loRange)
+ result = rtHi;
}
-
//------------------------------------------------------------------------------
+ return result;
+}
+//------------------------------------------------------------------------------
- bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
- const IntPoint pt3, bool UseFullInt64Range)
+bool Orientation(const Polygon &poly)
+{
+ int highI = (int)poly.size() -1;
+ if (highI < 2) return false;
+ bool UseFullInt64Range = false;
+
+ int j = 0, jplus, jminus;
+ for (int i = 0; i <= highI; ++i)
{
- if (pt1.Y == pt2.Y) return (pt2.Y == pt3.Y);
- else if (pt1.X == pt2.X) return (pt2.X == pt3.X);
- else if (UseFullInt64Range)
- return Int128(pt1.Y-pt2.Y) * Int128(pt2.X-pt3.X) ==
- Int128(pt1.X-pt2.X) * Int128(pt2.Y-pt3.Y);
- else return (pt1.Y-pt2.Y)*(pt2.X-pt3.X) == (pt1.X-pt2.X)*(pt2.Y-pt3.Y);
- }
-
//------------------------------------------------------------------------------
+ if (Abs(poly[i].X) > hiRange || Abs(poly[i].Y) > hiRange)
+ throw "Coordinate exceeds range bounds.";
+ if (Abs(poly[i].X) > loRange || Abs(poly[i].Y) > loRange)
+ UseFullInt64Range = true;
+ if (poly[i].Y < poly[j].Y) continue;
+ if ((poly[i].Y > poly[j].Y || poly[i].X < poly[j].X)) j = i;
+ };
+ if (j == highI) jplus = 0;
+ else jplus = j +1;
+ if (j == 0) jminus = highI;
+ else jminus = j -1;
- bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
- const IntPoint pt3, const IntPoint pt4, bool UseFullInt64Range)
+ IntPoint vec1, vec2;
+ //get cross product of vectors of the edges adjacent to highest point ...
+ vec1.X = poly[j].X - poly[jminus].X;
+ vec1.Y = poly[j].Y - poly[jminus].Y;
+ vec2.X = poly[jplus].X - poly[j].X;
+ vec2.Y = poly[jplus].Y - poly[j].Y;
+
+ if (UseFullInt64Range)
{
- if (pt1.Y == pt2.Y) return (pt3.Y == pt4.Y);
- else if (pt1.X == pt2.X) return (pt3.X == pt4.X);
- else if (UseFullInt64Range)
- return Int128(pt1.Y-pt2.Y) * Int128(pt3.X-pt4.X) ==
- Int128(pt1.X-pt2.X) * Int128(pt3.Y-pt4.Y);
- else return (pt1.Y-pt2.Y)*(pt3.X-pt4.X) == (pt1.X-pt2.X)*(pt3.Y-pt4.Y);
+ Int128 cross = Int128(vec1.X) * Int128(vec2.Y) -
+ Int128(vec2.X) * Int128(vec1.Y);
+ return cross > 0;
}
-
//------------------------------------------------------------------------------
-
- double GetDx(const IntPoint pt1, const IntPoint pt2)
+ else
{
- if (pt1.Y == pt2.Y) return HORIZONTAL;
- else return
- (double)(pt2.X - pt1.X) / (double)(pt2.Y - pt1.Y);
+ return (vec1.X * vec2.Y - vec2.X * vec1.Y) > 0;
}
-
//---------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- void SetDx(TEdge &e)
+bool Orientation(OutRec *outRec, bool UseFullInt64Range)
+{
+ OutPt *opBottom = outRec->pts, *op = outRec->pts->next;
+ while (op != outRec->pts)
{
- if (e.ybot == e.ytop) e.dx = HORIZONTAL;
- else e.dx =
- (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot);
+ if (op->pt.Y >= opBottom->pt.Y)
+ {
+ if (op->pt.Y > opBottom->pt.Y || op->pt.X < opBottom->pt.X)
+ opBottom = op;
+ }
+ op = op->next;
}
-
//---------------------------------------------------------------------------
- void SwapSides(TEdge &edge1, TEdge &edge2)
+ IntPoint vec1, vec2;
+ vec1.X = op->pt.X - op->prev->pt.X;
+ vec1.Y = op->pt.Y - op->prev->pt.Y;
+ vec2.X = op->next->pt.X - op->pt.X;
+ vec2.Y = op->next->pt.Y - op->pt.Y;
+
+ if (UseFullInt64Range)
{
- EdgeSide side = edge1.side;
- edge1.side = edge2.side;
- edge2.side = side;
+ Int128 cross = Int128(vec1.X) * Int128(vec2.Y) - Int128(vec2.X) *
Int128(vec1.Y);
+ return cross > 0;
}
-
//------------------------------------------------------------------------------
-
- void SwapPolyIndexes(TEdge &edge1, TEdge &edge2)
+ else
{
- int outIdx = edge1.outIdx;
- edge1.outIdx = edge2.outIdx;
- edge2.outIdx = outIdx;
+ return (vec1.X * vec2.Y - vec2.X * vec1.Y) > 0;
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- inline long64 Round(double val)
+inline bool PointsEqual( const IntPoint &pt1, const IntPoint &pt2)
+{
+ return ( pt1.X == pt2.X && pt1.Y == pt2.Y );
+}
+//------------------------------------------------------------------------------
+
+double Area(const Polygon &poly)
+{
+ int highI = (int)poly.size() -1;
+ if (highI < 2) return 0;
+ bool UseFullInt64Range;
+ RangeTest rt = TestRange(poly);
+ switch (rt) {
+ case rtLo:
+ UseFullInt64Range = false;
+ break;
+ case rtHi:
+ UseFullInt64Range = true;
+ break;
+ default:
+ throw "Coordinate exceeds range bounds.";
+ }
+
+ if (UseFullInt64Range) {
+ Int128 a(0);
+ a = (Int128(poly[highI].X) * Int128(poly[0].Y)) -
+ Int128(poly[0].X) * Int128(poly[highI].Y);
+ for (int i = 0; i < highI; ++i)
+ a += Int128(poly[i].X) * Int128(poly[i+1].Y) -
+ Int128(poly[i+1].X) * Int128(poly[i].Y);
+ return a.AsDouble() / 2;
+ }
+ else
{
- if ((val < 0)) return static_cast<long64>(val - 0.5);
- else return static_cast<long64>(val + 0.5);
+ double a;
+ a = (double)poly[highI].X * poly[0].Y - (double)poly[0].X *
poly[highI].Y;
+ for (int i = 0; i < highI; ++i)
+ a += (double)poly[i].X * poly[i+1].Y - (double)poly[i+1].X *
poly[i].Y;
+ return a/2;
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- long64 TopX(TEdge &edge, const long64 currentY)
+bool PointIsVertex(const IntPoint &pt, OutPt *pp)
+{
+ OutPt *pp2 = pp;
+ do
{
- if( currentY == edge.ytop ) return edge.xtop;
- return edge.xbot + Round(edge.dx *(currentY - edge.ybot));
+ if (PointsEqual(pp2->pt, pt)) return true;
+ pp2 = pp2->next;
}
-
//------------------------------------------------------------------------------
+ while (pp2 != pp);
+ return false;
+}
+//------------------------------------------------------------------------------
- long64 TopX(const IntPoint pt1, const IntPoint pt2, const long64 currentY)
- {
- //preconditions: pt1.Y <> pt2.Y and pt1.Y > pt2.Y
- if (currentY >= pt1.Y) return pt1.X;
- else if (currentY == pt2.Y) return pt2.X;
- else if (pt1.X == pt2.X) return pt1.X;
- else
+bool PointInPolygon(const IntPoint &pt, OutPt *pp, bool UseFullInt64Range)
+{
+ OutPt *pp2 = pp;
+ bool result = false;
+ if (UseFullInt64Range) {
+ do
{
- double q = (double)(pt1.X-pt2.X)/(double)(pt1.Y-pt2.Y);
- return Round(pt1.X + (currentY - pt1.Y) *q);
+ if ((((pp2->pt.Y <= pt.Y) && (pt.Y < pp2->prev->pt.Y)) ||
+ ((pp2->prev->pt.Y <= pt.Y) && (pt.Y < pp2->pt.Y))) &&
+ Int128(pt.X - pp2->pt.X) < (Int128(pp2->prev->pt.X - pp2->pt.X)
*
+ Int128(pt.Y - pp2->pt.Y)) /
Int128(pp2->prev->pt.Y - pp2->pt.Y))
+ result = !result;
+ pp2 = pp2->next;
}
+ while (pp2 != pp);
}
-
//------------------------------------------------------------------------------
-
- bool IntersectPoint(TEdge &edge1, TEdge &edge2,
- IntPoint &ip, bool UseFullInt64Range)
+ else
{
- double b1, b2;
- if (SlopesEqual(edge1, edge2, UseFullInt64Range)) return false;
- else if (NEAR_ZERO(edge1.dx))
+ do
{
- ip.X = edge1.xbot;
- if (NEAR_EQUAL(edge2.dx, HORIZONTAL))
- {
- ip.Y = edge2.ybot;
- } else
- {
- b2 = edge2.ybot - (edge2.xbot/edge2.dx);
- ip.Y = Round(ip.X/edge2.dx + b2);
- }
+ if ((((pp2->pt.Y <= pt.Y) && (pt.Y < pp2->prev->pt.Y)) ||
+ ((pp2->prev->pt.Y <= pt.Y) && (pt.Y < pp2->pt.Y))) &&
+ (pt.X < (pp2->prev->pt.X - pp2->pt.X) * (pt.Y - pp2->pt.Y) /
+ (pp2->prev->pt.Y - pp2->pt.Y) + pp2->pt.X )) result = !result;
+ pp2 = pp2->next;
}
- else if (NEAR_ZERO(edge2.dx))
- {
- ip.X = edge2.xbot;
- if (NEAR_EQUAL(edge1.dx, HORIZONTAL))
- {
- ip.Y = edge1.ybot;
- } else
- {
- b1 = edge1.ybot - (edge1.xbot/edge1.dx);
- ip.Y = Round(ip.X/edge1.dx + b1);
- }
- } else
- {
- b1 = edge1.xbot - edge1.ybot * edge1.dx;
- b2 = edge2.xbot - edge2.ybot * edge2.dx;
- b2 = (b2-b1)/(edge1.dx - edge2.dx);
- ip.Y = Round(b2);
- ip.X = Round(edge1.dx * b2 + b1);
- }
-
- return
- //can be *so close* to the top of one edge that the rounded Y
equals one ytop ...
- (ip.Y == edge1.ytop && ip.Y >= edge2.ytop && edge1.tmpX >
edge2.tmpX) ||
- (ip.Y == edge2.ytop && ip.Y >= edge1.ytop && edge1.tmpX >
edge2.tmpX) ||
- (ip.Y > edge1.ytop && ip.Y > edge2.ytop);
+ while (pp2 != pp);
}
-
//------------------------------------------------------------------------------
+ return result;
+}
+//------------------------------------------------------------------------------
- void ReversePolyPtLinks(OutPt &pp)
+bool SlopesEqual(TEdge &e1, TEdge &e2, bool UseFullInt64Range)
+{
+ if (e1.ybot == e1.ytop) return (e2.ybot == e2.ytop);
+ else if (e1.xbot == e1.xtop) return (e2.xbot == e2.xtop);
+ else if (UseFullInt64Range)
+ return Int128(e1.ytop - e1.ybot) * Int128(e2.xtop - e2.xbot) ==
+ Int128(e1.xtop - e1.xbot) * Int128(e2.ytop - e2.ybot);
+ else return (e1.ytop - e1.ybot)*(e2.xtop - e2.xbot) ==
+ (e1.xtop - e1.xbot)*(e2.ytop - e2.ybot);
+}
+//------------------------------------------------------------------------------
+
+bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
+ const IntPoint pt3, bool UseFullInt64Range)
+{
+ if (pt1.Y == pt2.Y) return (pt2.Y == pt3.Y);
+ else if (pt1.X == pt2.X) return (pt2.X == pt3.X);
+ else if (UseFullInt64Range)
+ return Int128(pt1.Y-pt2.Y) * Int128(pt2.X-pt3.X) ==
+ Int128(pt1.X-pt2.X) * Int128(pt2.Y-pt3.Y);
+ else return (pt1.Y-pt2.Y)*(pt2.X-pt3.X) == (pt1.X-pt2.X)*(pt2.Y-pt3.Y);
+}
+//------------------------------------------------------------------------------
+
+bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
+ const IntPoint pt3, const IntPoint pt4, bool UseFullInt64Range)
+{
+ if (pt1.Y == pt2.Y) return (pt3.Y == pt4.Y);
+ else if (pt1.X == pt2.X) return (pt3.X == pt4.X);
+ else if (UseFullInt64Range)
+ return Int128(pt1.Y-pt2.Y) * Int128(pt3.X-pt4.X) ==
+ Int128(pt1.X-pt2.X) * Int128(pt3.Y-pt4.Y);
+ else return (pt1.Y-pt2.Y)*(pt3.X-pt4.X) == (pt1.X-pt2.X)*(pt3.Y-pt4.Y);
+}
+//------------------------------------------------------------------------------
+
+double GetDx(const IntPoint pt1, const IntPoint pt2)
+{
+ if (pt1.Y == pt2.Y) return HORIZONTAL;
+ else return
+ (double)(pt2.X - pt1.X) / (double)(pt2.Y - pt1.Y);
+}
+//---------------------------------------------------------------------------
+
+void SetDx(TEdge &e)
+{
+ if (e.ybot == e.ytop) e.dx = HORIZONTAL;
+ else e.dx =
+ (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot);
+}
+//---------------------------------------------------------------------------
+
+void SwapSides(TEdge &edge1, TEdge &edge2)
+{
+ EdgeSide side = edge1.side;
+ edge1.side = edge2.side;
+ edge2.side = side;
+}
+//------------------------------------------------------------------------------
+
+void SwapPolyIndexes(TEdge &edge1, TEdge &edge2)
+{
+ int outIdx = edge1.outIdx;
+ edge1.outIdx = edge2.outIdx;
+ edge2.outIdx = outIdx;
+}
+//------------------------------------------------------------------------------
+
+inline long64 Round(double val)
+{
+ if ((val < 0)) return static_cast<long64>(val - 0.5);
+ else return static_cast<long64>(val + 0.5);
+}
+//------------------------------------------------------------------------------
+
+long64 TopX(TEdge &edge, const long64 currentY)
+{
+ if( currentY == edge.ytop ) return edge.xtop;
+ return edge.xbot + Round(edge.dx *(currentY - edge.ybot));
+}
+//------------------------------------------------------------------------------
+
+long64 TopX(const IntPoint pt1, const IntPoint pt2, const long64 currentY)
+{
+ //preconditions: pt1.Y <> pt2.Y and pt1.Y > pt2.Y
+ if (currentY >= pt1.Y) return pt1.X;
+ else if (currentY == pt2.Y) return pt2.X;
+ else if (pt1.X == pt2.X) return pt1.X;
+ else
{
- OutPt *pp1, *pp2;
- pp1 = &pp;
- do {
- pp2 = pp1->next;
- pp1->next = pp1->prev;
- pp1->prev = pp2;
- pp1 = pp2;
- } while( pp1 != &pp );
+ double q = (double)(pt1.X-pt2.X)/(double)(pt1.Y-pt2.Y);
+ return Round(pt1.X + (currentY - pt1.Y) *q);
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- void DisposeOutPts(OutPt*& pp)
+bool IntersectPoint(TEdge &edge1, TEdge &edge2,
+ IntPoint &ip, bool UseFullInt64Range)
+{
+ double b1, b2;
+ if (SlopesEqual(edge1, edge2, UseFullInt64Range)) return false;
+ else if (NEAR_ZERO(edge1.dx))
{
- if (pp == 0) return;
- pp->prev->next = 0;
- while( pp )
+ ip.X = edge1.xbot;
+ if (NEAR_EQUAL(edge2.dx, HORIZONTAL))
{
- OutPt *tmpPp = pp;
- pp = pp->next;
- delete tmpPp ;
+ ip.Y = edge2.ybot;
+ } else
+ {
+ b2 = edge2.ybot - (edge2.xbot/edge2.dx);
+ ip.Y = Round(ip.X/edge2.dx + b2);
}
}
-
//------------------------------------------------------------------------------
-
- void InitEdge(TEdge *e, TEdge *eNext,
- TEdge *ePrev, const IntPoint &pt, PolyType polyType)
+ else if (NEAR_ZERO(edge2.dx))
{
- std::memset( e, 0, sizeof( TEdge ));
-
- e->next = eNext;
- e->prev = ePrev;
- e->xcurr = pt.X;
- e->ycurr = pt.Y;
- if (e->ycurr >= e->next->ycurr)
+ ip.X = edge2.xbot;
+ if (NEAR_EQUAL(edge1.dx, HORIZONTAL))
{
- e->xbot = e->xcurr;
- e->ybot = e->ycurr;
- e->xtop = e->next->xcurr;
- e->ytop = e->next->ycurr;
- e->windDelta = 1;
+ ip.Y = edge1.ybot;
} else
{
- e->xtop = e->xcurr;
- e->ytop = e->ycurr;
- e->xbot = e->next->xcurr;
- e->ybot = e->next->ycurr;
- e->windDelta = -1;
+ b1 = edge1.ybot - (edge1.xbot/edge1.dx);
+ ip.Y = Round(ip.X/edge1.dx + b1);
}
- SetDx(*e);
- e->polyType = polyType;
- e->outIdx = -1;
- }
-
//------------------------------------------------------------------------------
-
- inline void SwapX(TEdge &e)
+ } else
{
- //swap horizontal edges' top and bottom x's so they follow the natural
- //progression of the bounds - ie so their xbots will align with the
- //adjoining lower edge. [Helpful in the ProcessHorizontal() method.]
- e.xcurr = e.xtop;
- e.xtop = e.xbot;
- e.xbot = e.xcurr;
+ b1 = edge1.xbot - edge1.ybot * edge1.dx;
+ b2 = edge2.xbot - edge2.ybot * edge2.dx;
+ b2 = (b2-b1)/(edge1.dx - edge2.dx);
+ ip.Y = Round(b2);
+ ip.X = Round(edge1.dx * b2 + b1);
}
-
//------------------------------------------------------------------------------
- void SwapPoints(IntPoint &pt1, IntPoint &pt2)
- {
- IntPoint tmp = pt1;
- pt1 = pt2;
- pt2 = tmp;
- }
-
//------------------------------------------------------------------------------
+ return
+ //can be *so close* to the top of one edge that the rounded Y equals
one ytop ...
+ (ip.Y == edge1.ytop && ip.Y >= edge2.ytop && edge1.tmpX > edge2.tmpX) ||
+ (ip.Y == edge2.ytop && ip.Y >= edge1.ytop && edge1.tmpX > edge2.tmpX) ||
+ (ip.Y > edge1.ytop && ip.Y > edge2.ytop);
+}
+//------------------------------------------------------------------------------
- bool GetOverlapSegment(IntPoint pt1a, IntPoint pt1b, IntPoint pt2a,
- IntPoint pt2b, IntPoint &pt1, IntPoint &pt2)
+void ReversePolyPtLinks(OutPt &pp)
+{
+ OutPt *pp1, *pp2;
+ pp1 = &pp;
+ do {
+ pp2 = pp1->next;
+ pp1->next = pp1->prev;
+ pp1->prev = pp2;
+ pp1 = pp2;
+ } while( pp1 != &pp );
+}
+//------------------------------------------------------------------------------
+
+void DisposeOutPts(OutPt*& pp)
+{
+ if (pp == 0) return;
+ pp->prev->next = 0;
+ while( pp )
{
- //precondition: segments are colinear.
- if ( pt1a.Y == pt1b.Y || Abs((pt1a.X - pt1b.X)/(pt1a.Y - pt1b.Y)) > 1 )
- {
- if (pt1a.X > pt1b.X) SwapPoints(pt1a, pt1b);
- if (pt2a.X > pt2b.X) SwapPoints(pt2a, pt2b);
- if (pt1a.X > pt2a.X) pt1 = pt1a; else pt1 = pt2a;
- if (pt1b.X < pt2b.X) pt2 = pt1b; else pt2 = pt2b;
- return pt1.X < pt2.X;
- } else
- {
- if (pt1a.Y < pt1b.Y) SwapPoints(pt1a, pt1b);
- if (pt2a.Y < pt2b.Y) SwapPoints(pt2a, pt2b);
- if (pt1a.Y < pt2a.Y) pt1 = pt1a; else pt1 = pt2a;
- if (pt1b.Y > pt2b.Y) pt2 = pt1b; else pt2 = pt2b;
- return pt1.Y > pt2.Y;
- }
+ OutPt *tmpPp = pp;
+ pp = pp->next;
+ delete tmpPp ;
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- OutPt* PolygonBottom(OutPt* pp)
+void InitEdge(TEdge *e, TEdge *eNext,
+ TEdge *ePrev, const IntPoint &pt, PolyType polyType)
+{
+ std::memset( e, 0, sizeof( TEdge ));
+
+ e->next = eNext;
+ e->prev = ePrev;
+ e->xcurr = pt.X;
+ e->ycurr = pt.Y;
+ if (e->ycurr >= e->next->ycurr)
{
- OutPt* p = pp->next;
- OutPt* result = pp;
- while (p != pp)
- {
- if (p->pt.Y > result->pt.Y) result = p;
- else if (p->pt.Y == result->pt.Y && p->pt.X < result->pt.X) result
= p;
- p = p->next;
- }
- return result;
+ e->xbot = e->xcurr;
+ e->ybot = e->ycurr;
+ e->xtop = e->next->xcurr;
+ e->ytop = e->next->ycurr;
+ e->windDelta = 1;
+ } else
+ {
+ e->xtop = e->xcurr;
+ e->ytop = e->ycurr;
+ e->xbot = e->next->xcurr;
+ e->ybot = e->next->ycurr;
+ e->windDelta = -1;
}
-
//------------------------------------------------------------------------------
+ SetDx(*e);
+ e->polyType = polyType;
+ e->outIdx = -1;
+}
+//------------------------------------------------------------------------------
- bool FindSegment(OutPt* &pp, IntPoint &pt1, IntPoint &pt2)
+inline void SwapX(TEdge &e)
+{
+ //swap horizontal edges' top and bottom x's so they follow the natural
+ //progression of the bounds - ie so their xbots will align with the
+ //adjoining lower edge. [Helpful in the ProcessHorizontal() method.]
+ e.xcurr = e.xtop;
+ e.xtop = e.xbot;
+ e.xbot = e.xcurr;
+}
+//------------------------------------------------------------------------------
+
+void SwapPoints(IntPoint &pt1, IntPoint &pt2)
+{
+ IntPoint tmp = pt1;
+ pt1 = pt2;
+ pt2 = tmp;
+}
+//------------------------------------------------------------------------------
+
+bool GetOverlapSegment(IntPoint pt1a, IntPoint pt1b, IntPoint pt2a,
+ IntPoint pt2b, IntPoint &pt1, IntPoint &pt2)
+{
+ //precondition: segments are colinear.
+ if ( pt1a.Y == pt1b.Y || Abs((pt1a.X - pt1b.X)/(pt1a.Y - pt1b.Y)) > 1 )
{
- //outPt1 & outPt2 => the overlap segment (if the function returns true)
- if (!pp) return false;
- OutPt* pp2 = pp;
- IntPoint pt1a = pt1, pt2a = pt2;
- do
- {
- if (SlopesEqual(pt1a, pt2a, pp->pt, pp->prev->pt, true) &&
- SlopesEqual(pt1a, pt2a, pp->pt, true) &&
- GetOverlapSegment(pt1a, pt2a, pp->pt, pp->prev->pt, pt1,
pt2))
- return true;
- pp = pp->next;
- }
- while (pp != pp2);
- return false;
+ if (pt1a.X > pt1b.X) SwapPoints(pt1a, pt1b);
+ if (pt2a.X > pt2b.X) SwapPoints(pt2a, pt2b);
+ if (pt1a.X > pt2a.X) pt1 = pt1a; else pt1 = pt2a;
+ if (pt1b.X < pt2b.X) pt2 = pt1b; else pt2 = pt2b;
+ return pt1.X < pt2.X;
+ } else
+ {
+ if (pt1a.Y < pt1b.Y) SwapPoints(pt1a, pt1b);
+ if (pt2a.Y < pt2b.Y) SwapPoints(pt2a, pt2b);
+ if (pt1a.Y < pt2a.Y) pt1 = pt1a; else pt1 = pt2a;
+ if (pt1b.Y > pt2b.Y) pt2 = pt1b; else pt2 = pt2b;
+ return pt1.Y > pt2.Y;
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- bool Pt3IsBetweenPt1AndPt2(const IntPoint pt1,
- const IntPoint pt2, const IntPoint pt3)
+OutPt* PolygonBottom(OutPt* pp)
+{
+ OutPt* p = pp->next;
+ OutPt* result = pp;
+ while (p != pp)
{
- if (PointsEqual(pt1, pt3) || PointsEqual(pt2, pt3)) return true;
- else if (pt1.X != pt2.X) return (pt1.X < pt3.X) == (pt3.X < pt2.X);
- else return (pt1.Y < pt3.Y) == (pt3.Y < pt2.Y);
+ if (p->pt.Y > result->pt.Y) result = p;
+ else if (p->pt.Y == result->pt.Y && p->pt.X < result->pt.X) result = p;
+ p = p->next;
}
-
//------------------------------------------------------------------------------
+ return result;
+}
+//------------------------------------------------------------------------------
- OutPt* InsertPolyPtBetween(OutPt* p1, OutPt* p2, const IntPoint pt)
+bool FindSegment(OutPt* &pp, IntPoint &pt1, IntPoint &pt2)
+{
+ //outPt1 & outPt2 => the overlap segment (if the function returns true)
+ if (!pp) return false;
+ OutPt* pp2 = pp;
+ IntPoint pt1a = pt1, pt2a = pt2;
+ do
{
- if (p1 == p2) throw "JoinError";
- OutPt* result = new OutPt;
- result->pt = pt;
- if (p2 == p1->next)
- {
- p1->next = result;
- p2->prev = result;
- result->next = p2;
- result->prev = p1;
- } else
- {
- p2->next = result;
- p1->prev = result;
- result->next = p1;
- result->prev = p2;
- }
- return result;
+ if (SlopesEqual(pt1a, pt2a, pp->pt, pp->prev->pt, true) &&
+ SlopesEqual(pt1a, pt2a, pp->pt, true) &&
+ GetOverlapSegment(pt1a, pt2a, pp->pt, pp->prev->pt, pt1, pt2))
+ return true;
+ pp = pp->next;
}
+ while (pp != pp2);
+ return false;
+}
+//------------------------------------------------------------------------------
-
//------------------------------------------------------------------------------
- // ClipperBase class methods ...
-
//------------------------------------------------------------------------------
+bool Pt3IsBetweenPt1AndPt2(const IntPoint pt1,
+ const IntPoint pt2, const IntPoint pt3)
+{
+ if (PointsEqual(pt1, pt3) || PointsEqual(pt2, pt3)) return true;
+ else if (pt1.X != pt2.X) return (pt1.X < pt3.X) == (pt3.X < pt2.X);
+ else return (pt1.Y < pt3.Y) == (pt3.Y < pt2.Y);
+}
+//------------------------------------------------------------------------------
- ClipperBase::ClipperBase() //constructor
+OutPt* InsertPolyPtBetween(OutPt* p1, OutPt* p2, const IntPoint pt)
+{
+ if (p1 == p2) throw "JoinError";
+ OutPt* result = new OutPt;
+ result->pt = pt;
+ if (p2 == p1->next)
{
- m_MinimaList = 0;
- m_CurrentLM = 0;
- m_UseFullRange = true;
- m_edges = new EdgeList();
- }
-
//------------------------------------------------------------------------------
-
- ClipperBase::~ClipperBase() //destructor
+ p1->next = result;
+ p2->prev = result;
+ result->next = p2;
+ result->prev = p1;
+ } else
{
- Clear();
- delete m_edges;
- m_edges = NULL;
+ p2->next = result;
+ p1->prev = result;
+ result->next = p1;
+ result->prev = p2;
}
-
//------------------------------------------------------------------------------
+ return result;
+}
- bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType)
- {
- int len = (int)pg.size();
- if (len < 3) return false;
- Polygon p(len);
- p[0] = pg[0];
- int j = 0;
+//------------------------------------------------------------------------------
+// ClipperBase class methods ...
+//------------------------------------------------------------------------------
- long64 maxVal;
- if (m_UseFullRange) maxVal = hiRange; else maxVal = loRange;
+ClipperBase::ClipperBase() //constructor
+{
+ m_MinimaList = 0;
+ m_CurrentLM = 0;
+ m_UseFullRange = true;
+ m_edges = new EdgeList();
+}
+//------------------------------------------------------------------------------
- for (int i = 0; i < len; ++i)
- {
- if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal)
- {
- if (m_UseFullRange)
- throw "Coordinate exceeds range bounds";
- maxVal = hiRange;
- if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal)
- throw "Coordinate exceeds range bounds";
- m_UseFullRange = true;
- }
+ClipperBase::~ClipperBase() //destructor
+{
+ Clear();
+ delete m_edges;
+ m_edges = NULL;
+}
+//------------------------------------------------------------------------------
- if (i == 0 || PointsEqual(p[j], pg[i])) continue;
- else if (j > 0 && SlopesEqual(p[j-1], p[j], pg[i], m_UseFullRange))
- {
- if (PointsEqual(p[j-1], pg[i])) j--;
- } else j++;
- p[j] = pg[i];
- }
- if (j < 2) return false;
+bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType)
+{
+ int len = (int)pg.size();
+ if (len < 3) return false;
+ Polygon p(len);
+ p[0] = pg[0];
+ int j = 0;
- len = j+1;
- for (;;)
+ long64 maxVal;
+ if (m_UseFullRange) maxVal = hiRange; else maxVal = loRange;
+
+ for (int i = 0; i < len; ++i)
+ {
+ if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal)
{
- //nb: test for point equality before testing slopes ...
- if (PointsEqual(p[j], p[0])) j--;
- else if (PointsEqual(p[0], p[1]) ||
- SlopesEqual(p[j], p[0], p[1], m_UseFullRange))
- p[0] = p[j--];
- else if (SlopesEqual(p[j-1], p[j], p[0], m_UseFullRange)) j--;
- else if (SlopesEqual(p[0], p[1], p[2], m_UseFullRange))
- {
- for (int i = 2; i <= j; ++i) p[i-1] = p[i];
- j--;
- }
- //exit loop if nothing is changed or there are too few vertices ...
- if (j == len-1 || j < 2) break;
- len = j +1;
+ if (m_UseFullRange)
+ throw "Coordinate exceeds range bounds";
+ maxVal = hiRange;
+ if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal)
+ throw "Coordinate exceeds range bounds";
+ m_UseFullRange = true;
}
- if (len < 3) return false;
- //create a new edge array ...
- TEdge *edges = new TEdge [len];
- m_edges->push_back(edges);
-
- //convert vertices to a double-linked-list of edges and initialize ...
- edges[0].xcurr = p[0].X;
- edges[0].ycurr = p[0].Y;
- InitEdge(&edges[len-1], &edges[0], &edges[len-2], p[len-1], polyType);
- for (int i = len-2; i > 0; --i)
- InitEdge(&edges[i], &edges[i+1], &edges[i-1], p[i], polyType);
- InitEdge(&edges[0], &edges[1], &edges[len-1], p[0], polyType);
-
- //reset xcurr & ycurr and find 'eHighest' (given the Y axis coordinates
- //increase downward so the 'highest' edge will have the smallest ytop)
...
- TEdge *e = &edges[0];
- TEdge *eHighest = e;
- do
+ if (i == 0 || PointsEqual(p[j], pg[i])) continue;
+ else if (j > 0 && SlopesEqual(p[j-1], p[j], pg[i], m_UseFullRange))
{
- e->xcurr = e->xbot;
- e->ycurr = e->ybot;
- if (e->ytop < eHighest->ytop) eHighest = e;
- e = e->next;
- }
- while ( e != &edges[0]);
-
- //make sure eHighest is positioned so the following loop works safely
...
- if (eHighest->windDelta > 0) eHighest = eHighest->next;
- if (NEAR_EQUAL(eHighest->dx, HORIZONTAL)) eHighest = eHighest->next;
-
- //finally insert each local minima ...
- e = eHighest;
- do {
- e = AddBoundsToLML(e);
- }
- while( e != eHighest );
- return true;
+ if (PointsEqual(p[j-1], pg[i])) j--;
+ } else j++;
+ p[j] = pg[i];
}
-
//------------------------------------------------------------------------------
+ if (j < 2) return false;
- void ClipperBase::InsertLocalMinima(LocalMinima *newLm)
+ len = j+1;
+ for (;;)
{
- if( ! m_MinimaList )
+ //nb: test for point equality before testing slopes ...
+ if (PointsEqual(p[j], p[0])) j--;
+ else if (PointsEqual(p[0], p[1]) ||
+ SlopesEqual(p[j], p[0], p[1], m_UseFullRange))
+ p[0] = p[j--];
+ else if (SlopesEqual(p[j-1], p[j], p[0], m_UseFullRange)) j--;
+ else if (SlopesEqual(p[0], p[1], p[2], m_UseFullRange))
{
- m_MinimaList = newLm;
+ for (int i = 2; i <= j; ++i) p[i-1] = p[i];
+ j--;
}
- else if( newLm->Y >= m_MinimaList->Y )
- {
- newLm->next = m_MinimaList;
- m_MinimaList = newLm;
- } else
- {
- LocalMinima* tmpLm = m_MinimaList;
- while( tmpLm->next && ( newLm->Y < tmpLm->next->Y ) )
- tmpLm = tmpLm->next;
- newLm->next = tmpLm->next;
- tmpLm->next = newLm;
- }
+ //exit loop if nothing is changed or there are too few vertices ...
+ if (j == len-1 || j < 2) break;
+ len = j +1;
}
-
//------------------------------------------------------------------------------
+ if (len < 3) return false;
- TEdge* ClipperBase::AddBoundsToLML(TEdge *e)
+ //create a new edge array ...
+ TEdge *edges = new TEdge [len];
+ m_edges->push_back(edges);
+
+ //convert vertices to a double-linked-list of edges and initialize ...
+ edges[0].xcurr = p[0].X;
+ edges[0].ycurr = p[0].Y;
+ InitEdge(&edges[len-1], &edges[0], &edges[len-2], p[len-1], polyType);
+ for (int i = len-2; i > 0; --i)
+ InitEdge(&edges[i], &edges[i+1], &edges[i-1], p[i], polyType);
+ InitEdge(&edges[0], &edges[1], &edges[len-1], p[0], polyType);
+
+ //reset xcurr & ycurr and find 'eHighest' (given the Y axis coordinates
+ //increase downward so the 'highest' edge will have the smallest ytop) ...
+ TEdge *e = &edges[0];
+ TEdge *eHighest = e;
+ do
{
- //Starting at the top of one bound we progress to the bottom where
there's
- //a local minima. We then go to the top of the next bound. These two
bounds
- //form the left and right (or right and left) bounds of the local
minima.
- e->nextInLML = 0;
+ e->xcurr = e->xbot;
+ e->ycurr = e->ybot;
+ if (e->ytop < eHighest->ytop) eHighest = e;
e = e->next;
- for (;;)
- {
- if (NEAR_EQUAL(e->dx, HORIZONTAL))
- {
- //nb: proceed through horizontals when approaching from their
right,
- // but break on horizontal minima if approaching from their
left.
- // This ensures 'local minima' are always on the left of
horizontals.
- if (e->next->ytop < e->ytop && e->next->xbot > e->prev->xbot)
break;
- if (e->xtop != e->prev->xbot) SwapX(*e);
- e->nextInLML = e->prev;
- }
- else if (e->ycurr == e->prev->ycurr) break;
- else e->nextInLML = e->prev;
- e = e->next;
- }
+ }
+ while ( e != &edges[0]);
- //e and e.prev are now at a local minima ...
- LocalMinima* newLm = new LocalMinima;
- newLm->next = 0;
- newLm->Y = e->prev->ybot;
+ //make sure eHighest is positioned so the following loop works safely ...
+ if (eHighest->windDelta > 0) eHighest = eHighest->next;
+ if (NEAR_EQUAL(eHighest->dx, HORIZONTAL)) eHighest = eHighest->next;
- if ( NEAR_EQUAL(e->dx, HORIZONTAL) ) //horizontal edges never start a
left bound
- {
- if (e->xbot != e->prev->xbot) SwapX(*e);
- newLm->leftBound = e->prev;
- newLm->rightBound = e;
- } else if (e->dx < e->prev->dx)
- {
- newLm->leftBound = e->prev;
- newLm->rightBound = e;
- } else
- {
- newLm->leftBound = e;
- newLm->rightBound = e->prev;
- }
- newLm->leftBound->side = esLeft;
- newLm->rightBound->side = esRight;
- InsertLocalMinima( newLm );
-
- for (;;)
- {
- if ( e->next->ytop == e->ytop && !NEAR_EQUAL(e->next->dx,
HORIZONTAL) ) break;
- e->nextInLML = e->next;
- e = e->next;
- if ( NEAR_EQUAL(e->dx, HORIZONTAL) && e->xbot != e->prev->xtop)
SwapX(*e);
- }
- return e->next;
+ //finally insert each local minima ...
+ e = eHighest;
+ do {
+ e = AddBoundsToLML(e);
}
-
//------------------------------------------------------------------------------
+ while( e != eHighest );
+ return true;
+}
+//------------------------------------------------------------------------------
- bool ClipperBase::AddPolygons(const Polygons &ppg, PolyType polyType)
+void ClipperBase::InsertLocalMinima(LocalMinima *newLm)
+{
+ if( ! m_MinimaList )
{
- bool result = true;
- for (Polygons::size_type i = 0; i < ppg.size(); ++i)
- if (AddPolygon(ppg[i], polyType)) result = false;
- return result;
+ m_MinimaList = newLm;
}
-
//------------------------------------------------------------------------------
-
- void ClipperBase::Clear()
+ else if( newLm->Y >= m_MinimaList->Y )
{
- DisposeLocalMinimaList();
- for (EdgeList::size_type i = 0; i < m_edges->size(); ++i) delete []
m_edges->at(i);
- m_edges->clear();
- m_UseFullRange = false;
+ newLm->next = m_MinimaList;
+ m_MinimaList = newLm;
+ } else
+ {
+ LocalMinima* tmpLm = m_MinimaList;
+ while( tmpLm->next && ( newLm->Y < tmpLm->next->Y ) )
+ tmpLm = tmpLm->next;
+ newLm->next = tmpLm->next;
+ tmpLm->next = newLm;
}
-
//------------------------------------------------------------------------------
+}
+//------------------------------------------------------------------------------
- void ClipperBase::Reset()
+TEdge* ClipperBase::AddBoundsToLML(TEdge *e)
+{
+ //Starting at the top of one bound we progress to the bottom where there's
+ //a local minima. We then go to the top of the next bound. These two bounds
+ //form the left and right (or right and left) bounds of the local minima.
+ e->nextInLML = 0;
+ e = e->next;
+ for (;;)
{
- m_CurrentLM = m_MinimaList;
- if( !m_CurrentLM ) return; //ie nothing to process
-
- //reset all edges ...
- LocalMinima* lm = m_MinimaList;
- while( lm )
+ if (NEAR_EQUAL(e->dx, HORIZONTAL))
{
- TEdge* e = lm->leftBound;
- while( e )
- {
- e->xcurr = e->xbot;
- e->ycurr = e->ybot;
- e->side = esLeft;
- e->outIdx = -1;
- e = e->nextInLML;
- }
- e = lm->rightBound;
- while( e )
- {
- e->xcurr = e->xbot;
- e->ycurr = e->ybot;
- e->side = esRight;
- e->outIdx = -1;
- e = e->nextInLML;
- }
- lm = lm->next;
+ //nb: proceed through horizontals when approaching from their right,
+ // but break on horizontal minima if approaching from their left.
+ // This ensures 'local minima' are always on the left of
horizontals.
+ if (e->next->ytop < e->ytop && e->next->xbot > e->prev->xbot) break;
+ if (e->xtop != e->prev->xbot) SwapX(*e);
+ e->nextInLML = e->prev;
}
+ else if (e->ycurr == e->prev->ycurr) break;
+ else e->nextInLML = e->prev;
+ e = e->next;
}
-
//------------------------------------------------------------------------------
- void ClipperBase::DisposeLocalMinimaList()
+ //e and e.prev are now at a local minima ...
+ LocalMinima* newLm = new LocalMinima;
+ newLm->next = 0;
+ newLm->Y = e->prev->ybot;
+
+ if ( NEAR_EQUAL(e->dx, HORIZONTAL) ) //horizontal edges never start a left
bound
{
- while( m_MinimaList )
- {
- LocalMinima* tmpLm = m_MinimaList->next;
- delete m_MinimaList;
- m_MinimaList = tmpLm;
- }
- m_CurrentLM = 0;
+ if (e->xbot != e->prev->xbot) SwapX(*e);
+ newLm->leftBound = e->prev;
+ newLm->rightBound = e;
+ } else if (e->dx < e->prev->dx)
+ {
+ newLm->leftBound = e->prev;
+ newLm->rightBound = e;
+ } else
+ {
+ newLm->leftBound = e;
+ newLm->rightBound = e->prev;
}
-
//------------------------------------------------------------------------------
+ newLm->leftBound->side = esLeft;
+ newLm->rightBound->side = esRight;
+ InsertLocalMinima( newLm );
- void ClipperBase::PopLocalMinima()
+ for (;;)
{
- if( ! m_CurrentLM ) return;
- m_CurrentLM = m_CurrentLM->next;
+ if ( e->next->ytop == e->ytop && !NEAR_EQUAL(e->next->dx, HORIZONTAL) )
break;
+ e->nextInLML = e->next;
+ e = e->next;
+ if ( NEAR_EQUAL(e->dx, HORIZONTAL) && e->xbot != e->prev->xtop)
SwapX(*e);
}
-
//------------------------------------------------------------------------------
+ return e->next;
+}
+//------------------------------------------------------------------------------
- IntRect ClipperBase::GetBounds()
+bool ClipperBase::AddPolygons(const Polygons &ppg, PolyType polyType)
+{
+ bool result = true;
+ for (Polygons::size_type i = 0; i < ppg.size(); ++i)
+ if (AddPolygon(ppg[i], polyType)) result = false;
+ return result;
+}
+//------------------------------------------------------------------------------
+
+void ClipperBase::Clear()
+{
+ DisposeLocalMinimaList();
+ for (EdgeList::size_type i = 0; i < m_edges->size(); ++i) delete []
m_edges->at(i);
+ m_edges->clear();
+ m_UseFullRange = false;
+}
+//------------------------------------------------------------------------------
+
+void ClipperBase::Reset()
+{
+ m_CurrentLM = m_MinimaList;
+ if( !m_CurrentLM ) return; //ie nothing to process
+
+ //reset all edges ...
+ LocalMinima* lm = m_MinimaList;
+ while( lm )
{
- IntRect result;
- LocalMinima* lm = m_MinimaList;
- if (!lm)
+ TEdge* e = lm->leftBound;
+ while( e )
{
- result.left = result.top = result.right = result.bottom = 0;
- return result;
+ e->xcurr = e->xbot;
+ e->ycurr = e->ybot;
+ e->side = esLeft;
+ e->outIdx = -1;
+ e = e->nextInLML;
}
- result.left = lm->leftBound->xbot;
- result.top = lm->leftBound->ybot;
- result.right = lm->leftBound->xbot;
- result.bottom = lm->leftBound->ybot;
- while (lm)
+ e = lm->rightBound;
+ while( e )
{
- if (lm->leftBound->ybot > result.bottom)
- result.bottom = lm->leftBound->ybot;
- TEdge* e = lm->leftBound;
- for (;;) {
- TEdge* bottomE = e;
- while (e->nextInLML)
- {
- if (e->xbot < result.left) result.left = e->xbot;
- if (e->xbot > result.right) result.right = e->xbot;
- e = e->nextInLML;
- }
- if (e->xbot < result.left) result.left = e->xbot;
- if (e->xbot > result.right) result.right = e->xbot;
- if (e->xtop < result.left) result.left = e->xtop;
- if (e->xtop > result.right) result.right = e->xtop;
- if (e->ytop < result.top) result.top = e->ytop;
-
- if (bottomE == lm->leftBound) e = lm->rightBound;
- else break;
- }
- lm = lm->next;
+ e->xcurr = e->xbot;
+ e->ycurr = e->ybot;
+ e->side = esRight;
+ e->outIdx = -1;
+ e = e->nextInLML;
}
- return result;
+ lm = lm->next;
}
+}
+//------------------------------------------------------------------------------
-
-
//------------------------------------------------------------------------------
- // TClipper methods ...
-
//------------------------------------------------------------------------------
-
- Clipper::Clipper() : ClipperBase() //constructor
+void ClipperBase::DisposeLocalMinimaList()
+{
+ while( m_MinimaList )
{
- m_Scanbeam = 0;
- m_ActiveEdges = 0;
- m_SortedEdges = 0;
- m_IntersectNodes = 0;
- m_ExecuteLocked = false;
- m_UseFullRange = false;
- m_ReverseOutput = false;
- m_PolyOuts = new PolyOutList();
- m_Joins = new JoinList();
- m_HorizJoins = new HorzJoinList();
+ LocalMinima* tmpLm = m_MinimaList->next;
+ delete m_MinimaList;
+ m_MinimaList = tmpLm;
}
-
//------------------------------------------------------------------------------
+ m_CurrentLM = 0;
+}
+//------------------------------------------------------------------------------
- Clipper::~Clipper() //destructor
- {
- Clear();
- DisposeScanbeamList();
- delete m_PolyOuts;
- delete m_Joins;
- delete m_HorizJoins;
- }
-
//------------------------------------------------------------------------------
+void ClipperBase::PopLocalMinima()
+{
+ if( ! m_CurrentLM ) return;
+ m_CurrentLM = m_CurrentLM->next;
+}
+//------------------------------------------------------------------------------
- void Clipper::Clear()
+IntRect ClipperBase::GetBounds()
+{
+ IntRect result;
+ LocalMinima* lm = m_MinimaList;
+ if (!lm)
{
- if (m_edges->size() == 0) return; //avoids problems with ClipperBase
destructor
- DisposeAllPolyPts();
- ClipperBase::Clear();
+ result.left = result.top = result.right = result.bottom = 0;
+ return result;
}
-
//------------------------------------------------------------------------------
-
- void Clipper::DisposeScanbeamList()
+ result.left = lm->leftBound->xbot;
+ result.top = lm->leftBound->ybot;
+ result.right = lm->leftBound->xbot;
+ result.bottom = lm->leftBound->ybot;
+ while (lm)
{
- while ( m_Scanbeam ) {
- Scanbeam* sb2 = m_Scanbeam->next;
- delete m_Scanbeam;
- m_Scanbeam = sb2;
- }
- }
-
//------------------------------------------------------------------------------
+ if (lm->leftBound->ybot > result.bottom)
+ result.bottom = lm->leftBound->ybot;
+ TEdge* e = lm->leftBound;
+ for (;;) {
+ TEdge* bottomE = e;
+ while (e->nextInLML)
+ {
+ if (e->xbot < result.left) result.left = e->xbot;
+ if (e->xbot > result.right) result.right = e->xbot;
+ e = e->nextInLML;
+ }
+ if (e->xbot < result.left) result.left = e->xbot;
+ if (e->xbot > result.right) result.right = e->xbot;
+ if (e->xtop < result.left) result.left = e->xtop;
+ if (e->xtop > result.right) result.right = e->xtop;
+ if (e->ytop < result.top) result.top = e->ytop;
- void Clipper::Reset()
- {
- ClipperBase::Reset();
- m_Scanbeam = 0;
- m_ActiveEdges = 0;
- m_SortedEdges = 0;
- LocalMinima* lm = m_MinimaList;
- while (lm)
- {
- InsertScanbeam(lm->Y);
- InsertScanbeam(lm->leftBound->ytop);
- lm = lm->next;
+ if (bottomE == lm->leftBound) e = lm->rightBound;
+ else break;
}
+ lm = lm->next;
}
-
//------------------------------------------------------------------------------
+ return result;
+}
- bool Clipper::Execute(ClipType clipType, Polygons &solution,
- PolyFillType subjFillType, PolyFillType clipFillType)
- {
- if( m_ExecuteLocked ) return false;
- m_ExecuteLocked = true;
- solution.resize(0);
- m_SubjFillType = subjFillType;
- m_ClipFillType = clipFillType;
- m_ClipType = clipType;
- bool succeeded = ExecuteInternal(false);
- if (succeeded) BuildResult(solution);
- m_ExecuteLocked = false;
- return succeeded;
+
+//------------------------------------------------------------------------------
+// TClipper methods ...
+//------------------------------------------------------------------------------
+
+Clipper::Clipper() : ClipperBase() //constructor
+{
+ m_Scanbeam = 0;
+ m_ActiveEdges = 0;
+ m_SortedEdges = 0;
+ m_IntersectNodes = 0;
+ m_ExecuteLocked = false;
+ m_UseFullRange = false;
+ m_ReverseOutput = false;
+ m_PolyOuts = new PolyOutList();
+ m_Joins = new JoinList();
+ m_HorizJoins = new HorzJoinList();
+}
+//------------------------------------------------------------------------------
+
+Clipper::~Clipper() //destructor
+{
+ Clear();
+ DisposeScanbeamList();
+ delete m_PolyOuts;
+ delete m_Joins;
+ delete m_HorizJoins;
+}
+//------------------------------------------------------------------------------
+
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
