From 1949940b9648121587baa531c2c6386e295ad28b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 09:44:28 +0100 Subject: [PATCH 1/6] Fix misleading parameter name --- ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h index a18a924218..3b165a1bfa 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h @@ -53,7 +53,7 @@ public: static std::vector> subtractPolygons( const std::vector& sourcePolygon, - const std::vector>& polygonToSubtract ); + const std::vector>& polygonsToSubtract ); enum ZInterpolationType { From cb664d9820b746522888e0f12f600c477cc0088e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 11:14:37 +0100 Subject: [PATCH 2/6] Improved EarClipTesselator. More resillient against zero area ears. Less prone to make invalid traingles due to vertices touching edge but might fail to produce triangles at all in that case. --- .../ReservoirDataModel/cvfGeometryTools.cpp | 41 +++++++++++++------ .../ReservoirDataModel/cvfGeometryTools.h | 13 ++++-- 2 files changed, 38 insertions(+), 16 deletions(-) diff --git a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp index 1f7bfa5175..e5188e1cf9 100644 --- a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp @@ -880,12 +880,17 @@ bool EarClipTesselator::calculateTriangles( std::vector* triangleIndices ++w; // v + 1; if ( w == m_polygonIndices.end() ) w = m_polygonIndices.begin(); // if (nv <= w) w = 0; - if ( isTriangleValid( u, v, w ) ) + EarClipTesselator::TriangleStatus triStatus = calculateTriangleStatus( u, v, w ); + + if ( triStatus != INVALID_TRIANGLE ) { // Indices of the vertices - triangleIndices->push_back( *u ); - triangleIndices->push_back( *v ); - triangleIndices->push_back( *w ); + if ( triStatus == VALID_TRIANGLE ) + { + triangleIndices->push_back( *u ); + triangleIndices->push_back( *v ); + triangleIndices->push_back( *w ); + } // Remove v from remaining polygon m_polygonIndices.erase( v ); @@ -904,9 +909,9 @@ bool EarClipTesselator::calculateTriangles( std::vector* triangleIndices /// Is this a valid triangle ? ( No points inside, and points not on a line. ) //-------------------------------------------------------------------------------------------------- -bool EarClipTesselator::isTriangleValid( std::list::const_iterator u, - std::list::const_iterator v, - std::list::const_iterator w ) const +EarClipTesselator::TriangleStatus EarClipTesselator::calculateTriangleStatus( std::list::const_iterator u, + std::list::const_iterator v, + std::list::const_iterator w ) const { CVF_ASSERT( m_X > -1 && m_Y > -1 ); @@ -914,9 +919,17 @@ bool EarClipTesselator::isTriangleValid( std::list::const_iterator u, cvf::Vec3d B = ( *m_nodeCoords )[*v]; cvf::Vec3d C = ( *m_nodeCoords )[*w]; - if ( m_areaTolerance > - ( ( ( B[m_X] - A[m_X] ) * ( C[m_Y] - A[m_Y] ) ) - ( ( B[m_Y] - A[m_Y] ) * ( C[m_X] - A[m_X] ) ) ) ) - return false; + double mainAxisProjectedArea = ( ( B[m_X] - A[m_X] ) * ( C[m_Y] - A[m_Y] ) ) - + ( ( B[m_Y] - A[m_Y] ) * ( C[m_X] - A[m_X] ) ); + if ( -m_areaTolerance > mainAxisProjectedArea ) + { + // Definite negative triangle + return INVALID_TRIANGLE; + } + else if ( fabs( mainAxisProjectedArea ) < m_areaTolerance ) + { + return NEAR_ZERO_AREA_TRIANGLE; + } std::list::const_iterator c; std::list::const_iterator outside; @@ -948,15 +961,16 @@ bool EarClipTesselator::isTriangleValid( std::list::const_iterator u, cvf::Vec3d P = ( *m_nodeCoords )[*c]; - if ( isPointInsideTriangle( A, B, C, P ) ) return false; + if ( isPointInsideTriangle( A, B, C, P ) ) return INVALID_TRIANGLE; } - return true; + return VALID_TRIANGLE; } //-------------------------------------------------------------------------------------------------- /// Decides if a point P is inside of the triangle defined by A, B, C. /// By calculating the "double area" (cross product) of Corner to corner x Corner to point vectors +/// If cross product is negative, the point is outside //-------------------------------------------------------------------------------------------------- bool EarClipTesselator::isPointInsideTriangle( const cvf::Vec3d& A, @@ -983,7 +997,8 @@ bool EarClipTesselator::isPointInsideTriangle( const cvf::Vec3d& A, double aCROSSbp = ax * bpy - ay * bpx; double cCROSSap = cx * apy - cy * apx; double bCROSScp = bx * cpy - by * cpx; - double tol = 0; + double tol = -m_areaTolerance; + return ( ( aCROSSbp >= tol ) && ( bCROSScp >= tol ) && ( cCROSSap >= tol ) ); }; diff --git a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.h b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.h index 2405689f3b..ee6b083095 100644 --- a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.h @@ -209,9 +209,16 @@ public: virtual bool calculateTriangles( std::vector* triangles ); protected: - bool isTriangleValid( std::list::const_iterator u, - std::list::const_iterator v, - std::list::const_iterator w ) const; + enum TriangleStatus + { + INVALID_TRIANGLE, + NEAR_ZERO_AREA_TRIANGLE, + VALID_TRIANGLE + }; + TriangleStatus calculateTriangleStatus( std::list::const_iterator u, + std::list::const_iterator v, + std::list::const_iterator w ) const; + bool isPointInsideTriangle( const cvf::Vec3d& A, const cvf::Vec3d& B, const cvf::Vec3d& C, const cvf::Vec3d& P ) const; double calculateProjectedPolygonArea() const; From 896a0f067eaef06790e6b73220041f7831972a5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 11:17:04 +0100 Subject: [PATCH 3/6] #5369 First sensible border geometry. Still not handling holes returned from clipper. Int-point size (for clipper) implications not yet understood/investigated --- ...ivSurfaceIntersectionGeometryGenerator.cpp | 256 +++++++++++++++++- 1 file changed, 253 insertions(+), 3 deletions(-) diff --git a/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp b/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp index 9ba1f43b02..b202442b47 100644 --- a/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp +++ b/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp @@ -47,6 +47,9 @@ #include "../cafHexInterpolator/cafHexInterpolator.h" #include "RivSectionFlattner.h" +#include "RiaLogging.h" +#include "clipper.hpp" + cvf::ref displayCoordTransform( const RimIntersection* intersection ) { Rim3dView* rimView = nullptr; @@ -129,6 +132,114 @@ private: cvf::cref m_hexGrid; }; +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +template +double closestAxisSignedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, const std::vector& polygon ) +{ + int Z = cvf::GeometryTools::findClosestAxis( planeNormal ); + int X = ( Z + 1 ) % 3; + int Y = ( Z + 2 ) % 3; + + // Use Shoelace formula to calculate signed area. + // https://en.wikipedia.org/wiki/Shoelace_formula + double signedArea = 0.0; + for ( size_t i = 0; i < polygon.size(); ++i ) + { + signedArea += ( polygon[( i + 1 ) % polygon.size()][X] + polygon[i][X] ) * + ( polygon[( i + 1 ) % polygon.size()][Y] - polygon[i][Y] ); + } + + return ( planeNormal[Z] > 0 ) ? signedArea : -signedArea; +} + +class ClipperInterface +{ +public: + static ClipperLib::IntPoint toClipperPoint( const cvf::Vec3d& cvfPoint ) + { + int xInt = cvfPoint.x() * clipperConversionFactor; + int yInt = cvfPoint.y() * clipperConversionFactor; + return ClipperLib::IntPoint( xInt, yInt, 0 ); + } + + static cvf::Vec3d fromClipperPoint( const ClipperLib::IntPoint& clipPoint ) + { + return cvf::Vec3d( clipPoint.X, clipPoint.Y, 0.0 ) / clipperConversionFactor; + } + //-------------------------------------------------------------------------------------------------- + /// + //-------------------------------------------------------------------------------------------------- + static std::vector> + subtractAndSimplifyPolygons( const std::vector& sourcePolygon, + const std::vector>& polygonsToSubtract ) + { + ClipperLib::Paths unionOfPolygonsToSubtract; + { + ClipperLib::Clipper unionator; + + for ( const auto& path : polygonsToSubtract ) + { + ClipperLib::Path polyToSubtractPath; + for ( const auto& v : path ) + { + polyToSubtractPath.push_back( toClipperPoint( v ) ); + } + + unionator.AddPath( polyToSubtractPath, ClipperLib::ptSubject, true ); + } + + unionator.Execute( ClipperLib::ctUnion, + unionOfPolygonsToSubtract, + ClipperLib::pftEvenOdd, + ClipperLib::pftEvenOdd ); + } + + ClipperLib::Path sourcePolygonPath; + + for ( const auto& v : sourcePolygon ) + { + sourcePolygonPath.push_back( toClipperPoint( v ) ); + } + + ClipperLib::Clipper subtractor; + + subtractor.AddPath( sourcePolygonPath, ClipperLib::ptSubject, true ); + subtractor.AddPaths( unionOfPolygonsToSubtract, ClipperLib::ptClip, true ); + + ClipperLib::Paths subtractionResultPaths; + subtractor.Execute( ClipperLib::ctDifference, + subtractionResultPaths, + ClipperLib::pftEvenOdd, + ClipperLib::pftEvenOdd ); + + ClipperLib::CleanPolygons( subtractionResultPaths, 3 ); + + std::vector> clippedPolygons; + + // Convert back to std::vector > + + for ( ClipperLib::Path pathInSol : subtractionResultPaths ) + { + std::vector clippedPolygon; + for ( ClipperLib::IntPoint IntPosition : pathInSol ) + { + clippedPolygon.push_back( fromClipperPoint( IntPosition ) ); + } + + clippedPolygons.push_back( clippedPolygon ); + } + + return clippedPolygons; + } + +private: + static double clipperConversionFactor; +}; + +double ClipperInterface::clipperConversionFactor = 100; // For transform to clipper int + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -177,7 +288,7 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() std::array cellCorners; std::array cornerIndices; - size_t startOfGeneratedTrianglesForNativeTriangles = outputTriangleVertices.size(); + size_t startOfGeneratedTrianglesForNativeTriangle = outputTriangleVertices.size(); for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx ) { @@ -253,7 +364,7 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() // Add triangles for the part of the native triangle outside any gridcells - if ( startOfGeneratedTrianglesForNativeTriangles == outputTriangleVertices.size() ) + if ( startOfGeneratedTrianglesForNativeTriangle == outputTriangleVertices.size() ) { // No triangles created, use the complete native triangle outputTriangleVertices.push_back( cvf::Vec3f( p0 - displayModelOffset ) ); @@ -268,9 +379,148 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() } else { - // Todo: + // Use area to check if the native triangle was completely covered by the intersection triangles + double nativeTriangleArea = 0.0; + double intersectionArea = 0.0; + double minSignificantTriangleArea = 0.0; + { + std::vector nativeTriangle = {cvf::Vec3f( p0 - displayModelOffset ), + cvf::Vec3f( p1 - displayModelOffset ), + cvf::Vec3f( p2 - displayModelOffset )}; + + nativeTriangleArea = closestAxisSignedAreaPlanarPolygon( plane.normal(), nativeTriangle ); + minSignificantTriangleArea = 1e-4 * nativeTriangleArea; + + std::vector intersectionTriangle( 3 ); + + for ( size_t tvxIdx = startOfGeneratedTrianglesForNativeTriangle; tvxIdx < outputTriangleVertices.size(); + tvxIdx += 3 ) + { + std::copy( outputTriangleVertices.begin() + tvxIdx, + outputTriangleVertices.begin() + tvxIdx + 3, + intersectionTriangle.begin() ); + intersectionArea += closestAxisSignedAreaPlanarPolygon( plane.normal(), intersectionTriangle ); + } + + // If we have covered enough, do not try to create triangles for the rest + + if ( ( nativeTriangleArea - intersectionArea ) < minSignificantTriangleArea ) continue; + } + // Subtract the created triangles from the native triangle + // We need to transform the triangles into x/y plane to do the polygon operation + + // Find the local CS for the native triangle in display coords + + cvf::Mat4d nativeTriangleCS; + { + cvf::Vec3d ez = plane.normal().getNormalized(); + cvf::Vec3d ex = ( p1 - p0 ).getNormalized(); + cvf::Vec3d ey = ez ^ ex; + nativeTriangleCS = cvf::Mat4d::fromCoordSystemAxes( &ex, &ey, &ez ); + nativeTriangleCS.setTranslation( p0 - displayModelOffset ); + } + + cvf::Mat4d invNativeTriangleCS = nativeTriangleCS.getInverted(); + + std::vector> polygonsToSubtract; + for ( size_t tvxIdx = startOfGeneratedTrianglesForNativeTriangle; tvxIdx < outputTriangleVertices.size(); + tvxIdx += 3 ) + { + std::vector triangle = + {cvf::Vec3d( outputTriangleVertices[tvxIdx + 0] ).getTransformedPoint( invNativeTriangleCS ), + cvf::Vec3d( outputTriangleVertices[tvxIdx + 1] ).getTransformedPoint( invNativeTriangleCS ), + cvf::Vec3d( outputTriangleVertices[tvxIdx + 2] ).getTransformedPoint( invNativeTriangleCS )}; + polygonsToSubtract.push_back( triangle ); + } + + std::vector nativeTrianglePoly = + {( cvf::Vec3d( cvf::Vec3f( p0 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS ), + ( cvf::Vec3d( cvf::Vec3f( p1 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS ), + ( cvf::Vec3d( cvf::Vec3f( p2 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS )}; + + std::vector> remainingPolygons; + + remainingPolygons = ClipperInterface::subtractAndSimplifyPolygons( nativeTrianglePoly, polygonsToSubtract ); + + // Check for holes in solution + bool hasHoles = false; + + for ( const auto& remainingPolygon : remainingPolygons ) + { + double area = closestAxisSignedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, remainingPolygon ); + if ( area < -minSignificantTriangleArea ) + { + hasHoles = true; + } + } + + if ( hasHoles ) continue; // Cant tesselate polygons with holes + // Add the remains + + for ( const auto& remainingPolygon : remainingPolygons ) + { + if ( remainingPolygon.empty() ) continue; + + cvf::EarClipTesselator tess; + tess.setNormal( plane.normal() ); + tess.setMinTriangleArea( minSignificantTriangleArea ); + cvf::Vec3dArray cvfNodes( remainingPolygon ); + tess.setGlobalNodeArray( cvfNodes ); + + std::vector polyIndexes; + for ( size_t idx = 0; idx < remainingPolygon.size(); ++idx ) + { + polyIndexes.push_back( idx ); + } + tess.setPolygonIndices( polyIndexes ); + + std::vector triangleIndices; + bool isTesselationOk = tess.calculateTriangles( &triangleIndices ); + + if ( !isTesselationOk ) + { + // continue; + // CVF_ASSERT( false ); + } + + double tesselatedArea = 0; + + for ( size_t idx = 0; idx < triangleIndices.size(); idx += 3 ) + { + cvf::Vec3f tp1( + ( remainingPolygon[triangleIndices[idx + 0]] ).getTransformedPoint( nativeTriangleCS ) ); + cvf::Vec3f tp2( + ( remainingPolygon[triangleIndices[idx + 1]] ).getTransformedPoint( nativeTriangleCS ) ); + cvf::Vec3f tp3( + ( remainingPolygon[triangleIndices[idx + 2]] ).getTransformedPoint( nativeTriangleCS ) ); + + outputTriangleVertices.push_back( tp1 ); + outputTriangleVertices.push_back( tp2 ); + outputTriangleVertices.push_back( tp3 ); + + std::vector nativeTriangle = {tp1, tp2, tp3}; + + tesselatedArea += closestAxisSignedAreaPlanarPolygon( plane.normal(), nativeTriangle ); + + m_triangleToCellIdxMap.push_back( cvf::UNDEFINED_SIZE_T ); + + m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); + m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); + m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); + } + + if ( ( tesselatedArea - 20 * minSignificantTriangleArea ) > ( nativeTriangleArea - intersectionArea ) ) + { + double overlapArea = tesselatedArea - ( nativeTriangleArea - intersectionArea ); + + RiaLogging::debug( "Surface intersection triangularization overlap detected : " + + QString::number( overlapArea ) ); + + // CVF_ASSERT( false ); + } + } } } From 091156dec1449e04a7e27ccd0ffe9aa97a10ac6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 15:47:48 +0100 Subject: [PATCH 4/6] Add a unit test start frame for EarClipTesselator --- .../UnitTests/cvfGeometryTools-Test.cpp | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/ApplicationCode/UnitTests/cvfGeometryTools-Test.cpp b/ApplicationCode/UnitTests/cvfGeometryTools-Test.cpp index 38b21f2f09..52b29f8146 100644 --- a/ApplicationCode/UnitTests/cvfGeometryTools-Test.cpp +++ b/ApplicationCode/UnitTests/cvfGeometryTools-Test.cpp @@ -531,3 +531,44 @@ TEST( CellFaceIntersectionTst, PolygonAreaNormal3D ) EXPECT_DOUBLE_EQ( 0.0, area.z() ); } } + +TEST( EarClipTesselator, ErrorTest ) +{ + std::vector remainingPolygon{ + cvf::Vec3d( 44.66, 20.17, 0 ), + cvf::Vec3d( 78.08, 35.26, 0 ), + cvf::Vec3d( 93.97, 35.83, 0 ), + cvf::Vec3d( 144.95, 44.42, 0 ), + cvf::Vec3d( 172.59, 39.73, 0 ), + cvf::Vec3d( 227.27, 24.01, 0 ), + cvf::Vec3d( 217.46, 45.72, 0 ), + cvf::Vec3d( 178.5, 57.61, 0 ), + cvf::Vec3d( 141.33, 63.82, 0 ), + cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 63.77, 0, 0 ), + }; + + double nativeTriangleArea = 21266; + + cvf::EarClipTesselator tess; + tess.setNormal( cvf::Vec3d::Z_AXIS ); + tess.setMinTriangleArea( 1e-3 * nativeTriangleArea ); + cvf::Vec3dArray cvfNodes( remainingPolygon ); + tess.setGlobalNodeArray( cvfNodes ); + + std::vector polyIndexes; + for ( size_t idx = 0; idx < remainingPolygon.size(); ++idx ) + { + polyIndexes.push_back( idx ); + } + tess.setPolygonIndices( polyIndexes ); + + std::vector triangleIndices; + bool isTesselationOk = tess.calculateTriangles( &triangleIndices ); + + // CVF_ASSERT( isTesselationOk ); + if ( !isTesselationOk ) + { + // continue; + } +} From e021c4635ce839607e7734bdd0a58864ab40da26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 15:57:14 +0100 Subject: [PATCH 5/6] Add a more correct version of signedAreaPlanarPolygon as reference for later. --- .../ReservoirDataModel/cvfGeometryTools.cpp | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp index e5188e1cf9..3161ab3e2e 100644 --- a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp @@ -175,6 +175,31 @@ double GeometryTools::signedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, co return signedArea; } +//-------------------------------------------------------------------------------------------------- +/// This method below is more correct than the one above, both in naming and behaviour. +/// Should be replaced, but is not done now to avoid possible sideeffects +/// The difference is the sign of the area. The one below retuns correct sign according to the plane normal +/// provided +//-------------------------------------------------------------------------------------------------- +template +double closestAxisSignedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, const std::vector& polygon ) +{ + int Z = cvf::GeometryTools::findClosestAxis( planeNormal ); + int X = ( Z + 1 ) % 3; + int Y = ( Z + 2 ) % 3; + + // Use Shoelace formula to calculate signed area. + // https://en.wikipedia.org/wiki/Shoelace_formula + double signedArea = 0.0; + for ( size_t i = 0; i < polygon.size(); ++i ) + { + signedArea += ( polygon[( i + 1 ) % polygon.size()][X] + polygon[i][X] ) * + ( polygon[( i + 1 ) % polygon.size()][Y] - polygon[i][Y] ); + } + + return ( planeNormal[Z] > 0 ) ? signedArea : -signedArea; +} + /* Determine the intersection point of two line segments From Paul Bourke, but modified to really handle coincident lines From de31cc32709977ee64986cab5cc7a73f6b63b9ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 24 Jan 2020 15:59:18 +0100 Subject: [PATCH 6/6] #5369 Final solution for surface intersection border : Use polygon offsett --- ...ivSurfaceIntersectionGeometryGenerator.cpp | 271 ------------------ .../Surfaces/RivSurfacePartMgr.cpp | 9 +- 2 files changed, 7 insertions(+), 273 deletions(-) diff --git a/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp b/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp index b202442b47..d3218a3d39 100644 --- a/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp +++ b/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp @@ -132,114 +132,6 @@ private: cvf::cref m_hexGrid; }; -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -template -double closestAxisSignedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, const std::vector& polygon ) -{ - int Z = cvf::GeometryTools::findClosestAxis( planeNormal ); - int X = ( Z + 1 ) % 3; - int Y = ( Z + 2 ) % 3; - - // Use Shoelace formula to calculate signed area. - // https://en.wikipedia.org/wiki/Shoelace_formula - double signedArea = 0.0; - for ( size_t i = 0; i < polygon.size(); ++i ) - { - signedArea += ( polygon[( i + 1 ) % polygon.size()][X] + polygon[i][X] ) * - ( polygon[( i + 1 ) % polygon.size()][Y] - polygon[i][Y] ); - } - - return ( planeNormal[Z] > 0 ) ? signedArea : -signedArea; -} - -class ClipperInterface -{ -public: - static ClipperLib::IntPoint toClipperPoint( const cvf::Vec3d& cvfPoint ) - { - int xInt = cvfPoint.x() * clipperConversionFactor; - int yInt = cvfPoint.y() * clipperConversionFactor; - return ClipperLib::IntPoint( xInt, yInt, 0 ); - } - - static cvf::Vec3d fromClipperPoint( const ClipperLib::IntPoint& clipPoint ) - { - return cvf::Vec3d( clipPoint.X, clipPoint.Y, 0.0 ) / clipperConversionFactor; - } - //-------------------------------------------------------------------------------------------------- - /// - //-------------------------------------------------------------------------------------------------- - static std::vector> - subtractAndSimplifyPolygons( const std::vector& sourcePolygon, - const std::vector>& polygonsToSubtract ) - { - ClipperLib::Paths unionOfPolygonsToSubtract; - { - ClipperLib::Clipper unionator; - - for ( const auto& path : polygonsToSubtract ) - { - ClipperLib::Path polyToSubtractPath; - for ( const auto& v : path ) - { - polyToSubtractPath.push_back( toClipperPoint( v ) ); - } - - unionator.AddPath( polyToSubtractPath, ClipperLib::ptSubject, true ); - } - - unionator.Execute( ClipperLib::ctUnion, - unionOfPolygonsToSubtract, - ClipperLib::pftEvenOdd, - ClipperLib::pftEvenOdd ); - } - - ClipperLib::Path sourcePolygonPath; - - for ( const auto& v : sourcePolygon ) - { - sourcePolygonPath.push_back( toClipperPoint( v ) ); - } - - ClipperLib::Clipper subtractor; - - subtractor.AddPath( sourcePolygonPath, ClipperLib::ptSubject, true ); - subtractor.AddPaths( unionOfPolygonsToSubtract, ClipperLib::ptClip, true ); - - ClipperLib::Paths subtractionResultPaths; - subtractor.Execute( ClipperLib::ctDifference, - subtractionResultPaths, - ClipperLib::pftEvenOdd, - ClipperLib::pftEvenOdd ); - - ClipperLib::CleanPolygons( subtractionResultPaths, 3 ); - - std::vector> clippedPolygons; - - // Convert back to std::vector > - - for ( ClipperLib::Path pathInSol : subtractionResultPaths ) - { - std::vector clippedPolygon; - for ( ClipperLib::IntPoint IntPosition : pathInSol ) - { - clippedPolygon.push_back( fromClipperPoint( IntPosition ) ); - } - - clippedPolygons.push_back( clippedPolygon ); - } - - return clippedPolygons; - } - -private: - static double clipperConversionFactor; -}; - -double ClipperInterface::clipperConversionFactor = 100; // For transform to clipper int - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -288,8 +180,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() std::array cellCorners; std::array cornerIndices; - size_t startOfGeneratedTrianglesForNativeTriangle = outputTriangleVertices.size(); - for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx ) { size_t globalCellIdx = triIntersectedCellCandidates[ticIdx]; @@ -361,167 +251,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() } } } - - // Add triangles for the part of the native triangle outside any gridcells - - if ( startOfGeneratedTrianglesForNativeTriangle == outputTriangleVertices.size() ) - { - // No triangles created, use the complete native triangle - outputTriangleVertices.push_back( cvf::Vec3f( p0 - displayModelOffset ) ); - outputTriangleVertices.push_back( cvf::Vec3f( p1 - displayModelOffset ) ); - outputTriangleVertices.push_back( cvf::Vec3f( p2 - displayModelOffset ) ); - - m_triangleToCellIdxMap.push_back( cvf::UNDEFINED_SIZE_T ); - - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - } - else - { - // Use area to check if the native triangle was completely covered by the intersection triangles - double nativeTriangleArea = 0.0; - double intersectionArea = 0.0; - double minSignificantTriangleArea = 0.0; - { - std::vector nativeTriangle = {cvf::Vec3f( p0 - displayModelOffset ), - cvf::Vec3f( p1 - displayModelOffset ), - cvf::Vec3f( p2 - displayModelOffset )}; - - nativeTriangleArea = closestAxisSignedAreaPlanarPolygon( plane.normal(), nativeTriangle ); - minSignificantTriangleArea = 1e-4 * nativeTriangleArea; - - std::vector intersectionTriangle( 3 ); - - for ( size_t tvxIdx = startOfGeneratedTrianglesForNativeTriangle; tvxIdx < outputTriangleVertices.size(); - tvxIdx += 3 ) - { - std::copy( outputTriangleVertices.begin() + tvxIdx, - outputTriangleVertices.begin() + tvxIdx + 3, - intersectionTriangle.begin() ); - intersectionArea += closestAxisSignedAreaPlanarPolygon( plane.normal(), intersectionTriangle ); - } - - // If we have covered enough, do not try to create triangles for the rest - - if ( ( nativeTriangleArea - intersectionArea ) < minSignificantTriangleArea ) continue; - } - - // Subtract the created triangles from the native triangle - // We need to transform the triangles into x/y plane to do the polygon operation - - // Find the local CS for the native triangle in display coords - - cvf::Mat4d nativeTriangleCS; - { - cvf::Vec3d ez = plane.normal().getNormalized(); - cvf::Vec3d ex = ( p1 - p0 ).getNormalized(); - cvf::Vec3d ey = ez ^ ex; - nativeTriangleCS = cvf::Mat4d::fromCoordSystemAxes( &ex, &ey, &ez ); - nativeTriangleCS.setTranslation( p0 - displayModelOffset ); - } - - cvf::Mat4d invNativeTriangleCS = nativeTriangleCS.getInverted(); - - std::vector> polygonsToSubtract; - for ( size_t tvxIdx = startOfGeneratedTrianglesForNativeTriangle; tvxIdx < outputTriangleVertices.size(); - tvxIdx += 3 ) - { - std::vector triangle = - {cvf::Vec3d( outputTriangleVertices[tvxIdx + 0] ).getTransformedPoint( invNativeTriangleCS ), - cvf::Vec3d( outputTriangleVertices[tvxIdx + 1] ).getTransformedPoint( invNativeTriangleCS ), - cvf::Vec3d( outputTriangleVertices[tvxIdx + 2] ).getTransformedPoint( invNativeTriangleCS )}; - polygonsToSubtract.push_back( triangle ); - } - - std::vector nativeTrianglePoly = - {( cvf::Vec3d( cvf::Vec3f( p0 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS ), - ( cvf::Vec3d( cvf::Vec3f( p1 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS ), - ( cvf::Vec3d( cvf::Vec3f( p2 - displayModelOffset ) ) ).getTransformedPoint( invNativeTriangleCS )}; - - std::vector> remainingPolygons; - - remainingPolygons = ClipperInterface::subtractAndSimplifyPolygons( nativeTrianglePoly, polygonsToSubtract ); - - // Check for holes in solution - bool hasHoles = false; - - for ( const auto& remainingPolygon : remainingPolygons ) - { - double area = closestAxisSignedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, remainingPolygon ); - if ( area < -minSignificantTriangleArea ) - { - hasHoles = true; - } - } - - if ( hasHoles ) continue; // Cant tesselate polygons with holes - - // Add the remains - - for ( const auto& remainingPolygon : remainingPolygons ) - { - if ( remainingPolygon.empty() ) continue; - - cvf::EarClipTesselator tess; - tess.setNormal( plane.normal() ); - tess.setMinTriangleArea( minSignificantTriangleArea ); - cvf::Vec3dArray cvfNodes( remainingPolygon ); - tess.setGlobalNodeArray( cvfNodes ); - - std::vector polyIndexes; - for ( size_t idx = 0; idx < remainingPolygon.size(); ++idx ) - { - polyIndexes.push_back( idx ); - } - tess.setPolygonIndices( polyIndexes ); - - std::vector triangleIndices; - bool isTesselationOk = tess.calculateTriangles( &triangleIndices ); - - if ( !isTesselationOk ) - { - // continue; - // CVF_ASSERT( false ); - } - - double tesselatedArea = 0; - - for ( size_t idx = 0; idx < triangleIndices.size(); idx += 3 ) - { - cvf::Vec3f tp1( - ( remainingPolygon[triangleIndices[idx + 0]] ).getTransformedPoint( nativeTriangleCS ) ); - cvf::Vec3f tp2( - ( remainingPolygon[triangleIndices[idx + 1]] ).getTransformedPoint( nativeTriangleCS ) ); - cvf::Vec3f tp3( - ( remainingPolygon[triangleIndices[idx + 2]] ).getTransformedPoint( nativeTriangleCS ) ); - - outputTriangleVertices.push_back( tp1 ); - outputTriangleVertices.push_back( tp2 ); - outputTriangleVertices.push_back( tp3 ); - - std::vector nativeTriangle = {tp1, tp2, tp3}; - - tesselatedArea += closestAxisSignedAreaPlanarPolygon( plane.normal(), nativeTriangle ); - - m_triangleToCellIdxMap.push_back( cvf::UNDEFINED_SIZE_T ); - - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() ); - } - - if ( ( tesselatedArea - 20 * minSignificantTriangleArea ) > ( nativeTriangleArea - intersectionArea ) ) - { - double overlapArea = tesselatedArea - ( nativeTriangleArea - intersectionArea ); - - RiaLogging::debug( "Surface intersection triangularization overlap detected : " + - QString::number( overlapArea ) ); - - // CVF_ASSERT( false ); - } - } - } } m_triangleVxes->assign( outputTriangleVertices ); diff --git a/ApplicationCode/ModelVisualization/Surfaces/RivSurfacePartMgr.cpp b/ApplicationCode/ModelVisualization/Surfaces/RivSurfacePartMgr.cpp index 89911c7b7d..b88ca2a5e0 100644 --- a/ApplicationCode/ModelVisualization/Surfaces/RivSurfacePartMgr.cpp +++ b/ApplicationCode/ModelVisualization/Surfaces/RivSurfacePartMgr.cpp @@ -111,7 +111,7 @@ void RivSurfacePartMgr::updateCellResultColor( size_t timeStepIndex ) m_intersectionFacesTextureCoords.p() ); } - if ( m_nativeTrianglesPart.notNull() ) + if ( false && m_nativeTrianglesPart.notNull() ) { if ( !m_nativeVertexToCellIndexMap.size() ) { @@ -230,6 +230,8 @@ void RivSurfacePartMgr::appendIntersectionGeometryPartsToModel( cvf::ModelBasicL m_intersectionFaultGridLines->setTransform( scaleTransform ); model->addPart( m_intersectionFaultGridLines.p() ); } + + appendNativeGeometryPartsToModel( model, scaleTransform ); } //-------------------------------------------------------------------------------------------------- @@ -241,9 +243,12 @@ void RivSurfacePartMgr::applySingleColor() caf::SurfaceEffectGenerator surfaceGen( cvf::Color4f( m_surfaceInView->surface()->color() ), caf::PO_1 ); cvf::ref eff = surfaceGen.generateCachedEffect(); + caf::SurfaceEffectGenerator surfaceGenBehind( cvf::Color4f( m_surfaceInView->surface()->color() ), caf::PO_2 ); + cvf::ref effBehind = surfaceGenBehind.generateCachedEffect(); + if ( m_nativeTrianglesPart.notNull() ) { - m_nativeTrianglesPart->setEffect( eff.p() ); + m_nativeTrianglesPart->setEffect( effBehind.p() ); } if ( m_intersectionFaces.notNull() )