OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Templatized some geometry tools

Browse Source
This commit is contained in:
Jacob Støren 2013-12-05 09:45:27 +01:00
parent 286bae9e85
commit f92ca7fba7
4 changed files with 956 additions and 876 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
ApplicationCode/ModelVisualization/ModelVisualization_UnitTests/RivCellFaceIntersection-Test.cpp
View File

@ -93,9 +93,6 @@ NodeType quadNormal (const ArrayWrapperConst<NodeArrayType, NodeType>& nodeCoor
( nodeCoords[cubeFaceIndices[3]] - nodeCoords[cubeFaceIndices[1]]);
}
#if 1
#endif
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -147,6 +144,29 @@ private:
QuadFaceIntersectorImplHandle * m_implementation;
};
std::vector<cvf::Vec3d> createVertices()
{
std::vector<cvf::Vec3d> vxs;
vxs.resize(14, cvf::Vec3d::ZERO);
vxs[ 0]= cvf::Vec3d( 0 , 0 , 0 );
vxs[ 1]= cvf::Vec3d( 1 , 0 , 0 );
vxs[ 2]= cvf::Vec3d( 1 , 1 , 0 );
vxs[ 3]= cvf::Vec3d( 0 , 1 , 0 );
vxs[ 4]= cvf::Vec3d(-0.4 ,-0.2 , 0.0 );
vxs[ 5]= cvf::Vec3d( 0.4 , 0.6 , 0.0 );
vxs[ 6]= cvf::Vec3d( 0.8 , 0.2 , 0.0 );
vxs[ 7]= cvf::Vec3d( 0.0 ,-0.6 , 0.0 );
vxs[ 8]= cvf::Vec3d( 1.0 , 1.2 , 0.0 );
vxs[ 9]= cvf::Vec3d( 1.4 , 0.8 , 0.0 );
vxs[10]= cvf::Vec3d( 0.4 ,-0.2 , 0.0 );
vxs[11]= cvf::Vec3d( 1.2 , 0.6 , 0.0 );
vxs[12]= cvf::Vec3d( 1.6 , 0.2 , 0.0 );
vxs[13]= cvf::Vec3d( 0.8 ,-0.6 , 0.0 );
return vxs;
}
//--------------------------------------------------------------------------------------------------
///
@ -165,7 +185,7 @@ TEST(CellFaceIntersectionTst, Intersection)
nodes.resize(8);
nodes.setAll(cvf::Vec3d(0, 0, 0));
EdgeIntersectStorage edgeIntersectionStorage;
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount(nodes.size());
// Face 1
@ -180,7 +200,8 @@ TEST(CellFaceIntersectionTst, Intersection)
nodes[7] = cvf::Vec3d(0, 1, 0);
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage, nodes, cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage,
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
EXPECT_EQ( 4, polygon.size());
EXPECT_EQ( (size_t)0, additionalVertices.size());
EXPECT_TRUE(isOk);
@ -197,7 +218,8 @@ TEST(CellFaceIntersectionTst, Intersection)
nodes[7] = cvf::Vec3d(-0.25, 0.5, 0);
polygon.clear();
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage, nodes, cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage,
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
EXPECT_EQ( 8, polygon.size());
EXPECT_EQ( (size_t)8, additionalVertices.size());
EXPECT_TRUE(isOk);
@ -224,7 +246,7 @@ TEST(CellFaceIntersectionTst, FreeFacePolygon)
nodes.resize(8);
nodes.setAll(cvf::Vec3d(0, 0, 0));
EdgeIntersectStorage edgeIntersectionStorage;
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount(nodes.size());
// Face 1
@ -239,13 +261,27 @@ TEST(CellFaceIntersectionTst, FreeFacePolygon)
nodes[7] = cvf::Vec3d(0, 1, 0);
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage, nodes, cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage,
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
EXPECT_EQ( 4, polygon.size());
EXPECT_EQ( (size_t)0, additionalVertices.size());
EXPECT_TRUE(isOk);
std::vector< bool > faceOverlapPolygonWinding;
std::vector< std::vector<size_t>* > faceOverlapPolygons;
faceOverlapPolygons.push_back(&polygon);
faceOverlapPolygonWinding.push_back(true);
//GeometryTools::calculatePartiallyFreeCubeFacePolygon(nodes, );
std::vector<size_t> partialFacePolygon;
bool hasHoles = false;
GeometryTools::calculatePartiallyFreeCubeFacePolygon(
wrapArrayConst(&nodes),
wrapArrayConst(cv1CubeFaceIndices, 4),
Vec3d(0,0,1),
faceOverlapPolygons,
faceOverlapPolygonWinding,
&partialFacePolygon,
&hasHoles);
// Face 1
nodes[0] = cvf::Vec3d(0, 0, 0);
@ -259,7 +295,8 @@ TEST(CellFaceIntersectionTst, FreeFacePolygon)
nodes[7] = cvf::Vec3d(-0.25, 0.5, 0);
polygon.clear();
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage, nodes, cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, edgeIntersectionStorage,
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
EXPECT_EQ( 8, polygon.size());
EXPECT_EQ( (size_t)8, additionalVertices.size());
EXPECT_TRUE(isOk);

843
ApplicationCode/ModelVisualization/cvfGeometryTools.cpp
View File

@ -307,106 +307,6 @@ GeometryTools::inPlaneLineIntersect3D( const cvf::Vec3d& planeNormal,
}
//--------------------------------------------------------------------------------------------------
/// \brief Test if a point touches a polygon within the specified tolerance
///
/// \param polygonNorm Polygon normal
/// \param pPolygonVerts Array of polygon vertice coordinates
/// \param piVertexIndices Array of integer node indices for this polygon
/// \param iNumVerts Number of vertices in polygon
/// \param point The point to be checked
/// \param tolerance Tolerance in length
/// \param touchedEdgeIndex returns -1 if point is inside, and edge index if point touches an edge.
/// \return true if point lies inside or on the border of the polygon.
///
/// \assumpt Assumes that the polygon is planar
/// \comment First check if point is on an edge, Then check if it is inside by
/// counting the number of times a ray from point along positive X axis
/// crosses an edge. Odd number says inside.
/// \author SP (really by Eric Haines) and JJS
//--------------------------------------------------------------------------------------------------
bool GeometryTools::isPointTouchingIndexedPolygon(const cvf::Vec3d& polygonNormal, const cvf::Vec3d* vertices, const size_t* indices, size_t numIndices, const cvf::Vec3d& point, int* touchedEdgeIndex, double tolerance )
{
int Z = findClosestAxis(polygonNormal);
int X = (Z + 1) % 3;
int Y = (Z + 2) % 3;
int crossings;
int xBelowVx0;
int yBelowVx0;
int yBelowVx1 = 0;
const double* vtx0;
const double* vtx1 = NULL;
double dv0;
cvf::uint j;
// Check if point is on an edge or vertex
size_t firstIdx;
size_t secondIdx;
CVF_TIGHT_ASSERT(touchedEdgeIndex);
*touchedEdgeIndex = -1;
for (firstIdx = 0, secondIdx = 1; firstIdx < numIndices; ++firstIdx, ++secondIdx)
{
if (secondIdx >= numIndices) secondIdx = 0;
const cvf::Vec3d& vx0 = vertices[indices[firstIdx]];
const cvf::Vec3d& vx1 = vertices[indices[secondIdx]];
double sqDist = GeometryTools::linePointSquareDist(vx0, vx1, point);
if (sqDist < tolerance*tolerance)
{
*touchedEdgeIndex = static_cast<int>(firstIdx);
return true;
}
}
vtx0 = vertices[indices[numIndices-1]].ptr();
// get test bit for above/below Y axis. Y of Point is under Y of vtx0
yBelowVx0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0;
crossings = 0;
for (j = 0; j < numIndices; j++)
{
// cleverness: bobble between filling endpoints of edges, so that the previous edge's shared endpoint is maintained.
if (j & 0x1)
{
vtx0 = vertices[indices[j]].ptr();
yBelowVx0 = (dv0 = vtx0[Y] - point[Y]) >= 0.0;
}
else
{
vtx1 = vertices[indices[j]].ptr();
yBelowVx1 = (vtx1[Y] >= point[Y]);
}
// check if Y of point is between Y of Vx0 and Vx1
if (yBelowVx0 != yBelowVx1)
{
// check if X of point is not between X of Vx0 and Vx1
if ( (xBelowVx0 = (vtx0[X] >= point[X])) == (vtx1[X] >= point[X]) )
{
if (xBelowVx0) crossings++;
}
else
{
// compute intersection of polygon segment with X ray, note if > point's X.
crossings += (vtx0[X] - dv0*(vtx1[X] - vtx0[X])/(vtx1[Y] - vtx0[Y])) >= point[X];
}
}
}
// test if crossings is odd. If we care about its winding number > 0, then just: inside_flag = crossings > 0;
if (crossings & 0x01) return true;
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -631,7 +531,7 @@ bool GeometryTools::insertVertexInPolygon(std::list<std::pair<cvf::uint, bool> >
//--------------------------------------------------------------------------------------------------
void GeometryTools::addMidEdgeNodes(std::list<std::pair<cvf::uint, bool> >* polygon, const cvf::Vec3dArray& nodes, EdgeSplitStorage& edgeSplitStorage, std::vector<cvf::Vec3d>* createdVertexes)
{
cvf::uint newVertexIndex = nodes.size() + createdVertexes->size();
size_t newVertexIndex = nodes.size() + createdVertexes->size();
std::list<std::pair<cvf::uint, bool> >::iterator it;
std::list<std::pair<cvf::uint, bool> >::iterator it2;
@ -670,743 +570,6 @@ void GeometryTools::addMidEdgeNodes(std::list<std::pair<cvf::uint, bool> >* poly
}
//--------------------------------------------------------------------------------------------------
/// Returns true if we get an actual polygon
//--------------------------------------------------------------------------------------------------
bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector<size_t> * polygon, std::vector<cvf::Vec3d>* createdVertexes,
EdgeIntersectStorage& edgeIntersectionStorage,
const cvf::Vec3dArray& nodes,
const size_t cv1CubeFaceIndices[4],
const size_t cv2CubeFaceIndices[4],
double tolerance)
{
// Topology analysis
if (createdVertexes == NULL) return false;
size_t newVertexIndex = nodes.size() + createdVertexes->size();
bool cv1VxTouchCv2[4] = { false, false, false, false };
bool cv2VxTouchCv1[4] = { false, false, false, false };
int cv1VxTouchCv2Edge[4] = { -1, -1, -1, -1 };
int cv2VxTouchCv1Edge[4] = { -1, -1, -1, -1 };
int cv1Idx, cv2Idx;
int numMatchedNodes = 0;
// First check for complete topological match.
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
{
bool found = false;
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx)
{
if (cv1CubeFaceIndices[cv1Idx] == cv2CubeFaceIndices[cv2Idx])
{
cv1VxTouchCv2[cv1Idx] = true;
cv2VxTouchCv1[cv2Idx] = true;
found = true;
++numMatchedNodes;
continue;
}
}
}
if (numMatchedNodes >= 4) // Todo: Handle collapsed cells
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv1CubeFaceIndices[k]);
}
return true;
}
cvf::Vec3d diag1 = nodes[cv1CubeFaceIndices[2]] - nodes[cv1CubeFaceIndices[0]];
cvf::Vec3d diag2 = nodes[cv1CubeFaceIndices[3]] - nodes[cv1CubeFaceIndices[1]];
cvf::Vec3d normal = diag1^diag2;
int numCv1VxesOnCv2 = numMatchedNodes;
int numCv2VxesOnCv1 = numMatchedNodes;
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
{
if (!cv1VxTouchCv2[cv1Idx])
{
cv1VxTouchCv2[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal, nodes.ptr(), &cv2CubeFaceIndices[0], 4, nodes[cv1CubeFaceIndices[cv1Idx]], &(cv1VxTouchCv2Edge[cv1Idx]), tolerance);
if (cv1VxTouchCv2[cv1Idx]) ++numCv1VxesOnCv2;
}
if (!cv2VxTouchCv1[cv1Idx])
{
cv2VxTouchCv1[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal, nodes.ptr(), &cv1CubeFaceIndices[0], 4, nodes[cv2CubeFaceIndices[cv1Idx]], &(cv2VxTouchCv1Edge[cv1Idx]), tolerance);
if (cv2VxTouchCv1[cv1Idx]) ++numCv2VxesOnCv1;
}
}
// Handle case where one of the faces are completely covered by the other
if (numCv1VxesOnCv2 >= 4)
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv1CubeFaceIndices[k]);
}
return true;
}
if (numCv2VxesOnCv1 >= 4)
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv2CubeFaceIndices[k]);
}
return true;
}
// Handle partial coverage
// Algorithm outline as follows:
// Loop over edges in the face of Cv1. Intersect each one with all the edges of the Cv2 face.
// Add first point of the cv1 edge to polygon if it really touches Cv2 ( touch of edge is considered as not touching)
// Add each intersection point along the Cv1 edge if present
// and finally: if the cv1 edge is going out of cv2, the add the cv2 vertexes from that intersection as long as they touch cv1.
int nextCv1Idx = 1;
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx, ++nextCv1Idx)
{
if (nextCv1Idx > 3) nextCv1Idx = 0;
if (cv1VxTouchCv2[cv1Idx] && cv1VxTouchCv2Edge[cv1Idx] == -1) // Start of cv1 edge is touching inside the cv2 polygon (not on an cv2 edge)
{
if (polygon->empty() || polygon->back() != cv1CubeFaceIndices[cv1Idx])
{
polygon->push_back(cv1CubeFaceIndices[cv1Idx]);
}
if (cv1VxTouchCv2[nextCv1Idx] && cv1VxTouchCv2Edge[nextCv1Idx] == -1)
{
// Both ends of this cv1 edge is touching inside(not on an edge) cv2 polygon, no intersections possible (assuming convex cube-face)
// Continue with next Cv1-edge.
continue;
}
}
// Find intersection(s) on this edge
std::vector<size_t> intersectionVxIndices;
std::vector<int> intersectedCv2EdgeIdxs;
std::vector<double> intersectionFractionsAlongEdge;
int nextCv2Idx = 1;
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx, ++nextCv2Idx)
{
if (nextCv2Idx > 3) nextCv2Idx = 0;
// Find a possible real intersection point.
cvf::Vec3d intersection(0,0,0);
double fractionAlongEdge1;
GeometryTools::IntersectionStatus intersectStatus = GeometryTools::NO_INTERSECTION;
size_t intersectionVxIndex = cvf::UNDEFINED_UINT;
// First handle some "trivial" ones to ease the burden for the real intersection calculation
// It could be tested whether it really is necessary to do
if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx && cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
{
intersectStatus = GeometryTools::LINES_OVERLAP;
fractionAlongEdge1 = 1;
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
}
else if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx )
{
// When this happens, the cv1 vertex will already have been added to the polygon
// by the statements in the top of the cv1 edge loop. Should it be treated specially ?
intersectStatus = GeometryTools::LINES_TOUCH;
fractionAlongEdge1 = 0;
intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
}
else if (cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
{
intersectStatus = GeometryTools::LINES_TOUCH;
fractionAlongEdge1 = 1;
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
}
else
{
double fractionAlongEdge2;
bool found = edgeIntersectionStorage.findIntersection( cv1CubeFaceIndices[cv1Idx],
cv1CubeFaceIndices[nextCv1Idx],
cv2CubeFaceIndices[cv2Idx],
cv2CubeFaceIndices[nextCv2Idx],
&intersectionVxIndex, &intersectStatus,
&fractionAlongEdge1, &fractionAlongEdge2);
if (!found)
{
intersectStatus = GeometryTools::inPlaneLineIntersect3D(normal,
nodes[cv1CubeFaceIndices[cv1Idx]],
nodes[cv1CubeFaceIndices[nextCv1Idx]],
nodes[cv2CubeFaceIndices[cv2Idx]],
nodes[cv2CubeFaceIndices[nextCv2Idx]],
&intersection, &fractionAlongEdge1, &fractionAlongEdge2,
tolerance);
switch (intersectStatus)
{
case GeometryTools::LINES_CROSSES:
{
intersectionVxIndex = newVertexIndex;
createdVertexes->push_back(intersection);
++newVertexIndex;
}
break;
case GeometryTools::LINES_TOUCH:
{
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
else CVF_ASSERT(false); // Tolerance trouble
}
break;
case GeometryTools::LINES_OVERLAP:
{
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
else CVF_ASSERT(false); // Tolerance trouble
}
break;
default:
break;
}
edgeIntersectionStorage.addIntersection( cv1CubeFaceIndices[cv1Idx],
cv1CubeFaceIndices[nextCv1Idx],
cv2CubeFaceIndices[cv2Idx],
cv2CubeFaceIndices[nextCv2Idx],
intersectionVxIndex, intersectStatus,
fractionAlongEdge1, fractionAlongEdge2);
}
}
// Store data for each intersection along the Cv1-edge
if ( (intersectStatus == GeometryTools::LINES_CROSSES)
|| (intersectStatus == GeometryTools::LINES_TOUCH)
|| (intersectStatus == GeometryTools::LINES_OVERLAP) )
{
CVF_ASSERT(intersectionVxIndex != cvf::UNDEFINED_UINT);
intersectionFractionsAlongEdge.push_back(fractionAlongEdge1);
intersectedCv2EdgeIdxs.push_back(cv2Idx);
intersectionVxIndices.push_back(intersectionVxIndex);
}
}
// Insert the intersections into the polygon in the correct order along the Cv1 edge.
// Find the last intersection in order to possibly continue the polygon along Cv2 into Cv1
size_t i;
size_t lastIntersection = std::numeric_limits<size_t>::max();
double largestFraction = -1;
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
{
if (intersectionFractionsAlongEdge[i] > largestFraction)
{
lastIntersection = i;
largestFraction = intersectionFractionsAlongEdge[i];
}
}
// Insert indices to the new intersection vertices into the polygon of
// this connection according to fraction along edge
std::map<double,size_t> sortingMap;
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
{
sortingMap[intersectionFractionsAlongEdge[i]] = intersectionVxIndices[i];
}
std::map<double, size_t>::iterator it;
for (it = sortingMap.begin(); it != sortingMap.end(); ++it)
{
if (polygon->empty() || polygon->back() != it->second)
{
polygon->push_back(it->second);
}
}
// Then if the Cv1 edge is going out of Cv2, add to the polygon, all the Cv2 face vertex-indices
// from the intersection that touches Cv1.
// if cv1 edge in any way touches cv2 and ends up outside, it went out.
/*
if cv1 edge is going out of cv2 then
if intersected cv2 edge has endpoint touching cv1 then
add endpoint to polygon. continue to add next endpoint until it does not touch Cv1
*/
if ( !cv1VxTouchCv2[nextCv1Idx]
&& ( cv1VxTouchCv2[cv1Idx] || ( intersectedCv2EdgeIdxs.size() ) ) ) // Two touches along edge also qualifies
{
if(lastIntersection < intersectedCv2EdgeIdxs.size())
{
cv2Idx = intersectedCv2EdgeIdxs[lastIntersection];
int count = 0;
// Continue the polygon along the Cv2 edges as long as they touch cv1.
// Depending on the faces having opposite winding, which is guaranteed as long as
// no intersecting CVs share a connection
while (cv2VxTouchCv1[cv2Idx] && count < 4 && (cv2VxTouchCv1Edge[cv2Idx] == -1)) // Touch of edge is regarded as being outside, so we must stop
{
if (polygon->empty() || polygon->back() != cv2CubeFaceIndices[cv2Idx])
{
polygon->push_back(cv2CubeFaceIndices[cv2Idx]);
}
--cv2Idx;
if (cv2Idx < 0 ) cv2Idx = 3;
++count;
}
}
else
{
CVF_ASSERT(lastIntersection < intersectedCv2EdgeIdxs.size());
}
}
}
if (polygon->size() > 2)
{
if (polygon->back() == polygon->front()) polygon->pop_back();
}
// Sanity checks
if (polygon->size() < 3)
{
// cvf::Trace::show(cvf::String("Degenerated connection polygon detected. (Less than 3 vertexes) Cv's probably not in contact: %1 , %2").arg(m_ownerCvId).arg(m_neighborCvId));
polygon->clear();
return false;
}
return true;
}
//--------------------------------------------------------------------------------------------------
/// This method assumes that all intersection and mid edge vertexes are created an are already
/// merged into all the polygons. We can also assume that all the connection polygons are completely
/// inside (or sharing edges with) the cube face polygon initially
//--------------------------------------------------------------------------------------------------
// Vertex Index to position in polygon
typedef std::map< size_t, std::vector<size_t >::const_iterator > VxIdxToPolygonPositionMap;
#define DEBUG_PRINT 0
//template <typename NodeArrayType, typename NodeType, typename IndicesArrayType, typename IndicesType>
//void setup( ArrayWrapper<NodeArrayType, NodeType> nodeArray, ArrayWrapper<IndicesArrayType, IndicesType> indices)
void GeometryTools::calculatePartiallyFreeCubeFacePolygon(const cvf::Vec3dArray& nodeCoords,
const std::vector<size_t>* completeFacePolygon,
const cvf::Vec3d& faceNormal,
const std::vector< std::vector<size_t>* >& faceOverlapPolygons,
const std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags,
std::vector<size_t>* partialFacePolygon,
bool* m_partiallyFreeCubeFaceHasHoles)
{
CVF_ASSERT(m_partiallyFreeCubeFaceHasHoles);
CVF_ASSERT(completeFacePolygon != NULL);
CVF_ASSERT(partialFacePolygon != NULL);
// Copy the start polygon
std::list<size_t> resultPolygon;
for (size_t pcIdx = 0; pcIdx < completeFacePolygon->size(); ++pcIdx)
{
resultPolygon.push_back((*completeFacePolygon)[pcIdx]);
}
// First build search maps to fast find whether and where an index is positioned in a polygon
// Map from Vertex-index to position in polygon
std::vector< VxIdxToPolygonPositionMap > polygonSearchMaps;
std::vector<bool> isConnectionPolygonMerged;
polygonSearchMaps.resize(faceOverlapPolygons.size());
isConnectionPolygonMerged.resize(faceOverlapPolygons.size(), false);
// Build search maps
{
size_t count;
for (size_t i = 0; i < faceOverlapPolygons.size(); ++i)
{
count = 0;
for (std::vector<size_t >::const_iterator pcIt = faceOverlapPolygons[i]->begin();
pcIt != faceOverlapPolygons[i]->end();
++pcIt)
{
polygonSearchMaps[i][*pcIt] = pcIt;
++count;
}
if (count < 3) isConnectionPolygonMerged[i] = true; // Ignore false polygons
}
}
#if DEBUG_PRINT
{
cvf::Trace::show("Circumference polygon: ");
std::list<size_t>::const_iterator polIt;
for ( polIt = resultPolygon.begin(); polIt != resultPolygon.end(); ++polIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg(*polIt)
.arg(nodeCoords[*polIt].x())
.arg(nodeCoords[*polIt].y())
.arg(nodeCoords[*polIt].z()));
}
}
#endif
#if DEBUG_PRINT
{
cvf::Trace::show("Connection polygons: ");
for (size_t cIdx = 0; cIdx < faceOverlapPolygons.size(); cIdx++)
{
std::vector<size_t >::const_iterator polIt;
cvf::Trace::show("Connection " + cvf::String((long long)cIdx));
for (polIt = faceOverlapPolygons[cIdx]->begin(); polIt != faceOverlapPolygons[cIdx]->end(); ++polIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg(*polIt)
.arg(nodeCoords[*polIt].x())
.arg(nodeCoords[*polIt].y())
.arg(nodeCoords[*polIt].z()));
}
}
}
#endif
// Merge connection polygons with the main polygon as long as one of them has something in common.
// For each vx in the m_freeFacePolygons[cubeFace] polygon .
// loop over all connections
// if it has the node in common and that the edge angle will decrease if inserting
// merge the connection polygon into the main polygon,
// and remove the connection polygon from the merge able connection polygons.
for (std::list<size_t>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
// Set iterator to previous node in polygon
std::list<size_t>::iterator prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
cvf::Vec3d pToPrev = nodeCoords[*prevPIt] - nodeCoords[*pIt];
// Set iterator to next node in polygon. Used to insert before and as pointer to the next point
std::list<size_t>::iterator nextPIt = pIt;
++nextPIt;
std::list<size_t>::iterator insertBeforePIt = nextPIt;
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
// Calculate existing edge to edge angle
cvf::Vec3d pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
double mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
// Find connections containing the pIt vertex index. Merge them into the main polygon
for (size_t opIdx = 0; opIdx < faceOverlapPolygons.size(); ++opIdx)
{
if (isConnectionPolygonMerged[opIdx]) continue; // Already merged
// Find position of pIt vertex index in the current connection polygon
VxIdxToPolygonPositionMap::iterator vxIndexPositionInPolygonIt = polygonSearchMaps[opIdx].find(*pIt);
if (vxIndexPositionInPolygonIt != polygonSearchMaps[opIdx].end())
{
// Merge the connection polygon into the main polygon
// if the angle prevPIt pIt nextPIt is larger than angle prevPIt pIt (startCPIt++)
std::vector<size_t>::const_iterator startCPIt;
startCPIt = vxIndexPositionInPolygonIt->second;
// First vx to insert is the one after the match
bool hasSameWinding = faceOverlapPolygonWindingSameAsCubeFaceFlags[opIdx];
if (hasSameWinding)
{
// Same winding as main polygon. We need to go the opposite way
if (startCPIt == faceOverlapPolygons[opIdx]->begin()) startCPIt = faceOverlapPolygons[opIdx]->end();
--startCPIt;
}
else
{
// Opposite winding. Go forward when merging
++startCPIt; if (startCPIt == faceOverlapPolygons[opIdx]->end()) startCPIt = faceOverlapPolygons[opIdx]->begin();
}
// Calculate possible new edge-to-edge angle and test against existing angle
cvf::Vec3d pToStart = nodeCoords[*startCPIt] - nodeCoords[*pIt];
double candidatePolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToStart , pToPrev);
if (candidatePolygonEdgeAngle < mainPolygonEdgeAngle )
{
// Merge ok
std::vector<size_t >::const_iterator pcIt = startCPIt;
if (hasSameWinding)
{
do
{
resultPolygon.insert(insertBeforePIt, (*pcIt));
if (pcIt == faceOverlapPolygons[opIdx]->begin()) pcIt = faceOverlapPolygons[opIdx]->end();
--pcIt;
} while (pcIt != startCPIt);
}
else
{
do
{
resultPolygon.insert(insertBeforePIt, (*pcIt));
++pcIt;
if (pcIt == faceOverlapPolygons[opIdx]->end()) pcIt = faceOverlapPolygons[opIdx]->begin();
} while (pcIt != startCPIt);
}
isConnectionPolygonMerged[opIdx] = true;
// Recalculate the next position to point into the new nodes
// Set iterator in the main polygon to insert before and to the next point
nextPIt = pIt;
++nextPIt;
insertBeforePIt = nextPIt;
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
// Recalculate the existing edge to edge angle
pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
}
}
}
}
// Now remove all double edges
bool goneAround = false;
for ( std::list<size_t>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end() && !goneAround; ++pIt)
{
// Set iterator to next node in polygon.
std::list<size_t>::iterator nextPIt = pIt;
++nextPIt;
if (nextPIt == resultPolygon.end())
{
nextPIt = resultPolygon.begin();
goneAround = true; // Gone around polygon. Stop even if pIt is jumping over end()
}
// Set iterator to previous node in polygon
std::list<size_t>::iterator prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
// If previous and next node are the same, erase
while(*nextPIt == *prevPIt)
{
resultPolygon.erase(pIt);
resultPolygon.erase(prevPIt);
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of while
pIt = nextPIt;
++nextPIt;
if (nextPIt == resultPolygon.end())
{
nextPIt = resultPolygon.begin();
goneAround = true; // Gone around polygon pIt is jumping over end()
}
prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
}
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of for loop
}
// Check for holes
bool hasHoles = false;
for (size_t i = 0; i < isConnectionPolygonMerged.size(); ++i)
{
hasHoles = !isConnectionPolygonMerged[i];
if(hasHoles)
{
*m_partiallyFreeCubeFaceHasHoles = true;
break;
}
}
#if DEBUG_PRINT
{
cvf::Trace::show("Polygon: ");
for (std::list<size_t>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg(*pIt)
.arg(nodeCoords[*pIt].x())
.arg(nodeCoords[*pIt].y())
.arg(nodeCoords[*pIt].z()));
}
}
#endif
// Copy the result polygon to the output variable
partialFacePolygon->clear();
for (std::list<size_t>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
partialFacePolygon->push_back(*pIt);
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void EdgeIntersectStorage::setVertexCount(size_t size)
{
m_edgeIntsectMap.resize(size);
}
void EdgeIntersectStorage::canonizeAddress(size_t& e1P1, size_t& e1P2, size_t& e2P1, size_t& e2P2, bool& flipE1, bool& flipE2, bool& flipE1E2)
{
flipE1 = e1P1 > e1P2;
flipE2 = e2P1 > e2P2;
flipE1E2 = (flipE1 ? e1P2: e1P1) > (flipE2 ? e2P2: e2P1);
static size_t temp;
if (flipE1)
{
temp = e1P1;
e1P1 = e1P2;
e1P2 = temp;
}
if (flipE2)
{
temp = e2P1;
e2P1 = e2P2;
e2P2 = temp;
}
if (flipE1E2)
{
temp = e1P1;
e1P1 = e2P1;
e2P1 = temp;
temp = e1P2;
e1P2 = e2P2;
e2P2 = temp;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void EdgeIntersectStorage::addIntersection(size_t e1P1, size_t e1P2, size_t e2P1, size_t e2P2,
size_t vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
double fractionAlongEdge1, double fractionAlongEdge2)
{
static bool flipE1 ;
static bool flipE2 ;
static bool flipE1E2;
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
static IntersectData iData;
iData.fractionAlongEdge1 = flipE1 ? 1 - fractionAlongEdge1 : fractionAlongEdge1;
iData.fractionAlongEdge2 = flipE2 ? 1 - fractionAlongEdge2 : fractionAlongEdge2;
iData.intersectionStatus = intersectionStatus;
if (flipE1E2)
{
double temp = iData.fractionAlongEdge1;
iData.fractionAlongEdge1 = iData.fractionAlongEdge2;
iData.fractionAlongEdge2 = temp;
}
iData.intersectionPointIndex = vxIndexIntersectionPoint;
CVF_ASSERT(e1P1 < m_edgeIntsectMap.size());
m_edgeIntsectMap[e1P1][e1P2][e2P1][e2P2] = iData;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool EdgeIntersectStorage::findIntersection(size_t e1P1, size_t e1P2, size_t e2P1, size_t e2P2,
size_t* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
double* fractionAlongEdge1, double* fractionAlongEdge2)
{
static bool flipE1 ;
static bool flipE2 ;
static bool flipE1E2;
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
if (!m_edgeIntsectMap[e1P1].size()) return false;
std::map<size_t, std::map<size_t, std::map<size_t, IntersectData > > >::iterator it;
it = m_edgeIntsectMap[e1P1].find(e1P2);
if (it == m_edgeIntsectMap[e1P1].end()) return false;
std::map<size_t, std::map<size_t, IntersectData > >::iterator it2;
it2 = it->second.find(e2P1);
if (it2 == it->second.end()) return false;
std::map<size_t, IntersectData >::iterator it3;
it3 = it2->second.find(e2P2);
if (it3 == it2->second.end()) return false;
*vxIndexIntersectionPoint = it3->second.intersectionPointIndex;
*intersectionStatus = it3->second.intersectionStatus;
if (flipE1E2)
{
*fractionAlongEdge1 = it3->second.fractionAlongEdge2;
*fractionAlongEdge2 = it3->second.fractionAlongEdge1;
}
else
{
*fractionAlongEdge1 = it3->second.fractionAlongEdge1;
*fractionAlongEdge2 = it3->second.fractionAlongEdge2;
}
if (flipE1) *fractionAlongEdge1 = 1 - *fractionAlongEdge1;
if (flipE2) *fractionAlongEdge2 = 1 - *fractionAlongEdge2;
return true;
}
//--------------------------------------------------------------------------------------------------
@ -1477,7 +640,7 @@ EarClipTesselator::EarClipTesselator():
/// \return true when a tesselation was successully created
//--------------------------------------------------------------------------------------------------
bool EarClipTesselator::calculateTriangles( std::vector<cvf::uint>* triangleIndices )
bool EarClipTesselator::calculateTriangles( std::vector<size_t>* triangleIndices )
{
CVF_ASSERT(m_nodeCoords != NULL);
CVF_ASSERT(m_X > -1 && m_Y > -1);
@ -1711,7 +874,7 @@ void FanEarClipTesselator::setCenterNode(size_t centerNodeIndex)
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool FanEarClipTesselator::calculateTriangles(std::vector<cvf::uint>* triangles)
bool FanEarClipTesselator::calculateTriangles(std::vector<size_t>* triangles)
{
CVF_ASSERT(m_centerNodeIndex != std::numeric_limits<size_t>::max());
CVF_ASSERT(m_nodeCoords != NULL);

60
ApplicationCode/ModelVisualization/cvfGeometryTools.h
View File

@ -4,13 +4,14 @@
#include <list>
#include <map>
#include <hash_map>
#include "cvfArrayWrapperConst.h"
namespace cvf
{
class EdgeSplitStorage;
class EdgeIntersectStorage;
template <typename IndexType> class EdgeIntersectStorage;
class GeometryTools
{
@ -44,56 +45,61 @@ public:
cvf::Vec3d* intersectionPoint, double* fractionAlongLine1, double* fractionAlongLine2,
double tolerance = 1e-6);
//template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
static bool isPointTouchingIndexedPolygon(const cvf::Vec3d& polygonNormal, const cvf::Vec3d* vertices, const size_t* indices, size_t numIndices,
const cvf::Vec3d& point, int* touchedEdgeIndex, double tolerance = 1e-6);
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
static bool isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNormal,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> vertices,
ArrayWrapperConst<PolygonArrayType, IndexType> indices,
const cvf::Vec3d& point,
int* touchedEdgeIndex,
double tolerance = 1e-6);
static bool calculateOverlapPolygonOfTwoQuads( std::vector<size_t> * polygon, std::vector<cvf::Vec3d>* createdVertexes,
EdgeIntersectStorage& edgeIntersectionStorage,
const cvf::Vec3dArray& nodes,
const size_t cv1CubeFaceIndices[4],
const size_t cv2CubeFaceIndices[4],
template<typename VerticeArrayType, typename IndexType>
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType> * polygon, std::vector<cvf::Vec3d>* createdVertexes,
EdgeIntersectStorage<IndexType>& edgeIntersectionStorage,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodes,
const IndexType cv1CubeFaceIndices[4],
const IndexType cv2CubeFaceIndices[4],
double tolerance);
static void calculatePartiallyFreeCubeFacePolygon(const cvf::Vec3dArray& nodeCoords,
const std::vector<size_t>* completeFacePolygon,
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
static void calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
ArrayWrapperConst<PolygonArrayType, IndexType> completeFacePolygon,
const cvf::Vec3d& faceNormal,
const std::vector< std::vector<size_t>* >& faceOverlapPolygons,
const std::vector< std::vector<IndexType>* >& faceOverlapPolygons,
const std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags,
std::vector<size_t>* partialFacePolygon,
std::vector<IndexType>* partialFacePolygon,
bool* m_partiallyFreeCubeFaceHasHoles);
};
template <typename IndexType>
class EdgeIntersectStorage
{
public:
void setVertexCount(size_t size);
bool findIntersection( size_t e1P1, size_t e1P2, size_t e2P1, size_t e2P2,
size_t* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
void setVertexCount(IndexType size);
bool findIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
double* fractionAlongEdge1, double* fractionAlongEdge2);
void addIntersection( size_t e1P1, size_t e1P2, size_t e2P1, size_t e2P2,
size_t vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
void addIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
IndexType vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
double fractionAlongEdge1, double fractionAlongEdge2);
private:
struct IntersectData
{
size_t intersectionPointIndex;
IndexType intersectionPointIndex;
GeometryTools::IntersectionStatus intersectionStatus;
double fractionAlongEdge1;
double fractionAlongEdge2;
};
void canonizeAddress(size_t& e1P1, size_t& e1P2, size_t& e2P1, size_t& e2P2, bool& flipE1, bool& flipE2, bool& flipE1E2);
void canonizeAddress(IndexType& e1P1, IndexType& e1P2, IndexType& e2P1, IndexType& e2P2, bool& flipE1, bool& flipE2, bool& flipE1E2);
// A map containing the intersection data. The addressing is :
// ( when leastVxIdxEdge1 < leastVxIdxEdge2 )
// leastVxIdxEdge1, largestVxIdxEdge1, leastVxIdxEdge2, largestVxIdxEdge2, { vxIdxIntersection, fractionAlongEdg1, fractionAlonEdge2 }
std::vector< std::map<size_t, std::map<size_t, std::map<size_t, IntersectData > > > > m_edgeIntsectMap;
std::vector< std::map<IndexType, std::map<IndexType, std::map<IndexType, IntersectData > > > > m_edgeIntsectMap;
};
class EdgeSplitStorage
@ -122,7 +128,7 @@ public:
void setPolygonIndices(const std::list<size_t>& polygon);
void setPolygonIndices(const std::vector<size_t>& polygon);
virtual bool calculateTriangles(std::vector<cvf::uint>* triangles);
virtual bool calculateTriangles(std::vector<size_t>* triangles);
protected:
bool isTriangleValid( std::list<size_t>::const_iterator u, std::list<size_t>::const_iterator v, std::list<size_t>::const_iterator w) const;
@ -145,10 +151,12 @@ public:
FanEarClipTesselator();
void setCenterNode(size_t centerNodeIndex );
virtual bool calculateTriangles(std::vector<cvf::uint>* triangles);
virtual bool calculateTriangles(std::vector<size_t>* triangles);
private:
bool isTriangleValid( size_t u, size_t v, size_t w);
size_t m_centerNodeIndex;
};
}
}
#include "cvfGeometryTools.inl"

872
ApplicationCode/ModelVisualization/cvfGeometryTools.inl Normal file
View File

@ -0,0 +1,872 @@
#pragma warning (disable : 4503)
namespace cvf
{
//--------------------------------------------------------------------------------------------------
/// \brief Test if a point touches a polygon within the specified tolerance
///
/// \param polygonNorm Polygon normal
/// \param pPolygonVerts Array of polygon vertice coordinates
/// \param piVertexIndices Array of integer node indices for this polygon
/// \param iNumVerts Number of vertices in polygon
/// \param point The point to be checked
/// \param tolerance Tolerance in length
/// \param touchedEdgeIndex returns -1 if point is inside, and edge index if point touches an edge.
/// \return true if point lies inside or on the border of the polygon.
///
/// \assumpt Assumes that the polygon is planar
/// \comment First check if point is on an edge, Then check if it is inside by
/// counting the number of times a ray from point along positive X axis
/// crosses an edge. Odd number says inside.
/// \author SP (really by Eric Haines) and JJS
//--------------------------------------------------------------------------------------------------
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNormal,
cvf::ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> vertices,
cvf::ArrayWrapperConst<PolygonArrayType, IndexType> indices,
const cvf::Vec3d& point,
int* touchedEdgeIndex,
double tolerance)
{
size_t numIndices = indices.size();
int Z = findClosestAxis(polygonNormal);
int X = (Z + 1) % 3;
int Y = (Z + 2) % 3;
int crossings;
int xBelowVx0;
int yBelowVx0;
int yBelowVx1 = 0;
const double* vtx0;
const double* vtx1 = NULL;
double dv0;
cvf::uint j;
// Check if point is on an edge or vertex
size_t firstIdx;
size_t secondIdx;
CVF_TIGHT_ASSERT(touchedEdgeIndex);
*touchedEdgeIndex = -1;
for (firstIdx = 0, secondIdx = 1; firstIdx < numIndices; ++firstIdx, ++secondIdx)
{
if (secondIdx >= numIndices) secondIdx = 0;
const cvf::Vec3d& vx0 = vertices[indices[firstIdx]];
const cvf::Vec3d& vx1 = vertices[indices[secondIdx]];
double sqDist = GeometryTools::linePointSquareDist(vx0, vx1, point);
if (sqDist < tolerance*tolerance)
{
*touchedEdgeIndex = static_cast<int>(firstIdx);
return true;
}
}
vtx0 = vertices[indices[numIndices-1]].ptr();
// get test bit for above/below Y axis. Y of Point is under Y of vtx0
yBelowVx0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0;
crossings = 0;
for (j = 0; j < numIndices; j++)
{
// cleverness: bobble between filling endpoints of edges, so that the previous edge's shared endpoint is maintained.
if (j & 0x1)
{
vtx0 = vertices[indices[j]].ptr();
yBelowVx0 = (dv0 = vtx0[Y] - point[Y]) >= 0.0;
}
else
{
vtx1 = vertices[indices[j]].ptr();
yBelowVx1 = (vtx1[Y] >= point[Y]);
}
// check if Y of point is between Y of Vx0 and Vx1
if (yBelowVx0 != yBelowVx1)
{
// check if X of point is not between X of Vx0 and Vx1
if ( (xBelowVx0 = (vtx0[X] >= point[X])) == (vtx1[X] >= point[X]) )
{
if (xBelowVx0) crossings++;
}
else
{
// compute intersection of polygon segment with X ray, note if > point's X.
crossings += (vtx0[X] - dv0*(vtx1[X] - vtx0[X])/(vtx1[Y] - vtx0[Y])) >= point[X];
}
}
}
// test if crossings is odd. If we care about its winding number > 0, then just: inside_flag = crossings > 0;
if (crossings & 0x01) return true;
return false;
}
//--------------------------------------------------------------------------------------------------
/// Returns true if we get an actual polygon
//--------------------------------------------------------------------------------------------------
template<typename VerticeArrayType, typename IndexType>
bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector<IndexType> * polygon,
std::vector<cvf::Vec3d>* createdVertexes,
EdgeIntersectStorage<IndexType>& edgeIntersectionStorage,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodes,
const IndexType cv1CubeFaceIndices[4],
const IndexType cv2CubeFaceIndices[4],
double tolerance)
{
// Topology analysis
if (createdVertexes == NULL) return false;
size_t newVertexIndex = nodes.size() + createdVertexes->size();
bool cv1VxTouchCv2[4] = { false, false, false, false };
bool cv2VxTouchCv1[4] = { false, false, false, false };
int cv1VxTouchCv2Edge[4] = { -1, -1, -1, -1 };
int cv2VxTouchCv1Edge[4] = { -1, -1, -1, -1 };
int cv1Idx, cv2Idx;
int numMatchedNodes = 0;
// First check for complete topological match.
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
{
bool found = false;
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx)
{
if (cv1CubeFaceIndices[cv1Idx] == cv2CubeFaceIndices[cv2Idx])
{
cv1VxTouchCv2[cv1Idx] = true;
cv2VxTouchCv1[cv2Idx] = true;
found = true;
++numMatchedNodes;
continue;
}
}
}
if (numMatchedNodes >= 4) // Todo: Handle collapsed cells
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv1CubeFaceIndices[k]);
}
return true;
}
cvf::Vec3d diag1 = nodes[cv1CubeFaceIndices[2]] - nodes[cv1CubeFaceIndices[0]];
cvf::Vec3d diag2 = nodes[cv1CubeFaceIndices[3]] - nodes[cv1CubeFaceIndices[1]];
cvf::Vec3d normal = diag1^diag2;
int numCv1VxesOnCv2 = numMatchedNodes;
int numCv2VxesOnCv1 = numMatchedNodes;
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
{
if (!cv1VxTouchCv2[cv1Idx])
{
cv1VxTouchCv2[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal,
nodes,
wrapArrayConst(cv2CubeFaceIndices, 4),
nodes[cv1CubeFaceIndices[cv1Idx]],
&(cv1VxTouchCv2Edge[cv1Idx]),
tolerance);
if (cv1VxTouchCv2[cv1Idx]) ++numCv1VxesOnCv2;
}
if (!cv2VxTouchCv1[cv1Idx])
{
cv2VxTouchCv1[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal,
nodes,
wrapArrayConst(cv1CubeFaceIndices, 4),
nodes[cv2CubeFaceIndices[cv1Idx]],
&(cv2VxTouchCv1Edge[cv1Idx]),
tolerance);
if (cv2VxTouchCv1[cv1Idx]) ++numCv2VxesOnCv1;
}
}
// Handle case where one of the faces are completely covered by the other
if (numCv1VxesOnCv2 >= 4)
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv1CubeFaceIndices[k]);
}
return true;
}
if (numCv2VxesOnCv1 >= 4)
{
int k;
for (k = 0; k < 4; ++k)
{
polygon->push_back(cv2CubeFaceIndices[k]);
}
return true;
}
// Handle partial coverage
// Algorithm outline as follows:
// Loop over edges in the face of Cv1. Intersect each one with all the edges of the Cv2 face.
// Add first point of the cv1 edge to polygon if it really touches Cv2 ( touch of edge is considered as not touching)
// Add each intersection point along the Cv1 edge if present
// and finally: if the cv1 edge is going out of cv2, the add the cv2 vertexes from that intersection as long as they touch cv1.
int nextCv1Idx = 1;
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx, ++nextCv1Idx)
{
if (nextCv1Idx > 3) nextCv1Idx = 0;
if (cv1VxTouchCv2[cv1Idx] && cv1VxTouchCv2Edge[cv1Idx] == -1) // Start of cv1 edge is touching inside the cv2 polygon (not on an cv2 edge)
{
if (polygon->empty() || polygon->back() != cv1CubeFaceIndices[cv1Idx])
{
polygon->push_back(cv1CubeFaceIndices[cv1Idx]);
}
if (cv1VxTouchCv2[nextCv1Idx] && cv1VxTouchCv2Edge[nextCv1Idx] == -1)
{
// Both ends of this cv1 edge is touching inside(not on an edge) cv2 polygon, no intersections possible (assuming convex cube-face)
// Continue with next Cv1-edge.
continue;
}
}
// Find intersection(s) on this edge
std::vector<IndexType> intersectionVxIndices;
std::vector<int> intersectedCv2EdgeIdxs;
std::vector<double> intersectionFractionsAlongEdge;
int nextCv2Idx = 1;
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx, ++nextCv2Idx)
{
if (nextCv2Idx > 3) nextCv2Idx = 0;
// Find a possible real intersection point.
cvf::Vec3d intersection(0,0,0);
double fractionAlongEdge1;
GeometryTools::IntersectionStatus intersectStatus = GeometryTools::NO_INTERSECTION;
IndexType intersectionVxIndex = cvf::UNDEFINED_UINT;
// First handle some "trivial" ones to ease the burden for the real intersection calculation
// It could be tested whether it really is necessary to do
if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx && cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
{
intersectStatus = GeometryTools::LINES_OVERLAP;
fractionAlongEdge1 = 1;
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
}
else if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx )
{
// When this happens, the cv1 vertex will already have been added to the polygon
// by the statements in the top of the cv1 edge loop. Should it be treated specially ?
intersectStatus = GeometryTools::LINES_TOUCH;
fractionAlongEdge1 = 0;
intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
}
else if (cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
{
intersectStatus = GeometryTools::LINES_TOUCH;
fractionAlongEdge1 = 1;
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
}
else
{
double fractionAlongEdge2;
bool found = edgeIntersectionStorage.findIntersection( cv1CubeFaceIndices[cv1Idx],
cv1CubeFaceIndices[nextCv1Idx],
cv2CubeFaceIndices[cv2Idx],
cv2CubeFaceIndices[nextCv2Idx],
&intersectionVxIndex, &intersectStatus,
&fractionAlongEdge1, &fractionAlongEdge2);
if (!found)
{
intersectStatus = GeometryTools::inPlaneLineIntersect3D(normal,
nodes[cv1CubeFaceIndices[cv1Idx]],
nodes[cv1CubeFaceIndices[nextCv1Idx]],
nodes[cv2CubeFaceIndices[cv2Idx]],
nodes[cv2CubeFaceIndices[nextCv2Idx]],
&intersection, &fractionAlongEdge1, &fractionAlongEdge2,
tolerance);
switch (intersectStatus)
{
case GeometryTools::LINES_CROSSES:
{
intersectionVxIndex = newVertexIndex;
createdVertexes->push_back(intersection);
++newVertexIndex;
}
break;
case GeometryTools::LINES_TOUCH:
{
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
else CVF_ASSERT(false); // Tolerance trouble
}
break;
case GeometryTools::LINES_OVERLAP:
{
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
else CVF_ASSERT(false); // Tolerance trouble
}
break;
default:
break;
}
edgeIntersectionStorage.addIntersection( cv1CubeFaceIndices[cv1Idx],
cv1CubeFaceIndices[nextCv1Idx],
cv2CubeFaceIndices[cv2Idx],
cv2CubeFaceIndices[nextCv2Idx],
intersectionVxIndex, intersectStatus,
fractionAlongEdge1, fractionAlongEdge2);
}
}
// Store data for each intersection along the Cv1-edge
if ( (intersectStatus == GeometryTools::LINES_CROSSES)
|| (intersectStatus == GeometryTools::LINES_TOUCH)
|| (intersectStatus == GeometryTools::LINES_OVERLAP) )
{
CVF_ASSERT(intersectionVxIndex != cvf::UNDEFINED_UINT);
intersectionFractionsAlongEdge.push_back(fractionAlongEdge1);
intersectedCv2EdgeIdxs.push_back(cv2Idx);
intersectionVxIndices.push_back(intersectionVxIndex);
}
}
// Insert the intersections into the polygon in the correct order along the Cv1 edge.
// Find the last intersection in order to possibly continue the polygon along Cv2 into Cv1
size_t i;
size_t lastIntersection = std::numeric_limits<size_t>::max();
double largestFraction = -1;
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
{
if (intersectionFractionsAlongEdge[i] > largestFraction)
{
lastIntersection = i;
largestFraction = intersectionFractionsAlongEdge[i];
}
}
// Insert indices to the new intersection vertices into the polygon of
// this connection according to fraction along edge
std::map<double,size_t> sortingMap;
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
{
sortingMap[intersectionFractionsAlongEdge[i]] = intersectionVxIndices[i];
}
std::map<double, size_t>::iterator it;
for (it = sortingMap.begin(); it != sortingMap.end(); ++it)
{
if (polygon->empty() || polygon->back() != it->second)
{
polygon->push_back(it->second);
}
}
// Then if the Cv1 edge is going out of Cv2, add to the polygon, all the Cv2 face vertex-indices
// from the intersection that touches Cv1.
// if cv1 edge in any way touches cv2 and ends up outside, it went out.
/*
if cv1 edge is going out of cv2 then
if intersected cv2 edge has endpoint touching cv1 then
add endpoint to polygon. continue to add next endpoint until it does not touch Cv1
*/
if ( !cv1VxTouchCv2[nextCv1Idx]
&& ( cv1VxTouchCv2[cv1Idx] || ( intersectedCv2EdgeIdxs.size() ) ) ) // Two touches along edge also qualifies
{
if(lastIntersection < intersectedCv2EdgeIdxs.size())
{
cv2Idx = intersectedCv2EdgeIdxs[lastIntersection];
int count = 0;
// Continue the polygon along the Cv2 edges as long as they touch cv1.
// Depending on the faces having opposite winding, which is guaranteed as long as
// no intersecting CVs share a connection
while (cv2VxTouchCv1[cv2Idx] && count < 4 && (cv2VxTouchCv1Edge[cv2Idx] == -1)) // Touch of edge is regarded as being outside, so we must stop
{
if (polygon->empty() || polygon->back() != cv2CubeFaceIndices[cv2Idx])
{
polygon->push_back(cv2CubeFaceIndices[cv2Idx]);
}
--cv2Idx;
if (cv2Idx < 0 ) cv2Idx = 3;
++count;
}
}
else
{
CVF_ASSERT(lastIntersection < intersectedCv2EdgeIdxs.size());
}
}
}
if (polygon->size() > 2)
{
if (polygon->back() == polygon->front()) polygon->pop_back();
}
// Sanity checks
if (polygon->size() < 3)
{
// cvf::Trace::show(cvf::String("Degenerated connection polygon detected. (Less than 3 vertexes) Cv's probably not in contact: %1 , %2").arg(m_ownerCvId).arg(m_neighborCvId));
polygon->clear();
return false;
}
return true;
}
//--------------------------------------------------------------------------------------------------
/// This method assumes that all intersection and mid edge vertexes are created an are already
/// merged into all the polygons. We can also assume that all the connection polygons are completely
/// inside (or sharing edges with) the cube face polygon initially
//--------------------------------------------------------------------------------------------------
#define DEBUG_PRINT 0
//template <typename NodeArrayType, typename NodeType, typename IndicesArrayType, typename IndicesType>
//void setup( ArrayWrapper<NodeArrayType, NodeType> nodeArray, ArrayWrapper<IndicesArrayType, IndicesType> indices)
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
ArrayWrapperConst<PolygonArrayType, IndexType> completeFacePolygon,
const cvf::Vec3d& faceNormal,
const std::vector< std::vector<IndexType>* >& faceOverlapPolygons,
const std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags,
std::vector<IndexType>* partialFacePolygon,
bool* m_partiallyFreeCubeFaceHasHoles)
{
// Vertex Index to position in polygon
typedef std::map< size_t, std::vector<size_t >::const_iterator > VxIdxToPolygonPositionMap;
CVF_ASSERT(m_partiallyFreeCubeFaceHasHoles);
CVF_ASSERT(partialFacePolygon != NULL);
// Copy the start polygon
std::list<IndexType> resultPolygon;
for (size_t pcIdx = 0; pcIdx < completeFacePolygon.size(); ++pcIdx)
{
resultPolygon.push_back(completeFacePolygon[pcIdx]);
}
// First build search maps to fast find whether and where an index is positioned in a polygon
// Map from Vertex-index to position in polygon
std::vector< VxIdxToPolygonPositionMap > polygonSearchMaps;
std::vector<bool> isConnectionPolygonMerged;
polygonSearchMaps.resize(faceOverlapPolygons.size());
isConnectionPolygonMerged.resize(faceOverlapPolygons.size(), false);
// Build search maps
{
size_t count;
for (size_t i = 0; i < faceOverlapPolygons.size(); ++i)
{
count = 0;
for (std::vector<IndexType >::const_iterator pcIt = faceOverlapPolygons[i]->begin();
pcIt != faceOverlapPolygons[i]->end();
++pcIt)
{
polygonSearchMaps[i][*pcIt] = pcIt;
++count;
}
if (count < 3) isConnectionPolygonMerged[i] = true; // Ignore false polygons
}
}
#if DEBUG_PRINT
{
cvf::Trace::show("Circumference polygon: ");
std::list<IndexType>::const_iterator polIt;
for ( polIt = resultPolygon.begin(); polIt != resultPolygon.end(); ++polIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*polIt))
.arg(nodeCoords[*polIt].x())
.arg(nodeCoords[*polIt].y())
.arg(nodeCoords[*polIt].z()));
}
}
#endif
#if DEBUG_PRINT
{
cvf::Trace::show("Connection polygons: ");
for (size_t cIdx = 0; cIdx < faceOverlapPolygons.size(); cIdx++)
{
std::vector<IndexType >::const_iterator polIt;
cvf::Trace::show("Connection " + cvf::String((long long)cIdx));
for (polIt = faceOverlapPolygons[cIdx]->begin(); polIt != faceOverlapPolygons[cIdx]->end(); ++polIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*polIt))
.arg(nodeCoords[*polIt].x())
.arg(nodeCoords[*polIt].y())
.arg(nodeCoords[*polIt].z()));
}
}
}
#endif
// Merge connection polygons with the main polygon as long as one of them has something in common.
// For each vx in the m_freeFacePolygons[cubeFace] polygon .
// loop over all connections
// if it has the node in common and that the edge angle will decrease if inserting
// merge the connection polygon into the main polygon,
// and remove the connection polygon from the merge able connection polygons.
for (std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
// Set iterator to previous node in polygon
std::list<IndexType>::iterator prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
cvf::Vec3d pToPrev = nodeCoords[*prevPIt] - nodeCoords[*pIt];
// Set iterator to next node in polygon. Used to insert before and as pointer to the next point
std::list<IndexType>::iterator nextPIt = pIt;
++nextPIt;
std::list<IndexType>::iterator insertBeforePIt = nextPIt;
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
// Calculate existing edge to edge angle
cvf::Vec3d pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
double mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
// Find connections containing the pIt vertex index. Merge them into the main polygon
for (size_t opIdx = 0; opIdx < faceOverlapPolygons.size(); ++opIdx)
{
if (isConnectionPolygonMerged[opIdx]) continue; // Already merged
// Find position of pIt vertex index in the current connection polygon
VxIdxToPolygonPositionMap::iterator vxIndexPositionInPolygonIt = polygonSearchMaps[opIdx].find(*pIt);
if (vxIndexPositionInPolygonIt != polygonSearchMaps[opIdx].end())
{
// Merge the connection polygon into the main polygon
// if the angle prevPIt pIt nextPIt is larger than angle prevPIt pIt (startCPIt++)
std::vector<IndexType>::const_iterator startCPIt;
startCPIt = vxIndexPositionInPolygonIt->second;
// First vx to insert is the one after the match
bool hasSameWinding = faceOverlapPolygonWindingSameAsCubeFaceFlags[opIdx];
if (hasSameWinding)
{
// Same winding as main polygon. We need to go the opposite way
if (startCPIt == faceOverlapPolygons[opIdx]->begin()) startCPIt = faceOverlapPolygons[opIdx]->end();
--startCPIt;
}
else
{
// Opposite winding. Go forward when merging
++startCPIt; if (startCPIt == faceOverlapPolygons[opIdx]->end()) startCPIt = faceOverlapPolygons[opIdx]->begin();
}
// Calculate possible new edge-to-edge angle and test against existing angle
cvf::Vec3d pToStart = nodeCoords[*startCPIt] - nodeCoords[*pIt];
double candidatePolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToStart , pToPrev);
if (candidatePolygonEdgeAngle < mainPolygonEdgeAngle )
{
// Merge ok
std::vector<IndexType >::const_iterator pcIt = startCPIt;
if (hasSameWinding)
{
do
{
resultPolygon.insert(insertBeforePIt, (*pcIt));
if (pcIt == faceOverlapPolygons[opIdx]->begin()) pcIt = faceOverlapPolygons[opIdx]->end();
--pcIt;
} while (pcIt != startCPIt);
}
else
{
do
{
resultPolygon.insert(insertBeforePIt, (*pcIt));
++pcIt;
if (pcIt == faceOverlapPolygons[opIdx]->end()) pcIt = faceOverlapPolygons[opIdx]->begin();
} while (pcIt != startCPIt);
}
isConnectionPolygonMerged[opIdx] = true;
// Recalculate the next position to point into the new nodes
// Set iterator in the main polygon to insert before and to the next point
nextPIt = pIt;
++nextPIt;
insertBeforePIt = nextPIt;
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
// Recalculate the existing edge to edge angle
pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
}
}
}
}
// Now remove all double edges
bool goneAround = false;
for ( std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end() && !goneAround; ++pIt)
{
// Set iterator to next node in polygon.
std::list<IndexType>::iterator nextPIt = pIt;
++nextPIt;
if (nextPIt == resultPolygon.end())
{
nextPIt = resultPolygon.begin();
goneAround = true; // Gone around polygon. Stop even if pIt is jumping over end()
}
// Set iterator to previous node in polygon
std::list<IndexType>::iterator prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
// If previous and next node are the same, erase
while(*nextPIt == *prevPIt)
{
resultPolygon.erase(pIt);
resultPolygon.erase(prevPIt);
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of while
pIt = nextPIt;
++nextPIt;
if (nextPIt == resultPolygon.end())
{
nextPIt = resultPolygon.begin();
goneAround = true; // Gone around polygon pIt is jumping over end()
}
prevPIt = pIt;
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
--prevPIt;
}
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of for loop
}
// Check for holes
bool hasHoles = false;
for (size_t i = 0; i < isConnectionPolygonMerged.size(); ++i)
{
hasHoles = !isConnectionPolygonMerged[i];
if(hasHoles)
{
*m_partiallyFreeCubeFaceHasHoles = true;
break;
}
}
#if DEBUG_PRINT
{
cvf::Trace::show("Polygon: ");
for (std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*pIt))
.arg(nodeCoords[*pIt].x())
.arg(nodeCoords[*pIt].y())
.arg(nodeCoords[*pIt].z()));
}
}
#endif
// Copy the result polygon to the output variable
partialFacePolygon->clear();
for (std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
{
partialFacePolygon->push_back(*pIt);
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template <typename IndexType>
void EdgeIntersectStorage<IndexType>::setVertexCount(IndexType size)
{
m_edgeIntsectMap.resize(size);
}
template <typename IndexType>
void EdgeIntersectStorage<IndexType>::canonizeAddress(IndexType& e1P1, IndexType& e1P2, IndexType& e2P1, IndexType& e2P2,
bool& flipE1, bool& flipE2, bool& flipE1E2)
{
flipE1 = e1P1 > e1P2;
flipE2 = e2P1 > e2P2;
flipE1E2 = (flipE1 ? e1P2: e1P1) > (flipE2 ? e2P2: e2P1);
static size_t temp;
if (flipE1)
{
temp = e1P1;
e1P1 = e1P2;
e1P2 = temp;
}
if (flipE2)
{
temp = e2P1;
e2P1 = e2P2;
e2P2 = temp;
}
if (flipE1E2)
{
temp = e1P1;
e1P1 = e2P1;
e2P1 = temp;
temp = e1P2;
e1P2 = e2P2;
e2P2 = temp;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template <typename IndexType>
void EdgeIntersectStorage<IndexType>::addIntersection(IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
IndexType vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
double fractionAlongEdge1, double fractionAlongEdge2)
{
static bool flipE1 ;
static bool flipE2 ;
static bool flipE1E2;
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
static IntersectData iData;
iData.fractionAlongEdge1 = flipE1 ? 1 - fractionAlongEdge1 : fractionAlongEdge1;
iData.fractionAlongEdge2 = flipE2 ? 1 - fractionAlongEdge2 : fractionAlongEdge2;
iData.intersectionStatus = intersectionStatus;
if (flipE1E2)
{
double temp = iData.fractionAlongEdge1;
iData.fractionAlongEdge1 = iData.fractionAlongEdge2;
iData.fractionAlongEdge2 = temp;
}
iData.intersectionPointIndex = vxIndexIntersectionPoint;
CVF_ASSERT(e1P1 < m_edgeIntsectMap.size());
m_edgeIntsectMap[e1P1][e1P2][e2P1][e2P2] = iData;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template <typename IndexType>
bool EdgeIntersectStorage<IndexType>::findIntersection(IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
double* fractionAlongEdge1, double* fractionAlongEdge2)
{
static bool flipE1 ;
static bool flipE2 ;
static bool flipE1E2;
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
if (!m_edgeIntsectMap[e1P1].size()) return false;
std::map<size_t, std::map<size_t, std::map<size_t, IntersectData > > >::iterator it;
it = m_edgeIntsectMap[e1P1].find(e1P2);
if (it == m_edgeIntsectMap[e1P1].end()) return false;
std::map<size_t, std::map<size_t, IntersectData > >::iterator it2;
it2 = it->second.find(e2P1);
if (it2 == it->second.end()) return false;
std::map<size_t, IntersectData >::iterator it3;
it3 = it2->second.find(e2P2);
if (it3 == it2->second.end()) return false;
*vxIndexIntersectionPoint = it3->second.intersectionPointIndex;
*intersectionStatus = it3->second.intersectionStatus;
if (flipE1E2)
{
*fractionAlongEdge1 = it3->second.fractionAlongEdge2;
*fractionAlongEdge2 = it3->second.fractionAlongEdge1;
}
else
{
*fractionAlongEdge1 = it3->second.fractionAlongEdge1;
*fractionAlongEdge2 = it3->second.fractionAlongEdge2;
}
if (flipE1) *fractionAlongEdge1 = 1 - *fractionAlongEdge1;
if (flipE2) *fractionAlongEdge2 = 1 - *fractionAlongEdge2;
return true;
}
}