#3876 Create Counter Clockwise closed polygons out of the contour line segments

ApplicationCode/ModelVisualization/RivContourMapProjectionPartMgr.cpp

@@ -174,29 +174,31 @@ cvf::ref<cvf::DrawableGeo> RivContourMapProjectionPartMgr::createProjectionMapDr
 //--------------------------------------------------------------------------------------------------
 std::vector<cvf::ref<cvf::DrawableGeo>> RivContourMapProjectionPartMgr::createContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform) const
 {
-    RimContourMapProjection::ContourPolygons contourPolygons = m_contourMapProjection->generateContourPolygons(displayCoordTransform);
+    RimContourMapProjection::ClosedContourPolygons contourPolygons = m_contourMapProjection->generateContourPolygons(displayCoordTransform);
 
     std::vector<cvf::ref<cvf::DrawableGeo>> contourDrawables;
-    contourDrawables.reserve(contourPolygons.size());
     for (size_t i = 0; i < contourPolygons.size(); ++i)
     {
-        cvf::ref<cvf::Vec3fArray> vertexArray = contourPolygons[i];
-        std::vector<cvf::uint> indices;
-        indices.reserve(contourPolygons[i]->size());
-        for (cvf::uint j = 0; j < contourPolygons[i]->size(); ++j)
+        for (size_t j = 0; j < contourPolygons[i].size(); ++j)
         {
-            indices.push_back(j);
+            cvf::ref<cvf::Vec3fArray> vertexArray = contourPolygons[i][j];
+            std::vector<cvf::uint> indices;
+            indices.reserve(contourPolygons[i][j]->size());
+            for (cvf::uint k = 0; k < contourPolygons[i][j]->size(); ++k)
+            {
+                indices.push_back(k);
+            }
+
+            cvf::ref<cvf::PrimitiveSetIndexedUInt> indexedUInt = new cvf::PrimitiveSetIndexedUInt(cvf::PrimitiveType::PT_LINES);
+            cvf::ref<cvf::UIntArray>               indexArray = new cvf::UIntArray(indices);
+            indexedUInt->setIndices(indexArray.p());
+
+            cvf::ref<cvf::DrawableGeo> geo = new cvf::DrawableGeo;
+
+            geo->addPrimitiveSet(indexedUInt.p());
+            geo->setVertexArray(vertexArray.p());
+            contourDrawables.push_back(geo);
         }
-
-        cvf::ref<cvf::PrimitiveSetIndexedUInt> indexedUInt = new cvf::PrimitiveSetIndexedUInt(cvf::PrimitiveType::PT_LINES);
-        cvf::ref<cvf::UIntArray>               indexArray = new cvf::UIntArray(indices);
-        indexedUInt->setIndices(indexArray.p());
-
-        cvf::ref<cvf::DrawableGeo> geo = new cvf::DrawableGeo;
-
-        geo->addPrimitiveSet(indexedUInt.p());
-        geo->setVertexArray(vertexArray.p());
-        contourDrawables.push_back(geo);
     }
     return contourDrawables;
 }

ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp

@@ -128,9 +128,9 @@ void RimContourMapProjection::generateVertices(cvf::Vec3fArray* vertices, const
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-RimContourMapProjection::ContourPolygons RimContourMapProjection::generateContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform)
+RimContourMapProjection::ClosedContourPolygons RimContourMapProjection::generateContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform)
 {    
-    std::vector<cvf::ref<cvf::Vec3fArray>> contourPolygons;
+    ClosedContourPolygons contourPolygons;
     if (minValue() != std::numeric_limits<double>::infinity() &&
         maxValue() != -std::numeric_limits<double>::infinity() &&
         std::fabs(maxValue() - minValue()) > 1.0e-8)
@@ -149,22 +149,22 @@ RimContourMapProjection::ContourPolygons RimContourMapProjection::generateContou
                     contourLevels[0] += (contourLevels[1] - contourLevels[0]) * 0.1;
                     contourLevels[nContourLevels - 1] -= (contourLevels[nContourLevels - 1] - contourLevels[nContourLevels - 2]) * 0.1;
                 }
-                std::vector<std::vector<cvf::Vec2d>> contourLines;
-                caf::ContourLines::create(m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels, &contourLines);
+                std::vector<caf::ContourLines::ClosedPolygons> closedContourLines =
+                    caf::ContourLines::create(m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels);
 
-                contourPolygons.reserve(contourLines.size());
-                for (size_t i = 0; i < contourLines.size(); ++i)
+                contourPolygons.resize(closedContourLines.size());
+                for (size_t i = 0; i < closedContourLines.size(); ++i)
                 {
-                    if (!contourLines[i].empty())
+                    for (size_t j = 0; j < closedContourLines[i].size(); ++j)
                     {
-                        cvf::ref<cvf::Vec3fArray> contourPolygon = new cvf::Vec3fArray(contourLines[i].size());
-                        for (size_t j = 0; j < contourLines[i].size(); ++j)
+                        cvf::ref<cvf::Vec3fArray> contourPolygon = new cvf::Vec3fArray(closedContourLines[i][j].size());
+                        for (size_t k = 0; k < closedContourLines[i][j].size(); ++k)
                         {
-                            cvf::Vec3d contourPoint3d = cvf::Vec3d(contourLines[i][j], m_fullBoundingBox.min().z());
+                            cvf::Vec3d contourPoint3d = cvf::Vec3d(closedContourLines[i][j][k], m_fullBoundingBox.min().z());
                             cvf::Vec3d displayPoint3d = displayCoordTransform->transformToDisplayCoord(contourPoint3d);
-                            (*contourPolygon)[j] = cvf::Vec3f(displayPoint3d);
+                            (*contourPolygon)[k] = cvf::Vec3f(displayPoint3d);
                         }
-                        contourPolygons.push_back(contourPolygon);
+                        contourPolygons[i].push_back(contourPolygon);
                     }
                 }
             }
@@ -1314,6 +1314,7 @@ void RimContourMapProjection::updateGridInformation()
     m_mainGrid        = eclipseCase()->eclipseCaseData()->mainGrid();
     m_sampleSpacing   = m_relativeSampleSpacing * m_mainGrid->characteristicIJCellSize();
     m_fullBoundingBox = eclipseCase()->activeCellsBoundingBox();
+    m_fullBoundingBox.expand(m_sampleSpacing * 2.0);
     m_mapSize         = calculateMapSize();
 
     // Re-jig max point to be an exact multiple of cell size

ApplicationCode/ProjectDataModel/RimContourMapProjection.h

@@ -60,13 +60,13 @@ public:
         RESULTS_HC_COLUMN
     };
     typedef caf::AppEnum<ResultAggregationEnum> ResultAggregation;
-    typedef std::vector<cvf::ref<cvf::Vec3fArray>> ContourPolygons;
+    typedef std::vector<std::vector<cvf::ref<cvf::Vec3fArray>>> ClosedContourPolygons;
 
     RimContourMapProjection();
     ~RimContourMapProjection() override;
 
     void                        generateVertices(cvf::Vec3fArray* vertices, const caf::DisplayCoordTransform* displayCoordTransform);    
-    ContourPolygons             generateContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform);
+    ClosedContourPolygons       generateContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform);
     cvf::ref<cvf::Vec3fArray>   generatePickPointPolygon(const caf::DisplayCoordTransform* displayCoordTransform);
     void                        generateResults();
 

Fwk/AppFwk/CommonCode/cafContourLines.cpp

@@ -22,7 +22,9 @@
 /////////////////////////////////////////////////////////////////////////////////
 
 #include "cafContourLines.h"
+
 #include <algorithm>
+#include <list>
 
 const int caf::ContourLines::s_castab[3][3][3] =
 {
@@ -206,6 +208,103 @@ void caf::ContourLines::create(const std::vector<double>& dataXY, const std::vec
     } /* j */
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<caf::ContourLines::ClosedPolygons> caf::ContourLines::create(const std::vector<double>& dataXY,
+    const std::vector<double>& xPositions,
+    const std::vector<double>& yPositions,
+    const std::vector<double>& contourLevels)
+{
+    const double eps = 1.0e-4;
+    std::vector<std::vector<cvf::Vec2d>> contourLineSegments;
+    caf::ContourLines::create(dataXY, xPositions, yPositions, contourLevels, &contourLineSegments);
+
+    std::vector<ClosedPolygons> closedPolygonsPerLevel(contourLevels.size());
+
+    for (size_t i = 0; i < contourLevels.size(); ++i)
+    {
+        size_t nPoints = contourLineSegments[i].size();
+        size_t nSegments = nPoints / 2;
+        if (nSegments >= 3u) // Need at least three segments for a closed polygon
+        {
+            std::list<std::pair<cvf::Vec2d, cvf::Vec2d>> unorderedSegments;
+            for (size_t j = 0; j < contourLineSegments[i].size(); j += 2)
+            {
+                unorderedSegments.push_back(std::make_pair(contourLineSegments[i][j], contourLineSegments[i][j + 1]));
+            }
+
+            std::deque<cvf::Vec2d> closedPolygonDeque;
+            while (!unorderedSegments.empty())
+            {                
+                bool expandedPolygon = false;
+                for (auto listIt = unorderedSegments.begin(); listIt != unorderedSegments.end(); ++listIt)
+                {
+                    if (closedPolygonDeque.empty() || listIt->first == closedPolygonDeque.back())
+                    {
+                        closedPolygonDeque.push_back(listIt->first);
+                        closedPolygonDeque.push_back(listIt->second);
+                        unorderedSegments.erase(listIt);
+                        expandedPolygon = true;
+                        break;
+                    }
+                    else if (listIt->second == closedPolygonDeque.back())
+                    {
+                        closedPolygonDeque.push_back(listIt->second);
+                        closedPolygonDeque.push_back(listIt->first);
+                        unorderedSegments.erase(listIt);
+                        expandedPolygon = true;
+                        break;
+                    }
+                    else if (listIt->first == closedPolygonDeque.front())
+                    {
+                        closedPolygonDeque.push_front(listIt->first);
+                        closedPolygonDeque.push_front(listIt->second);
+                        unorderedSegments.erase(listIt);
+                        expandedPolygon = true;
+                        break;
+                    }
+                    else if (listIt->second == closedPolygonDeque.front())
+                    {
+                        closedPolygonDeque.push_front(listIt->second);
+                        closedPolygonDeque.push_front(listIt->first);
+                        unorderedSegments.erase(listIt);
+                        expandedPolygon = true;
+                        break;
+                    }
+                }
+                if (!expandedPolygon || unorderedSegments.empty())
+                {
+                    if (closedPolygonDeque.back() != closedPolygonDeque.front())
+                    {
+                        closedPolygonDeque.push_back(closedPolygonDeque.back());
+                        closedPolygonDeque.push_back(closedPolygonDeque.front());
+                    }
+
+                    // Make sure it is counter clockwise
+                    double sum = 0.0;
+                    for (size_t j = 0; j < closedPolygonDeque.size() - 1; ++j)
+                    {
+                        sum += (closedPolygonDeque[j + 1].x() - closedPolygonDeque[j].x()) *
+                               (closedPolygonDeque[j + 1].y() + closedPolygonDeque[j].y());
+                    }
+                    if (sum < 0.0)
+                    {
+                        closedPolygonsPerLevel[i].emplace_back(closedPolygonDeque.rbegin(), closedPolygonDeque.rend());
+                    }
+                    else
+                    {
+                        closedPolygonsPerLevel[i].emplace_back(closedPolygonDeque.begin(), closedPolygonDeque.end());
+                    }
+
+                    closedPolygonDeque.clear();
+                }
+            }
+        }
+    }
+    return closedPolygonsPerLevel;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

Fwk/AppFwk/CommonCode/cafContourLines.h

@@ -25,6 +25,7 @@
 
 #include "cvfBase.h"
 #include "cvfVector2.h"
+#include <deque>
 #include <vector>
 
 namespace caf
@@ -32,12 +33,20 @@ namespace caf
 class ContourLines
 {
 public:
+    typedef std::vector<cvf::Vec2d> ClosedPolygon;
+    typedef std::vector<ClosedPolygon> ClosedPolygons;
+
+    static std::vector<ClosedPolygons> create(const std::vector<double>&            dataXY,
+                                              const std::vector<double>&            xPositions,
+                                              const std::vector<double>&            yPositions,
+                                              const std::vector<double>&            contourLevels);
+    
+private:
     static void create(const std::vector<double>& dataXY,
                        const std::vector<double>& xPositions,
                        const std::vector<double>& yPositions,
                        const std::vector<double>& contourLevels,
                        std::vector<std::vector<cvf::Vec2d>>* polygons);
-private:
     static double contourRange(const std::vector<double>& contourLevels);
     static double invalidValue(const std::vector<double>& contourLevels);
     static double saneValue(int index, const std::vector<double>& dataXY, const std::vector<double>& contourLevels);