diff --git a/ApplicationLibCode/Commands/CMakeLists_files.cmake b/ApplicationLibCode/Commands/CMakeLists_files.cmake index a9316b1a07..51a088b7ee 100644 --- a/ApplicationLibCode/Commands/CMakeLists_files.cmake +++ b/ApplicationLibCode/Commands/CMakeLists_files.cmake @@ -99,6 +99,7 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RicNewWellTargetCandidatesGeneratorFeature.h ${CMAKE_CURRENT_LIST_DIR}/RicNewStatisticsContourMapFeature.h ${CMAKE_CURRENT_LIST_DIR}/RicNewStatisticsContourMapViewFeature.h + ${CMAKE_CURRENT_LIST_DIR}/RicCreateWellTargetClusterPolygonsFeature.h ) set(SOURCE_GROUP_SOURCE_FILES @@ -202,6 +203,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RicNewWellTargetCandidatesGeneratorFeature.cpp ${CMAKE_CURRENT_LIST_DIR}/RicNewStatisticsContourMapFeature.cpp ${CMAKE_CURRENT_LIST_DIR}/RicNewStatisticsContourMapViewFeature.cpp + ${CMAKE_CURRENT_LIST_DIR}/RicCreateWellTargetClusterPolygonsFeature.cpp ) if(RESINSIGHT_USE_QT_CHARTS) diff --git a/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.cpp b/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.cpp new file mode 100644 index 0000000000..3a9a1a6d58 --- /dev/null +++ b/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.cpp @@ -0,0 +1,396 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2024- 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 "RicCreateWellTargetClusterPolygonsFeature.h" + +#include "RiaLogging.h" + +#include "RigActiveCellInfo.h" +#include "RigCaseCellResultsData.h" +#include "RigCell.h" +#include "RigConvexHull.h" +#include "RigMainGrid.h" +#include "RigStatisticsMath.h" + +#include "Polygons/RimPolygon.h" +#include "Polygons/RimPolygonCollection.h" +#include "Rim3dOverlayInfoConfig.h" +#include "RimEclipseCase.h" +#include "RimEclipseCellColors.h" +#include "RimEclipseView.h" +#include "RimOilField.h" +#include "RimProject.h" +#include "RimWellTargetCandidatesGenerator.h" + +#include "cafSelectionManager.h" +#include "cvfMath.h" + +#include + +CAF_CMD_SOURCE_INIT( RicCreateWellTargetClusterPolygonsFeature, "RicCreateWellTargetClusterPolygonsFeature" ); + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RicCreateWellTargetClusterPolygonsFeature::onActionTriggered( bool isChecked ) +{ + if ( auto generator = caf::SelectionManager::instance()->selectedItemOfType() ) + { + if ( auto ensembleStatisticsCase = generator->ensembleStatisticsCase() ) + { + if ( !ensembleStatisticsCase->reservoirViews().empty() ) + { + auto view = ensembleStatisticsCase->reservoirViews()[0]; + createWellTargetClusterPolygons( ensembleStatisticsCase, view ); + } + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// Function to smooth the histogram using a simple moving average +//-------------------------------------------------------------------------------------------------- +std::vector RicCreateWellTargetClusterPolygonsFeature::smoothHistogram( const std::vector& histogram, size_t windowSize ) +{ + std::vector smoothed( histogram.size(), 0.0 ); + size_t halfWindow = windowSize / 2; + + for ( size_t i = 0; i < histogram.size(); ++i ) + { + double sum = 0.0; + size_t count = 0; + for ( size_t j = ( i > halfWindow ? i - halfWindow : 0 ); j <= i + halfWindow && j < histogram.size(); ++j ) + { + sum += histogram[j]; + ++count; + } + smoothed[i] = sum / count; + } + + return smoothed; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector> + RicCreateWellTargetClusterPolygonsFeature::findPeaksWithDynamicProminence( const std::vector& smoothedHistogram ) +{ + std::vector> peaks; + std::vector prominences; + + // Find all peaks and calculate their prominences + for ( size_t i = 1; i < smoothedHistogram.size() - 1; ++i ) + { + if ( smoothedHistogram[i] > smoothedHistogram[i - 1] && smoothedHistogram[i] > smoothedHistogram[i + 1] ) + { + double leftMin = smoothedHistogram[i]; + for ( size_t j = i; j > 0; --j ) + { + if ( smoothedHistogram[j] < leftMin ) leftMin = smoothedHistogram[j]; + if ( smoothedHistogram[j] > smoothedHistogram[i] ) break; + } + + double rightMin = smoothedHistogram[i]; + for ( size_t j = i; j < smoothedHistogram.size(); ++j ) + { + if ( smoothedHistogram[j] < rightMin ) rightMin = smoothedHistogram[j]; + if ( smoothedHistogram[j] > smoothedHistogram[i] ) break; + } + + double prominence = smoothedHistogram[i] - std::max( leftMin, rightMin ); + prominences.push_back( prominence ); + peaks.emplace_back( i, prominence ); + } + } + + // Determine dynamic threshold + double min; + double max; + double sum; + double range; + double mean; + double dev; + RigStatisticsMath::calculateBasicStatistics( prominences, &min, &max, &sum, &range, &mean, &dev ); + double threshold = mean + dev; + + // Filter peaks based on the dynamic threshold + std::vector> filteredPeaks; + for ( const auto& [index, prominence] : peaks ) + { + if ( prominence >= threshold ) + { + filteredPeaks.emplace_back( index, prominence ); + } + } + + return filteredPeaks; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RicCreateWellTargetClusterPolygonsFeature::findCellCentersWithResultInRange( RimEclipseCase& eclipseCase, + const RigEclipseResultAddress& resultAddress, + size_t timeStepIndex, + double targetValue, + double rangeMinimum, + double rangeMaximum, + cvf::ref visibility, + std::vector& clusters, + int clusterId ) +{ + RigCaseCellResultsData* resultsData = eclipseCase.results( RiaDefines::PorosityModelType::MATRIX_MODEL ); + const RigMainGrid* mainGrid = eclipseCase.mainGrid(); + const RigActiveCellInfo* activeCellInfo = resultsData->activeCellInfo(); + + resultsData->ensureKnownResultLoaded( resultAddress ); + const std::vector& values = resultsData->cellScalarResults( resultAddress, timeStepIndex ); + + size_t startCell = findStartCell( *activeCellInfo, values, targetValue, rangeMinimum, rangeMaximum, visibility, clusters ); + + std::vector cellCenters; + + if ( startCell != std::numeric_limits::max() ) + { + growCluster( &eclipseCase, startCell, values, rangeMinimum, rangeMaximum, visibility, clusters, clusterId ); + + const size_t numReservoirCells = activeCellInfo->reservoirCellCount(); + + for ( size_t reservoirCellIndex = 0; reservoirCellIndex < numReservoirCells; reservoirCellIndex++ ) + { + size_t targetResultIndex = activeCellInfo->cellResultIndex( reservoirCellIndex ); + if ( reservoirCellIndex != cvf::UNDEFINED_SIZE_T && activeCellInfo->isActive( reservoirCellIndex ) && + targetResultIndex != cvf::UNDEFINED_SIZE_T && clusters[targetResultIndex] == clusterId ) + { + const RigCell& nativeCell = mainGrid->cell( reservoirCellIndex ); + cellCenters.push_back( nativeCell.center() ); + } + } + } + + return cellCenters; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RicCreateWellTargetClusterPolygonsFeature::findStartCell( const RigActiveCellInfo& activeCellInfo, + const std::vector& values, + double targetValue, + double minRange, + double maxRange, + cvf::ref visibility, + const std::vector& clusters ) +{ + size_t startCell = std::numeric_limits::max(); + double minDiff = std::numeric_limits::max(); + const size_t numReservoirCells = activeCellInfo.reservoirCellCount(); + for ( size_t reservoirCellIdx = 0; reservoirCellIdx < numReservoirCells; reservoirCellIdx++ ) + { + size_t resultIndex = activeCellInfo.cellResultIndex( reservoirCellIdx ); + if ( resultIndex != cvf::UNDEFINED_SIZE_T && clusters[resultIndex] == 0 && visibility->val( reservoirCellIdx ) ) + { + const double value = values[resultIndex]; + if ( !std::isinf( value ) && cvf::Math::valueInRange( value, minRange, maxRange ) ) + { + const double currentDiff = std::fabs( targetValue - value ); + if ( currentDiff < minDiff ) + { + minDiff = currentDiff; + startCell = reservoirCellIdx; + } + } + } + } + + return startCell; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RicCreateWellTargetClusterPolygonsFeature::assignClusterIdToCells( const RigActiveCellInfo& activeCellInfo, + const std::vector& cells, + std::vector& clusters, + int clusterId ) +{ + for ( size_t reservoirCellIdx : cells ) + { + size_t resultIndex = activeCellInfo.cellResultIndex( reservoirCellIdx ); + if ( resultIndex != cvf::UNDEFINED_SIZE_T ) clusters[resultIndex] = clusterId; + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RicCreateWellTargetClusterPolygonsFeature::growCluster( RimEclipseCase* eclipseCase, + size_t startCell, + const std::vector& volume, + double minRange, + double maxRange, + cvf::ref visibility, + std::vector& clusters, + int clusterId ) +{ + auto resultsData = eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL ); + + // Initially only the start cell is found + std::vector foundCells = { startCell }; + assignClusterIdToCells( *resultsData->activeCellInfo(), foundCells, clusters, clusterId ); + + while ( !foundCells.empty() ) + { + foundCells = findCandidates( *eclipseCase, foundCells, volume, minRange, maxRange, visibility, clusters ); + assignClusterIdToCells( *resultsData->activeCellInfo(), foundCells, clusters, clusterId ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RicCreateWellTargetClusterPolygonsFeature::findCandidates( const RimEclipseCase& eclipseCase, + const std::vector& previousCells, + const std::vector& volume, + double minRange, + double maxRange, + cvf::ref visibility, + std::vector& clusters ) +{ + std::vector candidates; + auto resultsData = eclipseCase.results( RiaDefines::PorosityModelType::MATRIX_MODEL ); + + const std::vector faces = { + cvf::StructGridInterface::FaceType::POS_I, + cvf::StructGridInterface::FaceType::NEG_I, + cvf::StructGridInterface::FaceType::POS_J, + cvf::StructGridInterface::FaceType::NEG_J, + cvf::StructGridInterface::FaceType::POS_K, + cvf::StructGridInterface::FaceType::NEG_K, + }; + + for ( size_t cellIdx : previousCells ) + { + const RigCell& cell = eclipseCase.mainGrid()->cell( cellIdx ); + if ( cell.isInvalid() ) continue; + + RigGridBase* grid = cell.hostGrid(); + size_t gridLocalCellIndex = cell.gridLocalCellIndex(); + size_t resultIndex = resultsData->activeCellInfo()->cellResultIndex( cellIdx ); + + size_t i, j, k; + + grid->ijkFromCellIndex( gridLocalCellIndex, &i, &j, &k ); + + for ( cvf::StructGridInterface::FaceType face : faces ) + { + size_t gridLocalNeighborCellIdx; + if ( grid->cellIJKNeighbor( i, j, k, face, &gridLocalNeighborCellIdx ) ) + { + size_t neighborResvCellIdx = grid->reservoirCellIndex( gridLocalNeighborCellIdx ); + size_t neighborResultIndex = resultsData->activeCellInfo()->cellResultIndex( neighborResvCellIdx ); + if ( neighborResultIndex != cvf::UNDEFINED_SIZE_T && clusters[neighborResultIndex] == 0 ) + { + double value = volume[neighborResultIndex]; + if ( !std::isinf( value ) && cvf::Math::valueInRange( value, minRange, maxRange ) && visibility->val( neighborResvCellIdx ) ) + { + candidates.push_back( neighborResvCellIdx ); + clusters[neighborResultIndex] = -1; + } + } + } + } + } + + return candidates; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RicCreateWellTargetClusterPolygonsFeature::createWellTargetClusterPolygons( RimEclipseCase* ensembleStatisticsCase, RimEclipseView* view ) +{ + auto overlay = view->overlayInfoConfig(); + RigHistogramData histogramData = overlay->histogramData(); + + auto computeBoundary = []( double histogramMin, double binSize, size_t index ) { return histogramMin + binSize * index; }; + + double binSize = ( histogramData.max - histogramData.min ) / histogramData.histogram.size(); + + // Smooth the histogram + size_t windowSize = 7; + std::vector smoothed = smoothHistogram( histogramData.histogram, windowSize ); + + // Find peaks with dynamic threshold + auto peaks = findPeaksWithDynamicProminence( smoothed ); + + RigCaseCellResultsData* resultsData = ensembleStatisticsCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL ); + const RigActiveCellInfo* activeCellInfo = resultsData->activeCellInfo(); + + std::vector clusterIds( activeCellInfo->reservoirActiveCellCount(), 0 ); + + int clusterId = 1; + for ( const auto& [index, prominence] : peaks ) + { + int minIndex = + cvf::Math::clamp( static_cast( index - ( windowSize / 2 ) ), 0, static_cast( histogramData.histogram.size() ) ); + int maxIndex = + cvf::Math::clamp( static_cast( index + ( windowSize / 2 ) ) + 1, 0, static_cast( histogramData.histogram.size() ) ); + + double minRange = computeBoundary( histogramData.min, binSize, minIndex ); + double maxRange = computeBoundary( histogramData.min, binSize, maxIndex ); + + auto resultAddress = view->cellResult()->eclipseResultAddress(); + auto cellCenters = findCellCentersWithResultInRange( *ensembleStatisticsCase, + resultAddress, + view->currentTimeStep(), + ( maxRange + minRange ) / 2.0, + minRange, + maxRange, + view->currentTotalCellVisibility(), + clusterIds, + clusterId ); + RiaLogging::info( QString( "Got %1 cell centers for range %2 %3" ).arg( cellCenters.size() ).arg( minRange ).arg( maxRange ) ); + clusterId++; + + // Need at least three points to make a polygon + if ( cellCenters.size() >= 3 ) + { + std::vector polygonBoundary = RigConvexHull::compute2d( cellCenters ); + polygonBoundary.push_back( polygonBoundary.front() ); + + auto proj = RimProject::current(); + auto polygonCollection = proj->activeOilField()->polygonCollection(); + + auto newPolygon = polygonCollection->appendUserDefinedPolygon(); + newPolygon->setPointsInDomainCoords( polygonBoundary ); + + polygonCollection->uiCapability()->updateAllRequiredEditors(); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RicCreateWellTargetClusterPolygonsFeature::setupActionLook( QAction* actionToSetup ) +{ + // actionToSetup->setIcon( QIcon( ":/GeoMechCasePropTable24x24.png" ) ); + actionToSetup->setText( "Create Well Target Cluster Polygons" ); +} diff --git a/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.h b/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.h new file mode 100644 index 0000000000..ddeb9500ee --- /dev/null +++ b/ApplicationLibCode/Commands/RicCreateWellTargetClusterPolygonsFeature.h @@ -0,0 +1,94 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2024- 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. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "cafCmdFeature.h" + +#include "cvfArray.h" +#include "cvfObject.h" + +class RimEclipseCase; +class RimEclipseView; +class RigEclipseResultAddress; +class RigActiveCellInfo; + +//================================================================================================== +/// +//================================================================================================== +class RicCreateWellTargetClusterPolygonsFeature : public caf::CmdFeature +{ + CAF_CMD_HEADER_INIT; + + static void createWellTargetClusterPolygons( RimEclipseCase* ensembleStatisticsCase, RimEclipseView* view ); + +protected: + static std::vector findCellCentersWithResultInRange( RimEclipseCase& eclipseCase, + RigEclipseResultAddress& resultAddress, + size_t timeStepIndex, + double rangeMinimum, + double rangeMaximum, + cvf::ref visibility ); + + static std::vector findCellCentersWithResultInRange( RimEclipseCase& eclipseCase, + const RigEclipseResultAddress& resultAddress, + size_t timeStepIndex, + double targetValue, + double rangeMinimum, + double rangeMaximum, + cvf::ref visibility, + std::vector& clusters, + int clusterId ); + + static size_t findStartCell( const RigActiveCellInfo& activeCellInfo, + const std::vector& values, + double targetValue, + double minRange, + double maxRange, + cvf::ref visibility, + const std::vector& clusters ); + + static std::vector findCandidates( const RimEclipseCase& eclipseCase, + const std::vector& previousCells, + const std::vector& volume, + double minRange, + double maxRange, + cvf::ref visibility, + std::vector& clusters ); + + static void assignClusterIdToCells( const RigActiveCellInfo& activeCellInfo, + const std::vector& cells, + std::vector& clusters, + int clusterId ); + + static void growCluster( RimEclipseCase* eclipseCase, + size_t startCell, + const std::vector& volume, + double minRange, + double maxRange, + cvf::ref visibility, + std::vector& clusters, + int clusterId ); + + static std::vector smoothHistogram( const std::vector& histogram, size_t windowSize ); + + static std::vector> findPeaksWithDynamicProminence( const std::vector& smoothedHistogram ); + + void onActionTriggered( bool isChecked ) override; + void setupActionLook( QAction* actionToSetup ) override; +}; diff --git a/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.cpp b/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.cpp index 10d8f6bbdc..bbc499b629 100644 --- a/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.cpp @@ -38,6 +38,7 @@ #include "RimRegularGridCase.h" #include "RimTools.h" +#include "cafCmdFeatureMenuBuilder.h" #include "cafPdmUiDoubleSliderEditor.h" #include "cafPdmUiPushButtonEditor.h" #include "cafPdmUiSliderTools.h" @@ -431,3 +432,11 @@ RimEclipseCase* RimWellTargetCandidatesGenerator::ensembleStatisticsCase() const { return m_ensembleStatisticsCase; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimWellTargetCandidatesGenerator::appendMenuItems( caf::CmdFeatureMenuBuilder& menuBuilder ) const +{ + menuBuilder << "RicCreateWellTargetClusterPolygonsFeature"; +} diff --git a/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.h b/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.h index c47fb01b5e..a38fc6bde5 100644 --- a/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.h +++ b/ApplicationLibCode/ProjectDataModel/RimWellTargetCandidatesGenerator.h @@ -50,6 +50,7 @@ class RimWellTargetCandidatesGenerator : public caf::PdmObject QList calculateValueOptions( const caf::PdmFieldHandle* fieldNeedingOptions ) override; void defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering ) override; void initAfterRead() override; + void appendMenuItems( caf::CmdFeatureMenuBuilder& menuBuilder ) const override; private: void generateCandidates( RimEclipseCase* eclipseCase );