mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
Refactor: Split out triangle generation from contour map
This commit is contained in:
Kristian Bendiksen
parent
9b6e441386
commit
a305464748
@ -18,12 +18,10 @@
|
||||
|
||||
#include "RimContourMapProjection.h"
|
||||
|
||||
#include "RiaOpenMPTools.h"
|
||||
|
||||
#include "RigCellGeometryTools.h"
|
||||
#include "RigContourMapCalculator.h"
|
||||
#include "RigContourMapGrid.h"
|
||||
#include "RigContourMapProjection.h"
|
||||
#include "RigContourMapTrianglesGenerator.h"
|
||||
|
||||
#include "RimCase.h"
|
||||
#include "RimGridView.h"
|
||||
@ -31,19 +29,15 @@
|
||||
#include "RimRegularLegendConfig.h"
|
||||
#include "RimTextAnnotation.h"
|
||||
|
||||
#include "cafContourLines.h"
|
||||
#include "cafPdmUiDoubleSliderEditor.h"
|
||||
#include "cafPdmUiTreeOrdering.h"
|
||||
#include "cafProgressInfo.h"
|
||||
|
||||
#include "cvfArray.h"
|
||||
#include "cvfGeometryUtils.h"
|
||||
#include "cvfScalarMapper.h"
|
||||
#include "cvfStructGridGeometryGenerator.h"
|
||||
#include "cvfVector2.h"
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
namespace caf
|
||||
{
|
||||
template <>
|
||||
@ -157,9 +151,39 @@ void RimContourMapProjection::generateGeometryIfNecessary()
|
||||
|
||||
if ( geometryNeedsUpdating() )
|
||||
{
|
||||
generateContourPolygons();
|
||||
progress.setProgress( 25 );
|
||||
generateTrianglesWithVertexValues();
|
||||
std::vector<double> contourLevels;
|
||||
|
||||
bool discrete = false;
|
||||
if ( legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER )
|
||||
{
|
||||
legendConfig()->scalarMapper()->majorTickValues( &contourLevels );
|
||||
|
||||
if ( resultRangeIsValid() )
|
||||
{
|
||||
std::tie( m_contourPolygons, m_contourLevelCumulativeAreas ) =
|
||||
RigContourMapTrianglesGenerator::generateContourPolygons( *m_contourMapGrid,
|
||||
*m_contourMapProjection,
|
||||
contourLevels,
|
||||
sampleSpacing(),
|
||||
sampleSpacingFactor(),
|
||||
m_smoothContourLines() );
|
||||
}
|
||||
progress.setProgress( 25 );
|
||||
|
||||
if ( legendConfig()->mappingMode() == RimRegularLegendConfig::MappingType::LINEAR_DISCRETE ||
|
||||
legendConfig()->mappingMode() == RimRegularLegendConfig::MappingType::LOG10_DISCRETE )
|
||||
{
|
||||
discrete = true;
|
||||
}
|
||||
}
|
||||
|
||||
m_trianglesWithVertexValues = RigContourMapTrianglesGenerator::generateTrianglesWithVertexValues( *m_contourMapGrid,
|
||||
*m_contourMapProjection,
|
||||
m_contourPolygons,
|
||||
contourLevels,
|
||||
m_contourLevelCumulativeAreas,
|
||||
discrete,
|
||||
sampleSpacing() );
|
||||
}
|
||||
progress.setProgress( 100 );
|
||||
}
|
||||
@ -507,333 +531,6 @@ void RimContourMapProjection::generateVertexResults()
|
||||
if ( m_contourMapProjection ) m_contourMapProjection->generateVertexResults();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RimContourMapProjection::generateTrianglesWithVertexValues()
|
||||
{
|
||||
std::vector<cvf::Vec3d> vertices = m_contourMapGrid->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 );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int numberOfThreads = RiaOpenMPTools::availableThreadCount();
|
||||
|
||||
std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( numberOfThreads );
|
||||
|
||||
const std::vector<double>& aggregatedVertexResults = m_contourMapProjection->aggregatedVertexResults();
|
||||
|
||||
#pragma omp parallel
|
||||
{
|
||||
int myThread = RiaOpenMPTools::currentThreadIndex();
|
||||
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 < aggregatedVertexResults.size() ? 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 = m_contourMapProjection->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 );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
// Add critical section here due to a weird bug when running in a single thread
|
||||
// Running multi threaded does not require this critical section, as we use a thread local data
|
||||
// structure
|
||||
#pragma omp critical
|
||||
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 = RigContourMapProjection::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;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RimContourMapProjection::generateContourPolygons()
|
||||
{
|
||||
std::vector<ContourPolygons> contourPolygons;
|
||||
|
||||
std::vector<double> contourLevels;
|
||||
if ( resultRangeIsValid() && legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER )
|
||||
{
|
||||
const std::vector<double>& aggregatedVertexResults = m_contourMapProjection->aggregatedVertexResults();
|
||||
|
||||
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( 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] = RigContourPolygonsTools::createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i],
|
||||
contourLevels[i],
|
||||
areaThreshold );
|
||||
|
||||
if ( m_smoothContourLines() )
|
||||
{
|
||||
RigContourPolygonsTools::smoothContourPolygons( contourPolygons[i], true, sampleSpacing() );
|
||||
}
|
||||
|
||||
for ( RigContourPolygonsTools::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 )
|
||||
{
|
||||
RigContourPolygonsTools::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 = RigContourPolygonsTools::sumPolygonArea( contourPolygons[i] );
|
||||
m_contourLevelCumulativeAreas[i] = levelOuterArea;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m_contourPolygons = contourPolygons;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
@ -132,8 +132,6 @@ protected:
|
||||
bool mapCellVisibilityNeedsUpdating( int timeStep );
|
||||
|
||||
void generateVertexResults();
|
||||
void generateTrianglesWithVertexValues();
|
||||
void generateContourPolygons();
|
||||
|
||||
double gridEdgeOffset() const;
|
||||
virtual void updateAfterResultGeneration( int timeStep ) = 0;
|
||||
|
||||
@ -85,6 +85,7 @@ set(SOURCE_GROUP_HEADER_FILES
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigWeightedMeanCalc.h
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapCalculator.h
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapProjection.h
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapTrianglesGenerator.h
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigEclipseContourMapProjection.h
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigGeoMechContourMapProjection.h
|
||||
)
|
||||
@ -174,6 +175,7 @@ set(SOURCE_GROUP_SOURCE_FILES
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourPolygonsTools.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapCalculator.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapProjection.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigContourMapTrianglesGenerator.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigEclipseContourMapProjection.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigGeoMechContourMapProjection.cpp
|
||||
)
|
||||
|
||||
@ -0,0 +1,354 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// 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 "RigContourMapTrianglesGenerator.h"
|
||||
|
||||
#include "RiaOpenMPTools.h"
|
||||
|
||||
#include "RigCellGeometryTools.h"
|
||||
#include "RigContourMapGrid.h"
|
||||
#include "RigContourMapProjection.h"
|
||||
|
||||
#include "cvfGeometryUtils.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<cvf::Vec4d>
|
||||
RigContourMapTrianglesGenerator::generateTrianglesWithVertexValues( const RigContourMapGrid& contourMapGrid,
|
||||
const RigContourMapProjection& contourMapProjection,
|
||||
const std::vector<ContourPolygons>& contourPolygons,
|
||||
const std::vector<double>& contourLevels,
|
||||
const std::vector<double>& contourLevelCumulativeAreas,
|
||||
bool discrete,
|
||||
double sampleSpacing )
|
||||
{
|
||||
std::vector<cvf::Vec3d> vertices = contourMapGrid.generateVertices();
|
||||
|
||||
cvf::Vec2ui patchSize = contourMapGrid.numberOfVerticesIJ();
|
||||
cvf::ref<cvf::UIntArray> faceList = new cvf::UIntArray;
|
||||
cvf::GeometryUtils::tesselatePatchAsTriangles( patchSize.x(), patchSize.y(), 0u, true, faceList.p() );
|
||||
|
||||
const double cellArea = sampleSpacing * sampleSpacing;
|
||||
const double areaThreshold = 1.0e-5 * 0.5 * cellArea;
|
||||
|
||||
std::vector<std::vector<std::vector<cvf::Vec3d>>> subtractPolygons;
|
||||
if ( !contourPolygons.empty() )
|
||||
{
|
||||
subtractPolygons.resize( contourPolygons.size() );
|
||||
for ( size_t i = 0; i < contourPolygons.size() - 1; ++i )
|
||||
{
|
||||
for ( size_t j = 0; j < contourPolygons[i + 1].size(); ++j )
|
||||
{
|
||||
subtractPolygons[i].push_back( contourPolygons[i + 1][j].vertices );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int numberOfThreads = RiaOpenMPTools::availableThreadCount();
|
||||
|
||||
std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( numberOfThreads );
|
||||
|
||||
const std::vector<double>& aggregatedVertexResults = contourMapProjection.aggregatedVertexResults();
|
||||
|
||||
#pragma omp parallel
|
||||
{
|
||||
int myThread = RiaOpenMPTools::currentThreadIndex();
|
||||
threadTriangles[myThread].resize( std::max( (size_t)1, 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 < aggregatedVertexResults.size() ? 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 ( contourPolygons.empty() )
|
||||
{
|
||||
threadTriangles[myThread][0].insert( threadTriangles[myThread][0].end(), triangleWithValues.begin(), triangleWithValues.end() );
|
||||
continue;
|
||||
}
|
||||
|
||||
bool outsideOuterLimit = false;
|
||||
for ( size_t c = 0; c < contourPolygons.size() && !outsideOuterLimit; ++c )
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d>> intersectPolygons;
|
||||
for ( size_t j = 0; j < contourPolygons[c].size(); ++j )
|
||||
{
|
||||
bool containsAtLeastOne = false;
|
||||
for ( size_t t = 0; t < 3; ++t )
|
||||
{
|
||||
if ( contourPolygons[c][j].bbox.contains( triangle[t] ) )
|
||||
{
|
||||
containsAtLeastOne = true;
|
||||
}
|
||||
}
|
||||
if ( containsAtLeastOne )
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d>> clippedPolygons =
|
||||
RigCellGeometryTools::intersectionWithPolygon( triangle, 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 = contourMapProjection.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 );
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
// Add critical section here due to a weird bug when running in a single thread
|
||||
// Running multi threaded does not require this critical section, as we use a thread local data
|
||||
// structure
|
||||
#pragma omp critical
|
||||
threadTriangles[myThread][c].insert( threadTriangles[myThread][c].end(),
|
||||
clippedTriangles.begin(),
|
||||
clippedTriangles.end() );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<std::vector<cvf::Vec4d>> trianglesPerLevel( std::max( (size_t)1, 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 = RigContourMapProjection::sumTriangleAreas( allTrianglesThisLevel );
|
||||
if ( c >= contourLevelCumulativeAreas.size() || triangleAreasThisLevel > 1.0e-3 * 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() );
|
||||
}
|
||||
|
||||
return finalTriangles;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::pair<std::vector<RigContourMapTrianglesGenerator::ContourPolygons>, std::vector<double>>
|
||||
RigContourMapTrianglesGenerator::generateContourPolygons( const RigContourMapGrid& contourMapGrid,
|
||||
const RigContourMapProjection& contourMapProjection,
|
||||
const std::vector<double>& initialContourLevels,
|
||||
double sampleSpacing,
|
||||
double sampleSpacingFactor,
|
||||
bool smoothContourLines )
|
||||
{
|
||||
if ( contourMapProjection.minValue() != std::numeric_limits<double>::infinity() &&
|
||||
contourMapProjection.maxValue() != -std::numeric_limits<double>::infinity() &&
|
||||
std::fabs( contourMapProjection.maxValue() - contourMapProjection.minValue() ) > 1.0e-8 )
|
||||
{
|
||||
std::vector<double> contourLevels = initialContourLevels;
|
||||
|
||||
int nContourLevels = static_cast<int>( contourLevels.size() );
|
||||
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 = smoothContourLines ? 5.0e-2 * sampleSpacing : 1.0e-3 * sampleSpacing;
|
||||
|
||||
if ( nContourLevels >= 10 )
|
||||
{
|
||||
simplifyEpsilon *= 2.0;
|
||||
}
|
||||
if ( contourMapGrid.numberOfCells() > 100000 )
|
||||
{
|
||||
simplifyEpsilon *= 2.0;
|
||||
}
|
||||
else if ( contourMapGrid.numberOfCells() > 1000000 )
|
||||
{
|
||||
simplifyEpsilon *= 4.0;
|
||||
}
|
||||
|
||||
const std::vector<double>& aggregatedVertexResults = contourMapProjection.aggregatedVertexResults();
|
||||
|
||||
std::vector<caf::ContourLines::ListOfLineSegments> unorderedLineSegmentsPerLevel =
|
||||
caf::ContourLines::create( aggregatedVertexResults,
|
||||
contourMapProjection.xVertexPositions(),
|
||||
contourMapProjection.yVertexPositions(),
|
||||
contourLevels );
|
||||
|
||||
std::vector<ContourPolygons> 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] = RigContourPolygonsTools::createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i],
|
||||
contourLevels[i],
|
||||
areaThreshold );
|
||||
|
||||
if ( smoothContourLines )
|
||||
{
|
||||
RigContourPolygonsTools::smoothContourPolygons( contourPolygons[i], true, sampleSpacing );
|
||||
}
|
||||
|
||||
for ( RigContourPolygonsTools::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 )
|
||||
{
|
||||
RigContourPolygonsTools::clipContourPolygons(&contourPolygons[i], &contourPolygons[i - 1] );
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
std::vector<double> contourLevelCumulativeAreas( contourPolygons.size(), 0.0 );
|
||||
for ( int64_t i = (int64_t)contourPolygons.size() - 1; i >= 0; --i )
|
||||
{
|
||||
double levelOuterArea = RigContourPolygonsTools::sumPolygonArea( contourPolygons[i] );
|
||||
contourLevelCumulativeAreas[i] = levelOuterArea;
|
||||
}
|
||||
|
||||
return { contourPolygons, contourLevelCumulativeAreas };
|
||||
}
|
||||
}
|
||||
|
||||
return {};
|
||||
}
|
||||
@ -0,0 +1,52 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// 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 <vector>
|
||||
|
||||
#include "RigContourPolygonsTools.h"
|
||||
|
||||
class RigContourMapGrid;
|
||||
class RigContourMapProjection;
|
||||
|
||||
//==================================================================================================
|
||||
///
|
||||
///
|
||||
//==================================================================================================
|
||||
class RigContourMapTrianglesGenerator
|
||||
{
|
||||
public:
|
||||
using ContourPolygons = std::vector<RigContourPolygonsTools::ContourPolygon>;
|
||||
|
||||
static std::vector<cvf::Vec4d> generateTrianglesWithVertexValues( const RigContourMapGrid& contourMapGrid,
|
||||
const RigContourMapProjection& contourMapProjection,
|
||||
const std::vector<ContourPolygons>& contourPolygons,
|
||||
const std::vector<double>& contourLevels,
|
||||
const std::vector<double>& contourLevelCumulativeAreas,
|
||||
bool discrete,
|
||||
double sampleSpacing );
|
||||
|
||||
static std::pair<std::vector<ContourPolygons>, std::vector<double>>
|
||||
generateContourPolygons( const RigContourMapGrid& contourMapGrid,
|
||||
const RigContourMapProjection& contourMapProjection,
|
||||
const std::vector<double>& initialContourLevels,
|
||||
double sampleSpacing,
|
||||
double sampleSpacingFactor,
|
||||
bool smoothContourLines );
|
||||
};
|
||||
Loading…
Reference in New Issue
Block a user