Skip to content

Commit

Permalink
Add polyline filter support (#11048)
Browse files Browse the repository at this point in the history 
* Support selecting cells below polyline, both eclipse and geomech cases
* Add some polyline intersection helpers, with tests
  • Loading branch information
@jonjenssen
jonjenssen authored Jan 12, 2024
1 parent ca119af commit 8db29e0
Show file tree
Hide file tree
Showing 5 changed files with 224 additions and 17 deletions.
79 changes: 63 additions & 16 deletions ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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" );
Expand Down Expand Up @@ -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 );
Expand Down Expand Up @@ -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 );
}
}

//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -527,6 +539,9 @@ bool RimPolygonFilter::cellInsidePolygon2D( cvf::Vec3d center, std::array<cvf::V
return bInside;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimPolygonFilter::updateCellsDepthEclipse( const std::vector<cvf::Vec3d>& points, const RigGridBase* grid )
{
// we should look in depth using Z coordinate
Expand Down Expand Up @@ -568,9 +583,10 @@ void RimPolygonFilter::updateCellsKIndexEclipse( const std::vector<cvf::Vec3d>&
const int gIdx = static_cast<int>( grid->gridIndex() );

std::list<size_t> 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++ )
Expand All @@ -584,11 +600,23 @@ void RimPolygonFilter::updateCellsKIndexEclipse( const std::vector<cvf::Vec3d>&
std::array<cvf::Vec3d, 8> 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 );
}
}
}
}
Expand Down Expand Up @@ -624,6 +652,8 @@ void RimPolygonFilter::updateCellsForEclipse( const std::vector<cvf::Vec3d>& 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 );
Expand Down Expand Up @@ -690,6 +720,8 @@ void RimPolygonFilter::updateCellsDepthGeoMech( const std::vector<cvf::Vec3d>& p
//--------------------------------------------------------------------------------------------------
void RimPolygonFilter::updateCellsKIndexGeoMech( const std::vector<cvf::Vec3d>& 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
Expand Down Expand Up @@ -753,11 +785,23 @@ void RimPolygonFilter::updateCellsKIndexGeoMech( const std::vector<cvf::Vec3d>&
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 );
}
}
else
{
// check if the polyline touches the top face of the cell
if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) )
{
#pragma omp critical
foundCells.push_back( cellIdx );
foundCells.push_back( cellIdx );
}
}
}
}
Expand Down Expand Up @@ -795,7 +839,10 @@ void RimPolygonFilter::updateCellsForGeoMech( const std::vector<cvf::Vec3d>& 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 )
{
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -855,7 +902,7 @@ cvf::ref<RigPolyLinesData> 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 );

Expand Down
1 change: 1 addition & 0 deletions ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ class RimPolygonFilter : public RimCellFilter, public RimPolylinePickerInterface
caf::PdmField<bool> m_lockPolygonToPlane;
caf::PdmField<bool> m_enableKFilter;
caf::PdmField<QString> m_kFilterStr;
caf::PdmField<bool> m_closePolygon;

std::shared_ptr<RicPolylineTargetsPickEventHandler> m_pickTargetsEventHandler;

Expand Down
105 changes: 104 additions & 1 deletion ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@

#include <algorithm>
#include <array>
#include <cmath>
#include <utility>
#include <vector>

//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -698,7 +700,7 @@ double RigCellGeometryTools::getLengthOfPolygonAlongLine( const std::pair<cvf::V
{
cvf::BoundingBox lineBoundingBox;

for ( cvf::Vec3d polygonPoint : polygon )
for ( const cvf::Vec3d& polygonPoint : polygon )
{
cvf::Vec3d pointOnLine = cvf::GeometryTools::projectPointOnLine( line.first, line.second, polygonPoint, nullptr );
lineBoundingBox.add( pointOnLine );
Expand Down Expand Up @@ -794,6 +796,7 @@ inline double RigCellGeometryTools::isLeftOfLine2D( const cvf::Vec3d& point1, co

//--------------------------------------------------------------------------------------------------
/// winding number test for a point in a polygon
/// Operates only in the XY plane
/// Input: point = the point to test,
/// polygon[] = vertex points of a closed polygon of size n, where polygon[n-1]=polygon[0]
///
Expand Down Expand Up @@ -831,3 +834,103 @@ bool RigCellGeometryTools::pointInsidePolygon2D( const cvf::Vec3d point, const s
}
return wn != 0;
}

//--------------------------------------------------------------------------------------------------
/// Returns the intersection of line 1 (a1 to b1) and line 2 (a2 to b2).
/// - operates only in the XY plane
/// - returns true and the x,y intersection if the lines intersect
/// - returns false if they do not intersect
/// ref. http://www.paulbourke.net/geometry/pointlineplane/pdb.c
//--------------------------------------------------------------------------------------------------
std::pair<bool, cvf::Vec2d>
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<cvf::Vec3d>& 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<cvf::Vec3d>& polyline, const std::array<cvf::Vec3d, 8>& 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;
}
7 changes: 7 additions & 0 deletions ApplicationLibCode/ReservoirDataModel/RigCellGeometryTools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,14 @@ class RigCellGeometryTools

static std::vector<cvf::Vec3d> unionOfPolygons( const std::vector<std::vector<cvf::Vec3d>>& 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<cvf::Vec3d>& polygon );
static std::pair<bool, cvf::Vec2d>
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<cvf::Vec3d>& polygon );
static bool polylineIntersectsCellNegK2D( const std::vector<cvf::Vec3d>& polyline, const std::array<cvf::Vec3d, 8>& cellCorners );

private:
static std::vector<cvf::Vec3d> ajustPolygonToAvoidIntersectionsAtVertex( const std::vector<cvf::Vec3d>& polyLine,
Expand Down
49 changes: 49 additions & 0 deletions ApplicationLibCode/UnitTests/RigCellGeometryTools-Test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<cvf::Vec3d> 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 ) );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
Expand Down

0 comments on commit 8db29e0

Please sign in to comment.