diff --git a/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.cpp b/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.cpp index dcfcc2de10..f255c2b640 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.cpp @@ -345,41 +345,52 @@ cvf::ref RigWellLogCurveData::calculateResampledCurveData( return reSampledData; } -void RigWellLogCurveData::interpolateSegment( RiaDefines::DepthTypeEnum resamplingDepthType, - double depthValue, - size_t firstIndex, - std::vector& xValues, - std::map>& resampledDepths, - const double eps ) const +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigWellLogCurveData::interpolateSegment( RiaDefines::DepthTypeEnum resamplingDepthType, + std::vector& resampledValues, + std::map>& resampledDepths, + double targetDepthValue, + size_t firstIndex, + const std::map>& originalDepths, + const std::vector& propertyValues, + double eps ) { - auto depthIt = m_depths.find( resamplingDepthType ); + if ( !originalDepths.contains( resamplingDepthType ) ) return; - size_t secondIndex = firstIndex + 1; + const size_t secondIndex = firstIndex + 1; + const auto& depthValues = originalDepths.find( resamplingDepthType )->second; + if ( secondIndex >= depthValues.size() ) return; - double depth0 = depthIt->second[firstIndex]; - double depth1 = depthIt->second[secondIndex]; - double x0 = m_propertyValues[firstIndex]; - double x1 = m_propertyValues[secondIndex]; - double slope = 0.0; + const double depth0 = depthValues[firstIndex]; + const double depth1 = depthValues[secondIndex]; + const double x0 = propertyValues[firstIndex]; + const double x1 = propertyValues[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 ); + const double resampledValue = slope * ( targetDepthValue - depth0 ) + x0; + resampledValues.push_back( resampledValue ); - for ( auto depthTypeValuesPair : m_depths ) + for ( const auto& [depthType, depthTypeValues] : originalDepths ) { - 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 ); - } + // Skip the depth type we are resampling and ensure depth values are available + if ( depthType == resamplingDepthType ) continue; + if ( depthTypeValues.size() < secondIndex - 1 ) continue; + + const double otherDepth0 = depthTypeValues[firstIndex]; + const double otherDepth1 = depthTypeValues[secondIndex]; + const double otherSlope = ( otherDepth1 - otherDepth0 ) / ( depth1 - depth0 ); + resampledDepths[depthType].push_back( otherSlope * ( targetDepthValue - depth0 ) + otherDepth0 ); } } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- bool isLeftOf( double x1, double x2, bool reverseOrder, double eps ) { if ( reverseOrder ) @@ -389,6 +400,9 @@ bool isLeftOf( double x1, double x2, bool reverseOrder, double eps ) return x2 - x1 > eps; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- bool isRightOf( double x1, double x2, bool reverseOrder, double eps ) { return isLeftOf( x2, x1, reverseOrder, eps ); @@ -397,53 +411,61 @@ bool isRightOf( double x1, double x2, bool reverseOrder, double eps ) //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -cvf::ref RigWellLogCurveData::calculateResampledCurveData( RiaDefines::DepthTypeEnum resamplingDepthType, - const std::vector& depths ) const +std::pair, std::map>> + RigWellLogCurveData::createResampledValuesAndDepths( RiaDefines::DepthTypeEnum resamplingDepthType, + const std::vector& targetDepths, + const std::map>& originalDepths, + const std::vector& propertyValues ) { const double eps = 1.0e-8; - std::vector xValues; + auto depthIt = originalDepths.find( resamplingDepthType ); + if ( depthIt == originalDepths.end() || depthIt->second.empty() ) return {}; + const auto& depthValues = depthIt->second; + std::vector resampledValues; std::map> resampledDepths; - resampledDepths.insert( std::make_pair( resamplingDepthType, depths ) ); + resampledDepths.insert( std::make_pair( resamplingDepthType, targetDepths ) ); - auto depthIt = m_depths.find( resamplingDepthType ); - - cvf::ref reSampledData = new RigWellLogCurveData; - - if ( depthIt == m_depths.end() || depthIt->second.empty() ) return reSampledData; + cvf::ref resampledCurveData = new RigWellLogCurveData; bool reverseOrder = resamplingDepthType == RiaDefines::DepthTypeEnum::CONNECTION_NUMBER; size_t segmentSearchStartIdx = 0; - for ( auto depth : depths ) + for ( const auto& targetDepth : targetDepths ) { bool foundPoint = false; - for ( size_t segmentStartIdx = segmentSearchStartIdx; segmentStartIdx < depthIt->second.size(); ++segmentStartIdx ) + for ( size_t segmentStartIdx = segmentSearchStartIdx; segmentStartIdx < depthValues.size(); ++segmentStartIdx ) { - if ( std::fabs( depthIt->second[segmentStartIdx] - depth ) < eps ) // already have this depth point, - // reuse it + if ( std::fabs( depthValues[segmentStartIdx] - targetDepth ) < eps ) // already have this depth point, reuse it { - xValues.push_back( m_propertyValues[segmentStartIdx] ); + resampledValues.push_back( propertyValues[segmentStartIdx] ); // Copy all depth types for this segment - for ( auto depthTypeValuesPair : m_depths ) + for ( const auto& depthTypeValuesPair : originalDepths ) { if ( depthTypeValuesPair.first != resamplingDepthType ) { resampledDepths[depthTypeValuesPair.first].push_back( depthTypeValuesPair.second[segmentStartIdx] ); } } - segmentSearchStartIdx = segmentStartIdx + 1; + segmentSearchStartIdx = segmentStartIdx; foundPoint = true; break; } - else if ( segmentStartIdx < depthIt->second.size() - 1 ) + else if ( segmentStartIdx < depthValues.size() - 1 ) { - double minDepthSegment = std::min( depthIt->second[segmentStartIdx], depthIt->second[segmentStartIdx + 1] ); - double maxDepthSegment = std::max( depthIt->second[segmentStartIdx], depthIt->second[segmentStartIdx + 1] ); - if ( cvf::Math::valueInRange( depth, minDepthSegment, maxDepthSegment ) ) + double minDepthSegment = std::min( depthValues[segmentStartIdx], depthValues[segmentStartIdx + 1] ); + double maxDepthSegment = std::max( depthValues[segmentStartIdx], depthValues[segmentStartIdx + 1] ); + if ( cvf::Math::valueInRange( targetDepth, minDepthSegment, maxDepthSegment ) ) { - interpolateSegment( resamplingDepthType, depth, segmentStartIdx, xValues, resampledDepths, eps ); + interpolateSegment( resamplingDepthType, + resampledValues, + resampledDepths, + targetDepth, + segmentStartIdx, + originalDepths, + propertyValues, + eps ); segmentSearchStartIdx = segmentStartIdx; foundPoint = true; break; @@ -452,17 +474,18 @@ cvf::ref RigWellLogCurveData::calculateResampledCurveData( } if ( !foundPoint ) { - if ( isLeftOf( depth, depthIt->second.front(), reverseOrder, eps ) ) + if ( isLeftOf( targetDepth, depthValues.front(), reverseOrder, eps ) ) { // Extrapolate from front two - interpolateSegment( resamplingDepthType, depth, 0, xValues, resampledDepths, eps ); + const size_t firstIndex = 0; + interpolateSegment( resamplingDepthType, resampledValues, resampledDepths, targetDepth, firstIndex, originalDepths, propertyValues, eps ); foundPoint = true; } - else if ( isRightOf( depth, depthIt->second.back(), reverseOrder, eps ) ) + else if ( isRightOf( targetDepth, depthValues.back(), reverseOrder, eps ) ) { // Extrapolate from end two - const size_t N = depthIt->second.size() - 1; - interpolateSegment( resamplingDepthType, depth, N - 1, xValues, resampledDepths, eps ); + const size_t N = depthValues.size() - 1; + interpolateSegment( resamplingDepthType, resampledValues, resampledDepths, targetDepth, N - 1, originalDepths, propertyValues, eps ); foundPoint = true; } } @@ -470,6 +493,18 @@ cvf::ref RigWellLogCurveData::calculateResampledCurveData( CAF_ASSERT( foundPoint ); } + return std::make_pair( resampledValues, resampledDepths ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::ref RigWellLogCurveData::calculateResampledCurveData( RiaDefines::DepthTypeEnum resamplingDepthType, + const std::vector& depths ) const +{ + const auto [xValues, resampledDepths] = createResampledValuesAndDepths( resamplingDepthType, depths, m_depths, m_propertyValues ); + + cvf::ref reSampledData = new RigWellLogCurveData; reSampledData->setValuesAndDepths( xValues, resampledDepths, m_rkbDiff, m_depthUnit, true, m_useLogarithmicScale ); return reSampledData; } diff --git a/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.h b/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.h index 1e1dbbc923..ee0f03d70b 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.h +++ b/ApplicationLibCode/ReservoirDataModel/RigWellLogCurveData.h @@ -82,12 +82,20 @@ public: cvf::ref calculateResampledCurveData( double newMeasuredDepthStepSize ) const; cvf::ref calculateResampledCurveData( RiaDefines::DepthTypeEnum resamplingDepthType, const std::vector& depths ) const; - void interpolateSegment( RiaDefines::DepthTypeEnum resamplingDepthType, - double depthValue, - size_t firstIndex, - std::vector& xValues, - std::map>& resampledDepths, - const double eps ) const; + static void interpolateSegment( RiaDefines::DepthTypeEnum resamplingDepthType, + std::vector& resampledValues, + std::map>& resampledDepths, + double targetDepthValue, + size_t firstIndex, + const std::map>& originalDepths, + const std::vector& propertyValues, + double eps ); + + static std::pair, std::map>> + createResampledValuesAndDepths( RiaDefines::DepthTypeEnum resamplingDepthType, + const std::vector& targetDepths, + const std::map>& originalDepths, + const std::vector& propertyValues ); private: void calculateIntervalsOfContinousValidValues(); diff --git a/ApplicationLibCode/UnitTests/CMakeLists_files.cmake b/ApplicationLibCode/UnitTests/CMakeLists_files.cmake index d4611833ba..66bd8d8a35 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists_files.cmake +++ b/ApplicationLibCode/UnitTests/CMakeLists_files.cmake @@ -91,6 +91,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RifStimPlanCsvSummaryReader-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RiaEnsembleNameTools-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigDeclineCurveCalculator-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigWellLogCurveData-Test.cpp ) if(RESINSIGHT_ENABLE_GRPC) diff --git a/ApplicationLibCode/UnitTests/RigWellLogCurveData-Test.cpp b/ApplicationLibCode/UnitTests/RigWellLogCurveData-Test.cpp new file mode 100644 index 0000000000..4ece723a29 --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigWellLogCurveData-Test.cpp @@ -0,0 +1,183 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2023- 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RiaDefines.h" + +#include "RigWellLogCurveData.h" + +#include "cvfVector3.h" + +#include + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigWellLogCurveData, interpolateSegment_first ) +{ + // Input data + const std::map> originalDepths = + { { RiaDefines::DepthTypeEnum::MEASURED_DEPTH, { 0.0, 20.0, 40.0 } }, + { RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH, { 0.0, 30.0, 60.0 } }, + { RiaDefines::DepthTypeEnum::PSEUDO_LENGTH, { 0.0, 40.0, 80.0 } } }; + const std::vector propertyValues = { 0.0, 100.0, 150.0 }; + const double eps = 1e-6; + + // Output data + const auto resamplingDepthType = RiaDefines::DepthTypeEnum::MEASURED_DEPTH; + std::vector resampledValues; + std::map> resampledDepths; + + // Target data (resampling with MEASURED_DEPTH) + const double targetDepthValue = 10.0; // Halfway between index 0 and 1 for MEASURED_DEPTH in originalDepths + const size_t firstIndex = 0; + + // Call the function under test + RigWellLogCurveData::interpolateSegment( resamplingDepthType, + resampledValues, + resampledDepths, + targetDepthValue, + firstIndex, + originalDepths, + propertyValues, + eps ); + + // Check the results + ASSERT_EQ( resampledValues.size(), size_t( 1 ) ); + ASSERT_DOUBLE_EQ( resampledValues[0], 50.0 ); + + ASSERT_EQ( resampledDepths.size(), size_t( 2 ) ); + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH].size(), size_t( 1 ) ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][0], 15.0 ); + + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH].size(), size_t( 1 ) ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][0], 20.0 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigWellLogCurveData, interpolateSegment_second ) +{ + // Input data + const std::map> originalDepths = + { { RiaDefines::DepthTypeEnum::MEASURED_DEPTH, { 0.0, 20.0, 40.0 } }, + { RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH, { 0.0, 30.0, 60.0 } }, + { RiaDefines::DepthTypeEnum::PSEUDO_LENGTH, { 0.0, 40.0, 80.0 } } }; + const std::vector propertyValues = { 0.0, 100.0, 150.0 }; + const double eps = 1e-6; + + // Output data + const auto resamplingDepthType = RiaDefines::DepthTypeEnum::MEASURED_DEPTH; + std::vector resampledValues; + std::map> resampledDepths; + + // Target data (resampling with MEASURED_DEPTH) + const double firstTargetDepthValue = 10.0; // Halfway between index 0 and 1 for MEASURED_DEPTH in originalDepths + const double secondTargetDepthValue = 30.0; // Halfway between index 1 and 2 for MEASURED_DEPTH in originalDepths + const size_t firstIndex = 0; + const size_t secondIndex = 1; + + // Call the function under test with first index and target value + RigWellLogCurveData::interpolateSegment( resamplingDepthType, + resampledValues, + resampledDepths, + firstTargetDepthValue, + firstIndex, + originalDepths, + propertyValues, + eps ); + + // Call the function under test with second index and target value + RigWellLogCurveData::interpolateSegment( resamplingDepthType, + resampledValues, + resampledDepths, + secondTargetDepthValue, + secondIndex, + originalDepths, + propertyValues, + eps ); + + // Check the results + ASSERT_EQ( resampledValues.size(), size_t( 2 ) ); + ASSERT_DOUBLE_EQ( resampledValues[0], 50.0 ); + ASSERT_DOUBLE_EQ( resampledValues[1], 125.0 ); + + ASSERT_EQ( resampledDepths.size(), size_t( 2 ) ); + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH].size(), size_t( 2 ) ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][0], 15.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][1], 45.0 ); + + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH].size(), size_t( 2 ) ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][0], 20.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][1], 60.0 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigWellLogCurveData, CreateResampledValuesAndDepthsTest ) +{ + // Input data + RiaDefines::DepthTypeEnum resamplingDepthType = RiaDefines::DepthTypeEnum::MEASURED_DEPTH; + const std::vector targetDepths = { 0.0, 5.0, 10.0, 15.0 }; + const std::map> originalDepths = + { { RiaDefines::DepthTypeEnum::MEASURED_DEPTH, { 0.0, 10.0, 20.0 } }, + { RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH, { 0.0, 20.0, 40.0 } }, + { RiaDefines::DepthTypeEnum::PSEUDO_LENGTH, { 0.0, 30.0, 60.0 } } }; + const std::vector propertyValues = { 0.0, 100.0, 200.0 }; + + // Call the function under test + auto result = RigWellLogCurveData::createResampledValuesAndDepths( resamplingDepthType, targetDepths, originalDepths, propertyValues ); + + // Check the results + std::vector& resampledPropertyValues = result.first; + std::map>& resampledDepths = result.second; + const auto expectedSize = targetDepths.size(); + + ASSERT_EQ( resampledDepths.size(), originalDepths.size() ); + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::MEASURED_DEPTH].size(), expectedSize ); + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH].size(), expectedSize ); + ASSERT_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH].size(), expectedSize ); + + ASSERT_EQ( resampledPropertyValues.size(), expectedSize ); + + // Example assertions for the specific values + ASSERT_DOUBLE_EQ( resampledPropertyValues[0], 0.0 ); + ASSERT_DOUBLE_EQ( resampledPropertyValues[1], 50.0 ); + ASSERT_DOUBLE_EQ( resampledPropertyValues[2], 100.0 ); + ASSERT_DOUBLE_EQ( resampledPropertyValues[3], 150.0 ); + + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::MEASURED_DEPTH][0], 0.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::MEASURED_DEPTH][1], 5.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::MEASURED_DEPTH][2], 10.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::MEASURED_DEPTH][3], 15.0 ); + + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][0], 0.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][1], 10.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][2], 20.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH][3], 30.0 ); + + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][0], 0.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][1], 15.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][2], 30.0 ); + ASSERT_DOUBLE_EQ( resampledDepths[RiaDefines::DepthTypeEnum::PSEUDO_LENGTH][3], 45.0 ); +}