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

#12170 Contour Map: Create multiple polygons based on area limit

Browse Source
This commit is contained in:
Magne Sjaastad
2025-02-19 15:40:58 +01:00
parent 346670f349
commit 7a82c2f145
6 changed files with 167 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ApplicationLibCode/Commands/RicCreateContourMapPolygonAdvancedFeature.cpp
View File

@@ -45,7 +45,7 @@ void RicCreateContourMapPolygonAdvancedFeature::onActionTriggered( bool isChecke
auto processedImage = dlg.processedImageData();
if ( processedImage.empty() ) return;
RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonFromImage( processedImage, rigContourMapProjection );
RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonsFromImage( processedImage, rigContourMapProjection );
}
//--------------------------------------------------------------------------------------------------

2
ApplicationLibCode/Commands/RicCreateContourMapPolygonFeature.cpp
View File

@@ -43,7 +43,7 @@ void RicCreateContourMapPolygonFeature::onActionTriggered( bool isChecked )
if ( dilated.empty() ) return;
RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonFromImage( dilated, rigContourMapProjection );
RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonsFromImage( dilated, rigContourMapProjection );
}
//--------------------------------------------------------------------------------------------------

100
ApplicationLibCode/Commands/RicCreateContourMapPolygonTools.cpp
View File

@@ -193,46 +193,86 @@ std::vector<std::vector<int>> RicCreateContourMapPolygonTools::convertToBinaryIm
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimPolygon* RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonFromImage( std::vector<std::vector<int>> image,
const RigContourMapProjection* contourMapProjection )
RimPolygon* RicCreateContourMapPolygonTools::createAndAddBoundaryPolygonsFromImage( std::vector<std::vector<int>> image,
const RigContourMapProjection* contourMapProjection )
{
if ( !contourMapProjection ) return nullptr;
if ( image.empty() ) return nullptr;
std::vector<cvf::Vec3d> polygonDomainCoords;
auto xVertexPositions = contourMapProjection->xVertexPositions();
auto yVertexPositions = contourMapProjection->yVertexPositions();
auto origin3d = contourMapProjection->origin3d();
auto depth = contourMapProjection->topDepthBoundingBox();
auto xVertexPositions = contourMapProjection->xVertexPositions();
auto yVertexPositions = contourMapProjection->yVertexPositions();
auto origin3d = contourMapProjection->origin3d();
auto depth = contourMapProjection->topDepthBoundingBox();
auto boundaryPoints = RigPolygonTools::boundary( image );
for ( auto [i, j] : boundaryPoints )
auto hasAnyValue = []( const std::vector<std::vector<int>>& image ) -> bool
{
double xDomain = xVertexPositions.at( i ) + origin3d.x();
double yDomain = yVertexPositions.at( j ) + origin3d.y();
for ( size_t i = 0; i < image.size(); ++i )
{
for ( size_t j = 0; j < image[i].size(); ++j )
{
if ( image[i][j] > 0 )
{
return true;
}
}
}
polygonDomainCoords.emplace_back( cvf::Vec3d( xDomain, yDomain, depth ) );
return false;
};
auto polygonCollection = RimTools::polygonCollection();
RimPolygon* currentPolygon = nullptr;
while ( hasAnyValue( image ) )
{
std::vector<cvf::Vec3d> polygonDomainCoords;
auto boundaryPoints = RigPolygonTools::boundary( image );
// The area threshold defines the minimum number of pixels to create a polygon for
const auto area = RigPolygonTools::area( boundaryPoints );
const double areaThreshold = 10;
if ( area > areaThreshold )
{
for ( const auto& [i, j] : boundaryPoints )
{
auto xDomain = xVertexPositions.at( i ) + origin3d.x();
auto yDomain = yVertexPositions.at( j ) + origin3d.y();
polygonDomainCoords.emplace_back( cvf::Vec3d( xDomain, yDomain, depth ) );
}
// Epsilon used to simplify polygon. Useful range typical value in [5..50]
const double defaultEpsilon = 40.0;
RigPolygonTools::simplifyPolygon( polygonDomainCoords, defaultEpsilon );
if ( polygonDomainCoords.size() >= 3 )
{
auto newPolygon = polygonCollection->appendUserDefinedPolygon();
newPolygon->setPointsInDomainCoords( polygonDomainCoords );
newPolygon->coordinatesChanged.send();
currentPolygon = newPolygon;
}
}
// Subtract all pixels inside the polygon by setting their value to 0
const int value = 0;
image = RigPolygonTools::assignValueInsidePolygon( image, boundaryPoints, value );
/*
// NB: Do not remove, debug code to export images to file
{
QString filename = "boundary_removed";
filename = QString( "f:/scratch/images/%1%2.png" ).arg( filename ).arg( QString::number( imageCounter ) );
RicCreateContourMapPolygonTools::exportVectorAsGrayscaleImage( currentImage, filename );
}
*/
}
// Epsilon used to simplify polygon. Useful range typical value in [5..50]
const double defaultEpsilon = 40.0;
RigPolygonTools::simplifyPolygon( polygonDomainCoords, defaultEpsilon );
polygonCollection->uiCapability()->updateAllRequiredEditors();
if ( polygonDomainCoords.size() >= 3 )
{
auto polygonCollection = RimTools::polygonCollection();
auto newPolygon = polygonCollection->appendUserDefinedPolygon();
newPolygon->setPointsInDomainCoords( polygonDomainCoords );
newPolygon->coordinatesChanged.send();
polygonCollection->uiCapability()->updateAllRequiredEditors();
return newPolygon;
}
return nullptr;
return currentPolygon;
}
//--------------------------------------------------------------------------------------------------

2
ApplicationLibCode/Commands/RicCreateContourMapPolygonTools.h
View File

@@ -40,7 +40,7 @@ void exportVectorAsGrayscaleImage( const std::vector<std::vector<int>>& data, co
std::vector<std::vector<int>> convertImageToBinary( QImage image );
std::vector<std::vector<int>> convertToBinaryImage( const RigContourMapProjection* contourMapProjection );
RimPolygon* createAndAddBoundaryPolygonFromImage( std::vector<std::vector<int>> image, const RigContourMapProjection* contourMapProjection );
RimPolygon* createAndAddBoundaryPolygonsFromImage( std::vector<std::vector<int>> image, const RigContourMapProjection* contourMapProjection );
const RigContourMapProjection* findCurrentContourMapProjection();

100
ApplicationLibCode/ReservoirDataModel/RigPolygonTools.cpp
View File

@@ -53,9 +53,9 @@ namespace internal
{
if ( !isValidImage( image ) ) return;
auto rows = static_cast<int>( image.size() );
auto cols = static_cast<int>( image[0].size() );
std::stack<std::pair<int, int>> stack;
auto rows = static_cast<int>( image.size() );
auto cols = static_cast<int>( image[0].size() );
std::stack<Point> stack;
stack.push( { x, y } );
while ( !stack.empty() )
@@ -74,6 +74,45 @@ namespace internal
}
}
// Function to check if a point is on a line segment (edge of the polygon)
bool isOnSegment( Point p, Point p1, Point p2 )
{
int x = p.first, y = p.second;
int x1 = p1.first, y1 = p1.second;
int x2 = p2.first, y2 = p2.second;
// Check if the point (x, y) lies between (x1, y1) and (x2, y2)
return ( ( x >= std::min( x1, x2 ) && x <= std::max( x1, x2 ) ) && ( y >= std::min( y1, y2 ) && y <= std::max( y1, y2 ) ) &&
( ( x2 - x1 ) * ( y - y1 ) == ( y2 - y1 ) * ( x - x1 ) ) ); // Collinearity check
}
// Check if a point is inside a polygon using the Ray-Casting Algorithm
bool isInsidePolygon( const Point& p, const std::vector<Point>& polygon )
{
int n = static_cast<int>( polygon.size() );
int count = 0;
for ( int i = 0; i < n; i++ )
{
auto p1 = polygon[i];
auto p2 = polygon[( i + 1 ) % n]; // Next vertex (looping back to first at end)
// Check if point is exactly on an edge
if ( isOnSegment( p, p1, p2 ) ) return true;
// Check if point is between y-bounds of edge
if ( ( p1.second > p.second ) != ( p2.second > p.second ) )
{
// Compute intersection point of the edge with horizontal line at p.y
double xIntersect = p1.first + (double)( p.second - p1.second ) * ( p2.first - p1.first ) / ( p2.second - p1.second );
if ( p.first < xIntersect )
{
count++;
}
}
}
return ( count % 2 ) == 1; // Odd count means inside, even means outside
}
}; // namespace internal
//--------------------------------------------------------------------------------------------------
@@ -197,18 +236,18 @@ IntegerImage fillInterior( IntegerImage sourceImage )
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<int, int>> boundary( const IntegerImage& image )
std::vector<Point> boundary( const IntegerImage& image )
{
if ( !internal::isValidImage( image ) ) return {};
std::vector<std::pair<int, int>> boundaries;
std::vector<Point> boundaries;
// Get dimensions of the image
int rows = static_cast<int>( image.size() );
int cols = static_cast<int>( image[0].size() );
// Direction vectors for clockwise search (8-connectivity)
const std::vector<std::pair<int, int>> directions = { { -1, 0 }, { -1, 1 }, { 0, 1 }, { 1, 1 }, { 1, 0 }, { 1, -1 }, { 0, -1 }, { -1, -1 } };
const std::vector<Point> directions = { { -1, 0 }, { -1, 1 }, { 0, 1 }, { 1, 1 }, { 1, 0 }, { 1, -1 }, { 0, -1 }, { -1, -1 } };
// Helper lambda to check if a pixel is a valid boundary pixel
auto isBoundaryPixel = [&]( int x, int y )
@@ -227,7 +266,7 @@ std::vector<std::pair<int, int>> boundary( const IntegerImage& image )
};
// Find the starting boundary pixel
std::pair<int, int> start( -1, -1 );
Point start( -1, -1 );
for ( int row = 0; row < rows; ++row )
{
for ( int col = 0; col < cols; ++col )
@@ -244,8 +283,8 @@ std::vector<std::pair<int, int>> boundary( const IntegerImage& image )
if ( start.first == -1 ) return boundaries; // No boundary found
// Contour following algorithm
std::pair<int, int> current = start;
int direction = 0; // Start search direction (arbitrary)
Point current = start;
int direction = 0; // Start search direction (arbitrary)
do
{
boundaries.push_back( current );
@@ -275,6 +314,49 @@ std::vector<std::pair<int, int>> boundary( const IntegerImage& image )
return boundaries;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigPolygonTools::IntegerImage assignValueInsidePolygon( IntegerImage image, const std::vector<Point>& polygon, int value )
{
if ( !internal::isValidImage( image ) ) return {};
auto rows = static_cast<int>( image.size() );
auto cols = static_cast<int>( image[0].size() );
for ( int i = 0; i < rows; ++i )
{
for ( int j = 0; j < cols; ++j )
{
if ( internal::isInsidePolygon( { i, j }, polygon ) )
{
image[i][j] = value;
}
}
}
return image;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double area( const std::vector<Point>& polygon )
{
int n = static_cast<int>( polygon.size() );
if ( n < 3 ) return 0.0; // A polygon must have at least 3 points
double area = 0.0;
for ( int i = 0; i < n; i++ )
{
int j = ( i + 1 ) % n; // Next vertex, wrapping around at the end
area += polygon[i].first * polygon[j].second;
area -= polygon[j].first * polygon[i].second;
}
return std::fabs( area ) / 2.0;
}
//--------------------------------------------------------------------------------------------------
/// Ramer-Douglas-Peucker simplification algorithm
///

3
ApplicationLibCode/ReservoirDataModel/RigPolygonTools.h
View File

@@ -26,12 +26,15 @@ namespace RigPolygonTools
{
// Integer images are assumed to be 2D arrays of integers, where 0 is background and 1 is foreground.
using IntegerImage = std::vector<std::vector<int>>;
using Point = std::pair<int, int>;
IntegerImage erode( IntegerImage image, int kernelSize );
IntegerImage dilate( IntegerImage image, int kernelSize );
IntegerImage fillInterior( IntegerImage sourceImage );
std::vector<std::pair<int, int>> boundary( const IntegerImage& image );
IntegerImage assignValueInsidePolygon( IntegerImage image, const std::vector<Point>& polygon, int value );
double area( const std::vector<Point>& polygon );
// Recursive function modifying the incoming vertices
void simplifyPolygon( std::vector<cvf::Vec3d>& vertices, double epsilon );