Better spline handling

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp

@@ -23,8 +23,6 @@
 #include "cvfMath.h"
 #include "cvfMatrix3.h"
 
-#include <QDebug>
-
 #include <algorithm>
 #include <cmath>
 
@@ -96,50 +94,41 @@ std::vector<double> RigWellPathGeometryTools::interpolateMdFromTvd(const std::ve
     std::vector<double> interpolatedMdValues;
     interpolatedMdValues.reserve(tvdValuesToInterpolateFrom.size());
 
-    QPolygonF originalPoints;
+   
-    for (size_t i = 0; i < originalMdValues.size(); ++i)
+    QwtSpline spline = createSpline(originalMdValues, originalTvdValues);
-    {
-        originalPoints << QPointF(originalMdValues[i], originalTvdValues[i]);
-    }
-    QwtSpline spline;
-    spline.setPoints(originalPoints);
 
-    double lastTVDValue = -1.0;
+    std::vector<int> segmentStartIndices = findSplineSegmentsContainingRoots(spline, originalMdValues, originalTvdValues, tvdValuesToInterpolateFrom);
-    std::vector<int> segmentStartIndices = findSegmentIndices(originalMdValues, originalTvdValues, tvdValuesToInterpolateFrom);
     for (size_t i = 0; i < segmentStartIndices.size(); ++i)
     {
         double currentTVDValue = tvdValuesToInterpolateFrom[i];
-        int startIndex = segmentStartIndices[i];
+        double startMD = spline.points().front().x();
-        int endIndex = startIndex + 1;
+        double endMD = spline.points().back().y();
+        if (segmentStartIndices[i] != -1)
+        {
+            int startIndex = segmentStartIndices[i];
+            int endIndex = startIndex + 1;
 
-        // Search interval for best MD value
+            // Search interval for best MD value
-        double startMD = originalMdValues[startIndex];
+            startMD = spline.points()[startIndex].x();
-        double endMD;
+            endMD = spline.points().back().y();
-        
-        double mdDiff = 0.0;
-        if (interpolatedMdValues.size() > 1)
-        {
-            mdDiff = interpolatedMdValues[i - 1] - interpolatedMdValues[i - 2];
-        }
 
-       
+            double mdDiff = 0.0;
-        if (endIndex == originalMdValues.size())
+            if (interpolatedMdValues.size() > 1)
-        {
-            endMD = originalMdValues[startIndex] + mdDiff;
-        }
-        else
-        {
-            if (!interpolatedMdValues.empty())
             {
-                startMD = std::max(startMD, interpolatedMdValues.back() + 0.25 * mdDiff);
+                mdDiff = interpolatedMdValues[i - 1] - interpolatedMdValues[i - 2];
             }
-            endMD = originalMdValues[endIndex] + mdDiff * 0.5;
-        }
 
+            if (endIndex < spline.points().size())
+            {
+                if (!interpolatedMdValues.empty())
+                {
+                    startMD = std::max(startMD, interpolatedMdValues.back() + 0.1 * mdDiff);
+                }
+                endMD = spline.points()[endIndex].x();
+            }
+        }
         double mdValue = solveForX(spline, startMD, endMD, currentTVDValue);
         interpolatedMdValues.push_back(mdValue);
-
-        lastTVDValue = currentTVDValue;
     }
     return interpolatedMdValues;
 }
@@ -147,65 +136,48 @@ std::vector<double> RigWellPathGeometryTools::interpolateMdFromTvd(const std::ve
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::vector<int> RigWellPathGeometryTools::findSegmentIndices(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom)
+std::vector<int> RigWellPathGeometryTools::findSplineSegmentsContainingRoots(const QwtSpline& spline, const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom)
 {
-    CVF_ASSERT(!originalMdValues.empty());
-    if (originalMdValues.size() < 2u)
-    {
-        return {0};
-    }
-
-    QPolygonF polygon;
-    for (size_t i = 0; i < originalMdValues.size(); ++i)
-    {
-        polygon << QPointF(originalMdValues[i], originalTvdValues[i]);
-    }
-    QwtSpline spline;
-    spline.setPoints(polygon);
-    for (double md = 0.0; md < 1000; md += 50)
-    {
-        qDebug() << md << ", " << spline.value(md);
-    }
-
     std::vector<int> segmentStartIndices;
     segmentStartIndices.reserve(tvdValuesToInterpolateFrom.size());
 
-    int lastStartIndex = 0;
+    int lastSplineStartIndex = 0;
     for (std::vector<double>::const_iterator it = tvdValuesToInterpolateFrom.begin(); it != tvdValuesToInterpolateFrom.end(); ++it)
     {
         double tvdValue = *it;
-        int currentStartIndex = lastStartIndex;
+
+        int currentSplineStartIndex = lastSplineStartIndex;
        
+        bool foundMatch = false;
         // Increment current_it until we find an interval containing our TVD
-        while (currentStartIndex < (int)(originalTvdValues.size() - 1))
+        while (currentSplineStartIndex < spline.points().size() - 2)
         {
-            double diffCurrent = originalTvdValues[currentStartIndex] - tvdValue;
+            double diffCurrent = spline.points()[currentSplineStartIndex].y() - tvdValue;
             if (std::abs(diffCurrent) < 1.0e-8) // Current is matching the point
             {
+                foundMatch = true;
                 break;
             }
 
-            int nextStartIndex = currentStartIndex + 1;
+            int nextStartIndex = currentSplineStartIndex + 1;
 
-            double diffNext = originalTvdValues[nextStartIndex] - tvdValue;
+            double diffNext = spline.points()[nextStartIndex].y() - tvdValue;
             if (diffCurrent * diffNext < 0.0) // One is above, the other is below
             {
+                foundMatch = true;
                 break;
             }
-            // Attempt to interpolate
+            currentSplineStartIndex = nextStartIndex;
-            double mdCurrent = originalMdValues[currentStartIndex];
+        }
-            double mdNext = originalMdValues[nextStartIndex];
+        if (foundMatch)
-            double tvdCurrent = spline.value(mdCurrent);
+        {
-            double tvdNext = spline.value(mdNext);
+            segmentStartIndices.push_back(currentSplineStartIndex);
-            if (std::abs(tvdCurrent - tvdValue) < std::abs(tvdNext - tvdValue))
+            lastSplineStartIndex = currentSplineStartIndex;
-            {
+        }
-                break;
+        else
-            }
+        {
-            currentStartIndex = nextStartIndex;
+            segmentStartIndices.push_back(-1);
         }
-        segmentStartIndices.push_back(currentStartIndex);
-
-        lastStartIndex = currentStartIndex;
     }
 
     return segmentStartIndices;
@@ -338,3 +310,46 @@ double RigWellPathGeometryTools::solveForX(const QwtSpline& spline, double minX,
     }
     return (a + b) / 2.0;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+QwtSpline RigWellPathGeometryTools::createSpline(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues)
+{
+    QPolygonF polygon;
+    for (size_t i = 0; i < originalMdValues.size(); ++i)
+    {
+        polygon << QPointF(originalMdValues[i], originalTvdValues[i]);
+    }
+    QwtSplineCurveFitter curveFitter;
+    QPolygonF splinePoints = curveFitter.fitCurve(polygon);
+
+    // Extend spline from 0.0 to a large value for MD
+    {
+        double x1 = splinePoints[0].x();
+        double x2 = splinePoints[1].x();
+        double y1 = splinePoints[0].y();
+        double y2 = splinePoints[1].y();
+        double M = (y2 - y1) / (x2 - x1);
+
+        QPointF startPoint(0.0f, M * (0.0f - x1) + y1);
+        splinePoints.push_front(startPoint);
+    }
+    {
+        int N = splinePoints.size() - 1;
+        double x1 = splinePoints[N - 1].x();
+        double x2 = splinePoints[N].x();
+        double y1 = splinePoints[N - 1].y();
+        double y2 = splinePoints[N].y();
+        double M = (y2 - y1) / (x2 - x1);
+        double endX = 2.0 * splinePoints[N].x();
+
+        QPointF endPoint(endX, M  * (endX - x1) + y1);
+        splinePoints.push_back(endPoint);
+    }
+
+    QwtSpline spline;
+    spline.setPoints(splinePoints);
+
+    return spline;
+}

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h

@@ -45,8 +45,6 @@ public:
                                                                double                         angle);
     static std::vector<double> interpolateMdFromTvd(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom);
 
-    static std::vector<int> findSegmentIndices(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom);
-
 private:
     static std::vector<cvf::Vec3d> interpolateUndefinedNormals(const cvf::Vec3d& planeNormal,
                                                                const std::vector<cvf::Vec3d>& normals,
@@ -54,4 +52,7 @@ private:
     static cvf::Vec3d estimateDominantDirectionInXYPlane(const std::vector<cvf::Vec3d>& vertices);
 
     static double solveForX(const QwtSpline& spline, double minX, double maxX, double y);
+
+    static QwtSpline createSpline(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues);
+    static std::vector<int> findSplineSegmentsContainingRoots(const QwtSpline& spline, const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom);
 };

ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp

@@ -20,9 +20,11 @@
 
 #include "RigWellPathGeometryTools.h"
 
+#include <algorithm>
 #include <vector>
 
 #define TOLERANCE 1.0e-7
+#define SOLVER_TOLERANCE 
 
 //--------------------------------------------------------------------------------------------------
 ///
@@ -32,14 +34,6 @@ TEST(RigWellPathGeometryTools, VerticalPath)
     std::vector<double> mdValues  = {100, 500, 1000};
     std::vector<double> tvdValues = {100, 500, 1000};
     std::vector<double> fullTVDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000};
-    std::vector<int> expectedSegmentIndices = {0, 0, 0, 0, 1, 1, 1, 1, 1, 2};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTVDValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
     std::vector<double> fullMDValues  = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTVDValues);
 
     EXPECT_EQ(fullTVDValues.size(), fullMDValues.size());
@@ -55,15 +49,6 @@ TEST(RigWellPathGeometryTools, LinearPath)
     std::vector<double> tvdValues = {50, 250, 500};
     std::vector<double> fullTVDValues = {50, 100, 150, 200, 250, 300, 350, 400, 450, 500};
 
-    std::vector<int> expectedSegmentIndices = {0, 0, 0, 0, 1, 1, 1, 1, 1, 2};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTVDValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
-
     std::vector<double> fullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTVDValues);
 
     EXPECT_EQ(fullTVDValues.size(), fullMDValues.size());
@@ -92,20 +77,15 @@ TEST(RigWellPathGeometryTools, QuadraticPath)
     {
         fullTvdValues.push_back(quadraticFunction(md));
     }
-
+    
-    std::vector<int> expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTvdValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
     std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
     EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
     {
-        EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
     }
 }
 
@@ -129,19 +109,14 @@ TEST(RigWellPathGeometryTools, CubicPath)
         fullTvdValues.push_back(cubicFunction(md));
     }
 
-    std::vector<int> expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTvdValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
     std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
     EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
     }
 }
 
@@ -160,49 +135,13 @@ TEST(RigWellPathGeometryTools, CubicPathPoorSampling)
         fullTvdValues.push_back(cubicFunction(md));
     }
 
-    std::vector<int> expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTvdValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
     std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
     EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
     }
 }
-
-TEST(RigWellPathGeometryTools, CubicPathVeryPoorSampling)
-{
-    std::vector<double> mdValues = {150, 400, 800, 900};
-    std::vector<double> tvdValues;
-    for (double md : mdValues)
-    {
-        tvdValues.push_back(cubicFunction(md));
-    }
-    std::vector<double> fullMDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000};
-    std::vector<double> fullTvdValues;
-    for (double md : fullMDValues)
-    {
-        fullTvdValues.push_back(cubicFunction(md));
-    }
-
-    std::vector<int> expectedSegmentIndices = {0, 0, 1, 1, 1, 1, 2, 2, 2, 3};
-    std::vector<int> segmentIndices = RigWellPathGeometryTools::findSegmentIndices(mdValues, tvdValues, fullTvdValues);
-    EXPECT_EQ(expectedSegmentIndices.size(), segmentIndices.size());
-    for (size_t i = 0; i < expectedSegmentIndices.size(); ++i)
-    {
-        EXPECT_EQ(expectedSegmentIndices[i], segmentIndices[i]);
-    }
-
-    std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
-    EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
-    for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
-    {
-        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
-    }
-}