#1514 Add polyline/polygon intersection method

To be used for wellpath stimplan intersection

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp

@@ -233,7 +233,7 @@ double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLen
         polygon.push_back(line2.first);
 
         //Use clipper to find overlap between bbpolygon and fracture
-        std::vector<std::vector<cvf::Vec3d> > clippedPolygons = clipPolygons(polygonToCalcLengthOf, polygon);
+        std::vector<std::vector<cvf::Vec3d> > clippedPolygons = intersectPolygons(polygonToCalcLengthOf, polygon);
 
         double area = 0;
         double length = 0;
@@ -264,29 +264,38 @@ double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLen
     return areaWeightedLength;
 }
 
+double clipperConversionFactor = 10000; //For transform to clipper int 
+
+ClipperLib::IntPoint toClipperPoint(const cvf::Vec3d& cvfPoint)
+{
+    int xInt = cvfPoint.x()*clipperConversionFactor;
+    int yInt = cvfPoint.y()*clipperConversionFactor;
+    int zInt = cvfPoint.z()*clipperConversionFactor;
+    return  ClipperLib::IntPoint(xInt, yInt, zInt); 
+}
+
+cvf::Vec3d fromClipperPoint(const ClipperLib::IntPoint& clipPoint)
+{
+    return cvf::Vec3d (clipPoint.X, clipPoint.Y, clipPoint.Z ) /clipperConversionFactor;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///  
 //--------------------------------------------------------------------------------------------------
-std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2)
+std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::intersectPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2)
 {
-    int polygonScaleFactor = 10000; //For transform to clipper int 
-    int xInt, yInt;
     
-    //Convert to int for clipper library and store as clipper "path"
+    // Convert to int for clipper library and store as clipper "path"
     ClipperLib::Path polygon1path;
     for (cvf::Vec3d& v : polygon1)
     {
-        xInt = v.x()*polygonScaleFactor;
+        polygon1path.push_back(toClipperPoint(v));
-        yInt = v.y()*polygonScaleFactor;
-        polygon1path.push_back(ClipperLib::IntPoint(xInt, yInt));
     }
 
     ClipperLib::Path polygon2path;
     for (cvf::Vec3d& v : polygon2)
     {
-        xInt = v.x()*polygonScaleFactor;
+        polygon2path.push_back(toClipperPoint(v));
-        yInt = v.y()*polygonScaleFactor;
-        polygon2path.push_back(ClipperLib::IntPoint(xInt, yInt));
     }
 
     ClipperLib::Clipper clpr;
@@ -296,17 +305,14 @@ std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolygons(std::ve
     ClipperLib::Paths solution;
     clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
 
-    //Convert back to std::vector<std::vector<cvf::Vec3d> >
+    // Convert back to std::vector<std::vector<cvf::Vec3d> >
     std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
     for (ClipperLib::Path pathInSol : solution)
     {
         std::vector<cvf::Vec3d> clippedPolygon;
         for (ClipperLib::IntPoint IntPosition : pathInSol)
         {
-            cvf::Vec3d v = cvf::Vec3d::ZERO;
+            clippedPolygon.push_back(fromClipperPoint(IntPosition));
-            v.x() = (float)IntPosition.X / (float)polygonScaleFactor;
-            v.y() = (float)IntPosition.Y / (float)polygonScaleFactor;
-            clippedPolygon.push_back(v);
         }
         clippedPolygons.push_back(clippedPolygon);
     }
@@ -315,6 +321,101 @@ std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolygons(std::ve
 
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void fillInterpolatedSubjectZ(ClipperLib::IntPoint& e1bot,
+                              ClipperLib::IntPoint& e1top,
+                              ClipperLib::IntPoint& e2bot,
+                              ClipperLib::IntPoint& e2top,
+                              ClipperLib::IntPoint& pt)
+{
+    int e2XRange = (e2top.X - e2bot.X);
+    int e2YRange = (e2top.Y - e2bot.Y);
+
+    double e2Length =  sqrt(e2XRange*e2XRange + e2YRange*e2YRange);
+    
+    if (e2Length <= 1)
+    {
+        pt.Z = e2bot.Z;
+        return;
+    }
+
+    int e2BotPtXRange = pt.X - e2bot.X;
+    int e2BotPtYRange = pt.Y - e2bot.Y;
+
+    double e2BotPtLength =  sqrt(e2BotPtXRange*e2BotPtXRange + e2BotPtYRange*e2BotPtYRange);
+
+    double fraction = e2BotPtLength/e2Length;
+
+    pt.Z = e2bot.Z + fraction*(e2top.Z - e2bot.Z);
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void fillUndefinedZ(ClipperLib::IntPoint& e1bot,
+                              ClipperLib::IntPoint& e1top,
+                              ClipperLib::IntPoint& e2bot,
+                              ClipperLib::IntPoint& e2top,
+                              ClipperLib::IntPoint& pt)
+{
+   pt.Z = HUGE_VAL;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygon(std::vector<cvf::Vec3d> polyLine, 
+                                                                                  std::vector<cvf::Vec3d> polygon, 
+                                                                                  ZInterpolationType interpolType)
+{
+    int polygonScaleFactor = 10000; //For transform to clipper int 
+
+    //Convert to int for clipper library and store as clipper "path"
+    ClipperLib::Path polyLinePath;
+    for (cvf::Vec3d& v : polyLine)
+    {
+        polyLinePath.push_back(toClipperPoint(v));
+    }
+
+    ClipperLib::Path polygonPath;
+    for (cvf::Vec3d& v : polygon)
+    {
+        polygonPath.push_back(toClipperPoint(v));
+    }
+
+    ClipperLib::Clipper clpr;
+    clpr.AddPath(polyLinePath, ClipperLib::ptSubject, false);
+    clpr.AddPath(polygonPath, ClipperLib::ptClip, true);
+
+    if ( interpolType == INTERPOLATE_LINE_Z )
+    {
+        clpr.ZFillFunction(&fillInterpolatedSubjectZ);
+    }
+    else if ( interpolType == USE_HUGEVAL )
+    {
+        clpr.ZFillFunction(&fillUndefinedZ);
+    }
+
+    ClipperLib::Paths solution;
+    clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
+
+    //Convert back to std::vector<std::vector<cvf::Vec3d> >
+    std::vector<std::vector<cvf::Vec3d> > clippedPolyline;
+    for (ClipperLib::Path pathInSol : solution)
+    {
+        std::vector<cvf::Vec3d> clippedPolygon;
+        for (ClipperLib::IntPoint IntPosition : pathInSol)
+        {
+            clippedPolygon.push_back(fromClipperPoint(IntPosition));
+        }
+        clippedPolyline.push_back(clippedPolygon);
+    }
+
+    return clippedPolyline;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///  
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h

@@ -43,7 +43,10 @@ public:
     
     static double polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygon2d);
    
-    static std::vector<std::vector<cvf::Vec3d> >  clipPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2);
+    static std::vector<std::vector<cvf::Vec3d> >  intersectPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2);
+
+    enum ZInterpolationType { INTERPOLATE_LINE_Z, USE_HUGEVAL, USE_ZERO};
+    static std::vector<std::vector<cvf::Vec3d> >  clipPolylineByPolygon(std::vector<cvf::Vec3d> polyLine, std::vector<cvf::Vec3d> polygon, ZInterpolationType interpolType = USE_ZERO);
 
     static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine);
 

ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.cpp

@@ -167,7 +167,7 @@ void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsM
 
         for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
         {
-            std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, stimPlanPolygon);
+            std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::intersectPolygons(planeCellPolygon, stimPlanPolygon);
             for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
             {
                 polygonsForStimPlanCellInEclipseCell.push_back(clippedPolygon);

ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp

@@ -176,7 +176,7 @@ std::vector<RigFracturedEclipseCellExportData>  RigFractureTransCalc::computeTra
 
         for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
         {
-            std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, fracPolygonDouble);
+            std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::intersectPolygons(planeCellPolygon, fracPolygonDouble);
             for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
             {
                 polygonsDescribingFractureInCell.push_back(clippedPolygon);

ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp

@@ -122,7 +122,7 @@ double RigStimPlanUpscalingCalc::computeHAupscale(RimStimPlanFractureTemplate* f
         {
             if (stimPlanCell->getConductivtyValue() > 1e-7)
             {
-                std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
+                std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::intersectPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
                 if (clippedStimPlanPolygons.size() > 0)
                 {
                     for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
@@ -191,7 +191,7 @@ double RigStimPlanUpscalingCalc::computeAHupscale(RimStimPlanFractureTemplate* f
         {
             if (stimPlanCell->getConductivtyValue() > 1e-7)
             {
-                std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
+                std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::intersectPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
                 if (clippedStimPlanPolygons.size() > 0)
                 {
                     for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)

ThirdParty/clipper/clipper.hpp

@@ -41,7 +41,7 @@
 //#define use_int32
 
 //use_xyz: adds a Z member to IntPoint. Adds a minor cost to perfomance.
-//#define use_xyz
+#define use_xyz
 
 //use_lines: Enables line clipping. Adds a very minor cost to performance.
 #define use_lines