Ensemble surface import and statistics

ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake

@@ -89,6 +89,8 @@ ${CMAKE_CURRENT_LIST_DIR}/RigTracer.h
 ${CMAKE_CURRENT_LIST_DIR}/RigStimPlanModelTools.h
 ${CMAKE_CURRENT_LIST_DIR}/RigSlice2D.h
 ${CMAKE_CURRENT_LIST_DIR}/RigEnsembleParameter.h
+${CMAKE_CURRENT_LIST_DIR}/RigSurfaceResampler.h
+${CMAKE_CURRENT_LIST_DIR}/RigSurfaceStatisticsCalculator.h
 )
 
 
@@ -175,6 +177,8 @@ ${CMAKE_CURRENT_LIST_DIR}/RigTracer.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigStimPlanModelTools.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigSlice2D.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigEnsembleParameter.cpp
+${CMAKE_CURRENT_LIST_DIR}/RigSurfaceResampler.cpp
+${CMAKE_CURRENT_LIST_DIR}/RigSurfaceStatisticsCalculator.cpp
 )
 
 list(APPEND CODE_HEADER_FILES

ApplicationLibCode/ReservoirDataModel/RigSurface.cpp

@@ -34,7 +34,7 @@ RigSurface::~RigSurface()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-const std::vector<unsigned>& RigSurface::triangleIndices()
+const std::vector<unsigned>& RigSurface::triangleIndices() const
 {
     return m_triangleIndices;
 }
@@ -42,7 +42,7 @@ const std::vector<unsigned>& RigSurface::triangleIndices()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-const std::vector<cvf::Vec3d>& RigSurface::vertices()
+const std::vector<cvf::Vec3d>& RigSurface::vertices() const
 {
     return m_vertices;
 }

ApplicationLibCode/ReservoirDataModel/RigSurface.h

@@ -31,8 +31,8 @@ public:
     RigSurface();
     ~RigSurface() override;
 
-    const std::vector<unsigned>&   triangleIndices();
-    const std::vector<cvf::Vec3d>& vertices();
+    const std::vector<unsigned>&   triangleIndices() const;
+    const std::vector<cvf::Vec3d>& vertices() const;
 
     void setTriangleData( const std::vector<unsigned>& tringleIndices, const std::vector<cvf::Vec3d>& vertices );
     void addVerticeResult( const QString resultName, const std::vector<float>& resultValues );

ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp

@@ -0,0 +1,119 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2021-     Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+#include "RigSurfaceResampler.h"
+
+#include "cvfGeometryTools.h"
+
+#include "cvfObject.h"
+#include <limits>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::ref<RigSurface> RigSurfaceResampler::resampleSurface( cvf::ref<RigSurface> targetSurface, cvf::ref<RigSurface> surface )
+{
+    cvf::ref<RigSurface> resampledSurface = cvf::make_ref<RigSurface>();
+
+    const std::vector<cvf::Vec3d>& targetVerts   = targetSurface->vertices();
+    const std::vector<unsigned>&   targetIndices = targetSurface->triangleIndices();
+
+    std::vector<cvf::Vec3d> resampledVerts;
+
+    for ( 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 );
+
+        cvf::Vec3d intersectionPoint;
+        bool       foundMatch =
+            resamplePoint( pointAbove, pointBelow, surface->triangleIndices(), surface->vertices(), intersectionPoint );
+        if ( !foundMatch )
+            intersectionPoint = cvf::Vec3d( targetVert.x(), targetVert.y(), std::numeric_limits<double>::infinity() );
+
+        resampledVerts.push_back( intersectionPoint );
+    }
+
+    resampledSurface->setTriangleData( targetIndices, resampledVerts );
+
+    return resampledSurface;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RigSurfaceResampler::resamplePoint( const cvf::Vec3d&                pointAbove,
+                                         const cvf::Vec3d&                pointBelow,
+                                         const std::vector<unsigned int>& indices,
+                                         const std::vector<cvf::Vec3d>&   vertices,
+                                         cvf::Vec3d&                      intersectionPoint )
+{
+    for ( size_t i = 0; i < indices.size(); i += 3 )
+    {
+        bool isLineDirDotNormalNegative = false;
+        if ( cvf::GeometryTools::intersectLineSegmentTriangle( pointAbove,
+                                                               pointBelow,
+                                                               vertices[indices[i]],
+                                                               vertices[indices[i + 1]],
+                                                               vertices[indices[i + 2]],
+                                                               &intersectionPoint,
+                                                               &isLineDirDotNormalNegative ) == 1 )
+            return true;
+    }
+
+    // Handle cases where no match is found due to floating point imprecision,
+    // or when falling off resulting grid slightly.
+    double maxDistance = 10.0;
+    return findClosestPointXY( pointAbove, vertices, maxDistance, intersectionPoint );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Find the closest vertex to targetPoint (must be closer than maxDistance) in XY plane.
+/// Unit maxDistance: meter.
+//--------------------------------------------------------------------------------------------------
+bool RigSurfaceResampler::findClosestPointXY( const cvf::Vec3d&              targetPoint,
+                                              const std::vector<cvf::Vec3d>& vertices,
+                                              double                         maxDistance,
+                                              cvf::Vec3d&                    intersectionPoint )
+{
+    // Find closest vertices
+    double shortestDistance = std::numeric_limits<double>::max();
+    double closestZ         = std::numeric_limits<double>::infinity();
+    for ( auto v : vertices )
+    {
+        // Ignore height (z) component when finding closest by
+        // moving point to same height as target point above
+        cvf::Vec3d p( v.x(), v.y(), targetPoint.z() );
+        double     distance = p.pointDistance( targetPoint );
+        if ( distance < shortestDistance )
+        {
+            shortestDistance = distance;
+            closestZ         = v.z();
+        }
+    }
+
+    // Check if the closest point is not to far away to be valid
+    if ( shortestDistance < maxDistance )
+    {
+        intersectionPoint = cvf::Vec3d( targetPoint.x(), targetPoint.y(), closestZ );
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}

ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h

@@ -0,0 +1,43 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2021-     Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+#pragma once
+
+#include "RigSurface.h"
+
+#include "cvfObject.h"
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+class RigSurfaceResampler
+{
+public:
+    static cvf::ref<RigSurface> resampleSurface( cvf::ref<RigSurface> targetSurface, cvf::ref<RigSurface> surface );
+
+    static bool resamplePoint( const cvf::Vec3d&                pointAbove,
+                               const cvf::Vec3d&                pointBelow,
+                               const std::vector<unsigned int>& indices,
+                               const std::vector<cvf::Vec3d>&   vertices,
+                               cvf::Vec3d&                      intersectionPoint );
+
+private:
+    static bool findClosestPointXY( const cvf::Vec3d&              targetPoint,
+                                    const std::vector<cvf::Vec3d>& vertices,
+                                    double                         maxDistance,
+                                    cvf::Vec3d&                    intersectionPoint );
+};

ApplicationLibCode/ReservoirDataModel/RigSurfaceStatisticsCalculator.cpp

@@ -0,0 +1,129 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2021 Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "RigSurfaceStatisticsCalculator.h"
+
+#include "RigStatisticsMath.h"
+
+#include "cafAppEnum.h"
+#include "cafAssert.h"
+
+#include "cvfObject.h"
+
+#include <limits>
+
+namespace caf
+{
+template <>
+void caf::AppEnum<RigSurfaceStatisticsCalculator::StatisticsType>::setUp()
+{
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::MEAN, "MEAN", "Mean" );
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::MIN, "MIN", "Minimum" );
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::MAX, "MAX", "Maximum" );
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::P10, "P10", "P10" );
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::P50, "P50", "P50" );
+    addItem( RigSurfaceStatisticsCalculator::StatisticsType::P90, "P90", "P90" );
+    setDefault( RigSurfaceStatisticsCalculator::StatisticsType::MEAN );
+}
+}; // namespace caf
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::ref<RigSurface> RigSurfaceStatisticsCalculator::computeStatistics( const std::vector<cvf::ref<RigSurface>>& surfaces )
+{
+    // Need at least one surface to generate statistics
+    if ( surfaces.empty() ) return nullptr;
+
+    // Check that the size of all the surfaces are the same
+    if ( !areSurfacesSameSize( surfaces ) ) return nullptr;
+
+    size_t vecSize = surfaces[0]->vertices().size();
+
+    std::vector<float> meanValues( vecSize, std::numeric_limits<double>::infinity() );
+    std::vector<float> minValues( vecSize, std::numeric_limits<double>::infinity() );
+    std::vector<float> maxValues( vecSize, std::numeric_limits<double>::infinity() );
+    std::vector<float> p10Values( vecSize, std::numeric_limits<double>::infinity() );
+    std::vector<float> p50Values( vecSize, std::numeric_limits<double>::infinity() );
+    std::vector<float> p90Values( vecSize, std::numeric_limits<double>::infinity() );
+
+    for ( size_t i = 0; i < vecSize; i++ )
+    {
+        std::vector<double> samples;
+        for ( auto surface : surfaces )
+        {
+            double z = surface->vertices()[i].z();
+            samples.push_back( z );
+        }
+
+        double min;
+        double max;
+        double sum;
+        double range;
+        double mean;
+        double dev;
+        RigStatisticsMath::calculateBasicStatistics( samples, &min, &max, &sum, &range, &mean, &dev );
+
+        double p10;
+        double p50;
+        double p90;
+        RigStatisticsMath::calculateStatisticsCurves( samples, &p10, &p50, &p90, &mean );
+
+        // TODO: improve handling of these cases
+        auto makeValid = []( double val ) {
+            if ( std::isinf( val ) || std::isnan( val ) ) return 0.0;
+            return val;
+        };
+
+        meanValues[i] = makeValid( mean );
+        minValues[i]  = makeValid( min );
+        maxValues[i]  = makeValid( max );
+        p10Values[i]  = makeValid( p10 );
+        p50Values[i]  = makeValid( p50 );
+        p90Values[i]  = makeValid( p90 );
+    }
+
+    cvf::ref<RigSurface> statSurface = cvf::make_ref<RigSurface>();
+    statSurface->setTriangleData( surfaces[0]->triangleIndices(), surfaces[0]->vertices() );
+
+    auto enumToText = []( auto statisticsType ) { return caf::AppEnum<StatisticsType>::text( statisticsType ); };
+
+    statSurface->addVerticeResult( enumToText( StatisticsType::MEAN ), meanValues );
+    statSurface->addVerticeResult( enumToText( StatisticsType::MIN ), minValues );
+    statSurface->addVerticeResult( enumToText( StatisticsType::MAX ), maxValues );
+    statSurface->addVerticeResult( enumToText( StatisticsType::P10 ), p10Values );
+    statSurface->addVerticeResult( enumToText( StatisticsType::P50 ), p50Values );
+    statSurface->addVerticeResult( enumToText( StatisticsType::P90 ), p90Values );
+
+    return statSurface;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RigSurfaceStatisticsCalculator::areSurfacesSameSize( const std::vector<cvf::ref<RigSurface>>& surfaces )
+{
+    CAF_ASSERT( !surfaces.empty() );
+
+    size_t vecSize = surfaces[0]->vertices().size();
+    for ( auto surface : surfaces )
+    {
+        if ( vecSize != surface->vertices().size() ) return false;
+    }
+    return true;
+}

ApplicationLibCode/ReservoirDataModel/RigSurfaceStatisticsCalculator.h

@@ -0,0 +1,46 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2021 Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#pragma once
+
+#include "RigSurface.h"
+
+#include "cvfObject.h"
+
+#include <vector>
+
+//==================================================================================================
+///
+///
+//==================================================================================================
+class RigSurfaceStatisticsCalculator
+{
+public:
+    enum class StatisticsType
+    {
+        MIN,
+        MAX,
+        P10,
+        P50,
+        P90,
+        MEAN
+    };
+
+    static cvf::ref<RigSurface> computeStatistics( const std::vector<cvf::ref<RigSurface>>& surfaces );
+    static bool                 areSurfacesSameSize( const std::vector<cvf::ref<RigSurface>>& surfaces );
+};

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -515,13 +515,13 @@ double GeometryTools::linePointSquareDist( const cvf::Vec3d& p1, const cvf::Vec3
 // dot product (3D) which allows vector operations in arguments
 #define dot( u, v ) ( ( u ).x() * ( v ).x() + ( u ).y() * ( v ).y() + ( u ).z() * ( v ).z() )
 
-int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d p0,
-                                                 const cvf::Vec3d p1,
-                                                 const cvf::Vec3d t0,
-                                                 const cvf::Vec3d t1,
-                                                 const cvf::Vec3d t2,
-                                                 cvf::Vec3d*      intersectionPoint,
-                                                 bool*            isLineDirDotNormalNegative )
+int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d& p0,
+                                                 const cvf::Vec3d& p1,
+                                                 const cvf::Vec3d& t0,
+                                                 const cvf::Vec3d& t1,
+                                                 const cvf::Vec3d& t2,
+                                                 cvf::Vec3d*       intersectionPoint,
+                                                 bool*             isLineDirDotNormalNegative )
 {
     CVF_TIGHT_ASSERT( intersectionPoint != nullptr );
     CVF_TIGHT_ASSERT( isLineDirDotNormalNegative != nullptr );

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h

@@ -47,13 +47,13 @@ public:
                                           double*           normalizedIntersection );
 
     static double linePointSquareDist( const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3 );
-    static int    intersectLineSegmentTriangle( const cvf::Vec3d p0,
-                                                const cvf::Vec3d p1,
-                                                const cvf::Vec3d t0,
-                                                const cvf::Vec3d t1,
-                                                const cvf::Vec3d t2,
-                                                cvf::Vec3d*      intersectionPoint,
-                                                bool*            isLineDirDotNormalNegative );
+    static int    intersectLineSegmentTriangle( const cvf::Vec3d& p0,
+                                                const cvf::Vec3d& p1,
+                                                const cvf::Vec3d& t0,
+                                                const cvf::Vec3d& t1,
+                                                const cvf::Vec3d& t2,
+                                                cvf::Vec3d*       intersectionPoint,
+                                                bool*             isLineDirDotNormalNegative );
     static cvf::Vec3d
                       barycentricCoords( const cvf::Vec3d& t0, const cvf::Vec3d& t1, const cvf::Vec3d& t2, const cvf::Vec3d& p );
     static cvf::Vec4d barycentricCoords( const cvf::Vec3d& v0,