From 8db29e0931ccd3dbb3ae7c815f5a8b55d31615b8 Mon Sep 17 00:00:00 2001 From: jonjenssen <69144954+jonjenssen@users.noreply.github.com> Date: Fri, 12 Jan 2024 09:52:29 +0100 Subject: [PATCH] Add polyline filter support (#11048) * Support selecting cells below polyline, both eclipse and geomech cases * Add some polyline intersection helpers, with tests --- .../CellFilters/RimPolygonFilter.cpp | 79 ++++++++++--- .../CellFilters/RimPolygonFilter.h | 1 + .../RigCellGeometryTools.cpp | 105 +++++++++++++++++- .../ReservoirDataModel/RigCellGeometryTools.h | 7 ++ .../UnitTests/RigCellGeometryTools-Test.cpp | 49 ++++++++ 5 files changed, 224 insertions(+), 17 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp index 1cc0af5fb7..7267dff235 100644 --- a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp +++ b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp @@ -141,6 +141,7 @@ RimPolygonFilter::RimPolygonFilter() CAF_PDM_InitField( &m_showLines, "ShowLines", true, "Show Lines" ); CAF_PDM_InitField( &m_showSpheres, "ShowSpheres", false, "Show Spheres" ); + CAF_PDM_InitField( &m_closePolygon, "ClosePolygon", true, "Closed Polygon" ); CAF_PDM_InitField( &m_lineThickness, "LineThickness", 3, "Line Thickness" ); CAF_PDM_InitField( &m_sphereRadiusFactor, "SphereRadiusFactor", 0.15, "Sphere Radius Factor" ); @@ -359,18 +360,19 @@ void RimPolygonFilter::defineUiOrdering( QString uiConfigName, caf::PdmUiOrderin auto group = uiOrdering.addNewGroup( "General" ); group->add( &m_filterMode ); group->add( &m_enableFiltering ); - group->add( &m_showLines ); - group->add( &m_showSpheres ); + group->add( &m_closePolygon ); auto group1 = uiOrdering.addNewGroup( "Polygon Selection" ); group1->add( &m_polyFilterMode ); - group1->add( &m_polyIncludeType ); + if ( m_closePolygon() ) group1->add( &m_polyIncludeType ); group1->add( &m_targets ); group1->add( &m_enablePicking ); m_polyIncludeType.uiCapability()->setUiName( "Cells to " + modeString() ); auto group2 = uiOrdering.addNewGroup( "Appearance" ); + group2->add( &m_showLines ); + group2->add( &m_showSpheres ); if ( m_showLines ) { group2->add( &m_lineThickness ); @@ -399,6 +401,16 @@ void RimPolygonFilter::defineUiOrdering( QString uiConfigName, caf::PdmUiOrderin { objField->uiCapability()->setUiReadOnly( readOnlyState ); } + + if ( !m_closePolygon() ) + { + m_polyFilterMode = RimPolygonFilter::PolygonFilterModeType::INDEX_K; + m_polyFilterMode.uiCapability()->setUiReadOnly( true ); + } + else + { + m_polyFilterMode.uiCapability()->setUiReadOnly( readOnlyState ); + } } //-------------------------------------------------------------------------------------------------- @@ -527,6 +539,9 @@ bool RimPolygonFilter::cellInsidePolygon2D( cvf::Vec3d center, std::array& points, const RigGridBase* grid ) { // we should look in depth using Z coordinate @@ -568,9 +583,10 @@ void RimPolygonFilter::updateCellsKIndexEclipse( const std::vector& const int gIdx = static_cast( grid->gridIndex() ); std::list foundCells; + const bool closedPolygon = m_closePolygon(); -#pragma omp parallel for // find all cells in the K layer that matches the polygon +#pragma omp parallel for for ( int i = 0; i < (int)grid->cellCountI(); i++ ) { for ( size_t j = 0; j < grid->cellCountJ(); j++ ) @@ -584,11 +600,23 @@ void RimPolygonFilter::updateCellsKIndexEclipse( const std::vector& std::array hexCorners; grid->cellCornerVertices( cellIdx, hexCorners.data() ); - // check if the polygon includes the cell - if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) + if ( closedPolygon ) { + // check if the polygon includes the cell + if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) + { #pragma omp critical - foundCells.push_back( cellIdx ); + foundCells.push_back( cellIdx ); + } + } + else + { + // check if the polyline touches the top face of the cell + if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + { +#pragma omp critical + foundCells.push_back( cellIdx ); + } } } } @@ -624,6 +652,8 @@ void RimPolygonFilter::updateCellsForEclipse( const std::vector& poi if ( m_polyFilterMode == PolygonFilterModeType::DEPTH_Z ) { + if ( !m_closePolygon() ) return; + for ( size_t gridIndex = 0; gridIndex < data->gridCount(); gridIndex++ ) { auto grid = data->grid( gridIndex ); @@ -690,6 +720,8 @@ void RimPolygonFilter::updateCellsDepthGeoMech( const std::vector& p //-------------------------------------------------------------------------------------------------- void RimPolygonFilter::updateCellsKIndexGeoMech( const std::vector& points, const RigFemPartGrid* grid, int partId ) { + const bool closedPolygon = m_closePolygon(); + // we need to find the K layer we hit with the first point size_t nk; // move the point a bit downwards to make sure it is inside something @@ -753,11 +785,23 @@ void RimPolygonFilter::updateCellsKIndexGeoMech( const std::vector& grid->cellCornerVertices( cellIdx, hexCorners.data() ); cvf::Vec3d center = grid->cellCentroid( cellIdx ); - // check if the polygon includes the cell - if ( cellInsidePolygon2D( center, hexCorners, points ) ) + if ( closedPolygon ) { + // check if the polygon includes the cell + if ( cellInsidePolygon2D( center, hexCorners, points ) ) + { #pragma omp critical - foundCells.push_back( cellIdx ); + foundCells.push_back( cellIdx ); + } + } + else + { + // check if the polyline touches the top face of the cell + if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + { +#pragma omp critical + foundCells.push_back( cellIdx ); + } } } } @@ -795,7 +839,10 @@ void RimPolygonFilter::updateCellsForGeoMech( const std::vector& poi if ( m_polyFilterMode == PolygonFilterModeType::DEPTH_Z ) { - updateCellsDepthGeoMech( points, grid, i ); + if ( m_closePolygon() ) + { + updateCellsDepthGeoMech( points, grid, i ); + } } else if ( m_polyFilterMode == PolygonFilterModeType::INDEX_K ) { @@ -823,11 +870,11 @@ void RimPolygonFilter::updateCells() if ( target->isEnabled() ) points.push_back( target->targetPointXYZ() ); } - // We need at least three points to make a sensible polygon - if ( points.size() < 3 ) return; + // We need at least three points to make a closed polygon, or just 2 for a polyline + if ( ( !m_closePolygon() && ( points.size() < 2 ) ) || ( m_closePolygon() && ( points.size() < 3 ) ) ) return; - // make sure first and last point is the same (req. by polygon methods used later) - points.push_back( points.front() ); + // make sure first and last point is the same (req. by closed polygon methods used later) + if ( m_closePolygon() ) points.push_back( points.front() ); RimEclipseCase* eCase = eclipseCase(); RimGeoMechCase* gCase = geoMechCase(); @@ -855,7 +902,7 @@ cvf::ref RimPolygonFilter::polyLinesData() const } pld->setPolyLine( line ); - pld->setLineAppearance( m_lineThickness, m_lineColor, true ); + pld->setLineAppearance( m_lineThickness, m_lineColor, m_closePolygon ); pld->setSphereAppearance( m_sphereRadiusFactor, m_sphereColor ); pld->setZPlaneLock( m_lockPolygonToPlane, -m_polygonPlaneDepth ); diff --git a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.h b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.h index 4ff7ca0ec5..94018cf843 100644 --- a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.h +++ b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.h @@ -133,6 +133,7 @@ private: caf::PdmField m_lockPolygonToPlane; caf::PdmField m_enableKFilter; caf::PdmField m_kFilterStr; + caf::PdmField m_closePolygon; std::shared_ptr m_pickTargetsEventHandler; diff --git a/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.cpp b/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.cpp index 3d51910552..fe30cdfbe6 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.cpp @@ -28,6 +28,8 @@ #include #include +#include +#include #include //-------------------------------------------------------------------------------------------------- @@ -698,7 +700,7 @@ double RigCellGeometryTools::getLengthOfPolygonAlongLine( const std::pair + RigCellGeometryTools::lineLineIntersection2D( const cvf::Vec3d a1, const cvf::Vec3d b1, const cvf::Vec3d a2, const cvf::Vec3d b2 ) +{ + constexpr double EPS = 0.000001; + double mua, mub; + double denom, numera, numerb; + const double x1 = a1.x(), x2 = b1.x(); + const double x3 = a2.x(), x4 = b2.x(); + const double y1 = a1.y(), y2 = b1.y(); + const double y3 = a2.y(), y4 = b2.y(); + + 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 ); + + // Are the lines coincident? + if ( ( std::abs( numera ) < EPS ) && ( std::abs( numerb ) < EPS ) && ( std::abs( denom ) < EPS ) ) + { + return std::make_pair( true, cvf::Vec2d( ( x1 + x2 ) / 2, ( y1 + y2 ) / 2 ) ); + } + + // Are the lines parallel? + if ( std::abs( denom ) < EPS ) + { + return std::make_pair( false, cvf::Vec2d() ); + } + + // Is the intersection along the segments? + mua = numera / denom; + mub = numerb / denom; + if ( mua < 0 || mua > 1 || mub < 0 || mub > 1 ) + { + return std::make_pair( false, cvf::Vec2d() ); + } + + return std::make_pair( true, cvf::Vec2d( x1 + mua * ( x2 - x1 ), y1 + mua * ( y2 - y1 ) ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// +/// Returns true if the line from a1 to b1 intersects the line from a2 to b2 +/// Operates only in the XY plane +/// +//-------------------------------------------------------------------------------------------------- +bool RigCellGeometryTools::lineIntersectsLine2D( const cvf::Vec3d a1, const cvf::Vec3d b1, const cvf::Vec3d a2, const cvf::Vec3d b2 ) +{ + return lineLineIntersection2D( a1, b1, a2, b2 ).first; +} + +//-------------------------------------------------------------------------------------------------- +/// +/// Returns true if the line from a to b intersects the closed, simple polygon defined by the corner +/// points in the input polygon vector, otherwise false +/// Operates only in the XY plane +/// +//-------------------------------------------------------------------------------------------------- +bool RigCellGeometryTools::lineIntersectsPolygon2D( const cvf::Vec3d a, const cvf::Vec3d b, const std::vector& polygon ) +{ + int nPolyLines = (int)polygon.size(); + + for ( int i = 1; i < nPolyLines; i++ ) + { + if ( lineIntersectsLine2D( a, b, polygon[i - 1], polygon[i] ) ) return true; + } + + if ( lineIntersectsLine2D( a, b, polygon[nPolyLines - 1], polygon[0] ) ) return true; + + return false; +} + +//-------------------------------------------------------------------------------------------------- +/// +/// Returns true if the polyline intersects the simple polygon defined by the NEGK face corners of the input cell +/// Operates only in the XY plane +/// +//-------------------------------------------------------------------------------------------------- +bool RigCellGeometryTools::polylineIntersectsCellNegK2D( const std::vector& polyline, const std::array& cellCorners ) +{ + const int nPoints = (int)polyline.size(); + const int nCorners = 4; + + for ( int i = 1; i < nPoints; i++ ) + { + for ( int j = 1; j < nCorners; j++ ) + { + if ( lineIntersectsLine2D( polyline[i - 1], polyline[i], cellCorners[j - 1], cellCorners[j] ) ) return true; + } + if ( lineIntersectsLine2D( polyline[i - 1], polyline[i], cellCorners[nCorners - 1], cellCorners[0] ) ) return true; + } + + return false; +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.h b/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.h index cbf2dbd7b6..28f3f91dc0 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.h +++ b/ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.h @@ -76,7 +76,14 @@ public: static std::vector unionOfPolygons( const std::vector>& polygons ); + // *** the 2D methods only looks at the X and Y coordinates of the input points *** + static bool pointInsidePolygon2D( const cvf::Vec3d point, const std::vector& polygon ); + static std::pair + lineLineIntersection2D( const cvf::Vec3d a1, const cvf::Vec3d b1, const cvf::Vec3d a2, const cvf::Vec3d b2 ); + static bool lineIntersectsLine2D( const cvf::Vec3d a1, const cvf::Vec3d b1, const cvf::Vec3d a2, const cvf::Vec3d b2 ); + static bool lineIntersectsPolygon2D( const cvf::Vec3d a, const cvf::Vec3d b, const std::vector& polygon ); + static bool polylineIntersectsCellNegK2D( const std::vector& polyline, const std::array& cellCorners ); private: static std::vector ajustPolygonToAvoidIntersectionsAtVertex( const std::vector& polyLine, diff --git a/ApplicationLibCode/UnitTests/RigCellGeometryTools-Test.cpp b/ApplicationLibCode/UnitTests/RigCellGeometryTools-Test.cpp index 714774fa6d..9a88034e67 100644 --- a/ApplicationLibCode/UnitTests/RigCellGeometryTools-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigCellGeometryTools-Test.cpp @@ -284,6 +284,55 @@ TEST( RigCellGeometryTools, lengthCalcTestTriangle ) EXPECT_LT( length, 1.8 ); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigCellGeometryTools, lineIntersectsLine2DTest ) +{ + cvf::Vec3d a1( 0, 0, 0 ); + cvf::Vec3d b1( 1, 1, 1 ); + + cvf::Vec3d a2( 0, 1, 500 ); + cvf::Vec3d b2( 1, 0, 7000 ); + + cvf::Vec3d a3( -10, -10, 0 ); + cvf::Vec3d b3( -4, -1, 0 ); + + EXPECT_TRUE( RigCellGeometryTools::lineIntersectsLine2D( a1, b1, a2, b2 ) ); + EXPECT_FALSE( RigCellGeometryTools::lineIntersectsLine2D( a1, b1, a3, b3 ) ); + + auto [intersect, point] = RigCellGeometryTools::lineLineIntersection2D( a1, b1, a2, b2 ); + EXPECT_TRUE( intersect ); + EXPECT_DOUBLE_EQ( 0.5, point.x() ); + EXPECT_DOUBLE_EQ( 0.5, point.y() ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigCellGeometryTools, lineIntersectsPolygon2DTest ) +{ + cvf::Vec3d a1( 0, 0, 0 ); + cvf::Vec3d b1( 1, 1, 1 ); + + cvf::Vec3d a2( 5, -5, 0 ); + cvf::Vec3d b2( 15, 25, 0 ); + + cvf::Vec3d p1( 10, 10, 0 ); + cvf::Vec3d p2( 11, 20, 1000 ); + cvf::Vec3d p3( 20, 7, -20 ); + cvf::Vec3d p4( 21, -1, 55 ); + + std::vector polygon; + polygon.push_back( p1 ); + polygon.push_back( p2 ); + polygon.push_back( p3 ); + polygon.push_back( p4 ); + + EXPECT_FALSE( RigCellGeometryTools::lineIntersectsPolygon2D( a1, b1, polygon ) ); + EXPECT_TRUE( RigCellGeometryTools::lineIntersectsPolygon2D( a2, b2, polygon ) ); +} + //-------------------------------------------------------------------------------------------------- /// //--------------------------------------------------------------------------------------------------