From b5964063a3bfcac35ce01c043a1a7d85b2762231 Mon Sep 17 00:00:00 2001 From: Gaute Lindkvist Date: Mon, 2 Sep 2019 15:11:00 +0200 Subject: [PATCH] Better spline handling --- .../RigWellPathGeometryTools.cpp | 161 ++++++++++-------- .../RigWellPathGeometryTools.h | 5 +- .../RigWellPathGeometryTools-Test.cpp | 91 ++-------- 3 files changed, 106 insertions(+), 151 deletions(-) diff --git a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp index 60b007e17e..501348a9df 100644 --- a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp @@ -23,8 +23,6 @@ #include "cvfMath.h" #include "cvfMatrix3.h" -#include - #include #include @@ -96,50 +94,41 @@ std::vector RigWellPathGeometryTools::interpolateMdFromTvd(const std::ve std::vector interpolatedMdValues; interpolatedMdValues.reserve(tvdValuesToInterpolateFrom.size()); - QPolygonF originalPoints; - for (size_t i = 0; i < originalMdValues.size(); ++i) - { - originalPoints << QPointF(originalMdValues[i], originalTvdValues[i]); - } - QwtSpline spline; - spline.setPoints(originalPoints); + + QwtSpline spline = createSpline(originalMdValues, originalTvdValues); - double lastTVDValue = -1.0; - std::vector segmentStartIndices = findSegmentIndices(originalMdValues, originalTvdValues, tvdValuesToInterpolateFrom); + std::vector segmentStartIndices = findSplineSegmentsContainingRoots(spline, originalMdValues, originalTvdValues, tvdValuesToInterpolateFrom); for (size_t i = 0; i < segmentStartIndices.size(); ++i) { double currentTVDValue = tvdValuesToInterpolateFrom[i]; - int startIndex = segmentStartIndices[i]; - int endIndex = startIndex + 1; + double startMD = spline.points().front().x(); + double endMD = spline.points().back().y(); + if (segmentStartIndices[i] != -1) + { + int startIndex = segmentStartIndices[i]; + int endIndex = startIndex + 1; - // Search interval for best MD value - double startMD = originalMdValues[startIndex]; - double endMD; - - double mdDiff = 0.0; - if (interpolatedMdValues.size() > 1) - { - mdDiff = interpolatedMdValues[i - 1] - interpolatedMdValues[i - 2]; - } + // Search interval for best MD value + startMD = spline.points()[startIndex].x(); + endMD = spline.points().back().y(); - - if (endIndex == originalMdValues.size()) - { - endMD = originalMdValues[startIndex] + mdDiff; - } - else - { - if (!interpolatedMdValues.empty()) + double mdDiff = 0.0; + if (interpolatedMdValues.size() > 1) { - 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 RigWellPathGeometryTools::interpolateMdFromTvd(const std::ve //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RigWellPathGeometryTools::findSegmentIndices(const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom) +std::vector RigWellPathGeometryTools::findSplineSegmentsContainingRoots(const QwtSpline& spline, const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& 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 segmentStartIndices; segmentStartIndices.reserve(tvdValuesToInterpolateFrom.size()); - int lastStartIndex = 0; + int lastSplineStartIndex = 0; for (std::vector::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 - double mdCurrent = originalMdValues[currentStartIndex]; - double mdNext = originalMdValues[nextStartIndex]; - double tvdCurrent = spline.value(mdCurrent); - double tvdNext = spline.value(mdNext); - if (std::abs(tvdCurrent - tvdValue) < std::abs(tvdNext - tvdValue)) - { - break; - } - currentStartIndex = nextStartIndex; + currentSplineStartIndex = nextStartIndex; + } + if (foundMatch) + { + segmentStartIndices.push_back(currentSplineStartIndex); + lastSplineStartIndex = currentSplineStartIndex; + } + else + { + 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& originalMdValues, const std::vector& 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; +} diff --git a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h index fcbb4fa3ee..5bb0971ef1 100644 --- a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h @@ -45,8 +45,6 @@ public: double angle); static std::vector interpolateMdFromTvd(const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom); - static std::vector findSegmentIndices(const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom); - private: static std::vector interpolateUndefinedNormals(const cvf::Vec3d& planeNormal, const std::vector& normals, @@ -54,4 +52,7 @@ private: static cvf::Vec3d estimateDominantDirectionInXYPlane(const std::vector& vertices); static double solveForX(const QwtSpline& spline, double minX, double maxX, double y); + + static QwtSpline createSpline(const std::vector& originalMdValues, const std::vector& originalTvdValues); + static std::vector findSplineSegmentsContainingRoots(const QwtSpline& spline, const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom); }; diff --git a/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp b/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp index 72a6c1c244..359e0d95a5 100644 --- a/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp +++ b/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp @@ -20,9 +20,11 @@ #include "RigWellPathGeometryTools.h" +#include #include #define TOLERANCE 1.0e-7 +#define SOLVER_TOLERANCE //-------------------------------------------------------------------------------------------------- /// @@ -32,14 +34,6 @@ TEST(RigWellPathGeometryTools, VerticalPath) std::vector mdValues = {100, 500, 1000}; std::vector tvdValues = {100, 500, 1000}; std::vector fullTVDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; - std::vector expectedSegmentIndices = {0, 0, 0, 0, 1, 1, 1, 1, 1, 2}; - std::vector 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 fullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTVDValues); EXPECT_EQ(fullTVDValues.size(), fullMDValues.size()); @@ -55,15 +49,6 @@ TEST(RigWellPathGeometryTools, LinearPath) std::vector tvdValues = {50, 250, 500}; std::vector fullTVDValues = {50, 100, 150, 200, 250, 300, 350, 400, 450, 500}; - std::vector expectedSegmentIndices = {0, 0, 0, 0, 1, 1, 1, 1, 1, 2}; - std::vector 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 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 expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3}; - std::vector 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 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 expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3}; - std::vector 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 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 expectedSegmentIndices = {0, 0, 1, 1, 1, 2, 2, 2, 2, 3}; - std::vector 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 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 mdValues = {150, 400, 800, 900}; - std::vector tvdValues; - for (double md : mdValues) - { - tvdValues.push_back(cubicFunction(md)); - } - std::vector fullMDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; - std::vector fullTvdValues; - for (double md : fullMDValues) - { - fullTvdValues.push_back(cubicFunction(md)); - } - - std::vector expectedSegmentIndices = {0, 0, 1, 1, 1, 1, 2, 2, 2, 3}; - std::vector 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 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]); - } -} \ No newline at end of file