OilfieldToolsNet/ResInsight
4
0 Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2026-10-24 22:27:55 -05:00
Code Issues Packages Projects Releases Wiki Activity

Improves geometry tools robustness and accuracy

Browse Source 
* Replace asserts with alternatives
* Refactors GeometryTools for improved robustness
This commit is contained in:
Magne Sjaastad authored and GitHub committed 2025-11-03 16:06:58 +01:00
1 parent fd5bdc96d2
commit 597cfd4d30
6 files changed
+402 -225
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

ApplicationLibCode/ReservoirDataModel/RigCaseCellResultsData.cpp
+2 -2
View File
@@ -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<size_t>*)nullptr,
cvf::wrapArrayConst( &nodes ),
face1.data(),
ApplicationLibCode/ReservoirDataModel/RigCellFaceGeometryTools.cpp
+4 -4
View File
@@ -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<size_t>*)nullptr,
cvf::wrapArrayConst( &mainGrid.nodes() ),
face1.data(),
@@ -335,8 +335,8 @@ void RigCellFaceGeometryTools::extractConnectionsForFace( const RigFault::FaultF
std::array<size_t, 4> candidateFaceIndices;
mainGrid->cell( candidateCellIndex ).faceIndices( candidateFace, &candidateFaceIndices );
bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&intersections,
bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon,
intersections,
(cvf::EdgeIntersectStorage<size_t>*)nullptr,
cvf::wrapArrayConst( &mainGridNodes ),
sourceFaceIndices.data(),
ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp
+55 -68
View File
@@ -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;
ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h
+12 -9
View File
@@ -23,6 +23,9 @@
#include "cvfArrayWrapperConst.h"
#include "cvfMatrix3.h"
#include "cvfMatrix4.h"
#include "cafPdmLogging.h"
#include <list>
#include <map>
@@ -93,7 +96,7 @@ public:
std::vector<cvf::Vec3d>* createdVertexes );
template <typename VerticeArrayType, typename IndexType>
static bool insertVertexInPolygon( std::vector<IndexType>* polygon,
static bool insertVertexInPolygon( std::vector<IndexType>& polygon,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> 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 <typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
@@ -113,12 +116,12 @@ public:
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> vertices,
ArrayWrapperConst<PolygonArrayType, IndexType> indices,
const cvf::Vec3d& point,
int* touchedEdgeIndex,
int& touchedEdgeIndex,
double tolerance = 1e-6 );
template <typename VerticeArrayType, typename IndexType>
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType>* polygon,
std::vector<cvf::Vec3d>* createdVertexes,
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType>& polygon,
std::vector<cvf::Vec3d>& createdVertexes,
EdgeIntersectStorage<IndexType>* edgeIntersectionStorage,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> 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,
ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.inl
+93 -74
View File
@@ -27,6 +27,11 @@ namespace cvf
template <typename Vec3Type>
Vec3Type GeometryTools::computePolygonCenter( const std::vector<Vec3Type>& 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 <typename VerticeArrayType, typename IndexType>
bool GeometryTools::insertVertexInPolygon( std::vector<IndexType>* polygon,
bool GeometryTools::insertVertexInPolygon( std::vector<IndexType>& polygon,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
IndexType vertexIndex,
double tolerance )
{
CVF_ASSERT( polygon );
// Check if vertex is directly included already
for ( typename std::vector<IndexType>::iterator it = polygon->begin(); it != polygon->end(); ++it )
for ( typename std::vector<IndexType>::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<IndexType>::iterator it = polygon->begin(); it != polygon->end(); ++it )
for ( typename std::vector<IndexType>::iterator it = polygon.begin(); it != polygon.end(); ++it )
{
if ( ( nodeCoords[*it] - nodeCoords[vertexIndex] ).length() < tolerance )
{
@@ -93,14 +95,13 @@ bool GeometryTools::insertVertexInPolygon( std::vector<IndexType>*
}
if ( existsOrInserted ) return true;
#endif
// Copy the start polygon to a list
std::list<IndexType> 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<IndexType>*
// Write polygon back into the vector
polygon->clear();
polygon.clear();
for ( typename std::list<IndexType>::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<VerticeArrayType, cvf::Vec3d> vertices,
cvf::ArrayWrapperConst<PolygonArrayType, IndexType> 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<int>( firstIdx );
touchedEdgeIndex = static_cast<int>( firstIdx );
return true;
}
}
@@ -247,20 +246,17 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d&
/// The second face must have opposite winding of the first
//--------------------------------------------------------------------------------------------------
template <typename VerticeArrayType, typename IndexType>
bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector<IndexType>* polygon,
std::vector<cvf::Vec3d>* createdVertexes,
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 )
{
CVF_ASSERT( polygon );
CVF_ASSERT( createdVertexes );
// Topology analysis
IndexType newVertexIndex = static_cast<IndexType>( nodes.size() + createdVertexes->size() );
IndexType newVertexIndex = static_cast<IndexType>( 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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
cv2CubeFaceIndices[nextCv2Idx],
&intersectionVxIndex,
&intersectStatus,
&fractionAlongEdge1,
&fractionAlongEdge2 );
fractionAlongEdge1,
fractionAlongEdge2 );
if ( !found )
{
@@ -441,9 +437,9 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads( std::vector<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
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<IndexType>*
typename std::map<double, IndexType>::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<IndexType>*
// 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<IndexType>*
}
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<IndexType>::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<IndexType>::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<IndexType>::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<IndexType>::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<IndexType, std::map<IndexType, std::map<IndexType, IntersectData>>>::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<IndexType, std::map<IndexType, IntersectData>>::iterator it2;
it2 = it->second.find( e2P1 );
auto it2 = it->second.find( e2P1 );
if ( it2 == it->second.end() ) return false;
typename std::map<IndexType, IntersectData>::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<IndexType>::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;
}
ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp
+236 -68
View File
@@ -31,61 +31,10 @@
#include "cvfGeometryTools.h"
#include <array>
#include <cmath>
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<Defines::CubeFace>(cubeFace), cubeFaceIndices);
std::vector<const brv::Connection*> conns;
connections(static_cast<Defines::CubeFace>(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<std::pair<cvf::uint, bool> >;
}
}
}
#endif
template <typename NodeArrayType, typename NodeType, typename IndexType>
NodeType quadNormal( ArrayWrapperConst<NodeArrayType, NodeType> nodeCoords, const IndexType cubeFaceIndices[4] )
{
@@ -159,8 +108,8 @@ TEST( CellFaceIntersectionTst, Intersection1 )
edgeIntersectionStorage.setVertexCount( nodes.size() );
{
std::vector<cvf::uint> 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<cvf::uint> 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<cvf::uint> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::uint> 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<cvf::Vec3d> 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<cvf::Vec3d> emptyPolygon;
cvf::Vec3d emptyCenter = GeometryTools::computePolygonCenter( emptyPolygon );
EXPECT_TRUE( emptyCenter.isUndefined() );
// Test with single point
std::vector<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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<cvf::Vec3d> 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 );
}