OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactor: Extractor contour polygon tools.

Browse Source
This commit is contained in:
Kristian Bendiksen
2024-10-17 09:19:52 +02:00
parent c636a80fe2
commit ebc805ab65
7 changed files with 268 additions and 182 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
ApplicationLibCode/ModelVisualization/RivContourMapProjectionPartMgr.cpp
View File

@@ -4,6 +4,7 @@
#include "RiaFontCache.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigCellGeometryTools.h"
#include "RigContourPolygonsTools.h"
#include "RivMeshLinesSourceInfo.h"
#include "RivPartPriority.h"
#include "RivScalarMapperUtils.h"
@@ -409,7 +410,7 @@ std::vector<cvf::ref<cvf::Drawable>>
std::vector<double> tickValues;
mapper->majorTickValues( &tickValues );
const RimContourMapProjection::ContourPolygons* previousLevel = nullptr;
const RigContourPolygonsTools::ContourPolygons* previousLevel = nullptr;
for ( int64_t i = (int64_t)m_contourLinePolygons.size() - 1; i > 0; --i )
{
cvf::Color3f backgroundColor( mapper->mapToColor( tickValues[i] ) );
@@ -435,7 +436,9 @@ std::vector<cvf::ref<cvf::Drawable>>
const cvf::Vec3d& localVertex2 = m_contourLinePolygons[i][j].vertices[nextVertex];
cvf::Vec3d lineCenter = ( localVertex1 + localVertex2 ) * 0.5;
if ( previousLevel && lineOverlapsWithPreviousContourLevel( lineCenter, previousLevel ) )
double tolerance = 1.0e-2 * m_contourMapProjection->sampleSpacing();
if ( previousLevel && lineOverlapsWithPreviousContourLevel( lineCenter, *previousLevel, tolerance ) )
{
continue;
}
@@ -557,41 +560,9 @@ cvf::ref<cvf::DrawableGeo>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RivContourMapProjectionPartMgr::lineOverlapsWithPreviousContourLevel( const cvf::Vec3d& lineCenter,
const RimContourMapProjection::ContourPolygons* previousLevel ) const
bool RivContourMapProjectionPartMgr::lineOverlapsWithPreviousContourLevel( const cvf::Vec3d& lineCenter,
const RigContourPolygonsTools::ContourPolygons& previousLevel,
double tolerance )
{
const int64_t jump = 50;
CVF_ASSERT( previousLevel );
double tolerance = 1.0e-2 * m_contourMapProjection->sampleSpacing();
for ( const RimContourMapProjection::ContourPolygon& edgePolygon : *previousLevel )
{
std::pair<int64_t, double> closestIndex( 0, std::numeric_limits<double>::infinity() );
for ( int64_t i = 0; i < (int64_t)edgePolygon.vertices.size(); i += jump )
{
const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[( i + 1 ) % edgePolygon.vertices.size()];
double dist1 = cvf::GeometryTools::linePointSquareDist( edgeVertex1, edgeVertex2, lineCenter );
if ( dist1 < tolerance )
{
return true;
}
if ( dist1 < closestIndex.second )
{
closestIndex = std::make_pair( i, dist1 );
}
}
for ( int64_t i = std::max( (int64_t)1, closestIndex.first - jump + 1 );
i < std::min( (int64_t)edgePolygon.vertices.size(), closestIndex.first + jump );
++i )
{
const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[( i + 1 ) % edgePolygon.vertices.size()];
double dist1 = cvf::GeometryTools::linePointSquareDist( edgeVertex1, edgeVertex2, lineCenter );
if ( dist1 < tolerance )
{
return true;
}
}
}
return false;
return RigContourPolygonsTools::lineOverlapsWithContourPolygons( lineCenter, previousLevel, tolerance );
}

7
ApplicationLibCode/ModelVisualization/RivContourMapProjectionPartMgr.h
View File

@@ -18,7 +18,7 @@
#pragma once
#include "RimEclipseContourMapProjection.h"
#include "RimContourMapProjection.h"
#include "cafDisplayCoordTransform.h"
#include "cafPdmPointer.h"
@@ -59,8 +59,9 @@ private:
const caf::DisplayCoordTransform* displayCoordTransform,
std::vector<std::vector<cvf::BoundingBox>>* labelBBoxes ) const;
cvf::ref<cvf::DrawableGeo> createPickPointVisDrawable( const caf::DisplayCoordTransform* displayCoordTransform ) const;
bool lineOverlapsWithPreviousContourLevel( const cvf::Vec3d& lineCenter,
const RimContourMapProjection::ContourPolygons* previousLevel ) const;
static bool lineOverlapsWithPreviousContourLevel( const cvf::Vec3d& lineCenter,
const RimContourMapProjection::ContourPolygons& previousLevel,
double tolerance );
private:
caf::PdmPointer<RimContourMapProjection> m_contourMapProjection;

135
ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp
View File

@@ -1061,18 +1061,21 @@ void RimContourMapProjection::generateContourPolygons()
caf::ContourLines::create( m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels );
contourPolygons = std::vector<ContourPolygons>( unorderedLineSegmentsPerLevel.size() );
const double areaThreshold = 1.5 * ( sampleSpacing() * sampleSpacing() ) / ( sampleSpacingFactor() * sampleSpacingFactor() );
#pragma omp parallel for
for ( int i = 0; i < (int)unorderedLineSegmentsPerLevel.size(); ++i )
{
contourPolygons[i] = createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i], contourLevels[i] );
contourPolygons[i] = RigContourPolygonsTools::createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i],
contourLevels[i],
areaThreshold );
if ( m_smoothContourLines() )
{
smoothContourPolygons( &contourPolygons[i], true );
RigContourPolygonsTools::smoothContourPolygons( contourPolygons[i], true, sampleSpacing() );
}
for ( ContourPolygon& polygon : contourPolygons[i] )
for ( RigContourPolygonsTools::ContourPolygon& polygon : contourPolygons[i] )
{
RigCellGeometryTools::simplifyPolygon( &polygon.vertices, simplifyEpsilon );
}
@@ -1086,7 +1089,7 @@ void RimContourMapProjection::generateContourPolygons()
{
for ( size_t i = 1; i < contourPolygons.size(); ++i )
{
clipContourPolygons(&contourPolygons[i], &contourPolygons[i - 1] );
RigContourPolygonsTools::clipContourPolygons(&contourPolygons[i], &contourPolygons[i - 1] );
}
}
*/
@@ -1094,7 +1097,7 @@ void RimContourMapProjection::generateContourPolygons()
m_contourLevelCumulativeAreas.resize( contourPolygons.size(), 0.0 );
for ( int64_t i = (int64_t)contourPolygons.size() - 1; i >= 0; --i )
{
double levelOuterArea = sumPolygonArea( contourPolygons[i] );
double levelOuterArea = RigContourPolygonsTools::sumPolygonArea( contourPolygons[i] );
m_contourLevelCumulativeAreas[i] = levelOuterArea;
}
}
@@ -1103,128 +1106,6 @@ void RimContourMapProjection::generateContourPolygons()
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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

18
ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h
View File

@@ -20,7 +20,8 @@
#include "RimCheckableNamedObject.h"
#include "cafContourLines.h"
#include "RigContourPolygonsTools.h"
#include "cafPdmField.h"
#include "cafPdmObject.h"
@@ -43,14 +44,6 @@ class RimContourMapProjection : public RimCheckableNamedObject
public:
using CellIndexAndResult = std::pair<size_t, double>;
struct ContourPolygon
{
std::vector<cvf::Vec3d> vertices;
double value;
double area;
cvf::BoundingBox bbox;
};
enum ResultAggregationEnum
{
RESULTS_TOP_VALUE,
@@ -66,7 +59,7 @@ public:
RESULTS_HC_COLUMN
};
using ResultAggregation = caf::AppEnum<ResultAggregationEnum>;
using ContourPolygons = std::vector<ContourPolygon>;
using ContourPolygons = std::vector<RigContourPolygonsTools::ContourPolygon>;
RimContourMapProjection();
~RimContourMapProjection() override;
@@ -164,10 +157,7 @@ protected:
void generateVertexResults();
void generateTrianglesWithVertexValues();
void generateContourPolygons();
ContourPolygons createContourPolygonsFromLineSegments( caf::ContourLines::ListOfLineSegments& unorderedLineSegments, double contourValue );
void smoothContourPolygons( ContourPolygons* contourPolygons, bool favourExpansion );
void clipContourPolygons( ContourPolygons* contourPolygons, const ContourPolygons* clipBy );
static double sumPolygonArea( const ContourPolygons& contourPolygons );
static double sumTriangleAreas( const std::vector<cvf::Vec4d>& triangles );
std::vector<CellIndexAndResult> cellOverlapVolumesAndResults( const cvf::Vec2d& globalPos2d,

2
ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake
View File

@@ -103,6 +103,7 @@ set(SOURCE_GROUP_HEADER_FILES
${CMAKE_CURRENT_LIST_DIR}/RigOsduWellLogData.h
${CMAKE_CURRENT_LIST_DIR}/RigWellTargetCandidatesGenerator.h
${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid.h
${CMAKE_CURRENT_LIST_DIR}/RigContourPolygonsTools.h
)
set(SOURCE_GROUP_SOURCE_FILES
@@ -204,6 +205,7 @@ set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RigOsduWellLogData.cpp
${CMAKE_CURRENT_LIST_DIR}/RigWellTargetCandidatesGenerator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid.cpp
${CMAKE_CURRENT_LIST_DIR}/RigContourPolygonsTools.cpp
)
list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES})

188
ApplicationLibCode/ReservoirDataModel/RigContourPolygonsTools.cpp Normal file
View File

@@ -0,0 +1,188 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RigContourPolygonsTools.h"
#include "RigCellGeometryTools.h"
#include "cafContourLines.h"
#include "cvfGeometryTools.h"
#include <algorithm>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigContourPolygonsTools::ContourPolygons
RigContourPolygonsTools::createContourPolygonsFromLineSegments( caf::ContourLines::ListOfLineSegments& unorderedLineSegments,
double contourValue,
double areaThreshold )
{
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 RigContourPolygonsTools::smoothContourPolygons( ContourPolygons& contourPolygons, bool favourExpansion, double sampleSpacing )
{
for ( RigContourPolygonsTools::ContourPolygon& polygon : contourPolygons )
{
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 RigContourPolygonsTools::clipContourPolygons( ContourPolygons& contourPolygons, const ContourPolygons& clipBy )
{
for ( RigContourPolygonsTools::ContourPolygon& polygon : contourPolygons )
{
for ( const RigContourPolygonsTools::ContourPolygon& clipPolygon : contourPolygons )
{
std::vector<std::vector<cvf::Vec3d>> intersections =
RigCellGeometryTools::intersectionWithPolygon( polygon.vertices, clipPolygon.vertices );
if ( !intersections.empty() )
{
polygon.vertices = intersections.front();
polygon.area = std::abs( cvf::GeometryTools::signedAreaPlanarPolygon( cvf::Vec3d::Z_AXIS, polygon.vertices ) );
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigContourPolygonsTools::sumPolygonArea( const ContourPolygons& contourPolygons )
{
double sumArea = 0.0;
for ( const ContourPolygon& polygon : contourPolygons )
{
sumArea += polygon.area;
}
return sumArea;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigContourPolygonsTools::lineOverlapsWithContourPolygons( const cvf::Vec3d& lineCenter,
const RigContourPolygonsTools::ContourPolygons& contourPolygons,
double tolerance )
{
const int64_t jump = 50;
for ( const RigContourPolygonsTools::ContourPolygon& edgePolygon : contourPolygons )
{
std::pair<int64_t, double> closestIndex( 0, std::numeric_limits<double>::infinity() );
for ( int64_t i = 0; i < (int64_t)edgePolygon.vertices.size(); i += jump )
{
const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[( i + 1 ) % edgePolygon.vertices.size()];
double dist1 = cvf::GeometryTools::linePointSquareDist( edgeVertex1, edgeVertex2, lineCenter );
if ( dist1 < tolerance )
{
return true;
}
if ( dist1 < closestIndex.second )
{
closestIndex = std::make_pair( i, dist1 );
}
}
for ( int64_t i = std::max( (int64_t)1, closestIndex.first - jump + 1 );
i < std::min( (int64_t)edgePolygon.vertices.size(), closestIndex.first + jump );
++i )
{
const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[( i + 1 ) % edgePolygon.vertices.size()];
double dist1 = cvf::GeometryTools::linePointSquareDist( edgeVertex1, edgeVertex2, lineCenter );
if ( dist1 < tolerance )
{
return true;
}
}
}
return false;
}

53
ApplicationLibCode/ReservoirDataModel/RigContourPolygonsTools.h Normal file
View File

@@ -0,0 +1,53 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "cafContourLines.h"
#include "cvfBoundingBox.h"
class RigContourMapGrid;
class RimGridView;
class RimRegularLegendConfig;
//==================================================================================================
///
///
//==================================================================================================
class RigContourPolygonsTools
{
public:
struct ContourPolygon
{
std::vector<cvf::Vec3d> vertices;
double value;
double area;
cvf::BoundingBox bbox;
};
using ContourPolygons = std::vector<ContourPolygon>;
static ContourPolygons createContourPolygonsFromLineSegments( caf::ContourLines::ListOfLineSegments& unorderedLineSegments,
double contourValue,
double areaThreshold );
static void smoothContourPolygons( ContourPolygons& contourPolygons, bool favourExpansion, double sampleSpacing );
static void clipContourPolygons( ContourPolygons& contourPolygons, const ContourPolygons& clipBy );
static double sumPolygonArea( const ContourPolygons& contourPolygons );
static bool lineOverlapsWithContourPolygons( const cvf::Vec3d& lineCenter, const ContourPolygons& contourPolygons, double tolerance );
};