Move line and surface intersection methods to RigSurfaceResampler

ApplicationLibCode/ModelVisualization/Intersections/RivExtrudedCurveIntersectionGeometryGenerator.cpp

@@ -23,7 +23,6 @@
 #include "RigResultAccessor.h"
 #include "RigSurface.h"
 #include "RigSurfaceResampler.h"
-#include "RigWellPath.h"
 
 #include "Rim3dView.h"
 #include "RimCase.h"
@@ -188,15 +187,15 @@ void RivExtrudedCurveIntersectionGeometryGenerator::calculateSurfaceIntersection
 
             // Resample polyline to required resolution
             const double maxLineSegmentLength = 1.0;
-            auto         resampledPolyline    = computeResampledPolyline( firstPolyLine, maxLineSegmentLength );
+            const auto resampledPolyline = RigSurfaceResampler::computeResampledPolylineWithSegmentInfo( firstPolyLine, maxLineSegmentLength );
 
-            for ( auto [point, segmentIndex] : resampledPolyline )
+            for ( const auto& [point, segmentIndex] : resampledPolyline )
             {
                 cvf::Vec3d pointAbove = cvf::Vec3d( point.x(), point.y(), 10000.0 );
                 cvf::Vec3d pointBelow = cvf::Vec3d( point.x(), point.y(), -10000.0 );
 
                 cvf::Vec3d intersectionPoint;
-                bool       foundMatch = RigSurfaceResampler::resamplePoint( surface, pointAbove, pointBelow, intersectionPoint );
+                bool foundMatch = RigSurfaceResampler::findClosestPointOnSurface( surface, pointAbove, pointBelow, intersectionPoint );
                 if ( foundMatch )
                 {
                     const size_t lineIndex = 0;
@@ -809,58 +808,6 @@ cvf::Vec3d RivExtrudedCurveIntersectionGeometryGenerator::transformPointByPolyli
     return domainCoord.getTransformedPoint( m_lineSegmentTransforms[lineIndex][segmentIndex] );
 }
 
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-std::vector<std::pair<cvf::Vec3d, size_t>>
-    RivExtrudedCurveIntersectionGeometryGenerator::computeResampledPolyline( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance )
-{
-    // Segments along the original polyline must be provided to be able to find the associated transform matrix
-    std::vector<std::pair<cvf::Vec3d, size_t>> resampledPolyline;
-
-    if ( polyline.size() > 1 )
-    {
-        std::vector<double> measuredDepth;
-        {
-            double aggregatedLength = 0.0;
-
-            cvf::Vec3d previousPoint = polyline.front();
-            measuredDepth.push_back( aggregatedLength );
-
-            for ( size_t i = 1; i < polyline.size(); i++ )
-            {
-                aggregatedLength += ( previousPoint - polyline[i] ).length();
-                previousPoint = polyline[i];
-                measuredDepth.push_back( aggregatedLength );
-            }
-        }
-
-        // Use RigWellPath to perform the interpolation along a line based on measured depth
-        RigWellPath dummyWellPath( polyline, measuredDepth );
-
-        for ( size_t i = 1; i < polyline.size(); i++ )
-        {
-            const auto& lineSegmentStart = polyline[i - 1];
-            const auto& lineSegmentEnd   = polyline[i];
-
-            auto startMD = measuredDepth[i - 1];
-            auto endMD   = measuredDepth[i];
-
-            const size_t segmentIndex = i - 1;
-            resampledPolyline.emplace_back( lineSegmentStart, segmentIndex );
-
-            for ( auto md = startMD + resamplingDistance; md < endMD; md += resamplingDistance )
-            {
-                resampledPolyline.emplace_back( dummyWellPath.interpolatedPointAlongWellPath( md ), segmentIndex );
-            }
-
-            resampledPolyline.emplace_back( lineSegmentEnd, segmentIndex );
-        }
-    }
-
-    return resampledPolyline;
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ModelVisualization/Intersections/RivExtrudedCurveIntersectionGeometryGenerator.h

@@ -97,9 +97,6 @@ private:
     const std::vector<std::vector<cvf::Vec3d>>& flattenedOrOffsettedPolyLines();
     cvf::Vec3d transformPointByPolylineSegmentIndex( const cvf::Vec3d& domainCoord, size_t lineIndex, size_t segmentIndex );
 
-    static std::vector<std::pair<cvf::Vec3d, size_t>> computeResampledPolyline( const std::vector<cvf::Vec3d>& polyline,
-                                                                                double                         resamplingDistance );
-
 private:
     RimExtrudedCurveIntersection*              m_intersection;
     cvf::ref<RivIntersectionHexGridInterface>  m_hexGrid;

ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp

@@ -17,6 +17,7 @@
 /////////////////////////////////////////////////////////////////////////////////
 
 #include "RigSurfaceResampler.h"
+#include "RigWellPath.h"
 
 #include "cvfGeometryTools.h"
 
@@ -26,7 +27,8 @@
 #include <limits>
 
 //--------------------------------------------------------------------------------------------------
-///
+/// Create a surface with same XY coordinates as targetSurface. Evaluate Z value at [X,Y] in the source surface. Search in XY plane to find
+/// closest point, limited by a max distance from [X,Y].
 //--------------------------------------------------------------------------------------------------
 cvf::ref<RigSurface> RigSurfaceResampler::resampleSurface( cvf::ref<RigSurface> targetSurface, cvf::ref<RigSurface> surface )
 {
@@ -37,13 +39,13 @@ cvf::ref<RigSurface> RigSurfaceResampler::resampleSurface( cvf::ref<RigSurface>
 
     std::vector<cvf::Vec3d> resampledVerts;
 
-    for ( auto targetVert : targetVerts )
+    for ( const auto& targetVert : targetVerts )
     {
-        cvf::Vec3d pointAbove = cvf::Vec3d( targetVert.x(), targetVert.y(), 10000.0 );
-        cvf::Vec3d pointBelow = cvf::Vec3d( targetVert.x(), targetVert.y(), -10000.0 );
+        const cvf::Vec3d pointAbove = cvf::Vec3d( targetVert.x(), targetVert.y(), 10000.0 );
+        const cvf::Vec3d pointBelow = cvf::Vec3d( targetVert.x(), targetVert.y(), -10000.0 );
 
         cvf::Vec3d intersectionPoint;
-        bool       foundMatch = resamplePoint( surface.p(), pointAbove, pointBelow, intersectionPoint );
+        bool       foundMatch = findClosestPointOnSurface( surface.p(), pointAbove, pointBelow, intersectionPoint );
         if ( !foundMatch ) intersectionPoint = cvf::Vec3d( targetVert.x(), targetVert.y(), std::numeric_limits<double>::infinity() );
 
         resampledVerts.push_back( intersectionPoint );
@@ -55,54 +57,106 @@ cvf::ref<RigSurface> RigSurfaceResampler::resampleSurface( cvf::ref<RigSurface>
 }
 
 //--------------------------------------------------------------------------------------------------
-///
+/// If an intersection point is found, return true and set intersectionPoint to the intersection point.
+/// If no intersection point is found, return false.
 //--------------------------------------------------------------------------------------------------
-bool RigSurfaceResampler::resamplePoint( RigSurface* surface, const cvf::Vec3d& pointAbove, const cvf::Vec3d& pointBelow, cvf::Vec3d& intersectionPoint )
+bool RigSurfaceResampler::computeIntersectionWithLine( RigSurface* surface, const cvf::Vec3d& p1, const cvf::Vec3d& p2, cvf::Vec3d& intersectionPoint )
 {
+    if ( !surface ) return false;
+
     surface->ensureIntersectionSearchTreeIsBuilt();
 
     cvf::BoundingBox bb;
-    bb.add( pointAbove );
-    bb.add( pointBelow );
+    bb.add( p1 );
+    bb.add( p2 );
 
     const std::vector<unsigned int>& triIndices = surface->triangleIndices();
     const std::vector<cvf::Vec3d>&   vertices   = surface->vertices();
 
-    {
-        std::vector<size_t> triangleStartIndices;
-        surface->findIntersectingTriangles( bb, &triangleStartIndices );
+    bool dummy = false;
 
-        if ( !triangleStartIndices.empty() )
-        {
-            for ( auto triangleStartIndex : triangleStartIndices )
-            {
-                bool isLineDirDotNormalNegative = false;
-                if ( cvf::GeometryTools::intersectLineSegmentTriangle( pointAbove,
-                                                                       pointBelow,
-                                                                       vertices[triIndices[triangleStartIndex + 0]],
-                                                                       vertices[triIndices[triangleStartIndex + 1]],
-                                                                       vertices[triIndices[triangleStartIndex + 2]],
-                                                                       &intersectionPoint,
-                                                                       &isLineDirDotNormalNegative ) == 1 )
-                    return true;
-            }
-        }
+    std::vector<size_t> triangleStartIndices;
+    surface->findIntersectingTriangles( bb, &triangleStartIndices );
+
+    for ( auto startIndex : triangleStartIndices )
+    {
+        if ( cvf::GeometryTools::intersectLineSegmentTriangle( p1,
+                                                               p2,
+                                                               vertices[triIndices[startIndex + 0]],
+                                                               vertices[triIndices[startIndex + 1]],
+                                                               vertices[triIndices[startIndex + 2]],
+                                                               &intersectionPoint,
+                                                               &dummy ) == 1 )
+            return true;
     }
 
+    return false;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Find the closest point on the surface to the line. The search is limited to a max distance from the line based on the resolution
+/// of the surface.
+//--------------------------------------------------------------------------------------------------
+bool RigSurfaceResampler::findClosestPointOnSurface( RigSurface* surface, const cvf::Vec3d& p1, const cvf::Vec3d& p2, cvf::Vec3d& intersectionPoint )
+{
+    if ( computeIntersectionWithLine( surface, p1, p2, intersectionPoint ) ) return true;
+
+    cvf::BoundingBox bb;
+    bb.add( p1 );
+    bb.add( p2 );
+
+    const std::vector<unsigned int>& triIndices = surface->triangleIndices();
+    const std::vector<cvf::Vec3d>&   vertices   = surface->vertices();
+
     double maxDistance = computeMaxDistance( surface );
 
-    // Expand the bounding box to cover a larger volume. Use this volume to find intersections.
+    // Expand the bounding box to cover a larger volume around bounding box of the line segment.
     bb.expand( maxDistance );
 
     std::vector<size_t> triangleStartIndices;
     surface->findIntersectingTriangles( bb, &triangleStartIndices );
 
-    return findClosestPointXY( pointAbove, vertices, triIndices, triangleStartIndices, maxDistance, intersectionPoint );
+    return findClosestPointXY( p1, vertices, triIndices, triangleStartIndices, maxDistance, intersectionPoint );
 }
 
 //--------------------------------------------------------------------------------------------------
-/// Find the closest vertex to targetPoint (must be closer than maxDistance) in XY plane.
-/// Unit maxDistance: meter.
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<cvf::Vec3d> RigSurfaceResampler::computeResampledPolyline( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance )
+{
+    auto polylineAndSegmentData = computeResampledPolylineWithSegmentInfoImpl( polyline, resamplingDistance );
+
+    return polylineAndSegmentData.first;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<std::pair<cvf::Vec3d, size_t>>
+    RigSurfaceResampler::computeResampledPolylineWithSegmentInfo( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance )
+{
+    // Segments along the original polyline must be provided to be able to find the associated transform matrix
+
+    std::vector<cvf::Vec3d> resampledPolyline;
+    std::vector<size_t>     segmentIndices;
+
+    std::tie( resampledPolyline, segmentIndices ) = computeResampledPolylineWithSegmentInfoImpl( polyline, resamplingDistance );
+
+    if ( resampledPolyline.size() != segmentIndices.size() ) return {};
+
+    // Create vector of pairs based on pair of two vectors
+
+    std::vector<std::pair<cvf::Vec3d, size_t>> polyLineAndSegment;
+    for ( size_t i = 0; i < polyLineAndSegment.size(); i++ )
+    {
+        polyLineAndSegment.push_back( std::make_pair( resampledPolyline[i], segmentIndices[i] ) );
+    }
+
+    return polyLineAndSegment;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Find the closest point in XY plane closer than maxDistance
 //--------------------------------------------------------------------------------------------------
 bool RigSurfaceResampler::findClosestPointXY( const cvf::Vec3d&                targetPoint,
                                               const std::vector<cvf::Vec3d>&   vertices,
@@ -174,3 +228,59 @@ double RigSurfaceResampler::computeMaxDistance( RigSurface* surface )
 
     return distance;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<std::vector<cvf::Vec3d>, std::vector<size_t>>
+    RigSurfaceResampler::computeResampledPolylineWithSegmentInfoImpl( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance )
+{
+    // Segments along the original polyline must be provided to be able to find the associated transform matrix
+    std::vector<cvf::Vec3d> resampledPolyline;
+    std::vector<size_t>     segmentIndices;
+
+    if ( polyline.size() > 1 )
+    {
+        std::vector<double> measuredDepth;
+        {
+            double aggregatedLength = 0.0;
+
+            cvf::Vec3d previousPoint = polyline.front();
+            measuredDepth.push_back( aggregatedLength );
+
+            for ( size_t i = 1; i < polyline.size(); i++ )
+            {
+                aggregatedLength += ( previousPoint - polyline[i] ).length();
+                previousPoint = polyline[i];
+                measuredDepth.push_back( aggregatedLength );
+            }
+        }
+
+        // Use RigWellPath to perform the interpolation along a line based on measured depth
+        RigWellPath dummyWellPath( polyline, measuredDepth );
+
+        for ( size_t i = 1; i < polyline.size(); i++ )
+        {
+            const auto& lineSegmentStart = polyline[i - 1];
+            const auto& lineSegmentEnd   = polyline[i];
+
+            auto startMD = measuredDepth[i - 1];
+            auto endMD   = measuredDepth[i];
+
+            const size_t segmentIndex = i - 1;
+            resampledPolyline.emplace_back( lineSegmentStart );
+            segmentIndices.emplace_back( segmentIndex );
+
+            for ( auto md = startMD + resamplingDistance; md < endMD; md += resamplingDistance )
+            {
+                resampledPolyline.emplace_back( dummyWellPath.interpolatedPointAlongWellPath( md ) );
+                segmentIndices.emplace_back( segmentIndex );
+            }
+
+            resampledPolyline.emplace_back( lineSegmentEnd );
+            segmentIndices.emplace_back( segmentIndex );
+        }
+    }
+
+    return std::make_pair( resampledPolyline, segmentIndices );
+}

ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h

@@ -29,7 +29,12 @@ class RigSurfaceResampler
 public:
     static cvf::ref<RigSurface> resampleSurface( cvf::ref<RigSurface> targetSurface, cvf::ref<RigSurface> surface );
 
-    static bool resamplePoint( RigSurface* surface, const cvf::Vec3d& pointAbove, const cvf::Vec3d& pointBelow, cvf::Vec3d& intersectionPoint );
+    static bool computeIntersectionWithLine( RigSurface* surface, const cvf::Vec3d& p1, const cvf::Vec3d& p2, cvf::Vec3d& intersectionPoint );
+    static bool findClosestPointOnSurface( RigSurface* surface, const cvf::Vec3d& p1, const cvf::Vec3d& p2, cvf::Vec3d& intersectionPoint );
+
+    static std::vector<cvf::Vec3d> computeResampledPolyline( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance );
+    static std::vector<std::pair<cvf::Vec3d, size_t>> computeResampledPolylineWithSegmentInfo( const std::vector<cvf::Vec3d>& polyline,
+                                                                                               double resamplingDistance );
 
 private:
     static bool findClosestPointXY( const cvf::Vec3d&                targetPoint,
@@ -40,4 +45,7 @@ private:
                                     cvf::Vec3d&                      intersectionPoint );
 
     static double computeMaxDistance( RigSurface* surface );
+
+    static std::pair<std::vector<cvf::Vec3d>, std::vector<size_t>>
+        computeResampledPolylineWithSegmentInfoImpl( const std::vector<cvf::Vec3d>& polyline, double resamplingDistance );
 };