diff --git a/ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp b/ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp index b2cee4d60d..a7e86ca762 100644 --- a/ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp +++ b/ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp @@ -44,6 +44,8 @@ #include #include +#include + namespace caf { template<> @@ -732,13 +734,12 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() int myThread = omp_get_thread_num(); threadTriangles[myThread].resize(std::max((size_t) 1, m_contourPolygons.size())); - std::set excludedFaceIndices; - #pragma omp for schedule(dynamic) for (int64_t i = 0; i < (int64_t)faceList->size(); i += 3) { std::vector triangle(3); std::vector triangleWithValues(3); + bool anyValidVertex = false; for (size_t n = 0; n < 3; ++n) { uint vn = (*faceList)[i + n]; @@ -746,6 +747,15 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() : std::numeric_limits::infinity(); triangle[n] = vertices[vn]; triangleWithValues[n] = cvf::Vec4d(vertices[vn], value); + if (value != std::numeric_limits::infinity()) + { + anyValidVertex = true; + } + } + + if (!anyValidVertex) + { + continue; } if (m_contourPolygons.empty()) @@ -755,7 +765,8 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() continue; } - for (size_t c = 0; c < m_contourPolygons.size() && excludedFaceIndices.count(i) == 0u; ++c) + bool outsideOuterLimit = false; + for (size_t c = 0; c < m_contourPolygons.size() && !outsideOuterLimit; ++c) { std::vector> intersectPolygons; for (size_t j = 0; j < m_contourPolygons[c].size(); ++j) @@ -778,7 +789,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() if (intersectPolygons.empty()) { - excludedFaceIndices.insert(i); + outsideOuterLimit = true; continue; } @@ -1001,6 +1012,14 @@ void RimContourMapProjection::generateContourPolygons() } } + double simplifyEpsilon = m_smoothContourLines() ? 5.0e-2 * m_sampleSpacing : 1.0e-3 * m_sampleSpacing; + for (size_t i = 0; i < contourPolygons.size(); ++i) + { + for (ContourPolygon& polygon : contourPolygons[i]) + { + RigCellGeometryTools::simplifyPolygon(&polygon.vertices, simplifyEpsilon); + } + } m_contourLevelCumulativeAreas.resize(contourPolygons.size(), 0.0); for (int64_t i = (int64_t) contourPolygons.size() - 1; i >= 0; --i) { @@ -1045,15 +1064,16 @@ void RimContourMapProjection::smoothContourPolygons(ContourPolygons* conto } // Only expand. cvf::Vec3d modifiedVertex = 0.5 * (v + 0.5 * (vm1 + vp1)); - cvf::Vec3d delta = (modifiedVertex - v).getNormalized(); + 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 half the amount - modifiedVertex = v + 0.5 * delta; + // 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()); @@ -1572,7 +1592,7 @@ void RimContourMapProjection::defineEditorAttribute(const caf::PdmFieldHandle* f caf::PdmUiDoubleSliderEditorAttribute* myAttr = dynamic_cast(attribute); if (myAttr) { - myAttr->m_minimum = 0.25; + myAttr->m_minimum = 0.2; myAttr->m_maximum = 2.0; } } diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp index 0037ba924e..d7d3d691f0 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp @@ -194,6 +194,50 @@ void RigCellGeometryTools::createPolygonFromLineSegments(std::list* vertices, double epsilon) +{ + CVF_ASSERT(vertices); + if (vertices->size() < 3) return; + + std::pair maxDistPoint(0u, 0.0); + + for (size_t i = 1; i < vertices->size() - 1; ++i) + { + cvf::Vec3d v = vertices->at(i); + double u; + cvf::Vec3d v_proj = cvf::GeometryTools::projectPointOnLine(vertices->front(), vertices->back(), v, &u); + double distance = (v_proj - v).length(); + if (distance > maxDistPoint.second) + { + maxDistPoint = std::make_pair(i, distance); + } + } + + if (maxDistPoint.second > epsilon) + { + std::vector newVertices1(vertices->begin(), vertices->begin() + maxDistPoint.first + 1); + std::vector newVertices2(vertices->begin() + maxDistPoint.first, vertices->end()); + + // Recurse + simplifyPolygon(&newVertices1, epsilon); + simplifyPolygon(&newVertices2, epsilon); + + std::vector newVertices(newVertices1.begin(), newVertices1.end() - 1); + newVertices.insert(newVertices.end(), newVertices2.begin(), newVertices2.end()); + *vertices = newVertices; + } + else + { + std::vector newVertices = {vertices->front(), vertices->back()}; + *vertices = newVertices; + } +} + //================================================================================================== /// //================================================================================================== diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h index fcc5b264fb..49d2f40b8b 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h @@ -39,6 +39,7 @@ public: static void createPolygonFromLineSegments(std::list>& intersectionLineSegments, std::vector>& polygons, double tolerance = 1.0e-4); + static void simplifyPolygon(std::vector* vertices, double epsilon); static void findCellLocalXYZ(const std::array& hexCorners, cvf::Vec3d& localXdirection,