Refactor PLT stacked curve code to be more general

ApplicationCode/ProjectDataModel/Flow/RimWellFlowRateCurve.cpp

@@ -165,7 +165,9 @@ void RimWellFlowRateCurve::onLoadDataAndUpdate( bool updateParentPlot )
     {
         m_qwtPlotCurve->setTitle( createCurveAutoName() );
 
-        updateStackedPlotData();
+        RimWellLogTrack* track = nullptr;
+        this->firstAncestorOrThisOfTypeAsserted( track );
+        track->updateStackedCurveData();
 
         updateZoomInParentPlot();
 
@@ -184,9 +186,9 @@ void RimWellFlowRateCurve::updateCurveAppearance()
     {
         RimWellLogTrack* wellLogTrack;
         firstAncestorOrThisOfTypeAsserted( wellLogTrack );
-        std::map<int, std::vector<RimWellFlowRateCurve*>> stackedCurveGroups = wellLogTrack->visibleStackedCurves();
+        std::map<int, std::vector<RimWellLogCurve*>> stackedCurveGroups = wellLogTrack->visibleStackedCurves();
-        const std::vector<RimWellFlowRateCurve*>&         curveGroup         = stackedCurveGroups[this->m_groupId];
+        const std::vector<RimWellLogCurve*>&         curveGroup         = stackedCurveGroups[this->m_groupId];
-        isLastCurveInGroup                                                   = ( curveGroup.back() == this );
+        isLastCurveInGroup                                              = ( curveGroup.back() == this );
     }
 
     if ( isUsingConnectionNumberDepthType() )
@@ -241,81 +243,6 @@ void RimWellFlowRateCurve::defineUiOrdering( QString uiConfigName, caf::PdmUiOrd
     uiOrdering.skipRemainingFields();
 }
 
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-void RimWellFlowRateCurve::updateStackedPlotData()
-{
-    RimWellLogPlot* wellLogPlot;
-    firstAncestorOrThisOfTypeAsserted( wellLogPlot );
-
-    RimWellLogTrack* wellLogTrack;
-    firstAncestorOrThisOfTypeAsserted( wellLogTrack );
-
-    RimWellLogPlot::DepthTypeEnum depthType   = wellLogPlot->depthType();
-    RiaDefines::DepthUnitType     displayUnit = wellLogPlot->depthUnit();
-    if ( depthType == RiaDefines::DepthTypeEnum::CONNECTION_NUMBER )
-    {
-        displayUnit = RiaDefines::DepthUnitType::UNIT_NONE;
-    }
-
-    std::vector<double>                    depthValues;
-    std::vector<double>                    stackedValues;
-    std::vector<std::pair<size_t, size_t>> polyLineStartStopIndices;
-
-    // Z-position of curve, to draw them in correct order
-    double zPos = -10000.0 + 100.0 * static_cast<double>( groupId() );
-
-    // Starting way behind the grid (z == 0) at -10000 giving room for 100 groups with 100 curves each before getting
-    // above the grid
-    {
-        std::map<int, std::vector<RimWellFlowRateCurve*>> stackedCurveGroups = wellLogTrack->visibleStackedCurves();
-
-        std::vector<RimWellFlowRateCurve*> stackedCurves;
-
-        if ( stackedCurveGroups.count( groupId() ) > 0 )
-        {
-            stackedCurves = stackedCurveGroups[groupId()];
-        }
-
-        std::vector<double> allDepthValues = curveData()->depths( depthType );
-        if ( allDepthValues.empty() ) return;
-
-        std::vector<double> allStackedValues( allDepthValues.size() );
-
-        for ( RimWellFlowRateCurve* stCurve : stackedCurves )
-        {
-            std::vector<double> allValues = stCurve->curveData()->xValues();
-
-            for ( size_t i = 0; i < allValues.size(); ++i )
-            {
-                if ( allValues[i] != HUGE_VAL )
-                {
-                    allStackedValues[i] += allValues[i];
-                }
-            }
-
-            if ( stCurve == this ) break;
-            zPos -= 1.0;
-        }
-
-        RigWellLogCurveData tempCurveData;
-        tempCurveData.setValuesAndDepths( allStackedValues, allDepthValues, depthType, 0.0, displayUnit, false );
-
-        depthValues              = tempCurveData.depthPlotValues( depthType, displayUnit );
-        stackedValues            = tempCurveData.xPlotValues();
-        polyLineStartStopIndices = tempCurveData.polylineStartStopIndices();
-    }
-
-    auto minmax_it = std::minmax_element( stackedValues.begin(), stackedValues.end() );
-
-    this->setOverrideCurveDataXRange( *( minmax_it.first ), *( minmax_it.second ) );
-    m_qwtPlotCurve->setSamples( stackedValues.data(), depthValues.data(), static_cast<int>( depthValues.size() ) );
-    m_qwtPlotCurve->setLineSegmentStartStopIndices( polyLineStartStopIndices );
-
-    m_qwtPlotCurve->setZ( zPos );
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/Flow/RimWellFlowRateCurve.h

@@ -42,7 +42,6 @@ public:
                                     RiaDefines::DepthTypeEnum  depthType,
                                     const std::vector<double>& depthValues,
                                     const std::vector<double>& flowRates );
-    void updateStackedPlotData();
 
     RimEclipseResultCase* rimCase();
     int                   timeStep();

ApplicationCode/ProjectDataModel/RimWellLogCurve.cpp

@@ -195,6 +195,22 @@ QString RimWellLogCurve::wellLogCurveIconName()
     return ":/WellLogCurve16x16.png";
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimWellLogCurve::setOverrideCurveData( const std::vector<double>&               xValues,
+                                            const std::vector<double>&               depthValues,
+                                            const RiaCurveDataTools::CurveIntervals& curveIntervals )
+{
+    auto minmax_it = std::minmax_element( xValues.begin(), xValues.end() );
+    this->setOverrideCurveDataXRange( *( minmax_it.first ), *( minmax_it.second ) );
+    if ( m_qwtPlotCurve )
+    {
+        m_qwtPlotCurve->setSamples( xValues.data(), depthValues.data(), static_cast<int>( depthValues.size() ) );
+        m_qwtPlotCurve->setLineSegmentStartStopIndices( curveIntervals );
+    }
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimWellLogCurve.h

@@ -75,6 +75,10 @@ public:
 
     static QString wellLogCurveIconName();
 
+    void setOverrideCurveData( const std::vector<double>&               xValues,
+                               const std::vector<double>&               depthValues,
+                               const RiaCurveDataTools::CurveIntervals& curveIntervals );
+
 protected:
     void updateZoomInParentPlot() override;
     void updateLegendsInPlot() override;

ApplicationCode/ProjectDataModel/RimWellLogTrack.cpp

@@ -89,6 +89,7 @@
 #include <QWheelEvent>
 
 #include <algorithm>
+#include <set>
 
 #define RI_LOGPLOTTRACK_MINX_DEFAULT -10.0
 #define RI_LOGPLOTTRACK_MAXX_DEFAULT 100.0
@@ -185,6 +186,7 @@ RimWellLogTrack::RimWellLogTrack()
     m_visibleDepthRangeMax.uiCapability()->setUiHidden( true );
     m_visibleDepthRangeMax.xmlCapability()->disableIO();
 
+    CAF_PDM_InitField( &m_stackCurves, "StackCurves", false, "Stack Curves", "", "", "" );
     CAF_PDM_InitField( &m_isAutoScaleXEnabled, "AutoScaleX", true, "Auto Scale", "", "", "" );
     m_isAutoScaleXEnabled.uiCapability()->setUiHidden( true );
 
@@ -357,12 +359,7 @@ void RimWellLogTrack::detachAllPlotItems()
 //--------------------------------------------------------------------------------------------------
 void RimWellLogTrack::calculateXZoomRange()
 {
-    std::map<int, std::vector<RimWellFlowRateCurve*>> stackCurveGroups = visibleStackedCurves();
+    updateStackedCurveData();
-    for ( const std::pair<int, std::vector<RimWellFlowRateCurve*>>& curveGroup : stackCurveGroups )
-    {
-        for ( RimWellFlowRateCurve* stCurve : curveGroup.second )
-            stCurve->updateStackedPlotData();
-    }
 
     double minValue = HUGE_VAL;
     double maxValue = -HUGE_VAL;
@@ -2002,18 +1999,22 @@ void RimWellLogTrack::setLogarithmicScale( bool enable )
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::map<int, std::vector<RimWellFlowRateCurve*>> RimWellLogTrack::visibleStackedCurves()
+std::map<int, std::vector<RimWellLogCurve*>> RimWellLogTrack::visibleStackedCurves()
 {
-    std::map<int, std::vector<RimWellFlowRateCurve*>> stackedCurves;
+    std::map<int, std::vector<RimWellLogCurve*>> stackedCurves;
     for ( RimWellLogCurve* curve : m_curves )
     {
         if ( curve && curve->isCurveVisible() )
         {
             RimWellFlowRateCurve* wfrCurve = dynamic_cast<RimWellFlowRateCurve*>( curve );
-            if ( wfrCurve != nullptr )
+            if ( wfrCurve != nullptr ) // Flow rate curves are always stacked
             {
                 stackedCurves[wfrCurve->groupId()].push_back( wfrCurve );
             }
+            else if ( m_stackCurves() )
+            {
+                stackedCurves[-1].push_back( curve );
+            }
         }
     }
 
@@ -2252,6 +2253,91 @@ std::vector<QString> RimWellLogTrack::formationNamesVector( RimCase* rimCase )
     return std::vector<QString>();
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimWellLogTrack::updateStackedCurveData()
+{
+    RimWellLogPlot* wellLogPlot;
+    firstAncestorOrThisOfTypeAsserted( wellLogPlot );
+
+    RimWellLogPlot::DepthTypeEnum depthType    = wellLogPlot->depthType();
+    RiaDefines::DepthUnitType     displayUnit  = wellLogPlot->depthUnit();
+    bool                          reverseOrder = false;
+    if ( depthType == RiaDefines::DepthTypeEnum::CONNECTION_NUMBER )
+    {
+        displayUnit  = RiaDefines::DepthUnitType::UNIT_NONE;
+        reverseOrder = true;
+    }
+
+    std::map<int, std::vector<RimWellLogCurve*>> stackedCurves = visibleStackedCurves();
+
+    for ( auto groupCurvePair : stackedCurves )
+    {
+        int                                  groupId              = groupCurvePair.first;
+        const std::vector<RimWellLogCurve*>& stackedCurvesInGroup = groupCurvePair.second;
+        if ( stackedCurvesInGroup.empty() ) continue;
+
+        // Z-position of curve, to draw them in correct order
+        double zPos = -10000.0 + 100.0 * static_cast<double>( groupId );
+
+        // Find common depths. We retain all depths from the first curve and insert ones that aren't already added.
+        std::vector<double> allDepthValues;
+
+        for ( auto curve : stackedCurvesInGroup )
+        {
+            auto depths = curve->curveData()->depths( depthType );
+            if ( allDepthValues.empty() )
+            {
+                allDepthValues.insert( allDepthValues.end(), depths.begin(), depths.end() );
+            }
+            else
+            {
+                for ( double depth : depths )
+                {
+                    // Finds the first larger or equal depth.
+                    auto it = std::lower_bound( allDepthValues.begin(),
+                                                allDepthValues.end(),
+                                                depth,
+                                                [reverseOrder]( double lhs, double rhs ) {
+                                                    return reverseOrder ? rhs < lhs : lhs < rhs;
+                                                } );
+
+                    // Insert if there is no larger or equal depths or if the first equal or if it actually larger.
+                    if ( it == allDepthValues.end() || std::fabs( depth - *it ) > 1.0e-8 )
+                    {
+                        allDepthValues.insert( it, depth );
+                    }
+                }
+            }
+        }
+        if ( allDepthValues.empty() ) continue;
+
+        std::vector<double> allStackedValues( allDepthValues.size(), 0.0 );
+        for ( auto curve : stackedCurvesInGroup )
+        {
+            auto interpolatedCurveValues = curve->curveData()->calculateResampledCurveData( depthType, allDepthValues );
+            auto xValues                 = interpolatedCurveValues->xValues();
+            for ( size_t i = 0; i < xValues.size(); ++i )
+            {
+                if ( xValues[i] != HUGE_VAL )
+                {
+                    allStackedValues[i] += xValues[i];
+                }
+            }
+
+            RigWellLogCurveData tempCurveData;
+            tempCurveData.setValuesAndDepths( allStackedValues, allDepthValues, depthType, 0.0, displayUnit, false );
+            auto plotDepthValues          = tempCurveData.depthPlotValues( depthType, displayUnit );
+            auto polyLineStartStopIndices = tempCurveData.polylineStartStopIndices();
+
+            curve->setOverrideCurveData( allStackedValues, plotDepthValues, polyLineStartStopIndices );
+            curve->setZOrder( zPos );
+            zPos -= 1.0;
+        }
+    }
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimWellLogTrack.h

@@ -191,7 +191,7 @@ public:
 
     void setLogarithmicScale( bool enable );
 
-    std::map<int, std::vector<RimWellFlowRateCurve*>> visibleStackedCurves();
+    std::map<int, std::vector<RimWellLogCurve*>> visibleStackedCurves();
 
     std::vector<RimWellLogCurve*> curves() const;
     std::vector<RimWellLogCurve*> visibleCurves() const;
@@ -224,6 +224,7 @@ public:
                                        std::vector<std::pair<double, double>>* yValues );
 
     static std::vector<QString> formationNamesVector( RimCase* rimCase );
+    void                        updateStackedCurveData();
 
     static void addOverburden( std::vector<QString>& namesVector, CurveSamplingPointData& curveData, double height );
     static void addUnderburden( std::vector<QString>& namesVector, CurveSamplingPointData& curveData, double height );
@@ -304,6 +305,7 @@ private:
     caf::PdmField<double>                     m_visibleDepthRangeMin;
     caf::PdmField<double>                     m_visibleDepthRangeMax;
 
+    caf::PdmField<bool>                         m_stackCurves;
     caf::PdmField<bool>                         m_isAutoScaleXEnabled;
     caf::PdmField<bool>                         m_isLogarithmicScaleEnabled;
     caf::PdmField<RimWellLogPlot::AxisGridEnum> m_xAxisGridVisibility;

ApplicationCode/ReservoirDataModel/RigWellLogCurveData.cpp

@@ -325,6 +325,124 @@ cvf::ref<RigWellLogCurveData> RigWellLogCurveData::calculateResampledCurveData(
     return reSampledData;
 }
 
+void RigWellLogCurveData::interpolateSegment( RiaDefines::DepthTypeEnum resamplingDepthType,
+                                              double                    depthValue,
+                                              size_t                    firstIndex,
+                                              std::vector<double>&      xValues,
+                                              std::map<RiaDefines::DepthTypeEnum, std::vector<double>>& resampledDepths,
+                                              const double                                              eps ) const
+{
+    auto depthIt = m_depths.find( resamplingDepthType );
+
+    size_t secondIndex = firstIndex + 1;
+
+    if ( secondIndex >= depthIt->second.size() ) return;
+
+    double depth0 = depthIt->second[firstIndex];
+    double depth1 = depthIt->second[secondIndex];
+    double x0     = m_xValues[firstIndex];
+    double x1     = m_xValues[secondIndex];
+    double slope  = 0.0;
+    if ( std::fabs( depth1 - depth0 ) > eps )
+    {
+        slope = ( x1 - x0 ) / ( depth1 - depth0 );
+    }
+    double xValue = slope * ( depthValue - depth0 ) + x0;
+    xValues.push_back( xValue );
+
+    for ( auto depthTypeValuesPair : m_depths )
+    {
+        if ( depthTypeValuesPair.first != resamplingDepthType )
+        {
+            double otherDepth0 = depthTypeValuesPair.second[0];
+            double otherDepth1 = depthTypeValuesPair.second[1];
+            double otherSlope  = ( otherDepth1 - otherDepth0 ) / ( depth1 - depth0 );
+            resampledDepths[depthTypeValuesPair.first].push_back( otherSlope * ( depthValue - depth0 ) + otherDepth0 );
+        }
+    }
+}
+
+bool isLeftOf( double x1, double x2, bool reverseOrder, double eps )
+{
+    if ( reverseOrder )
+    {
+        return x1 - x2 > eps;
+    }
+    return x2 - x1 > eps;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Assumes the data is well ordered
+//--------------------------------------------------------------------------------------------------
+cvf::ref<RigWellLogCurveData> RigWellLogCurveData::calculateResampledCurveData( RiaDefines::DepthTypeEnum resamplingDepthType,
+                                                                                const std::vector<double>& depths ) const
+{
+    const double eps = 1.0e-8;
+
+    std::vector<double> xValues;
+
+    bool                                                     isTVDAvailable = false;
+    std::map<RiaDefines::DepthTypeEnum, std::vector<double>> resampledDepths;
+
+    auto depthIt = m_depths.find( resamplingDepthType );
+
+    cvf::ref<RigWellLogCurveData> reSampledData = new RigWellLogCurveData;
+
+    if ( depthIt->second.empty() ) return reSampledData;
+
+    bool reverseOrder = resamplingDepthType == RiaDefines::DepthTypeEnum::CONNECTION_NUMBER;
+
+    size_t segmentSearchStartIdx = 0;
+    for ( double depth : depths )
+    {
+        if ( isLeftOf( depth, depthIt->second.front(), reverseOrder, eps ) )
+        {
+            // Extrapolate from front two
+            interpolateSegment( resamplingDepthType, depth, 0, xValues, resampledDepths, eps );
+        }
+        else if ( isLeftOf( depthIt->second.back(), depth, reverseOrder, eps ) )
+        {
+            // Extrapolate from end two
+            const size_t N = depthIt->second.size() - 1;
+            interpolateSegment( resamplingDepthType, depth, N - 1, xValues, resampledDepths, eps );
+        }
+        else
+        {
+            for ( size_t segmentStartIdx = segmentSearchStartIdx; segmentStartIdx < depthIt->second.size();
+                  ++segmentStartIdx )
+            {
+                if ( std::fabs( depthIt->second[segmentStartIdx] - depth ) < eps ) // already have this depth point,
+                                                                                   // reuse it
+                {
+                    xValues.push_back( m_xValues[segmentStartIdx] );
+                    for ( auto depthTypeValuesPair : m_depths )
+                    {
+                        if ( depthTypeValuesPair.first != resamplingDepthType )
+                        {
+                            resampledDepths[depthTypeValuesPair.first].push_back(
+                                depthTypeValuesPair.second[segmentStartIdx] );
+                        }
+                    }
+                    segmentSearchStartIdx = segmentStartIdx + 1;
+                }
+                else if ( segmentStartIdx < depthIt->second.size() - 1 &&
+                          isLeftOf( depthIt->second[segmentStartIdx], depth, reverseOrder, eps ) &&
+                          isLeftOf( depth, depthIt->second[segmentStartIdx + 1], reverseOrder, eps ) ) // interpolate
+                                                                                                       // between two
+                                                                                                       // closest points
+                {
+                    interpolateSegment( resamplingDepthType, depth, segmentStartIdx, xValues, resampledDepths, eps );
+                    segmentSearchStartIdx = segmentStartIdx;
+                }
+            }
+        }
+    }
+
+    resampledDepths.insert( std::make_pair( resamplingDepthType, depths ) );
+    reSampledData->setValuesAndDepths( xValues, resampledDepths, m_rkbDiff, m_depthUnit, true );
+    return reSampledData;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -335,7 +453,7 @@ void RigWellLogCurveData::calculateIntervalsOfContinousValidValues()
 
     m_intervalsOfContinousValidValues.clear();
 
-    if ( !m_isExtractionCurve )
+    if ( !m_isExtractionCurve || !m_depths.count( RiaDefines::DepthTypeEnum::MEASURED_DEPTH ) )
     {
         m_intervalsOfContinousValidValues = intervalsOfValidValues;
     }

ApplicationCode/ReservoirDataModel/RigWellLogCurveData.h

@@ -75,6 +75,14 @@ public:
     std::vector<std::pair<size_t, size_t>> polylineStartStopIndices() const;
 
     cvf::ref<RigWellLogCurveData> calculateResampledCurveData( double newMeasuredDepthStepSize ) const;
+    cvf::ref<RigWellLogCurveData> calculateResampledCurveData( RiaDefines::DepthTypeEnum  resamplingDepthType,
+                                                               const std::vector<double>& depths ) const;
+    void                          interpolateSegment( RiaDefines::DepthTypeEnum                                 resamplingDepthType,
+                                                      double                                                    depthValue,
+                                                      size_t                                                    firstIndex,
+                                                      std::vector<double>&                                      xValues,
+                                                      std::map<RiaDefines::DepthTypeEnum, std::vector<double>>& resampledDepths,
+                                                      const double                                              eps ) const;
 
 private:
     void calculateIntervalsOfContinousValidValues();