OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-01-01 03:37:15 -06:00
Code Issues Packages Projects Releases Wiki Activity

Prevent duplicated nodes in vertex array of intersection response

Browse Source
This commit is contained in:
Jørgen Herje 2024-04-04 09:10:38 +02:00
parent 4d87fc5927
commit 614a9f1213
7 changed files with 249 additions and 448 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
ApplicationLibCode/ModelVisualization/Intersections/RivEnclosingPolygonGenerator.cpp
View File

@ -1,6 +1,6 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2018- Equinor ASA
// 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
@ -21,110 +21,7 @@
#include <map>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
PolygonVertexWelder::PolygonVertexWelder( double weldEpsilon )
: m_epsilonSquared( weldEpsilon * weldEpsilon )
, m_first( cvf::UNDEFINED_UINT )
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void PolygonVertexWelder::reserveVertices( cvf::uint vertexCount )
{
m_vertex.reserve( vertexCount );
m_next.reserve( vertexCount );
}
//--------------------------------------------------------------------------------------------------
/// Add a vertex to the welder. If the vertex is within the tolerance of an existing vertex, the existing
// vertex index is returned
//--------------------------------------------------------------------------------------------------
cvf::uint PolygonVertexWelder::weldVertexAndGetIndex( const cvf::Vec3d& vertex )
{
// Call function to step through linked list of bucket, testing
// if vertex is within the epsilon of one of the vertices in the bucket
cvf::uint indexOfLocatedVertex = locateVertexInPolygon( vertex );
if ( indexOfLocatedVertex != cvf::UNDEFINED_UINT )
{
// if ( wasWelded ) *wasWelded = true;
return indexOfLocatedVertex;
}
// Vertex not found in epsilon neighborhood, add it to the list
cvf::uint indexOfAddedVertex = addVertexToPolygon( vertex );
return indexOfAddedVertex;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const cvf::Vec3d& PolygonVertexWelder::vertex( cvf::uint index ) const
{
return m_vertex[index];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::Vec3dArray> PolygonVertexWelder::createVertexArray() const
{
cvf::ref<cvf::Vec3dArray> vertexArray = new cvf::Vec3dArray( m_vertex );
return vertexArray;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::uint PolygonVertexWelder::locateVertexInPolygon( const cvf::Vec3d& vertex ) const
{
cvf::uint currentIndex = m_first;
while ( currentIndex != cvf::UNDEFINED_UINT )
{
// Weld point within tolerance
const auto distanceSquared = ( m_vertex[currentIndex] - vertex ).lengthSquared();
if ( distanceSquared < m_epsilonSquared )
{
return currentIndex;
}
currentIndex = m_next[currentIndex];
}
// No vertex found to weld to
return cvf::UNDEFINED_UINT;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::uint PolygonVertexWelder::addVertexToPolygon( const cvf::Vec3d& vertex )
{
// Add vertex and update linked list
m_vertex.push_back( vertex );
m_next.push_back( m_first );
CVF_TIGHT_ASSERT( m_vertex.size() == m_next.size() );
// Update index of first vertex
cvf::uint indexOfAddedVertex = static_cast<cvf::uint>( m_vertex.size() - 1 );
m_first = indexOfAddedVertex;
return indexOfAddedVertex;
}
//--------------------------------------------------------------------------------------------------
///
/// ************************************************************************************************
/// ************************************************************************************************
/// ************************************************************************************************
/// ************************************************************************************************
///
//--------------------------------------------------------------------------------------------------
RivEnclosingPolygonGenerator::RivEnclosingPolygonGenerator()
: m_polygonVertexWelder( 1e-3 )
{
}
@ -221,12 +118,7 @@ void RivEnclosingPolygonGenerator::constructEnclosingPolygon()
}
}
// Convert vertex indices to vertices
m_polygonVertices.clear();
for ( cvf::uint vertexIndex : enclosingPolygonVertexIndices )
{
m_polygonVertices.push_back( m_polygonVertexWelder.vertex( vertexIndex ) );
}
m_polygonVertexIndices = enclosingPolygonVertexIndices;
}
//--------------------------------------------------------------------------------------------------
@ -249,9 +141,9 @@ cvf::EdgeKey RivEnclosingPolygonGenerator::findNextEdge( cvf::uint vertexIndex,
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RivEnclosingPolygonGenerator::getPolygonVertices() const
std::vector<cvf::uint> RivEnclosingPolygonGenerator::getPolygonVertexIndices() const
{
return m_polygonVertices;
return m_polygonVertexIndices;
}
//--------------------------------------------------------------------------------------------------
@ -263,24 +155,21 @@ bool RivEnclosingPolygonGenerator::isValidPolygon() const
}
//--------------------------------------------------------------------------------------------------
/// Add triangle vertices to the polygon. The vertices are welded to prevent duplicated vertices
/// Add triangle vertices to the polygon. The vertex indices are welded and controlled outside
/// of this class.
///
/// Assumes that the vertices are given in counter-clockwise order
/// Assumes the vertices are given in counter-clockwise order
//--------------------------------------------------------------------------------------------------
void RivEnclosingPolygonGenerator::addTriangleVertices( const cvf::Vec3d& p0, const cvf::Vec3d& p1, const cvf::Vec3d& p2 )
void RivEnclosingPolygonGenerator::addTriangleVertexIndices( cvf::uint idxVx0, cvf::uint idxVx1, cvf::uint idxVx2 )
{
cvf::uint i0 = m_polygonVertexWelder.weldVertexAndGetIndex( p0 );
cvf::uint i1 = m_polygonVertexWelder.weldVertexAndGetIndex( p1 );
cvf::uint i2 = m_polygonVertexWelder.weldVertexAndGetIndex( p2 );
// Verify three unique vertices - i.e. no degenerate triangle
if ( i0 == i1 || i0 == i2 || i1 == i2 )
if ( idxVx0 == idxVx1 || idxVx0 == idxVx2 || idxVx1 == idxVx2 )
{
return;
}
// Add edges keys to list of all edges
m_allEdgeKeys.emplace_back( cvf::EdgeKey( i0, i1 ).toKeyVal() );
m_allEdgeKeys.emplace_back( cvf::EdgeKey( i1, i2 ).toKeyVal() );
m_allEdgeKeys.emplace_back( cvf::EdgeKey( i2, i0 ).toKeyVal() );
m_allEdgeKeys.emplace_back( cvf::EdgeKey( idxVx0, idxVx1 ).toKeyVal() );
m_allEdgeKeys.emplace_back( cvf::EdgeKey( idxVx1, idxVx2 ).toKeyVal() );
m_allEdgeKeys.emplace_back( cvf::EdgeKey( idxVx2, idxVx0 ).toKeyVal() );
}

49
ApplicationLibCode/ModelVisualization/Intersections/RivEnclosingPolygonGenerator.h
View File

@ -1,6 +1,6 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2018- Equinor ASA
// 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
@ -18,67 +18,36 @@
#pragma once
#include "cafLine.h"
#include "cvfArray.h"
#include "cvfBase.h"
#include "cvfEdgeKey.h"
#include "cvfObject.h"
#include "cvfVector3.h"
#include <algorithm>
#include <set>
#include <vector>
/*
* Class for handling welding of vertices internally in a polygon to prevent duplicated vertex 3D points within a tolerance margin
*/
class PolygonVertexWelder
{
public:
PolygonVertexWelder( double weldEpsilon );
void reserveVertices( cvf::uint vertexCount );
cvf::uint weldVertexAndGetIndex( const cvf::Vec3d& vertex );
const cvf::Vec3d& vertex( cvf::uint index ) const;
cvf::ref<cvf::Vec3dArray> createVertexArray() const;
private:
cvf::uint locateVertexInPolygon( const cvf::Vec3d& vertex ) const;
cvf::uint addVertexToPolygon( const cvf::Vec3d& vertex );
private:
const double m_epsilonSquared; // Tolerance for vertex welding, radius around vertex defining welding neighborhood
cvf::uint m_first; // Start of linked list
std::vector<cvf::uint> m_next; // Links each vertex to next in linked list. Always numVertices long, will grow as vertices are added
std::vector<cvf::Vec3d> m_vertex; // Unique vertices within tolerance
};
/*
* Class for generating an enclosing polygon from a set of vertices.
*
* The class will weld triangle vertices close to each other and provide a vertex index for
* the resulting set of vertices. These indices are used for algorithms constructing the enclosing polygon.
* The class add triangle vertex indices and construct constructing the enclosing polygon from the indices.
* The vertices must be provided in counter clock-wise order.
*
* The welding is done using a tolerance value to handle floating point errors.
* The vertex welding, and mapping of vertex indices to vertex 3D points, must be handled externally.
*/
class RivEnclosingPolygonGenerator
{
public:
RivEnclosingPolygonGenerator();
std::vector<cvf::Vec3d> getPolygonVertices() const;
std::vector<cvf::uint> getPolygonVertexIndices() const;
bool isValidPolygon() const;
void addTriangleVertices( const cvf::Vec3d& p0, const cvf::Vec3d& p1, const cvf::Vec3d& p2 );
void addTriangleVertexIndices( cvf::uint idxVx0, cvf::uint idxVx1, cvf::uint idxVx2 );
void constructEnclosingPolygon();
private:
static cvf::EdgeKey findNextEdge( cvf::uint vertextIndex, const std::set<cvf::EdgeKey>& boundaryEdges );
private:
PolygonVertexWelder m_polygonVertexWelder; // Add and weld vertices for a polygon, provides vertex index
std::vector<cvf::int64> m_allEdgeKeys; // Create edge defined by vertex indices when adding triangle. Using cvf::EdgeKey::toKeyVal()
std::vector<cvf::Vec3d> m_polygonVertices; // List polygon vertices counter clock-wise
std::vector<cvf::uint> m_polygonVertexIndices; // List polygon vertex indices counter clock-wise
};

401
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.cpp
View File

@ -36,7 +36,6 @@ RivPolylineIntersectionGeometryGenerator::RivPolylineIntersectionGeometryGenerat
: m_polylineUtm( initializePolylineUtmFromPolylineUtmXy( polylineUtmXy ) )
, m_hexGrid( grid )
{
m_polygonVertices = new cvf::Vec3fArray;
}
//--------------------------------------------------------------------------------------------------
@ -63,41 +62,6 @@ std::vector<cvf::Vec3d>
return polylineUtm;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RivPolylineIntersectionGeometryGenerator::isAnyGeometryPresent() const
{
return m_polylineSegmentsMeshData.size() > 0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const cvf::Vec3fArray* RivPolylineIntersectionGeometryGenerator::polygonVxes() const
{
CVF_ASSERT( m_polygonVertices->size() > 0 );
return m_polygonVertices.p();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::vertiesPerPolygon() const
{
CVF_ASSERT( m_verticesPerPolygon.size() > 0 );
return m_verticesPerPolygon;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::polygonToCellIndex() const
{
CVF_ASSERT( m_polygonToCellIdxMap.size() > 0 );
return m_polygonToCellIdxMap;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -107,6 +71,14 @@ const std::vector<PolylineSegmentMeshData>& RivPolylineIntersectionGeometryGener
return m_polylineSegmentsMeshData;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RivPolylineIntersectionGeometryGenerator::isAnyGeometryPresent() const
{
return m_polylineSegmentsMeshData.size() > 0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -121,10 +93,7 @@ void RivPolylineIntersectionGeometryGenerator::generateIntersectionGeometry( cvf
void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray* visibleCells )
{
if ( m_hexGrid == nullptr || m_polylineSegmentsMeshData.size() != 0 ) return;
// Mesh data per polyline segment
std::vector<PolylineSegmentMeshData> polylineSegmentMeshData = {};
std::vector<cvf::Vec3f> calculatedPolygonVertices = {};
if ( m_polylineUtm.size() < 2 ) return;
cvf::BoundingBox gridBBox = m_hexGrid->boundingBox();
const double topDepth = gridBBox.max().z();
@ -132,181 +101,191 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
const auto zAxisDirection = -cvf::Vec3d::Z_AXIS; // NOTE: Negative or positive direction?
const cvf::Vec3d maxHeightVec = zAxisDirection * gridBBox.radius();
if ( m_polylineUtm.size() > 1 )
// Weld vertices per polyline segment
// - Low welding distance, as the goal is to weld duplicate vertices
// - Number of buckets is set per segment, utilizing number of cells intersecting the segment
const double weldingDistance = 1.0e-3;
const double weldingCellSize = 4.0 * weldingDistance;
const size_t numPoints = m_polylineUtm.size();
size_t pointIdx = 0;
// Loop over polyline segments
//
// Create intersection geometry for each polyline segment and clip triangles outside the
// polyline segment, using UTM-coordinates.
//
// Afterwards convert to "local" coordinates (p1 as origin), where the local uz-coordinates
// for each vertex is calculated using Pythagoras theorem.
//
// As the plane is parallel to the z-axis, the local uz-coordinates per polyline segment
// can be calculated using Pythagoras theorem. Where a segment is defined between p1 and p2,
// which implies p1 is origin of the local coordinate system uz.
//
// For a segment a UTM coordinate (x_utm,y_utm,z_utm) converts to u = sqrt(x^2 + y^2) and z = z.
// Where x, y and z are local vertex coordinates, after subtracting the (x_utm, y_utm, z_utm)-values
// for p1 from the vertex UTM-coordinates.
while ( pointIdx < numPoints - 1 )
{
const size_t numPoints = m_polylineUtm.size();
size_t pointIdx = 0;
// Loop over polyline segments
//
// Create intersection geometry for each polyline segment and clip triangles outside the
// polyline segment, using UTM-coordinates.
//
// Afterwards convert to "local" coordinates (p1 as origin), where the local uz-coordinates
// for each vertex is calculated using Pythagoras theorem.
//
// As the plane is parallel to the z-axis, the local uz-coordinates per polyline segment
// can be calculated using Pythagoras theorem. Where a segment is defined between p1 and p2,
// which implies p1 is origin of the local coordinate system uz.
//
// For a segment a UTM coordinate (x_utm,y_utm,z_utm) converts to u = sqrt(x^2 + y^2) and z = z.
// Where x, y and z are local vertex coordinates, after subtracting the (x_utm, y_utm, z_utm)-values
// for p1 from the vertex UTM-coordinates.
while ( pointIdx < numPoints - 1 )
size_t nextPointIdx = RivSectionFlattener::indexToNextValidPoint( m_polylineUtm, zAxisDirection, pointIdx );
if ( nextPointIdx == size_t( -1 ) || nextPointIdx >= m_polylineUtm.size() )
{
size_t nextPointIdx = RivSectionFlattener::indexToNextValidPoint( m_polylineUtm, zAxisDirection, pointIdx );
if ( nextPointIdx == size_t( -1 ) || nextPointIdx >= m_polylineUtm.size() )
{
break;
}
// Start and end point of polyline segment in UTM-coordinates
const cvf::Vec3d p1 = m_polylineUtm[pointIdx];
const cvf::Vec3d p2 = m_polylineUtm[nextPointIdx];
// Get cell candidates for the polyline segment (subset of global cells)
std::vector<size_t> columnCellCandidates =
createPolylineSegmentCellCandidates( *m_hexGrid, p1, p2, maxHeightVec, topDepth, bottomDepth );
// Plane for the polyline segment
cvf::Plane plane;
plane.setFromPoints( p1, p2, p2 + maxHeightVec );
// Planes parallel to z-axis for p1 and p2, to prevent triangles outside the polyline segment
cvf::Plane p1Plane;
p1Plane.setFromPoints( p1, p1 + maxHeightVec, p1 + plane.normal() );
cvf::Plane p2Plane;
p2Plane.setFromPoints( p2, p2 + maxHeightVec, p2 - plane.normal() );
// Placeholder for triangle vertices per cell
std::vector<caf::HexGridIntersectionTools::ClipVx> hexPlaneCutTriangleVxes;
hexPlaneCutTriangleVxes.reserve( 5 * 3 );
std::vector<int> cellFaceForEachTriangleEdge;
cellFaceForEachTriangleEdge.reserve( 5 * 3 );
cvf::Vec3d cellCorners[8];
size_t cornerIndices[8];
// Mesh data for polyline segment
std::vector<float> polygonVerticesUz = {};
std::vector<cvf::uint> verticesPerPolygon = {};
std::vector<cvf::uint> polygonToCellIndexMap = {};
// Handle triangles per cell
for ( const auto globalCellIdx : columnCellCandidates )
{
if ( ( visibleCells != nullptr ) && ( ( *visibleCells )[globalCellIdx] == 0 ) ) continue;
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
// Perform intersection and clipping of triangles using UTM-coordinates
hexPlaneCutTriangleVxes.clear();
m_hexGrid->cellCornerVertices( globalCellIdx, cellCorners );
m_hexGrid->cellCornerIndices( globalCellIdx, cornerIndices );
// Triangle vertices for polyline segment
caf::HexGridIntersectionTools::planeHexIntersectionMC( plane,
cellCorners,
cornerIndices,
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
// Clip triangles outside the polyline segment using the planes for point p1 and p2
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
std::vector<int> cellFaceForEachClippedTriangleEdge;
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes( hexPlaneCutTriangleVxes,
cellFaceForEachTriangleEdge,
p1Plane,
p2Plane,
&clippedTriangleVxes,
&cellFaceForEachClippedTriangleEdge );
for ( caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes )
{
if ( !clvx.isVxIdsNative ) clvx.derivedVxLevel = 0;
}
// Object to for adding triangle vertices, well vertices and generate polygon vertices
RivEnclosingPolygonGenerator enclosingPolygonGenerator;
// Add clipped triangle vertices to the polygon generator using local coordinates
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
for ( size_t triangleIdx = 0; triangleIdx < clippedTriangleCount; ++triangleIdx )
{
// Get triangle vertices
const size_t vxIdx0 = triangleIdx * 3;
const auto& vx0 = clippedTriangleVxes[vxIdx0 + 0].vx;
const auto& vx1 = clippedTriangleVxes[vxIdx0 + 1].vx;
const auto& vx2 = clippedTriangleVxes[vxIdx0 + 2].vx;
// Convert to local coordinates, where p1 is origin.
// The z-values are global z-values in the uz-coordinates.
const cvf::Vec3d point0( vx0.x() - p1.x(), vx0.y() - p1.y(), vx0.z() );
const cvf::Vec3d point1( vx1.x() - p1.x(), vx1.y() - p1.y(), vx1.z() );
const cvf::Vec3d point2( vx2.x() - p1.x(), vx2.y() - p1.y(), vx2.z() );
// TODO: Ensure counter clockwise order of vertices point0, point1, point2?
// Add triangle to enclosing polygon line handler
enclosingPolygonGenerator.addTriangleVertices( point0, point1, point2 );
}
// Must be a valid polygon to continue
if ( !enclosingPolygonGenerator.isValidPolygon() )
{
continue;
}
// Construct enclosing polygon after adding each triangle
enclosingPolygonGenerator.constructEnclosingPolygon();
const auto& vertices = enclosingPolygonGenerator.getPolygonVertices();
// Construct local uz-coordinates using Pythagoras theorem
for ( const auto& vertex : vertices )
{
// NOTE: Can welding provide drifting of vertex positions?
// TODO: Project (x,y) into plane instead?
//
// Convert to local uz-coordinates, u is the distance along the normalized U-axis
const auto u = std::sqrt( vertex.x() * vertex.x() + vertex.y() * vertex.y() );
const auto z = vertex.z();
polygonVerticesUz.push_back( u );
polygonVerticesUz.push_back( z );
// Keep old code for debugging purposes
calculatedPolygonVertices.push_back( cvf::Vec3f( vertex + p1 ) );
}
verticesPerPolygon.push_back( static_cast<cvf::uint>( vertices.size() ) );
polygonToCellIndexMap.push_back( static_cast<cvf::uint>( globalCellIdx ) );
// Keep old code for debugging purposes
m_verticesPerPolygon.push_back( vertices.size() ); // TODO: Remove when not needed for debug
m_polygonToCellIdxMap.push_back( globalCellIdx ); // TODO: Remove when not needed for debug
}
// Create polygon indices array
const auto numVertices = static_cast<size_t>( polygonVerticesUz.size() / 2 );
std::vector<cvf::uint> polygonIndices( numVertices );
std::iota( polygonIndices.begin(), polygonIndices.end(), 0 );
// Construct polyline segment mesh data
PolylineSegmentMeshData polylineSegmentData;
polylineSegmentData.startUtmXY = cvf::Vec2d( p1.x(), p1.y() );
polylineSegmentData.endUtmXY = cvf::Vec2d( p2.x(), p2.y() );
polylineSegmentData.vertexArrayUZ = polygonVerticesUz;
polylineSegmentData.verticesPerPolygon = verticesPerPolygon;
polylineSegmentData.polygonIndices = polygonIndices;
polylineSegmentData.polygonToCellIndexMap = polygonToCellIndexMap;
// Add polyline segment mesh data to list
m_polylineSegmentsMeshData.push_back( polylineSegmentData );
// Set next polyline segment start index
pointIdx = nextPointIdx;
break;
}
}
m_polygonVertices->assign( calculatedPolygonVertices ); // TODO: Remove when not needed for debug
// Start and end point of polyline segment in UTM-coordinates
const cvf::Vec3d p1 = m_polylineUtm[pointIdx];
const cvf::Vec3d p2 = m_polylineUtm[nextPointIdx];
// Get cell candidates for the polyline segment (subset of global cells)
std::vector<size_t> columnCellCandidates =
createPolylineSegmentCellCandidates( *m_hexGrid, p1, p2, maxHeightVec, topDepth, bottomDepth );
// Plane for the polyline segment
cvf::Plane plane;
plane.setFromPoints( p1, p2, p2 + maxHeightVec );
// Planes parallel to z-axis for p1 and p2, to prevent triangles outside the polyline segment
cvf::Plane p1Plane;
p1Plane.setFromPoints( p1, p1 + maxHeightVec, p1 + plane.normal() );
cvf::Plane p2Plane;
p2Plane.setFromPoints( p2, p2 + maxHeightVec, p2 - plane.normal() );
// Placeholder for triangle vertices per cell
std::vector<caf::HexGridIntersectionTools::ClipVx> hexPlaneCutTriangleVxes;
hexPlaneCutTriangleVxes.reserve( 5 * 3 );
std::vector<int> cellFaceForEachTriangleEdge;
cellFaceForEachTriangleEdge.reserve( 5 * 3 );
cvf::Vec3d cellCorners[8];
size_t cornerIndices[8];
// Polyline segment data
std::vector<cvf::uint> polygonIndices = {};
std::vector<cvf::uint> verticesPerPolygon = {};
std::vector<cvf::uint> polygonToCellIndexMap = {};
// Welder for segment vertices
// - Number of buckets is size of columnCellCandidates divided by 8 to avoid too many buckets (Random selected value).
// Usage of columnCellCandidates is to get a dynamic number of buckets usable for respective segment.
const cvf::uint numWelderBuckets = static_cast<cvf::uint>( columnCellCandidates.size() / size_t( 8 ) );
cvf::VertexWelder segmentVertexWelder;
segmentVertexWelder.initialize( weldingDistance, weldingCellSize, numWelderBuckets );
// Intersection per grid cell - transform from set of triangles to polygon for cell
for ( const auto globalCellIdx : columnCellCandidates )
{
if ( ( visibleCells != nullptr ) && ( ( *visibleCells )[globalCellIdx] == 0 ) ) continue;
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
// Perform intersection and clipping of triangles using UTM-coordinates
hexPlaneCutTriangleVxes.clear();
m_hexGrid->cellCornerVertices( globalCellIdx, cellCorners );
m_hexGrid->cellCornerIndices( globalCellIdx, cornerIndices );
// Triangle vertices for polyline segment
caf::HexGridIntersectionTools::planeHexIntersectionMC( plane,
cellCorners,
cornerIndices,
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
// Clip triangles outside the polyline segment using the planes for point p1 and p2
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
std::vector<int> cellFaceForEachClippedTriangleEdge;
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes( hexPlaneCutTriangleVxes,
cellFaceForEachTriangleEdge,
p1Plane,
p2Plane,
&clippedTriangleVxes,
&cellFaceForEachClippedTriangleEdge );
for ( caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes )
{
if ( !clvx.isVxIdsNative ) clvx.derivedVxLevel = 0;
}
// Object to for adding triangle vertices, well vertices and generate polygon vertices
RivEnclosingPolygonGenerator enclosingPolygonGenerator;
// Add clipped triangle vertices to the polygon generator using local coordinates
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
for ( size_t triangleIdx = 0; triangleIdx < clippedTriangleCount; ++triangleIdx )
{
// Get triangle vertices
const size_t tVxIdx0 = triangleIdx * 3;
const auto& tVx0 = clippedTriangleVxes[tVxIdx0 + 0].vx;
const auto& tVx1 = clippedTriangleVxes[tVxIdx0 + 1].vx;
const auto& tVx2 = clippedTriangleVxes[tVxIdx0 + 2].vx;
// TODO: Ensure counter clockwise order of vertices point0, point1, point2?
// Convert to local coordinates, where p1 is origin.
// The z-values are global z-values in the uz-coordinates.
const cvf::Vec3f vx0( tVx0.x() - p1.x(), tVx0.y() - p1.y(), tVx0.z() );
const cvf::Vec3f vx1( tVx1.x() - p1.x(), tVx1.y() - p1.y(), tVx1.z() );
const cvf::Vec3f vx2( tVx2.x() - p1.x(), tVx2.y() - p1.y(), tVx2.z() );
// Weld vertices and get vertex index
bool isWelded = false;
const auto& wVxIdx0 = segmentVertexWelder.weldVertex( vx0, &isWelded );
const auto& wVxIdx1 = segmentVertexWelder.weldVertex( vx1, &isWelded );
const auto& wVxIdx2 = segmentVertexWelder.weldVertex( vx2, &isWelded );
// Add triangle to enclosing polygon line handler
enclosingPolygonGenerator.addTriangleVertexIndices( wVxIdx0, wVxIdx1, wVxIdx2 );
}
// Must be a valid polygon to continue
if ( !enclosingPolygonGenerator.isValidPolygon() )
{
continue;
}
// Construct enclosing polygon after adding all triangles for cell
enclosingPolygonGenerator.constructEnclosingPolygon();
// Add vertex index to polygon indices
const auto& vertexIndices = enclosingPolygonGenerator.getPolygonVertexIndices();
polygonIndices.insert( polygonIndices.end(), vertexIndices.begin(), vertexIndices.end() );
verticesPerPolygon.push_back( static_cast<cvf::uint>( vertexIndices.size() ) );
polygonToCellIndexMap.push_back( static_cast<cvf::uint>( globalCellIdx ) );
}
// Build vertex array for polyline segment
std::vector<float> polygonVerticesUz;
for ( cvf::uint i = 0; i < segmentVertexWelder.vertexCount(); i++ )
{
const auto& vertex = segmentVertexWelder.vertex( i );
// NOTE: Can welding provide drifting of vertex positions?
// TODO: Project (x,y) into plane instead?
//
// auto projectedVertex = plane.projectPoint( vertex );
// Convert to local uz-coordinates, u is the distance along the normalized U-axis.
// Construct local uz-coordinates using Pythagoras theorem
const auto u = std::sqrt( vertex.x() * vertex.x() + vertex.y() * vertex.y() );
const auto z = vertex.z();
polygonVerticesUz.push_back( u );
polygonVerticesUz.push_back( z );
}
// Construct polyline segment mesh data
PolylineSegmentMeshData polylineSegmentData;
polylineSegmentData.startUtmXY = cvf::Vec2d( p1.x(), p1.y() );
polylineSegmentData.endUtmXY = cvf::Vec2d( p2.x(), p2.y() );
polylineSegmentData.vertexArrayUZ = polygonVerticesUz;
polylineSegmentData.verticesPerPolygon = verticesPerPolygon;
polylineSegmentData.polygonIndices = polygonIndices;
polylineSegmentData.polygonToCellIndexMap = polygonToCellIndexMap;
// Add polyline segment mesh data to list
m_polylineSegmentsMeshData.push_back( polylineSegmentData );
// Set next polyline segment start index
pointIdx = nextPointIdx;
}
}
//--------------------------------------------------------------------------------------------------

12
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.h
View File

@ -26,8 +26,6 @@
#include <vector>
#include <QString>
class RigMainGrid;
class RigActiveCellInfo;
class RigResultAccessor;
@ -60,11 +58,6 @@ public:
void generateIntersectionGeometry( cvf::UByteArray* visibleCells );
bool isAnyGeometryPresent() const;
// TODO: Remove after testing?
const cvf::Vec3fArray* polygonVxes() const;
const std::vector<size_t>& vertiesPerPolygon() const;
const std::vector<size_t>& polygonToCellIndex() const;
const std::vector<PolylineSegmentMeshData>& polylineSegmentsMeshData() const;
private:
@ -84,9 +77,4 @@ private:
// Output
std::vector<PolylineSegmentMeshData> m_polylineSegmentsMeshData;
// TMP Output arrays for debug
std::vector<size_t> m_polygonToCellIdxMap;
cvf::ref<cvf::Vec3fArray> m_polygonVertices;
std::vector<size_t> m_verticesPerPolygon;
};

17
GrpcInterface/GrpcProtos/GridGeometryExtraction.proto
View File

@ -68,27 +68,18 @@ message CutAlongPolylineRequest
message FenceMeshSection
{
// U-axis defined by the unit length vector from start to end, Z is global Z
repeated float vertexArrayUZ = 1; // U coordinate is length along U-axis
repeated fixed32 polyIndicesArr = 2; // Note: Redundant if no shared nodes?
repeated float vertexArrayUZ = 1; // U coordinate is length along U-axis. Unique vertex coordinates (u,z). Stored [u0, z0, u1, z1, ...]
repeated fixed32 polyIndicesArr = 2; // Polygon vertex indices. Index of vertex (u,z) in the unique vertex coordinate array.
repeated fixed32 verticesPerPolygonArr = 3; // Number of vertices per polygon, numPolygons long
repeated fixed32 sourceCellIndicesArr = 4; // The originating cell index per polygon, numPolygons long
Vec2d startUtmXY = 6;
Vec2d endUtmXY = 7;
}
message PolylineTestResponse
{
// Test response returning all vertices without "local" section coordinates
repeated float polygonVertexArray = 1;
repeated fixed32 verticesPerPolygonArr = 2; // Number of vertices per polygon, numPolygons long
repeated fixed32 sourceCellIndicesArr = 3; // The originating cell index per polygon, numPolygons long
}
message CutAlongPolylineResponse
{
repeated FenceMeshSection fenceMeshSections = 1;
PolylineTestResponse polylineTestResponse = 2;
TimeElapsedInfo timeElapsedInfo = 3;
Vec3i gridDimensions = 4;
TimeElapsedInfo timeElapsedInfo = 2;
Vec3i gridDimensions = 3;
}

45
GrpcInterface/Python/rips/WebvizPythonExamples/grid_geometry_extraction_cut_along_polyline.py
View File

@ -44,7 +44,7 @@ norne_case_single_segment_poly_line_gap_utm_xy = [460877, 7.3236e06, 459279, 7.3
# fence_poly_line_utm_xy = norne_case_fence_poly_line_utm_xy
fence_poly_line_utm_xy = norne_case_single_segment_poly_line_gap_utm_xy
fence_poly_line_utm_xy = norne_case_fence_poly_line_utm_xy
cut_along_polyline_request = GridGeometryExtraction__pb2.CutAlongPolylineRequest(
gridFilename=grid_file_name,
@ -68,10 +68,20 @@ section_polygon_edges_3d = []
# Use first segment start (x,y) as origin
x_origin = fence_mesh_sections[0].startUtmXY.x if len(fence_mesh_sections) > 0 else 0
y_origin = fence_mesh_sections[0].startUtmXY.y if len(fence_mesh_sections) > 0 else 0
section_number = 1
for section in fence_mesh_sections:
# Continue to next section
polygon_vertex_array_uz = section.vertexArrayUZ
vertices_per_polygon = section.verticesPerPolygonArr
polygon_indices_array = section.polyIndicesArr
num_vertices = sum(vertices_per_polygon)
print(f"**** Section number: {section_number} ****")
print(f"Number of vertices in vertex uz array: {len(polygon_vertex_array_uz)/2}")
print(f"Number of polygon vertices: {len(polygon_indices_array)}")
print(f"Number of vertices: {num_vertices}")
section_number += 1
# Get start and end coordinates (local coordinates)
start_x = section.startUtmXY.x - x_origin
@ -94,18 +104,31 @@ for section in fence_mesh_sections:
x_array = []
y_array = []
z_array = []
for i in range(0, len(polygon_vertex_array_uz), vertex_step):
u = polygon_vertex_array_uz[i]
z = polygon_vertex_array_uz[i + 1]
# Calculate x, y from u and directional vector,
# where u is the length along the direction vector
x = start_x + u * direction_vector_norm[0]
y = start_y + u * direction_vector_norm[1]
vertex_offset = 0
for _, value in enumerate(vertices_per_polygon, 0):
num_polygon_vertices = value
polygon_indices = polygon_indices_array[
vertex_offset : vertex_offset + num_polygon_vertices
]
x_array.append(x)
y_array.append(y)
z_array.append(z)
for vertex_idx in polygon_indices:
idx = vertex_idx * vertex_step
# Extract u, z values for the polygon
u = polygon_vertex_array_uz[idx]
z = polygon_vertex_array_uz[idx + 1]
# Calculate x, y from u and directional vector,
# where u is the length along the direction vector
x = start_x + u * direction_vector_norm[0]
y = start_y + u * direction_vector_norm[1]
x_array.append(x)
y_array.append(y)
z_array.append(z)
vertex_offset = vertex_offset + num_polygon_vertices
i = []
j = []

38
GrpcInterface/RiaGrpcGridGeometryExtractionService.cpp
View File

@ -366,44 +366,6 @@ grpc::Status RiaGrpcGridGeometryExtractionService::CutAlongPolyline( grpc::Serve
m_elapsedTimeInfo.elapsedTimePerEventMs["FillResponse"] =
static_cast<std::uint32_t>( fillResponseTimeCount.elapsedMsCount() );
// Add temporary test response
{
auto fillTestResponseTimeCount = ElapsedTimeCount();
rips::PolylineTestResponse* polylineTestResponse = new rips::PolylineTestResponse;
// Polygon vertices
const auto& polygonVertices = polylineIntersectionGenerator.polygonVxes();
if ( polygonVertices->size() == 0 )
{
return grpc::Status( grpc::StatusCode::NOT_FOUND, "No polygon vertices found for polyline" );
}
for ( size_t i = 0; i < polygonVertices->size(); ++i )
{
const auto& vertex = polygonVertices->get( i );
polylineTestResponse->add_polygonvertexarray( vertex.x() );
polylineTestResponse->add_polygonvertexarray( vertex.y() );
polylineTestResponse->add_polygonvertexarray( vertex.z() );
}
// Vertices per polygon
const auto& verticesPerPolygon = polylineIntersectionGenerator.vertiesPerPolygon();
for ( const auto& elm : verticesPerPolygon )
{
polylineTestResponse->add_verticesperpolygonarr( static_cast<google::protobuf::uint32>( elm ) );
}
// Polygon to cell indices
const auto& polygonCellIndices = polylineIntersectionGenerator.polygonToCellIndex();
for ( const auto& elm : polygonCellIndices )
{
polylineTestResponse->add_sourcecellindicesarr( static_cast<google::protobuf::uint32>( elm ) );
}
response->set_allocated_polylinetestresponse( polylineTestResponse );
m_elapsedTimeInfo.elapsedTimePerEventMs["FillResponse"] =
static_cast<std::uint32_t>( fillTestResponseTimeCount.elapsedMsCount() );
}
// Clear existing view
tearDownExistingViewsInEclipseCase();