ResInsight/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp
Magne Sjaastad e1e7201b9b #7742 Contour map : Fix mix of global grid and result index
Always use globalGridCellIndex for grid cells
Always use resultGridIndex for result property vectors
2021-06-04 08:30:15 +02:00

1798 lines
72 KiB
C++

/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2018- 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 "RimContourMapProjection.h"
#include "RiaWeightedGeometricMeanCalculator.h"
#include "RiaWeightedHarmonicMeanCalculator.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigCellGeometryTools.h"
#include "RigHexIntersectionTools.h"
#include "RimCase.h"
#include "RimGridView.h"
#include "RimProject.h"
#include "RimRegularLegendConfig.h"
#include "RimTextAnnotation.h"
#include "cafContourLines.h"
#include "cafPdmUiDoubleSliderEditor.h"
#include "cafPdmUiTreeOrdering.h"
#include "cafProgressInfo.h"
#include "cvfArray.h"
#include "cvfCellRange.h"
#include "cvfGeometryTools.h"
#include "cvfGeometryUtils.h"
#include "cvfScalarMapper.h"
#include "cvfStructGridGeometryGenerator.h"
#include <algorithm>
#ifdef USE_OPENMP
#include <omp.h>
#endif
namespace caf
{
template <>
void RimContourMapProjection::ResultAggregation::setUp()
{
addItem( RimContourMapProjection::RESULTS_OIL_COLUMN, "OIL_COLUMN", "Oil Column" );
addItem( RimContourMapProjection::RESULTS_GAS_COLUMN, "GAS_COLUMN", "Gas Column" );
addItem( RimContourMapProjection::RESULTS_HC_COLUMN, "HC_COLUMN", "Hydrocarbon Column" );
addItem( RimContourMapProjection::RESULTS_MEAN_VALUE, "MEAN_VALUE", "Arithmetic Mean" );
addItem( RimContourMapProjection::RESULTS_HARM_VALUE, "HARM_VALUE", "Harmonic Mean" );
addItem( RimContourMapProjection::RESULTS_GEOM_VALUE, "GEOM_VALUE", "Geometric Mean" );
addItem( RimContourMapProjection::RESULTS_VOLUME_SUM, "VOLUME_SUM", "Volume Weighted Sum" );
addItem( RimContourMapProjection::RESULTS_SUM, "SUM", "Sum" );
addItem( RimContourMapProjection::RESULTS_TOP_VALUE, "TOP_VALUE", "Top Value" );
addItem( RimContourMapProjection::RESULTS_MIN_VALUE, "MIN_VALUE", "Min Value" );
addItem( RimContourMapProjection::RESULTS_MAX_VALUE, "MAX_VALUE", "Max Value" );
setDefault( RimContourMapProjection::RESULTS_MEAN_VALUE );
}
} // namespace caf
CAF_PDM_ABSTRACT_SOURCE_INIT( RimContourMapProjection, "RimContourMapProjection" );
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::RimContourMapProjection()
: m_pickPoint( cvf::Vec2d::UNDEFINED )
, m_mapSize( cvf::Vec2ui( 0u, 0u ) )
, m_currentResultTimestep( -1 )
, m_minResultAllTimeSteps( std::numeric_limits<double>::infinity() )
, m_maxResultAllTimeSteps( -std::numeric_limits<double>::infinity() )
{
CAF_PDM_InitObject( "RimContourMapProjection", ":/2DMapProjection16x16.png", "", "" );
CAF_PDM_InitField( &m_relativeSampleSpacing, "SampleSpacing", 0.9, "Sample Spacing Factor", "", "", "" );
m_relativeSampleSpacing.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleSliderEditor::uiEditorTypeName() );
CAF_PDM_InitFieldNoDefault( &m_resultAggregation, "ResultAggregation", "Result Aggregation", "", "", "" );
CAF_PDM_InitField( &m_showContourLines, "ContourLines", true, "Show Contour Lines", "", "", "" );
CAF_PDM_InitField( &m_showContourLabels, "ContourLabels", true, "Show Contour Labels", "", "", "" );
CAF_PDM_InitField( &m_smoothContourLines, "SmoothContourLines", true, "Smooth Contour Lines", "", "", "" );
setName( "Map Projection" );
nameField()->uiCapability()->setUiReadOnly( true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::~RimContourMapProjection()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateResultsIfNecessary( int timeStep )
{
caf::ProgressInfo progress( 100, "Generate Results", true );
updateGridInformation();
progress.setProgress( 10 );
if ( gridMappingNeedsUpdating() || mapCellVisibilityNeedsUpdating() || resultVariableChanged() )
{
clearResults();
clearTimeStepRange();
if ( gridMappingNeedsUpdating() ) m_projected3dGridIndices = generateGridMapping();
progress.setProgress( 20 );
m_mapCellVisibility = getMapCellVisibility();
progress.setProgress( 30 );
}
else
{
progress.setProgress( 30 );
}
if ( resultsNeedsUpdating( timeStep ) )
{
clearGeometry();
m_aggregatedResults = generateResults( timeStep );
progress.setProgress( 80 );
generateVertexResults();
}
progress.setProgress( 100 );
updateAfterResultGeneration( timeStep );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateGeometryIfNecessary()
{
caf::ProgressInfo progress( 100, "Generate Geometry", true );
if ( geometryNeedsUpdating() )
{
generateContourPolygons();
progress.setProgress( 25 );
generateTrianglesWithVertexValues();
}
progress.setProgress( 100 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RimContourMapProjection::generatePickPointPolygon()
{
std::vector<cvf::Vec3d> points;
if ( !m_pickPoint.isUndefined() )
{
{
#ifndef NDEBUG
cvf::Vec2d cellDiagonal( sampleSpacing() * 0.5, sampleSpacing() * 0.5 );
cvf::Vec2ui pickedCell = ijFromLocalPos( m_pickPoint );
cvf::Vec2d cellCenter = cellCenterPosition( pickedCell.x(), pickedCell.y() );
cvf::Vec2d cellCorner = cellCenter - cellDiagonal;
points.push_back( cvf::Vec3d( cellCorner, 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( sampleSpacing(), 0.0 ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( sampleSpacing(), 0.0 ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( sampleSpacing(), sampleSpacing() ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( sampleSpacing(), sampleSpacing() ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( 0.0, sampleSpacing() ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( 0.0, sampleSpacing() ), 0.0 ) );
points.push_back( cvf::Vec3d( cellCorner, 0.0 ) );
#endif
points.push_back( cvf::Vec3d( m_pickPoint - cvf::Vec2d( 0.5 * sampleSpacing(), 0.0 ), 0.0 ) );
points.push_back( cvf::Vec3d( m_pickPoint + cvf::Vec2d( 0.5 * sampleSpacing(), 0.0 ), 0.0 ) );
points.push_back( cvf::Vec3d( m_pickPoint - cvf::Vec2d( 0.0, 0.5 * sampleSpacing() ), 0.0 ) );
points.push_back( cvf::Vec3d( m_pickPoint + cvf::Vec2d( 0.0, 0.5 * sampleSpacing() ), 0.0 ) );
}
}
return points;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::clearGeometry()
{
m_contourPolygons.clear();
m_trianglesWithVertexValues.clear();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<RimContourMapProjection::ContourPolygons>& RimContourMapProjection::contourPolygons() const
{
return m_contourPolygons;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<cvf::Vec4d>& RimContourMapProjection::trianglesWithVertexValues()
{
return m_trianglesWithVertexValues;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sampleSpacingFactor() const
{
return m_relativeSampleSpacing();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::setSampleSpacingFactor( double spacingFactor )
{
m_relativeSampleSpacing = spacingFactor;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::showContourLines() const
{
return m_showContourLines();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::showContourLabels() const
{
return m_showContourLabels();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::resultAggregationText() const
{
return m_resultAggregation().uiText();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::caseName() const
{
RimCase* rimCase = baseView()->ownerCase();
if ( !rimCase )
{
return QString();
}
return rimCase->caseUserDescription.value();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::currentTimeStepName() const
{
RimCase* rimCase = baseView()->ownerCase();
if ( !rimCase || m_currentResultTimestep == -1 )
{
return QString();
}
return rimCase->timeStepName( m_currentResultTimestep );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::maxValue() const
{
return maxValue( m_aggregatedResults );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::minValue() const
{
return minValue( m_aggregatedResults );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::meanValue() const
{
return sumAllValues() / numberOfValidCells();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sumAllValues() const
{
double sum = 0.0;
for ( size_t index = 0; index < m_aggregatedResults.size(); ++index )
{
if ( m_aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
sum += m_aggregatedResults[index];
}
}
return sum;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::numberOfElementsIJ() const
{
return m_mapSize;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::numberOfVerticesIJ() const
{
cvf::Vec2ui mapSize = this->numberOfElementsIJ();
mapSize.x() += 1u;
mapSize.y() += 1u;
return mapSize;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isColumnResult() const
{
return m_resultAggregation() == RESULTS_OIL_COLUMN || m_resultAggregation() == RESULTS_GAS_COLUMN ||
m_resultAggregation() == RESULTS_HC_COLUMN;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::valueAtVertex( uint i, uint j ) const
{
size_t index = vertexIndexFromIJ( i, j );
if ( index < numberOfVertices() )
{
return m_aggregatedVertexResults.at( index );
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
uint RimContourMapProjection::numberOfCells() const
{
return m_mapSize.x() * m_mapSize.y();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
uint RimContourMapProjection::numberOfValidCells() const
{
uint validCount = 0u;
for ( uint i = 0; i < numberOfCells(); ++i )
{
cvf::Vec2ui ij = ijFromCellIndex( i );
if ( hasResultInCell( ij.x(), ij.y() ) )
{
validCount++;
}
}
return validCount;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimContourMapProjection::numberOfVertices() const
{
cvf::Vec2ui gridSize = numberOfVerticesIJ();
return static_cast<size_t>( gridSize.x() ) * static_cast<size_t>( gridSize.y() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::checkForMapIntersection( const cvf::Vec3d& localPoint3d,
cvf::Vec2d* contourMapPoint,
double* valueAtPoint ) const
{
CVF_TIGHT_ASSERT( contourMapPoint );
CVF_TIGHT_ASSERT( valueAtPoint );
cvf::Vec3d mapPos3d = localPoint3d - m_expandedBoundingBox.min() + m_gridBoundingBox.min();
cvf::Vec2d mapPos2d( mapPos3d.x(), mapPos3d.y() );
cvf::Vec2d gridorigin( m_expandedBoundingBox.min().x(), m_expandedBoundingBox.min().y() );
double value = interpolateValue( mapPos2d );
if ( value != std::numeric_limits<double>::infinity() )
{
*valueAtPoint = value;
*contourMapPoint = mapPos2d + gridorigin;
return true;
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::setPickPoint( cvf::Vec2d globalPickPoint )
{
m_pickPoint = globalPickPoint - origin2d();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RimContourMapProjection::origin3d() const
{
return m_expandedBoundingBox.min();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimContourMapProjection::gridResultIndex( size_t globalCellIdx ) const
{
return globalCellIdx;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::calculateValueInMapCell( uint i, uint j, const std::vector<double>& gridCellValues ) const
{
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtIJ( i, j );
if ( !matchingCells.empty() )
{
switch ( m_resultAggregation() )
{
case RESULTS_TOP_VALUE:
{
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
return cellValue;
}
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_MEAN_VALUE:
{
RiaWeightedMeanCalculator<double> calculator;
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
calculator.addValueAndWeight( cellValue, weight );
}
}
if ( calculator.validAggregatedWeight() )
{
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_GEOM_VALUE:
{
RiaWeightedGeometricMeanCalculator calculator;
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
if ( cellValue < 1.0e-8 )
{
return 0.0;
}
calculator.addValueAndWeight( cellValue, weight );
}
}
if ( calculator.validAggregatedWeight() )
{
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_HARM_VALUE:
{
RiaWeightedHarmonicMeanCalculator calculator;
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::fabs( cellValue ) < 1.0e-8 )
{
return 0.0;
}
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
calculator.addValueAndWeight( cellValue, weight );
}
}
if ( calculator.validAggregatedWeight() )
{
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_MAX_VALUE:
{
double maxValue = -std::numeric_limits<double>::infinity();
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
maxValue = std::max( maxValue, cellValue );
}
}
if ( maxValue == -std::numeric_limits<double>::infinity() )
{
maxValue = std::numeric_limits<double>::infinity();
}
return maxValue;
}
case RESULTS_MIN_VALUE:
{
double minValue = std::numeric_limits<double>::infinity();
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
minValue = std::min( minValue, cellValue );
}
return minValue;
}
case RESULTS_VOLUME_SUM:
case RESULTS_SUM:
case RESULTS_OIL_COLUMN:
case RESULTS_GAS_COLUMN:
case RESULTS_HC_COLUMN:
{
double sum = 0.0;
for ( auto [cellIdx, weight] : matchingCells )
{
double cellValue = gridCellValues[gridResultIndex( cellIdx )];
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
{
sum += cellValue * weight;
}
}
return sum;
}
default:
CVF_TIGHT_ASSERT( false );
}
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::gridMappingNeedsUpdating() const
{
if ( m_projected3dGridIndices.size() != numberOfCells() )
{
return true;
}
if ( m_cellGridIdxVisibility.isNull() )
{
return true;
}
cvf::ref<cvf::UByteArray> currentVisibility = getCellVisibility();
CVF_ASSERT( currentVisibility->size() == m_cellGridIdxVisibility->size() );
for ( size_t i = 0; i < currentVisibility->size(); ++i )
{
if ( ( *currentVisibility )[i] != ( *m_cellGridIdxVisibility )[i] ) return true;
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::resultsNeedsUpdating( int timeStep ) const
{
if ( m_aggregatedResults.size() != numberOfCells() )
{
return true;
}
if ( m_aggregatedVertexResults.size() != numberOfVertices() )
{
return true;
}
if ( timeStep != m_currentResultTimestep )
{
return true;
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::geometryNeedsUpdating() const
{
return m_contourPolygons.empty() || m_trianglesWithVertexValues.empty();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::resultRangeIsValid() const
{
if ( m_minResultAllTimeSteps == std::numeric_limits<double>::infinity() ||
m_maxResultAllTimeSteps == -std::numeric_limits<double>::infinity() )
{
return false;
}
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::clearGridMapping()
{
clearResults();
clearTimeStepRange();
m_projected3dGridIndices.clear();
m_mapCellVisibility.clear();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::clearResults()
{
clearGeometry();
m_aggregatedResults.clear();
m_aggregatedVertexResults.clear();
m_currentResultTimestep = -1;
clearResultVariable();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::clearTimeStepRange()
{
m_minResultAllTimeSteps = std::numeric_limits<double>::infinity();
m_maxResultAllTimeSteps = -std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::maxValue( const std::vector<double>& aggregatedResults ) const
{
double maxV = -std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
maxV = std::max( maxV, aggregatedResults[index] );
}
}
return maxV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::minValue( const std::vector<double>& aggregatedResults ) const
{
double minV = std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
minV = std::min( minV, aggregatedResults[index] );
}
}
return minV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, double> RimContourMapProjection::minmaxValuesAllTimeSteps()
{
if ( !resultRangeIsValid() )
{
clearTimeStepRange();
m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( m_aggregatedResults ) );
m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( m_aggregatedResults ) );
for ( int i = 0; i < (int)baseView()->ownerCase()->timeStepStrings().size() - 1; ++i )
{
if ( i != m_currentResultTimestep )
{
std::vector<double> aggregatedResults = generateResults( i );
m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( aggregatedResults ) );
m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( aggregatedResults ) );
}
}
}
return std::make_pair( m_minResultAllTimeSteps, m_maxResultAllTimeSteps );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::UByteArray> RimContourMapProjection::getCellVisibility() const
{
return baseView()->currentTotalCellVisibility();
}
//--------------------------------------------------------------------------------------------------
/// Empty default implementation
//--------------------------------------------------------------------------------------------------
std::vector<bool> RimContourMapProjection::getMapCellVisibility()
{
return std::vector<bool>( numberOfCells(), true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::mapCellVisibilityNeedsUpdating()
{
std::vector<bool> mapCellVisiblity = getMapCellVisibility();
return !( mapCellVisiblity == m_mapCellVisibility );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<std::pair<size_t, double>>> RimContourMapProjection::generateGridMapping()
{
m_cellGridIdxVisibility = getCellVisibility();
int nCells = numberOfCells();
std::vector<std::vector<std::pair<size_t, double>>> projected3dGridIndices( nCells );
std::vector<double> weightingResultValues = retrieveParameterWeights();
if ( isStraightSummationResult() )
{
#pragma omp parallel for
for ( int index = 0; index < nCells; ++index )
{
cvf::Vec2ui ij = ijFromCellIndex( index );
cvf::Vec2d globalPos = cellCenterPosition( ij.x(), ij.y() ) + origin2d();
projected3dGridIndices[index] = cellRayIntersectionAndResults( globalPos, weightingResultValues );
}
}
else
{
#pragma omp parallel for
for ( int index = 0; index < nCells; ++index )
{
cvf::Vec2ui ij = ijFromCellIndex( index );
cvf::Vec2d globalPos = cellCenterPosition( ij.x(), ij.y() ) + origin2d();
projected3dGridIndices[index] = cellOverlapVolumesAndResults( globalPos, weightingResultValues );
}
}
return projected3dGridIndices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateVertexResults()
{
size_t nCells = numberOfCells();
if ( nCells != m_aggregatedResults.size() ) return;
size_t nVertices = numberOfVertices();
m_aggregatedVertexResults = std::vector<double>( nVertices, std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( nVertices ); ++index )
{
cvf::Vec2ui ij = ijFromVertexIndex( index );
m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() );
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateTrianglesWithVertexValues()
{
std::vector<cvf::Vec3d> vertices = generateVertices();
cvf::Vec2ui patchSize = numberOfVerticesIJ();
cvf::ref<cvf::UIntArray> faceList = new cvf::UIntArray;
cvf::GeometryUtils::tesselatePatchAsTriangles( patchSize.x(), patchSize.y(), 0u, true, faceList.p() );
bool discrete = false;
std::vector<double> contourLevels;
if ( legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER )
{
legendConfig()->scalarMapper()->majorTickValues( &contourLevels );
if ( legendConfig()->mappingMode() == RimRegularLegendConfig::MappingType::LINEAR_DISCRETE ||
legendConfig()->mappingMode() == RimRegularLegendConfig::MappingType::LOG10_DISCRETE )
{
discrete = true;
}
}
const double cellArea = sampleSpacing() * sampleSpacing();
const double areaThreshold = 1.0e-5 * 0.5 * cellArea;
std::vector<std::vector<std::vector<cvf::Vec3d>>> subtractPolygons;
if ( !m_contourPolygons.empty() )
{
subtractPolygons.resize( m_contourPolygons.size() );
for ( size_t i = 0; i < m_contourPolygons.size() - 1; ++i )
{
for ( size_t j = 0; j < m_contourPolygons[i + 1].size(); ++j )
{
subtractPolygons[i].push_back( m_contourPolygons[i + 1][j].vertices );
}
}
}
#ifdef USE_OPENMP
std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( omp_get_max_threads() );
#else
std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( 1 );
#endif
#pragma omp parallel
{
#ifdef USE_OPENMP
int myThread = omp_get_thread_num();
#else
int myThread = 0;
#endif
threadTriangles[myThread].resize( std::max( (size_t)1, m_contourPolygons.size() ) );
#pragma omp for schedule( dynamic )
for ( int64_t i = 0; i < (int64_t)faceList->size(); i += 3 )
{
std::vector<cvf::Vec3d> triangle( 3 );
std::vector<cvf::Vec4d> triangleWithValues( 3 );
bool anyValidVertex = false;
for ( size_t n = 0; n < 3; ++n )
{
uint vn = ( *faceList )[i + n];
double value = vn < m_aggregatedVertexResults.size() ? m_aggregatedVertexResults[vn]
: std::numeric_limits<double>::infinity();
triangle[n] = vertices[vn];
triangleWithValues[n] = cvf::Vec4d( vertices[vn], value );
if ( value != std::numeric_limits<double>::infinity() )
{
anyValidVertex = true;
}
}
if ( !anyValidVertex )
{
continue;
}
if ( m_contourPolygons.empty() )
{
threadTriangles[myThread][0].insert( threadTriangles[myThread][0].end(),
triangleWithValues.begin(),
triangleWithValues.end() );
continue;
}
bool outsideOuterLimit = false;
for ( size_t c = 0; c < m_contourPolygons.size() && !outsideOuterLimit; ++c )
{
std::vector<std::vector<cvf::Vec3d>> intersectPolygons;
for ( size_t j = 0; j < m_contourPolygons[c].size(); ++j )
{
bool containsAtLeastOne = false;
for ( size_t t = 0; t < 3; ++t )
{
if ( m_contourPolygons[c][j].bbox.contains( triangle[t] ) )
{
containsAtLeastOne = true;
}
}
if ( containsAtLeastOne )
{
std::vector<std::vector<cvf::Vec3d>> clippedPolygons =
RigCellGeometryTools::intersectionWithPolygon( triangle, m_contourPolygons[c][j].vertices );
intersectPolygons.insert( intersectPolygons.end(), clippedPolygons.begin(), clippedPolygons.end() );
}
}
if ( intersectPolygons.empty() )
{
outsideOuterLimit = true;
continue;
}
std::vector<std::vector<cvf::Vec3d>> clippedPolygons;
if ( !subtractPolygons[c].empty() )
{
for ( const std::vector<cvf::Vec3d>& polygon : intersectPolygons )
{
std::vector<std::vector<cvf::Vec3d>> fullyClippedPolygons =
RigCellGeometryTools::subtractPolygons( polygon, subtractPolygons[c] );
clippedPolygons.insert( clippedPolygons.end(),
fullyClippedPolygons.begin(),
fullyClippedPolygons.end() );
}
}
else
{
clippedPolygons.swap( intersectPolygons );
}
{
std::vector<cvf::Vec4d> clippedTriangles;
for ( std::vector<cvf::Vec3d>& clippedPolygon : clippedPolygons )
{
std::vector<std::vector<cvf::Vec3d>> polygonTriangles;
if ( clippedPolygon.size() == 3u )
{
polygonTriangles.push_back( clippedPolygon );
}
else
{
cvf::Vec3d baryCenter = cvf::Vec3d::ZERO;
for ( size_t v = 0; v < clippedPolygon.size(); ++v )
{
cvf::Vec3d& clippedVertex = clippedPolygon[v];
baryCenter += clippedVertex;
}
baryCenter /= clippedPolygon.size();
for ( size_t v = 0; v < clippedPolygon.size(); ++v )
{
std::vector<cvf::Vec3d> clippedTriangle;
if ( v == clippedPolygon.size() - 1 )
{
clippedTriangle = { clippedPolygon[v], clippedPolygon[0], baryCenter };
}
else
{
clippedTriangle = { clippedPolygon[v], clippedPolygon[v + 1], baryCenter };
}
polygonTriangles.push_back( clippedTriangle );
}
}
for ( const std::vector<cvf::Vec3d>& polygonTriangle : polygonTriangles )
{
// Check triangle area
double area = 0.5 * ( ( polygonTriangle[1] - polygonTriangle[0] ) ^
( polygonTriangle[2] - polygonTriangle[0] ) )
.length();
if ( area < areaThreshold ) continue;
for ( const cvf::Vec3d& localVertex : polygonTriangle )
{
double value = std::numeric_limits<double>::infinity();
if ( discrete )
{
value = contourLevels[c] + 0.01 * ( contourLevels.back() - contourLevels.front() ) /
contourLevels.size();
}
else
{
for ( size_t n = 0; n < 3; ++n )
{
if ( ( triangle[n] - localVertex ).length() < sampleSpacing() * 0.01 &&
triangleWithValues[n].w() != std::numeric_limits<double>::infinity() )
{
value = triangleWithValues[n].w();
break;
}
}
if ( value == std::numeric_limits<double>::infinity() )
{
value = interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) );
if ( value == std::numeric_limits<double>::infinity() )
{
value = contourLevels[c];
}
}
}
cvf::Vec4d globalVertex( localVertex, value );
clippedTriangles.push_back( globalVertex );
}
}
}
threadTriangles[myThread][c].insert( threadTriangles[myThread][c].end(),
clippedTriangles.begin(),
clippedTriangles.end() );
}
}
}
}
std::vector<std::vector<cvf::Vec4d>> trianglesPerLevel( std::max( (size_t)1, m_contourPolygons.size() ) );
for ( size_t c = 0; c < trianglesPerLevel.size(); ++c )
{
std::vector<cvf::Vec4d> allTrianglesThisLevel;
for ( size_t i = 0; i < threadTriangles.size(); ++i )
{
allTrianglesThisLevel.insert( allTrianglesThisLevel.end(),
threadTriangles[i][c].begin(),
threadTriangles[i][c].end() );
}
double triangleAreasThisLevel = sumTriangleAreas( allTrianglesThisLevel );
if ( c >= m_contourLevelCumulativeAreas.size() ||
triangleAreasThisLevel > 1.0e-3 * m_contourLevelCumulativeAreas[c] )
{
trianglesPerLevel[c] = allTrianglesThisLevel;
}
}
std::vector<cvf::Vec4d> finalTriangles;
for ( size_t i = 0; i < trianglesPerLevel.size(); ++i )
{
finalTriangles.insert( finalTriangles.end(), trianglesPerLevel[i].begin(), trianglesPerLevel[i].end() );
}
m_trianglesWithVertexValues = finalTriangles;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RimContourMapProjection::generateVertices() const
{
size_t nVertices = numberOfVertices();
std::vector<cvf::Vec3d> vertices( nVertices, cvf::Vec3d::ZERO );
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( nVertices ); ++index )
{
cvf::Vec2ui ij = ijFromVertexIndex( index );
cvf::Vec2d mapPos = cellCenterPosition( ij.x(), ij.y() );
// Shift away from sample point to vertex
mapPos.x() -= sampleSpacing() * 0.5;
mapPos.y() -= sampleSpacing() * 0.5;
cvf::Vec3d vertexPos( mapPos, 0.0 );
vertices[index] = vertexPos;
}
return vertices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateContourPolygons()
{
std::vector<ContourPolygons> contourPolygons;
std::vector<double> contourLevels;
if ( resultRangeIsValid() && legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER )
{
legendConfig()->scalarMapper()->majorTickValues( &contourLevels );
int nContourLevels = static_cast<int>( contourLevels.size() );
if ( minValue() != std::numeric_limits<double>::infinity() &&
maxValue() != -std::numeric_limits<double>::infinity() && std::fabs( maxValue() - minValue() ) > 1.0e-8 )
{
if ( nContourLevels > 2 )
{
const size_t N = contourLevels.size();
// Adjust contour levels slightly to avoid weird visual artifacts due to numerical error.
double fudgeFactor = 1.0e-3;
double fudgeAmountMin = fudgeFactor * ( contourLevels[1] - contourLevels[0] );
double fudgeAmountMax = fudgeFactor * ( contourLevels[N - 1u] - contourLevels[N - 2u] );
contourLevels.front() += fudgeAmountMin;
contourLevels.back() -= fudgeAmountMax;
double simplifyEpsilon = m_smoothContourLines() ? 5.0e-2 * sampleSpacing() : 1.0e-3 * sampleSpacing();
if ( nContourLevels >= 10 )
{
simplifyEpsilon *= 2.0;
}
if ( numberOfCells() > 100000 )
{
simplifyEpsilon *= 2.0;
}
else if ( numberOfCells() > 1000000 )
{
simplifyEpsilon *= 4.0;
}
std::vector<caf::ContourLines::ListOfLineSegments> unorderedLineSegmentsPerLevel =
caf::ContourLines::create( m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels );
contourPolygons = std::vector<ContourPolygons>( unorderedLineSegmentsPerLevel.size() );
#pragma omp parallel for
for ( int i = 0; i < (int)unorderedLineSegmentsPerLevel.size(); ++i )
{
contourPolygons[i] =
createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i], contourLevels[i] );
if ( m_smoothContourLines() )
{
smoothContourPolygons( &contourPolygons[i], true );
}
for ( ContourPolygon& polygon : contourPolygons[i] )
{
RigCellGeometryTools::simplifyPolygon( &polygon.vertices, simplifyEpsilon );
}
}
// The clipping of contour polygons is intended to detect and fix a smoothed contour polygons
// crossing into an outer contour line. The current implementation has some side effects causing
// several contour lines to disappear. Disable this clipping for now
/*
if ( m_smoothContourLines() )
{
for ( size_t i = 1; i < contourPolygons.size(); ++i )
{
clipContourPolygons(&contourPolygons[i], &contourPolygons[i - 1] );
}
}
*/
m_contourLevelCumulativeAreas.resize( contourPolygons.size(), 0.0 );
for ( int64_t i = (int64_t)contourPolygons.size() - 1; i >= 0; --i )
{
double levelOuterArea = sumPolygonArea( contourPolygons[i] );
m_contourLevelCumulativeAreas[i] = levelOuterArea;
}
}
}
}
m_contourPolygons = contourPolygons;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::ContourPolygons
RimContourMapProjection::createContourPolygonsFromLineSegments( caf::ContourLines::ListOfLineSegments& unorderedLineSegments,
double contourValue )
{
const double areaThreshold =
1.5 * ( sampleSpacing() * sampleSpacing() ) / ( sampleSpacingFactor() * sampleSpacingFactor() );
ContourPolygons contourPolygons;
std::vector<std::vector<cvf::Vec3d>> polygons;
RigCellGeometryTools::createPolygonFromLineSegments( unorderedLineSegments, polygons, 1.0e-8 );
for ( size_t j = 0; j < polygons.size(); ++j )
{
double signedArea = cvf::GeometryTools::signedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, polygons[j] );
ContourPolygon contourPolygon;
contourPolygon.value = contourValue;
if ( signedArea < 0.0 )
{
contourPolygon.vertices.insert( contourPolygon.vertices.end(), polygons[j].rbegin(), polygons[j].rend() );
}
else
{
contourPolygon.vertices = polygons[j];
}
contourPolygon.area = cvf::GeometryTools::signedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, contourPolygon.vertices );
if ( contourPolygon.area > areaThreshold )
{
for ( const cvf::Vec3d& vertex : contourPolygon.vertices )
{
contourPolygon.bbox.add( vertex );
}
contourPolygons.push_back( contourPolygon );
}
}
return contourPolygons;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::smoothContourPolygons( ContourPolygons* contourPolygons, bool favourExpansion )
{
CVF_ASSERT( contourPolygons );
for ( size_t i = 0; i < contourPolygons->size(); ++i )
{
ContourPolygon& polygon = contourPolygons->at( i );
for ( size_t n = 0; n < 20; ++n )
{
std::vector<cvf::Vec3d> newVertices;
newVertices.resize( polygon.vertices.size() );
double maxChange = 0.0;
for ( size_t j = 0; j < polygon.vertices.size(); ++j )
{
cvf::Vec3d vm1 = polygon.vertices.back();
cvf::Vec3d v = polygon.vertices[j];
cvf::Vec3d vp1 = polygon.vertices.front();
if ( j > 0u )
{
vm1 = polygon.vertices[j - 1];
}
if ( j < polygon.vertices.size() - 1 )
{
vp1 = polygon.vertices[j + 1];
}
// Only expand.
cvf::Vec3d modifiedVertex = 0.5 * ( v + 0.5 * ( vm1 + vp1 ) );
cvf::Vec3d delta = modifiedVertex - v;
cvf::Vec3d tangent3d = vp1 - vm1;
cvf::Vec2d tangent2d( tangent3d.x(), tangent3d.y() );
cvf::Vec3d norm3d( tangent2d.getNormalized().perpendicularVector() );
if ( delta * norm3d > 0 && favourExpansion )
{
// Normal is always inwards facing so a positive dot product means inward movement
// Favour expansion rather than contraction by only contracting by a fraction.
// The fraction is empirically found to give a decent result.
modifiedVertex = v + 0.2 * delta;
}
newVertices[j] = modifiedVertex;
maxChange = std::max( maxChange, ( modifiedVertex - v ).length() );
}
polygon.vertices.swap( newVertices );
if ( maxChange < sampleSpacing() * 1.0e-2 ) break;
}
}
}
void RimContourMapProjection::clipContourPolygons( ContourPolygons* contourPolygons, const ContourPolygons* clipBy )
{
CVF_ASSERT( clipBy );
for ( size_t i = 0; i < contourPolygons->size(); ++i )
{
ContourPolygon& polygon = contourPolygons->at( i );
for ( size_t j = 0; j < clipBy->size(); ++j )
{
std::vector<std::vector<cvf::Vec3d>> intersections =
RigCellGeometryTools::intersectionWithPolygon( polygon.vertices, clipBy->at( j ).vertices );
if ( !intersections.empty() )
{
polygon.vertices = intersections.front();
polygon.area =
std::abs( cvf::GeometryTools::signedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, polygon.vertices ) );
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sumPolygonArea( const ContourPolygons& contourPolygons )
{
double sumArea = 0.0;
for ( const ContourPolygon& polygon : contourPolygons )
{
sumArea += polygon.area;
}
return sumArea;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sumTriangleAreas( const std::vector<cvf::Vec4d>& triangles )
{
double sumArea = 0.0;
for ( size_t i = 0; i < triangles.size(); i += 3 )
{
cvf::Vec3d v1( triangles[i].x(), triangles[i].y(), triangles[i].z() );
cvf::Vec3d v2( triangles[i + 1].x(), triangles[i + 1].y(), triangles[i + 1].z() );
cvf::Vec3d v3( triangles[i + 2].x(), triangles[i + 2].y(), triangles[i + 2].z() );
double area = 0.5 * ( ( v3 - v1 ) ^ ( v2 - v1 ) ).length();
sumArea += area;
}
return sumArea;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RimContourMapProjection::CellIndexAndResult>
RimContourMapProjection::cellOverlapVolumesAndResults( const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues ) const
{
cvf::Vec3d top2dElementCentroid( globalPos2d, m_expandedBoundingBox.max().z() );
cvf::Vec3d bottom2dElementCentroid( globalPos2d, m_expandedBoundingBox.min().z() );
cvf::Vec3d planarDiagonalVector( 0.5 * sampleSpacing(), 0.5 * sampleSpacing(), 0.0 );
cvf::Vec3d topNECorner = top2dElementCentroid + planarDiagonalVector;
cvf::Vec3d bottomSWCorner = bottom2dElementCentroid - planarDiagonalVector;
cvf::BoundingBox bbox2dElement( bottomSWCorner, topNECorner );
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
// Bounding box has been expanded, so 2d element may be outside actual 3d grid
if ( !bbox2dElement.intersects( m_gridBoundingBox ) )
{
return matchingVisibleCellsAndWeight;
}
std::vector<size_t> allCellIndices = findIntersectingCells( bbox2dElement );
std::vector<std::vector<size_t>> kLayerCellIndexVector;
kLayerCellIndexVector.resize( kLayers() );
for ( size_t globalCellIdx : allCellIndices )
{
if ( ( *m_cellGridIdxVisibility )[globalCellIdx] )
{
kLayerCellIndexVector[kLayer( globalCellIdx )].push_back( globalCellIdx );
}
}
for ( const auto& kLayerIndices : kLayerCellIndexVector )
{
for ( size_t globalCellIdx : kLayerIndices )
{
double overlapVolume = calculateOverlapVolume( globalCellIdx, bbox2dElement );
if ( overlapVolume > 0.0 )
{
double weight = overlapVolume *
getParameterWeightForCell( gridResultIndex( globalCellIdx ), weightingResultValues );
if ( weight > 0.0 )
{
matchingVisibleCellsAndWeight.push_back( std::make_pair( globalCellIdx, weight ) );
}
}
}
}
return matchingVisibleCellsAndWeight;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RimContourMapProjection::CellIndexAndResult>
RimContourMapProjection::cellRayIntersectionAndResults( const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues ) const
{
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
cvf::Vec3d highestPoint( globalPos2d, m_expandedBoundingBox.max().z() );
cvf::Vec3d lowestPoint( globalPos2d, m_expandedBoundingBox.min().z() );
// Bounding box has been expanded, so ray may be outside actual 3d grid
if ( !m_gridBoundingBox.contains( highestPoint ) )
{
return matchingVisibleCellsAndWeight;
}
cvf::BoundingBox rayBBox;
rayBBox.add( highestPoint );
rayBBox.add( lowestPoint );
std::vector<size_t> allCellIndices = findIntersectingCells( rayBBox );
std::map<size_t, std::vector<size_t>> kLayerIndexMap;
for ( size_t globalCellIdx : allCellIndices )
{
if ( ( *m_cellGridIdxVisibility )[globalCellIdx] )
{
kLayerIndexMap[kLayer( globalCellIdx )].push_back( globalCellIdx );
}
}
for ( const auto& kLayerIndexPair : kLayerIndexMap )
{
double weightSumThisKLayer = 0.0;
std::vector<std::pair<size_t, double>> cellsAndWeightsThisLayer;
for ( size_t globalCellIdx : kLayerIndexPair.second )
{
double lengthInCell = calculateRayLengthInCell( globalCellIdx, highestPoint, lowestPoint );
if ( lengthInCell > 0.0 )
{
cellsAndWeightsThisLayer.push_back( std::make_pair( globalCellIdx, lengthInCell ) );
weightSumThisKLayer += lengthInCell;
}
}
for ( auto& cellWeightPair : cellsAndWeightsThisLayer )
{
cellWeightPair.second /= weightSumThisKLayer;
matchingVisibleCellsAndWeight.push_back( cellWeightPair );
}
}
return matchingVisibleCellsAndWeight;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isMeanResult() const
{
return m_resultAggregation() == RESULTS_MEAN_VALUE || m_resultAggregation() == RESULTS_HARM_VALUE ||
m_resultAggregation() == RESULTS_GEOM_VALUE;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isStraightSummationResult() const
{
return isStraightSummationResult( m_resultAggregation() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isStraightSummationResult( ResultAggregationEnum aggregationType )
{
return aggregationType == RESULTS_OIL_COLUMN || aggregationType == RESULTS_GAS_COLUMN ||
aggregationType == RESULTS_HC_COLUMN || aggregationType == RESULTS_SUM;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) const
{
cvf::Vec2ui cellContainingPoint = ijFromLocalPos( gridPos2d );
cvf::Vec2d cellCenter = cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() );
std::array<cvf::Vec3d, 4> x;
x[0] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing() * 0.5, -sampleSpacing() * 0.5 ), 0.0 );
x[1] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing() * 0.5, -sampleSpacing() * 0.5 ), 0.0 );
x[2] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing() * 0.5, sampleSpacing() * 0.5 ), 0.0 );
x[3] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing() * 0.5, sampleSpacing() * 0.5 ), 0.0 );
cvf::Vec4d baryCentricCoords =
cvf::GeometryTools::barycentricCoords( x[0], x[1], x[2], x[3], cvf::Vec3d( gridPos2d, 0.0 ) );
std::array<cvf::Vec2ui, 4> v;
v[0] = cellContainingPoint;
v[1] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() );
v[2] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() + 1u );
v[3] = cvf::Vec2ui( cellContainingPoint.x(), cellContainingPoint.y() + 1u );
std::array<double, 4> vertexValues;
double validBarycentricCoordsSum = 0.0;
for ( int i = 0; i < 4; ++i )
{
double vertexValue = valueAtVertex( v[i].x(), v[i].y() );
if ( vertexValue == std::numeric_limits<double>::infinity() )
{
return std::numeric_limits<double>::infinity();
}
else
{
vertexValues[i] = vertexValue;
validBarycentricCoordsSum += baryCentricCoords[i];
}
}
if ( validBarycentricCoordsSum < 1.0e-8 )
{
return std::numeric_limits<double>::infinity();
}
// Calculate final value
double value = 0.0;
for ( int i = 0; i < 4; ++i )
{
value += baryCentricCoords[i] / validBarycentricCoordsSum * vertexValues[i];
}
return value;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::valueInCell( uint i, uint j ) const
{
size_t index = cellIndexFromIJ( i, j );
if ( index < numberOfCells() )
{
return m_aggregatedResults.at( index );
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::hasResultInCell( uint i, uint j ) const
{
return !cellsAtIJ( i, j ).empty();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::calculateValueAtVertex( uint vi, uint vj ) const
{
std::vector<uint> averageIs;
std::vector<uint> averageJs;
if ( vi > 0u ) averageIs.push_back( vi - 1 );
if ( vj > 0u ) averageJs.push_back( vj - 1 );
if ( vi < m_mapSize.x() ) averageIs.push_back( vi );
if ( vj < m_mapSize.y() ) averageJs.push_back( vj );
RiaWeightedMeanCalculator<double> calc;
for ( uint j : averageJs )
{
for ( uint i : averageIs )
{
if ( hasResultInCell( i, j ) )
{
calc.addValueAndWeight( valueInCell( i, j ), 1.0 );
}
}
}
if ( calc.validAggregatedWeight() )
{
return calc.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>> RimContourMapProjection::cellsAtIJ( uint i, uint j ) const
{
size_t cellIndex = this->cellIndexFromIJ( i, j );
if ( cellIndex < m_projected3dGridIndices.size() )
{
return m_projected3dGridIndices[cellIndex];
}
return std::vector<std::pair<size_t, double>>();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimContourMapProjection::cellIndexFromIJ( uint i, uint j ) const
{
CVF_ASSERT( i < m_mapSize.x() );
CVF_ASSERT( j < m_mapSize.y() );
return i + j * m_mapSize.x();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimContourMapProjection::vertexIndexFromIJ( uint i, uint j ) const
{
return i + j * ( m_mapSize.x() + 1 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::ijFromVertexIndex( size_t gridIndex ) const
{
cvf::Vec2ui gridSize = numberOfVerticesIJ();
uint quotientX = static_cast<uint>( gridIndex ) / gridSize.x();
uint remainderX = static_cast<uint>( gridIndex ) % gridSize.x();
return cvf::Vec2ui( remainderX, quotientX );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::ijFromCellIndex( size_t cellIndex ) const
{
CVF_TIGHT_ASSERT( cellIndex < numberOfCells() );
uint quotientX = static_cast<uint>( cellIndex ) / m_mapSize.x();
uint remainderX = static_cast<uint>( cellIndex ) % m_mapSize.x();
return cvf::Vec2ui( remainderX, quotientX );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::ijFromLocalPos( const cvf::Vec2d& localPos2d ) const
{
uint i = localPos2d.x() / sampleSpacing();
uint j = localPos2d.y() / sampleSpacing();
return cvf::Vec2ui( i, j );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2d RimContourMapProjection::cellCenterPosition( uint i, uint j ) const
{
cvf::Vec3d gridExtent = m_expandedBoundingBox.extent();
cvf::Vec2d cellCorner =
cvf::Vec2d( ( i * gridExtent.x() ) / ( m_mapSize.x() ), ( j * gridExtent.y() ) / ( m_mapSize.y() ) );
return cellCorner + cvf::Vec2d( sampleSpacing() * 0.5, sampleSpacing() * 0.5 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2d RimContourMapProjection::origin2d() const
{
return cvf::Vec2d( m_expandedBoundingBox.min().x(), m_expandedBoundingBox.min().y() );
}
//--------------------------------------------------------------------------------------------------
/// Vertex positions in local coordinates (add origin2d.x() for UTM x)
//--------------------------------------------------------------------------------------------------
std::vector<double> RimContourMapProjection::xVertexPositions() const
{
double gridExtent = m_expandedBoundingBox.extent().x();
cvf::Vec2ui gridSize = numberOfVerticesIJ();
std::vector<double> positions;
positions.reserve( gridSize.x() );
for ( uint i = 0; i < gridSize.x(); ++i )
{
positions.push_back( ( i * gridExtent ) / ( gridSize.x() - 1 ) );
}
return positions;
}
//--------------------------------------------------------------------------------------------------
/// Vertex positions in local coordinates (add origin2d.y() for UTM y)
//--------------------------------------------------------------------------------------------------
std::vector<double> RimContourMapProjection::yVertexPositions() const
{
double gridExtent = m_expandedBoundingBox.extent().y();
cvf::Vec2ui gridSize = numberOfVerticesIJ();
std::vector<double> positions;
positions.reserve( gridSize.y() );
for ( uint j = 0; j < gridSize.y(); ++j )
{
positions.push_back( ( j * gridExtent ) / ( gridSize.y() - 1 ) );
}
return positions;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::calculateMapSize() const
{
cvf::Vec3d gridExtent = m_expandedBoundingBox.extent();
uint projectionSizeX = static_cast<uint>( std::ceil( gridExtent.x() / sampleSpacing() ) );
uint projectionSizeY = static_cast<uint>( std::ceil( gridExtent.y() / sampleSpacing() ) );
return cvf::Vec2ui( projectionSizeX, projectionSizeY );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::gridEdgeOffset() const
{
return sampleSpacing() * 2.0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::fieldChangedByUi( const caf::PdmFieldHandle* changedField,
const QVariant& oldValue,
const QVariant& newValue )
{
if ( changedField == &m_resultAggregation )
{
ResultAggregation previousAggregation = static_cast<ResultAggregationEnum>( oldValue.toInt() );
if ( isStraightSummationResult( previousAggregation ) != isStraightSummationResult() )
{
clearGridMapping();
}
else
{
clearResults();
}
clearTimeStepRange();
}
else if ( changedField == &m_smoothContourLines )
{
clearGeometry();
}
else if ( changedField == &m_relativeSampleSpacing )
{
clearGridMapping();
clearResults();
clearTimeStepRange();
}
baseView()->updateConnectedEditors();
RimProject* proj;
firstAncestorOrThisOfTypeAsserted( proj );
proj->scheduleCreateDisplayModelAndRedrawAllViews();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineEditorAttribute( const caf::PdmFieldHandle* field,
QString uiConfigName,
caf::PdmUiEditorAttribute* attribute )
{
if ( &m_relativeSampleSpacing == field )
{
caf::PdmUiDoubleSliderEditorAttribute* myAttr = dynamic_cast<caf::PdmUiDoubleSliderEditorAttribute*>( attribute );
if ( myAttr )
{
myAttr->m_minimum = 0.2;
myAttr->m_maximum = 2.0;
myAttr->m_sliderTickCount = 9;
myAttr->m_delaySliderUpdateUntilRelease = true;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering )
{
caf::PdmUiGroup* mainGroup = uiOrdering.addNewGroup( "Projection Settings" );
mainGroup->add( &m_resultAggregation );
legendConfig()->uiOrdering( "NumLevelsOnly", *mainGroup );
mainGroup->add( &m_relativeSampleSpacing );
mainGroup->add( &m_showContourLines );
mainGroup->add( &m_showContourLabels );
m_showContourLabels.uiCapability()->setUiReadOnly( !m_showContourLines() );
mainGroup->add( &m_smoothContourLines );
m_smoothContourLines.uiCapability()->setUiReadOnly( !m_showContourLines() );
uiOrdering.skipRemainingFields( true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineUiTreeOrdering( caf::PdmUiTreeOrdering& uiTreeOrdering, QString uiConfigName /*= ""*/ )
{
uiTreeOrdering.skipRemainingChildren( true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::initAfterRead()
{
}