diff --git a/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp b/ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp index 9ba1f43b02..d3218a3d39 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; @@ -177,8 +180,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() std::array cellCorners; std::array cornerIndices; - size_t startOfGeneratedTrianglesForNativeTriangles = outputTriangleVertices.size(); - for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx ) { size_t globalCellIdx = triIntersectedCellCandidates[ticIdx]; @@ -250,28 +251,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays() } } } - - // Add triangles for the part of the native triangle outside any gridcells - - if ( startOfGeneratedTrianglesForNativeTriangles == 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 - { - // Todo: - // Subtract the created triangles from the native triangle - // Add the remains - } } 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() ) 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 { diff --git a/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp index 1f7bfa5175..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 @@ -880,12 +905,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 +934,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 +944,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 +986,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 +1022,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; 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; + } +}