From 597cfd4d30144513bb9be3bad69bc5c2dd5e8dcc Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Mon, 3 Nov 2025 16:06:58 +0100 Subject: [PATCH] Improves geometry tools robustness and accuracy * Replace asserts with alternatives * Refactors GeometryTools for improved robustness --- .../RigCaseCellResultsData.cpp | 4 +- .../RigCellFaceGeometryTools.cpp | 8 +- .../ReservoirDataModel/cvfGeometryTools.cpp | 123 ++++--- .../ReservoirDataModel/cvfGeometryTools.h | 21 +- .../ReservoirDataModel/cvfGeometryTools.inl | 167 +++++----- .../UnitTests/cvfGeometryTools-Test.cpp | 304 ++++++++++++++---- 6 files changed, 402 insertions(+), 225 deletions(-) diff --git a/ApplicationLibCode/ReservoirDataModel/RigCaseCellResultsData.cpp b/ApplicationLibCode/ReservoirDataModel/RigCaseCellResultsData.cpp index 12b0980f72..ca787329ff 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigCaseCellResultsData.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigCaseCellResultsData.cpp @@ -1994,8 +1994,8 @@ void calculateConnectionGeometry( const RigCell& c1, c1.faceIndices( faceId, &face1 ); c2.faceIndices( cvf::StructGridInterface::oppositeFace( faceId ), &face2 ); - bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &intersections, + bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + intersections, (cvf::EdgeIntersectStorage*)nullptr, cvf::wrapArrayConst( &nodes ), face1.data(), diff --git a/ApplicationLibCode/ReservoirDataModel/RigCellFaceGeometryTools.cpp b/ApplicationLibCode/ReservoirDataModel/RigCellFaceGeometryTools.cpp index 4fc5450221..17b4e78412 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigCellFaceGeometryTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigCellFaceGeometryTools.cpp @@ -94,8 +94,8 @@ cvf::StructGridInterface::FaceType RigCellFaceGeometryTools::calculateCellFaceOv c1.faceIndices( ( cvf::StructGridInterface::FaceType )( fIdx ), &face1 ); c2.faceIndices( cvf::StructGridInterface::oppositeFace( ( cvf::StructGridInterface::FaceType )( fIdx ) ), &face2 ); - bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &intersections, + bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + intersections, (cvf::EdgeIntersectStorage*)nullptr, cvf::wrapArrayConst( &mainGrid.nodes() ), face1.data(), @@ -335,8 +335,8 @@ void RigCellFaceGeometryTools::extractConnectionsForFace( const RigFault::FaultF std::array candidateFaceIndices; mainGrid->cell( candidateCellIndex ).faceIndices( candidateFace, &candidateFaceIndices ); - bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &intersections, + bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + intersections, (cvf::EdgeIntersectStorage*)nullptr, cvf::wrapArrayConst( &mainGridNodes ), sourceFaceIndices.data(), diff --git a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp index 78620eb1f3..ac92c8cb57 100644 --- a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp @@ -239,30 +239,20 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect( double x1, double y4, double l1NormalizedTolerance, double l2NormalizedTolerance, - double* x, - double* y, - double* fractionAlongLine1, - double* fractionAlongLine2 ) + double& x, + double& y, + double& fractionAlongLine1, + double& fractionAlongLine2 ) { - double mua, mub; - double denom, numera, numerb; - - denom = ( y4 - y3 ) * ( x2 - x1 ) - ( x4 - x3 ) * ( y2 - y1 ); - numera = ( x4 - x3 ) * ( y1 - y3 ) - ( y4 - y3 ) * ( x1 - x3 ); - numerb = ( x2 - x1 ) * ( y1 - y3 ) - ( y2 - y1 ) * ( x1 - x3 ); + auto denom = ( y4 - y3 ) * ( x2 - x1 ) - ( x4 - x3 ) * ( y2 - y1 ); + auto numera = ( x4 - x3 ) * ( y1 - y3 ) - ( y4 - y3 ) * ( x1 - x3 ); + auto numerb = ( x2 - x1 ) * ( y1 - y3 ) - ( y2 - y1 ) * ( x1 - x3 ); double EPS = 1e-40; // Are the line coincident? if ( fabs( numera ) < EPS && fabs( numerb ) < EPS && fabs( denom ) < EPS ) { -#if 0 - *x = 0; - *y = 0; - *fractionAlongLine1 = 0; - *fractionAlongLine2 = 0; - return GeometryTools::LINES_OVERLAP; -#else cvf::Vec3d p12( x2 - x1, y2 - y1, 0 ); cvf::Vec3d p13( x3 - x1, y3 - y1, 0 ); cvf::Vec3d p34( x4 - x3, y4 - y3, 0 ); @@ -274,10 +264,10 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect( double x1, if ( length12 < EPS ) { - *x = x1; - *y = y1; - *fractionAlongLine1 = 1; - *fractionAlongLine2 = p13.length() / p34.length(); + x = x1; + y = y1; + fractionAlongLine1 = 1; + fractionAlongLine2 = p13.length() / p34.length(); return GeometryTools::LINES_OVERLAP; } @@ -292,10 +282,10 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect( double x1, if ( ( normDist13 < 0 - l1NormalizedTolerance && normDist14 < 0 - l1NormalizedTolerance ) || ( normDist13 > 1 + l1NormalizedTolerance && normDist14 > 1 + l1NormalizedTolerance ) ) { - *x = 0; - *y = 0; - *fractionAlongLine1 = 0; - *fractionAlongLine2 = 0; + x = 0; + y = 0; + fractionAlongLine1 = 0; + fractionAlongLine2 = 0; return GeometryTools::NO_INTERSECTION; } @@ -311,68 +301,67 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect( double x1, if ( pt3IsInside && !pt4IsInside ) { - *fractionAlongLine1 = normDist13; - *fractionAlongLine2 = 0.0; - *x = x3; - *y = y3; + fractionAlongLine1 = normDist13; + fractionAlongLine2 = 0.0; + x = x3; + y = y3; } else if ( pt4IsInside && !pt3IsInside ) { - *fractionAlongLine1 = normDist14; - *fractionAlongLine2 = 1.0; - *x = x4; - *y = y4; + fractionAlongLine1 = normDist14; + fractionAlongLine2 = 1.0; + x = x4; + y = y4; } else if ( pt3IsInside && pt4IsInside ) { // Return edge 2 vertex furthest along edge 1 if ( normDist13 <= normDist14 ) { - *fractionAlongLine1 = normDist14; - *fractionAlongLine2 = 1.0; - *x = x4; - *y = y4; + fractionAlongLine1 = normDist14; + fractionAlongLine2 = 1.0; + x = x4; + y = y4; } else { - *fractionAlongLine1 = normDist13; - *fractionAlongLine2 = 0.0; - *x = x3; - *y = y3; + fractionAlongLine1 = normDist13; + fractionAlongLine2 = 0.0; + x = x3; + y = y3; } } else // both outside on each side { // Return End of edge 1 - *fractionAlongLine1 = 1.0; - *fractionAlongLine2 = normDist32; - *x = x2; - *y = y2; + fractionAlongLine1 = 1.0; + fractionAlongLine2 = normDist32; + x = x2; + y = y2; } return GeometryTools::LINES_OVERLAP; -#endif } /* Are the line parallel */ if ( fabs( denom ) < EPS ) { - *x = 0; - *y = 0; - *fractionAlongLine1 = 0; - *fractionAlongLine2 = 0; + x = 0; + y = 0; + fractionAlongLine1 = 0; + fractionAlongLine2 = 0; return GeometryTools::NO_INTERSECTION; } /* Is the intersection along the the segments */ - mua = numera / denom; - mub = numerb / denom; + auto mua = numera / denom; + auto mub = numerb / denom; - *x = x1 + mua * ( x2 - x1 ); - *y = y1 + mua * ( y2 - y1 ); - *fractionAlongLine1 = mua; - *fractionAlongLine2 = mub; + x = x1 + mua * ( x2 - x1 ); + y = y1 + mua * ( y2 - y1 ); + fractionAlongLine1 = mua; + fractionAlongLine2 = mub; if ( mua < 0 - l1NormalizedTolerance || 1 + l1NormalizedTolerance < mua || mub < 0 - l2NormalizedTolerance || 1 + l2NormalizedTolerance < mub ) @@ -382,10 +371,10 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect( double x1, else if ( fabs( mua ) < l1NormalizedTolerance || fabs( 1 - mua ) < l1NormalizedTolerance || fabs( mub ) < l2NormalizedTolerance || fabs( 1 - mub ) < l2NormalizedTolerance ) { - if ( fabs( mua ) < l1NormalizedTolerance ) *fractionAlongLine1 = 0; - if ( fabs( 1 - mua ) < l1NormalizedTolerance ) *fractionAlongLine1 = 1; - if ( fabs( mub ) < l2NormalizedTolerance ) *fractionAlongLine2 = 0; - if ( fabs( 1 - mub ) < l2NormalizedTolerance ) *fractionAlongLine2 = 1; + if ( fabs( mua ) < l1NormalizedTolerance ) fractionAlongLine1 = 0; + if ( fabs( 1 - mua ) < l1NormalizedTolerance ) fractionAlongLine1 = 1; + if ( fabs( mub ) < l2NormalizedTolerance ) fractionAlongLine2 = 0; + if ( fabs( 1 - mub ) < l2NormalizedTolerance ) fractionAlongLine2 = 1; return GeometryTools::LINES_TOUCH; } else @@ -403,13 +392,11 @@ GeometryTools::IntersectionStatus GeometryTools::inPlaneLineIntersect3D( const c const cvf::Vec3d& p2, const cvf::Vec3d& p3, const cvf::Vec3d& p4, - cvf::Vec3d* intersectionPoint, - double* fractionAlongLine1, - double* fractionAlongLine2, + cvf::Vec3d& intersectionPoint, + double& fractionAlongLine1, + double& fractionAlongLine2, double tolerance ) { - CVF_ASSERT( intersectionPoint != nullptr ); - int Z = findClosestAxis( planeNormal ); int X = ( Z + 1 ) % 3; int Y = ( Z + 2 ) % 3; @@ -420,16 +407,16 @@ GeometryTools::IntersectionStatus GeometryTools::inPlaneLineIntersect3D( const c double l2NormTol = tolerance / ( p4 - p3 ).length(); IntersectionStatus intersectionStatus = - inPlaneLineIntersect( p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], p4[X], p4[Y], l1NormTol, l2NormTol, &x, &y, fractionAlongLine1, fractionAlongLine2 ); + inPlaneLineIntersect( p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], p4[X], p4[Y], l1NormTol, l2NormTol, x, y, fractionAlongLine1, fractionAlongLine2 ); // Check if we have a valid intersection point if ( intersectionStatus == NO_INTERSECTION || intersectionStatus == LINES_OVERLAP ) { - intersectionPoint->setZero(); + intersectionPoint.setZero(); } else { - *intersectionPoint = p1 + ( *fractionAlongLine1 ) * ( p2 - p1 ); + intersectionPoint = p1 + ( fractionAlongLine1 ) * ( p2 - p1 ); } return intersectionStatus; diff --git a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h index c00bc417a5..e539c0e534 100644 --- a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h +++ b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h @@ -23,6 +23,9 @@ #include "cvfArrayWrapperConst.h" #include "cvfMatrix3.h" #include "cvfMatrix4.h" + +#include "cafPdmLogging.h" + #include #include @@ -93,7 +96,7 @@ public: std::vector* createdVertexes ); template - static bool insertVertexInPolygon( std::vector* polygon, + static bool insertVertexInPolygon( std::vector& polygon, ArrayWrapperConst nodeCoords, IndexType vertexIndex, double tolerance ); @@ -103,9 +106,9 @@ public: const cvf::Vec3d& p2, const cvf::Vec3d& p3, const cvf::Vec3d& p4, - cvf::Vec3d* intersectionPoint, - double* fractionAlongLine1, - double* fractionAlongLine2, + cvf::Vec3d& intersectionPoint, + double& fractionAlongLine1, + double& fractionAlongLine2, double tolerance = 1e-6 ); template @@ -113,12 +116,12 @@ public: ArrayWrapperConst vertices, ArrayWrapperConst indices, const cvf::Vec3d& point, - int* touchedEdgeIndex, + int& touchedEdgeIndex, double tolerance = 1e-6 ); template - static bool calculateOverlapPolygonOfTwoQuads( std::vector* polygon, - std::vector* createdVertexes, + static bool calculateOverlapPolygonOfTwoQuads( std::vector& polygon, + std::vector& createdVertexes, EdgeIntersectStorage* edgeIntersectionStorage, ArrayWrapperConst nodes, const IndexType cv1CubeFaceIndices[4], @@ -148,8 +151,8 @@ public: IndexType e2P2, IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus, - double* fractionAlongEdge1, - double* fractionAlongEdge2 ); + double& fractionAlongEdge1, + double& fractionAlongEdge2 ); void addIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, diff --git a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.inl b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.inl index b0d0d27036..f5468527f5 100644 --- a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.inl +++ b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.inl @@ -27,6 +27,11 @@ namespace cvf template Vec3Type GeometryTools::computePolygonCenter( const std::vector& polygon ) { + if ( polygon.empty() ) + { + return Vec3Type::UNDEFINED; + } + Vec3Type s; for ( size_t i = 0; i < polygon.size(); i++ ) @@ -65,25 +70,22 @@ DataType GeometryTools::interpolateQuad( const cvf::Vec3d& v1, /// Todo: If a vertex is replaced, the VxToCv map in TimeStepGeometry should be updated //-------------------------------------------------------------------------------------------------- template -bool GeometryTools::insertVertexInPolygon( std::vector* polygon, +bool GeometryTools::insertVertexInPolygon( std::vector& polygon, ArrayWrapperConst nodeCoords, IndexType vertexIndex, double tolerance ) { - CVF_ASSERT( polygon ); - // Check if vertex is directly included already - for ( typename std::vector::iterator it = polygon->begin(); it != polygon->end(); ++it ) + for ( typename std::vector::iterator it = polygon.begin(); it != polygon.end(); ++it ) { if ( *it == vertexIndex ) return true; } -#if 1 // Check if the new point is within tolerance of one of the polygon vertices bool existsOrInserted = false; - for ( typename std::vector::iterator it = polygon->begin(); it != polygon->end(); ++it ) + for ( typename std::vector::iterator it = polygon.begin(); it != polygon.end(); ++it ) { if ( ( nodeCoords[*it] - nodeCoords[vertexIndex] ).length() < tolerance ) { @@ -93,14 +95,13 @@ bool GeometryTools::insertVertexInPolygon( std::vector* } if ( existsOrInserted ) return true; -#endif // Copy the start polygon to a list std::list listPolygon; - for ( size_t pcIdx = 0; pcIdx < polygon->size(); ++pcIdx ) + for ( size_t pcIdx = 0; pcIdx < polygon.size(); ++pcIdx ) { - listPolygon.push_back( ( *polygon )[pcIdx] ); + listPolygon.push_back( polygon[pcIdx] ); } // Insert vertex in polygon if the distance to one of the edges is small enough @@ -125,10 +126,10 @@ bool GeometryTools::insertVertexInPolygon( std::vector* // Write polygon back into the vector - polygon->clear(); + polygon.clear(); for ( typename std::list::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it ) { - polygon->push_back( *it ); + polygon.push_back( *it ); } return existsOrInserted; @@ -157,7 +158,7 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& cvf::ArrayWrapperConst vertices, cvf::ArrayWrapperConst indices, const cvf::Vec3d& point, - int* touchedEdgeIndex, + int& touchedEdgeIndex, double tolerance ) { size_t numIndices = indices.size(); @@ -183,9 +184,7 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& size_t firstIdx; size_t secondIdx; - CVF_TIGHT_ASSERT( touchedEdgeIndex ); - - *touchedEdgeIndex = -1; + touchedEdgeIndex = -1; for ( firstIdx = 0, secondIdx = 1; firstIdx < numIndices; ++firstIdx, ++secondIdx ) { if ( secondIdx >= numIndices ) secondIdx = 0; @@ -195,7 +194,7 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& double sqDist = GeometryTools::linePointSquareDist( vx0, vx1, point ); if ( sqDist < tolerance * tolerance ) { - *touchedEdgeIndex = static_cast( firstIdx ); + touchedEdgeIndex = static_cast( firstIdx ); return true; } } @@ -247,20 +246,17 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& /// The second face must have opposite winding of the first //-------------------------------------------------------------------------------------------------- template -bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* polygon, - std::vector* createdVertexes, +bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector& polygon, + std::vector& createdVertexes, EdgeIntersectStorage* edgeIntersectionStorage, ArrayWrapperConst nodes, const IndexType cv1CubeFaceIndices[4], const IndexType cv2CubeFaceIndices[4], double tolerance ) { - CVF_ASSERT( polygon ); - CVF_ASSERT( createdVertexes ); - // Topology analysis - IndexType newVertexIndex = static_cast( nodes.size() + createdVertexes->size() ); + IndexType newVertexIndex = static_cast( nodes.size() + createdVertexes.size() ); bool cv1VxTouchCv2[4] = { false, false, false, false }; bool cv2VxTouchCv1[4] = { false, false, false, false }; @@ -291,7 +287,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* int k; for ( k = 0; k < 4; ++k ) { - polygon->push_back( cv1CubeFaceIndices[k] ); + polygon.push_back( cv1CubeFaceIndices[k] ); } return true; @@ -311,7 +307,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* nodes, wrapArrayConst( cv2CubeFaceIndices, 4 ), nodes[cv1CubeFaceIndices[cv1Idx]], - &( cv1VxTouchCv2Edge[cv1Idx] ), + cv1VxTouchCv2Edge[cv1Idx], tolerance ); if ( cv1VxTouchCv2[cv1Idx] ) ++numCv1VxesOnCv2; } @@ -322,7 +318,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* nodes, wrapArrayConst( cv1CubeFaceIndices, 4 ), nodes[cv2CubeFaceIndices[cv1Idx]], - &( cv2VxTouchCv1Edge[cv1Idx] ), + cv2VxTouchCv1Edge[cv1Idx], tolerance ); if ( cv2VxTouchCv1[cv1Idx] ) ++numCv2VxesOnCv1; } @@ -335,7 +331,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* int k; for ( k = 0; k < 4; ++k ) { - polygon->push_back( cv1CubeFaceIndices[k] ); + polygon.push_back( cv1CubeFaceIndices[k] ); } return true; @@ -346,7 +342,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* int k; for ( k = 3; k >= 0; --k ) // Return opposite winding, to match winding of face 1 { - polygon->push_back( cv2CubeFaceIndices[k] ); + polygon.push_back( cv2CubeFaceIndices[k] ); } return true; } @@ -367,9 +363,9 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* 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] ) + if ( polygon.empty() || polygon.back() != cv1CubeFaceIndices[cv1Idx] ) { - polygon->push_back( cv1CubeFaceIndices[cv1Idx] ); + polygon.push_back( cv1CubeFaceIndices[cv1Idx] ); } if ( cv1VxTouchCv2[nextCv1Idx] && cv1VxTouchCv2Edge[nextCv1Idx] == -1 ) @@ -431,8 +427,8 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* cv2CubeFaceIndices[nextCv2Idx], &intersectionVxIndex, &intersectStatus, - &fractionAlongEdge1, - &fractionAlongEdge2 ); + fractionAlongEdge1, + fractionAlongEdge2 ); if ( !found ) { @@ -441,9 +437,9 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* nodes[cv1CubeFaceIndices[nextCv1Idx]], nodes[cv2CubeFaceIndices[cv2Idx]], nodes[cv2CubeFaceIndices[nextCv2Idx]], - &intersection, - &fractionAlongEdge1, - &fractionAlongEdge2, + intersection, + fractionAlongEdge1, + fractionAlongEdge2, tolerance ); switch ( intersectStatus ) @@ -451,7 +447,7 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* case GeometryTools::LINES_CROSSES: { intersectionVxIndex = newVertexIndex; - createdVertexes->push_back( intersection ); + createdVertexes.push_back( intersection ); ++newVertexIndex; } break; @@ -466,7 +462,13 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* else if ( fractionAlongEdge2 >= 1.0 ) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx]; else - CVF_ASSERT( false ); // Tolerance trouble + { + // Numerical precision issue at edge boundary + // Create intersection vertex at the computed intersection point as fallback + intersectionVxIndex = newVertexIndex; + createdVertexes.push_back( intersection ); + ++newVertexIndex; + } } break; case GeometryTools::LINES_OVERLAP: @@ -480,7 +482,13 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* else if ( fractionAlongEdge2 >= 1.0 ) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx]; else - CVF_ASSERT( false ); // Tolerance trouble + { + // Numerical precision issue at edge boundary + // Create intersection vertex at the computed intersection point as fallback + intersectionVxIndex = newVertexIndex; + createdVertexes.push_back( intersection ); + ++newVertexIndex; + } } break; default: @@ -504,7 +512,11 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* if ( ( intersectStatus == GeometryTools::LINES_CROSSES ) || ( intersectStatus == GeometryTools::LINES_TOUCH ) || ( intersectStatus == GeometryTools::LINES_OVERLAP ) ) { - CVF_ASSERT( intersectionVxIndex != cvf::UNDEFINED_UINT ); + if ( intersectionVxIndex == cvf::UNDEFINED_UINT ) + { + // Invalid intersection - skip and continue processing + continue; + } intersectionFractionsAlongEdge.push_back( fractionAlongEdge1 ); intersectedCv2EdgeIdxs.push_back( cv2Idx ); @@ -539,9 +551,9 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* typename std::map::iterator it; for ( it = sortingMap.begin(); it != sortingMap.end(); ++it ) { - if ( polygon->empty() || polygon->back() != it->second ) + if ( polygon.empty() || polygon.back() != it->second ) { - polygon->push_back( it->second ); + polygon.push_back( it->second ); } } @@ -570,9 +582,9 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* // outside, so we // must stop { - if ( polygon->empty() || polygon->back() != cv2CubeFaceIndices[cv2Idx] ) + if ( polygon.empty() || polygon.back() != cv2CubeFaceIndices[cv2Idx] ) { - polygon->push_back( cv2CubeFaceIndices[cv2Idx] ); + polygon.push_back( cv2CubeFaceIndices[cv2Idx] ); } --cv2Idx; if ( cv2Idx < 0 ) cv2Idx = 3; @@ -581,22 +593,28 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector* } else { - CVF_ASSERT( lastIntersection < intersectedCv2EdgeIdxs.size() ); + if ( lastIntersection >= intersectedCv2EdgeIdxs.size() ) + { + // Array bounds exceeded - exit processing gracefully + break; + } } } } - if ( polygon->size() > 2 ) + if ( polygon.size() > 2 ) { - if ( polygon->back() == polygon->front() ) polygon->pop_back(); + if ( polygon.back() == polygon.front() ) polygon.pop_back(); } - // Sanity checks - if ( polygon->size() < 3 ) + // Sanity checks with logging + 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(); + caf::PdmLogging::warning( QString( "GeometryTools: Generated degenerate polygon with %1 vertices. " + "Input faces may not overlap significantly. Tolerance=%2" ) + .arg( polygon.size() ) + .arg( tolerance ) ); + polygon.clear(); return false; } @@ -910,7 +928,7 @@ void EdgeIntersectStorage::canonizeAddress( IndexType& e1P1, flipE1E2 = ( flipE1 ? e1P2 : e1P1 ) > ( flipE2 ? e2P2 : e2P1 ); - static IndexType temp; + IndexType temp; if ( flipE1 ) { temp = e1P1; @@ -949,13 +967,13 @@ void EdgeIntersectStorage::addIntersection( IndexType double fractionAlongEdge1, double fractionAlongEdge2 ) { - static bool flipE1; - static bool flipE2; - static bool flipE1E2; + bool flipE1; + bool flipE2; + bool flipE1E2; canonizeAddress( e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2 ); - static IntersectData iData; + IntersectData iData; iData.fractionAlongEdge1 = flipE1 ? 1 - fractionAlongEdge1 : fractionAlongEdge1; iData.fractionAlongEdge2 = flipE2 ? 1 - fractionAlongEdge2 : fractionAlongEdge2; @@ -969,7 +987,9 @@ void EdgeIntersectStorage::addIntersection( IndexType } iData.intersectionPointIndex = vxIndexIntersectionPoint; - CVF_ASSERT( e1P1 < m_edgeIntsectMap.size() ); + if ( e1P1 >= m_edgeIntsectMap.size() || e1P2 >= m_edgeIntsectMap.size() || e2P1 >= m_edgeIntsectMap.size() || + e2P2 >= m_edgeIntsectMap.size() ) + return; m_edgeIntsectMap[e1P1][e1P2][e2P1][e2P2] = iData; } @@ -983,27 +1003,26 @@ bool EdgeIntersectStorage::findIntersection( IndexType IndexType e2P2, IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus, - double* fractionAlongEdge1, - double* fractionAlongEdge2 ) + double& fractionAlongEdge1, + double& fractionAlongEdge2 ) { - static bool flipE1; - static bool flipE2; - static bool flipE1E2; + bool flipE1; + bool flipE2; + bool flipE1E2; canonizeAddress( e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2 ); - if ( !m_edgeIntsectMap[e1P1].size() ) return false; + if ( e1P1 >= m_edgeIntsectMap.size() || e1P2 >= m_edgeIntsectMap.size() || e2P1 >= m_edgeIntsectMap.size() || + e2P2 >= m_edgeIntsectMap.size() || !m_edgeIntsectMap[e1P1].size() ) + return false; - typename std::map>>::iterator it; - it = m_edgeIntsectMap[e1P1].find( e1P2 ); + auto it = m_edgeIntsectMap[e1P1].find( e1P2 ); if ( it == m_edgeIntsectMap[e1P1].end() ) return false; - typename std::map>::iterator it2; - it2 = it->second.find( e2P1 ); + auto it2 = it->second.find( e2P1 ); if ( it2 == it->second.end() ) return false; - typename std::map::iterator it3; - it3 = it2->second.find( e2P2 ); + auto it3 = it2->second.find( e2P2 ); if ( it3 == it2->second.end() ) return false; *vxIndexIntersectionPoint = it3->second.intersectionPointIndex; @@ -1011,17 +1030,17 @@ bool EdgeIntersectStorage::findIntersection( IndexType if ( flipE1E2 ) { - *fractionAlongEdge1 = it3->second.fractionAlongEdge2; - *fractionAlongEdge2 = it3->second.fractionAlongEdge1; + fractionAlongEdge1 = it3->second.fractionAlongEdge2; + fractionAlongEdge2 = it3->second.fractionAlongEdge1; } else { - *fractionAlongEdge1 = it3->second.fractionAlongEdge1; - *fractionAlongEdge2 = it3->second.fractionAlongEdge2; + fractionAlongEdge1 = it3->second.fractionAlongEdge1; + fractionAlongEdge2 = it3->second.fractionAlongEdge2; } - if ( flipE1 ) *fractionAlongEdge1 = 1 - *fractionAlongEdge1; - if ( flipE2 ) *fractionAlongEdge2 = 1 - *fractionAlongEdge2; + if ( flipE1 ) fractionAlongEdge1 = 1 - fractionAlongEdge1; + if ( flipE2 ) fractionAlongEdge2 = 1 - fractionAlongEdge2; return true; } diff --git a/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp b/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp index 7146e3bff3..7c7b3abc93 100644 --- a/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp +++ b/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp @@ -31,61 +31,10 @@ #include "cvfGeometryTools.h" #include +#include using namespace cvf; -#if 0 -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void ControlVolume::calculateCubeFaceStatus(const cvf::Vec3dArray& nodeCoords, double areaTolerance) -{ - int cubeFace; - cvf::uint cubeFaceIndices[4]; - for (cubeFace = 0; cubeFace < 6; ++cubeFace) - { - surfaceNodeIndices(static_cast(cubeFace), cubeFaceIndices); - - std::vector conns; - connections(static_cast(cubeFace), &conns); - - if (!conns.size()) - { - m_cubeFaceStatus[cubeFace] = FREE_FACE; - } - else - { - double area = 0.5 * (nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[0]] ^ nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[0]]).length(); - area += 0.5 * (nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[2]] ^ nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[2]]).length(); - double totConnectionArea = 0; - size_t i; - for (i = 0; i < conns.size(); ++i) - { - totConnectionArea += conns[i]->brfArea(); - } - - if ( totConnectionArea < area - areaTolerance ) - { - m_cubeFaceStatus[cubeFace] = PARTIALLY_COVERED; - } - else - { - m_cubeFaceStatus[cubeFace] = COMPLETELY_COVERED; - } - } - - // Create a polygon to store the complete polygon of the faces - // not completely covered by connections - // This polygon will be filled with nodes later - - if (m_cubeFaceStatus[cubeFace] != COMPLETELY_COVERED ) - { - m_freeFacePolygons[cubeFace] = new std::list >; - } - } -} -#endif - template NodeType quadNormal( ArrayWrapperConst nodeCoords, const IndexType cubeFaceIndices[4] ) { @@ -159,8 +108,8 @@ TEST( CellFaceIntersectionTst, Intersection1 ) edgeIntersectionStorage.setVertexCount( nodes.size() ); { std::vector polygon; - bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), faces[0].data(), @@ -176,8 +125,8 @@ TEST( CellFaceIntersectionTst, Intersection1 ) { std::vector polygon; - bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), faces[0].data(), @@ -193,8 +142,8 @@ TEST( CellFaceIntersectionTst, Intersection1 ) { std::vector polygon; - bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), faces[0].data(), @@ -214,7 +163,7 @@ TEST( CellFaceIntersectionTst, Intersection1 ) for ( cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx ) { - GeometryTools::insertVertexInPolygon( &basePolygon, wrapArrayConst( &nodes ), vxIdx, 1e-6 ); + GeometryTools::insertVertexInPolygon( basePolygon, wrapArrayConst( &nodes ), vxIdx, 1e-6 ); } EXPECT_EQ( (size_t)8, basePolygon.size() ); @@ -224,7 +173,7 @@ TEST( CellFaceIntersectionTst, Intersection1 ) { for ( cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx ) { - GeometryTools::insertVertexInPolygon( &overlapPolygons[pIdx], wrapArrayConst( &nodes ), vxIdx, 1e-6 ); + GeometryTools::insertVertexInPolygon( overlapPolygons[pIdx], wrapArrayConst( &nodes ), vxIdx, 1e-6 ); } if ( pIdx == 0 ) @@ -325,8 +274,8 @@ TEST( CellFaceIntersectionTst, Intersection ) nodes[6] = cvf::Vec3d( 1, 1, 0 ); nodes[7] = cvf::Vec3d( 0, 1, 0 ); - bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), cv1CubeFaceIndices, @@ -348,8 +297,8 @@ TEST( CellFaceIntersectionTst, Intersection ) nodes[7] = cvf::Vec3d( -0.25, 0.5, 0 ); polygon.clear(); - isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), cv1CubeFaceIndices, @@ -389,8 +338,8 @@ TEST( CellFaceIntersectionTst, FreeFacePolygon ) nodes[6] = cvf::Vec3d( 1, 1, 0 ); nodes[7] = cvf::Vec3d( 0, 1, 0 ); - bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), cv1CubeFaceIndices, @@ -427,8 +376,8 @@ TEST( CellFaceIntersectionTst, FreeFacePolygon ) nodes[7] = cvf::Vec3d( -0.25, 0.5, 0 ); polygon.clear(); - isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon, - &additionalVertices, + isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon, + additionalVertices, &edgeIntersectionStorage, wrapArrayConst( &nodes ), cv1CubeFaceIndices, @@ -596,3 +545,222 @@ TEST( EarClipTesselator, ErrorTest ) // continue; } } + +//-------------------------------------------------------------------------------------------------- +/// Test robust geometry calculation functions +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsErrorHandling, RobustCalculations ) +{ + // Test with degenerate polygons + std::vector degeneratePolygon = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 0, 0, 0 ) }; + + cvf::Vec3d center = GeometryTools::computePolygonCenter( degeneratePolygon ); + EXPECT_TRUE( center.isZero() ); + + double area = GeometryTools::polygonArea( degeneratePolygon ); + EXPECT_DOUBLE_EQ( 0.0, area ); + + // Test with collinear points + std::vector collinearPolygon = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 1, 0, 0 ), cvf::Vec3d( 2, 0, 0 ) }; + + area = GeometryTools::polygonArea( collinearPolygon ); + EXPECT_DOUBLE_EQ( 0.0, area ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test simple error handling cases +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsErrorHandling, BasicErrorHandling ) +{ + // Test empty polygon center calculation + std::vector emptyPolygon; + cvf::Vec3d center = GeometryTools::computePolygonCenter( emptyPolygon ); + EXPECT_TRUE( center.isUndefined() ); + + // Test empty polygon area calculation + double area = GeometryTools::polygonArea( emptyPolygon ); + EXPECT_DOUBLE_EQ( 0.0, area ); + + // Test simple geometric functions with valid inputs + std::vector triangle = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 1, 0, 0 ), cvf::Vec3d( 0, 1, 0 ) }; + + center = GeometryTools::computePolygonCenter( triangle ); + EXPECT_NEAR( 1.0 / 3.0, center.x(), 1e-6 ); + EXPECT_NEAR( 1.0 / 3.0, center.y(), 1e-6 ); + + area = GeometryTools::polygonArea( triangle ); + EXPECT_NEAR( 0.5, area, 1e-6 ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test bounds checking in EdgeIntersectStorage +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsErrorHandling, EdgeIntersectStorageBoundsCheck ) +{ + EdgeIntersectStorage storage; + storage.setVertexCount( 4 ); + + // Test adding intersection with out-of-bounds index + storage.addIntersection( 10, 11, 0, 1, 0, GeometryTools::LINES_CROSSES, 0.5, 0.5 ); + // Should return early without crashing + + cvf::uint intersectionIndex; + GeometryTools::IntersectionStatus status; + double fraction1, fraction2; + + // Test finding intersection with out-of-bounds index + bool found = storage.findIntersection( 10, 11, 0, 1, &intersectionIndex, &status, fraction1, fraction2 ); + EXPECT_FALSE( found ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test computePolygonCenter template function +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsTemplate, ComputePolygonCenter ) +{ + // Test with Vec3d + std::vector polygon = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 2, 0, 0 ), cvf::Vec3d( 2, 2, 0 ), cvf::Vec3d( 0, 2, 0 ) }; + + cvf::Vec3d center = GeometryTools::computePolygonCenter( polygon ); + EXPECT_DOUBLE_EQ( 1.0, center.x() ); + EXPECT_DOUBLE_EQ( 1.0, center.y() ); + EXPECT_DOUBLE_EQ( 0.0, center.z() ); + + // Test with empty polygon + std::vector emptyPolygon; + cvf::Vec3d emptyCenter = GeometryTools::computePolygonCenter( emptyPolygon ); + EXPECT_TRUE( emptyCenter.isUndefined() ); + + // Test with single point + std::vector singlePoint = { cvf::Vec3d( 5, 3, 1 ) }; + cvf::Vec3d singleCenter = GeometryTools::computePolygonCenter( singlePoint ); + EXPECT_DOUBLE_EQ( 5.0, singleCenter.x() ); + EXPECT_DOUBLE_EQ( 3.0, singleCenter.y() ); + EXPECT_DOUBLE_EQ( 1.0, singleCenter.z() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test interpolateQuad template function +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsTemplate, InterpolateQuad ) +{ + cvf::Vec3d v1( 0, 0, 0 ); + cvf::Vec3d v2( 1, 0, 0 ); + cvf::Vec3d v3( 1, 1, 0 ); + cvf::Vec3d v4( 0, 1, 0 ); + + // Test interpolation at corners + cvf::Vec3d center( 0.5, 0.5, 0 ); + double result = GeometryTools::interpolateQuad( v1, 10.0, v2, 20.0, v3, 30.0, v4, 40.0, center ); + EXPECT_DOUBLE_EQ( 25.0, result ); // Average of corner values + + // Test interpolation at v1 position + result = GeometryTools::interpolateQuad( v1, 10.0, v2, 20.0, v3, 30.0, v4, 40.0, v1 ); + EXPECT_DOUBLE_EQ( 10.0, result ); + + // Test with float values + float floatResult = GeometryTools::interpolateQuad( v1, 1.0f, v2, 2.0f, v3, 3.0f, v4, 4.0f, center ); + EXPECT_FLOAT_EQ( 2.5f, floatResult ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test geometric edge cases with simple functions +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsEdgeCases, GeometricEdgeCases ) +{ + // Test polygon with repeated vertices + std::vector repeatedVertices = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 1, 0, 0 ), + cvf::Vec3d( 1, 0, 0 ), // Repeated + cvf::Vec3d( 1, 1, 0 ), + cvf::Vec3d( 0, 1, 0 ) }; + + cvf::Vec3d center = GeometryTools::computePolygonCenter( repeatedVertices ); + EXPECT_GT( center.x(), 0.0 ); + EXPECT_GT( center.y(), 0.0 ); + + // Test very small polygon + std::vector smallPolygon = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 1e-10, 0, 0 ), + cvf::Vec3d( 1e-10, 1e-10, 0 ), + cvf::Vec3d( 0, 1e-10, 0 ) }; + + double area = GeometryTools::polygonArea( smallPolygon ); + EXPECT_GT( area, 0.0 ); + EXPECT_LT( area, 1e-15 ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test additional geometric calculations +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsEdgeCases, AdditionalGeometricTests ) +{ + // Test polygon normal calculation + std::vector square = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 1, 0, 0 ), cvf::Vec3d( 1, 1, 0 ), cvf::Vec3d( 0, 1, 0 ) }; + + cvf::Vec3d areaNormal = GeometryTools::polygonAreaNormal3D( square ); + EXPECT_NEAR( 0.0, areaNormal.x(), 1e-6 ); + EXPECT_NEAR( 0.0, areaNormal.y(), 1e-6 ); + EXPECT_NEAR( 1.0, areaNormal.z(), 1e-6 ); // Positive due to opposite winding + + // Test with different winding + std::vector squareReversed = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 0, 1, 0 ), cvf::Vec3d( 1, 1, 0 ), cvf::Vec3d( 1, 0, 0 ) }; + + areaNormal = GeometryTools::polygonAreaNormal3D( squareReversed ); + EXPECT_NEAR( 0.0, areaNormal.x(), 1e-6 ); + EXPECT_NEAR( 0.0, areaNormal.y(), 1e-6 ); + EXPECT_NEAR( -1.0, areaNormal.z(), 1e-6 ); // Negative due to winding +} + +//-------------------------------------------------------------------------------------------------- +/// Test numerical precision handling +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsEdgeCases, NumericalPrecision ) +{ + // Test with very small differences + std::vector nearIdenticalPoints = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 1, 0, 0 ), + cvf::Vec3d( 1.0000001, 1e-10, 0 ), // Very close to (1, 0, 0) + cvf::Vec3d( 0, 1, 0 ) }; + + double area = GeometryTools::polygonArea( nearIdenticalPoints ); + EXPECT_GT( area, 0.0 ); + EXPECT_LT( area, 1.0 ); // Should be less than unit square + + // Test with extreme coordinates + std::vector extremeCoords = { cvf::Vec3d( 1e6, 1e6, 0 ), + cvf::Vec3d( 1e6 + 1, 1e6, 0 ), + cvf::Vec3d( 1e6 + 1, 1e6 + 1, 0 ), + cvf::Vec3d( 1e6, 1e6 + 1, 0 ) }; + + area = GeometryTools::polygonArea( extremeCoords ); + EXPECT_NEAR( 1.0, area, 1e-6 ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test with various polygon sizes and configurations +//-------------------------------------------------------------------------------------------------- +TEST( GeometryToolsEdgeCases, PolygonVariations ) +{ + // Test with triangle + std::vector triangleNodes = { cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 1, 0, 0 ), cvf::Vec3d( 0.5, 1, 0 ) }; + + cvf::Vec3d triangleCenter = GeometryTools::computePolygonCenter( triangleNodes ); + EXPECT_NEAR( 0.5, triangleCenter.x(), 1e-6 ); + EXPECT_NEAR( 1.0 / 3.0, triangleCenter.y(), 1e-6 ); + + // Test with hexagon + std::vector hexagonNodes; + for ( int i = 0; i < 6; ++i ) + { + double angle = i * 2.0 * M_PI / 6.0; + hexagonNodes.push_back( cvf::Vec3d( cos( angle ), sin( angle ), 0 ) ); + } + + cvf::Vec3d hexagonCenter = GeometryTools::computePolygonCenter( hexagonNodes ); + EXPECT_NEAR( 0.0, hexagonCenter.x(), 1e-6 ); + EXPECT_NEAR( 0.0, hexagonCenter.y(), 1e-6 ); + + double hexagonArea = GeometryTools::polygonArea( hexagonNodes ); + EXPECT_GT( hexagonArea, 0.0 ); +}