#3970 Contour Maps: improve edge look

* Remove excess tiny triangles around edges
* Remove labels from contour lines that overlap inner contour lines

ApplicationCode/ModelVisualization/RivContourMapProjectionPartMgr.cpp

@@ -16,6 +16,7 @@
 #include "cvfCamera.h"
 #include "cvfDrawableText.h"
 #include "cvfGeometryBuilderFaceList.h"
+#include "cvfGeometryTools.h"
 #include "cvfGeometryUtils.h"
 #include "cvfMeshEdgeExtractor.h"
 #include "cvfPart.h"
@@ -25,6 +26,7 @@
 
 #include <cmath>
 
+#include <QDebug>
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -369,6 +371,7 @@ std::vector<cvf::ref<cvf::Drawable>>
     std::vector<cvf::ref<cvf::Drawable>>              labelDrawables;
     labelBBoxes->clear();
     labelBBoxes->resize(m_contourLinePolygons.size());
+    const RimContourMapProjection::ContourPolygons* previousLevel = nullptr;
     for (int64_t i = (int64_t)m_contourLinePolygons.size() - 1; i > 0; --i)
     {
         cvf::Color3f                backgroundColor(mapper->mapToColor(tickValues[i]));
@@ -388,9 +391,19 @@ std::vector<cvf::ref<cvf::Drawable>>
             {
                 size_t     nVertex = (nVertices * l) / nLabels;
                 size_t     nextVertex = (nVertex + 1) % nVertices;
-                cvf::Vec3d globalVertex1 = m_contourLinePolygons[i][j].vertices[nVertex] + m_contourMapProjection->origin3d();
+
-                cvf::Vec3d globalVertex2 =
+                const cvf::Vec3d& localVertex1 = m_contourLinePolygons[i][j].vertices[nVertex];
-                    m_contourLinePolygons[i][j].vertices[nextVertex] + m_contourMapProjection->origin3d();
+                const cvf::Vec3d& localVertex2 = m_contourLinePolygons[i][j].vertices[nextVertex];
+
+                cvf::Vec3d lineCenter = (localVertex1 + localVertex2) * 0.5;
+                if (false && previousLevel && lineOverlapsWithPreviousContourLevel(lineCenter, previousLevel))
+                {
+                    continue;
+                }
+
+                cvf::Vec3d globalVertex1 = localVertex1 + m_contourMapProjection->origin3d();
+                cvf::Vec3d globalVertex2 = localVertex2 + m_contourMapProjection->origin3d();
+
                 cvf::Vec3d globalVertex = 0.5 * (globalVertex1 + globalVertex2);
 
                 cvf::Vec3d segment = globalVertex2 - globalVertex1;
@@ -451,6 +464,8 @@ std::vector<cvf::ref<cvf::Drawable>>
         {
             labelDrawables.push_back(label);
         }
+
+        previousLevel = &m_contourLinePolygons[i];
     }
     return labelDrawables;
 }
@@ -497,3 +512,44 @@ cvf::ref<cvf::DrawableGeo>
     }
     return geo;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RivContourMapProjectionPartMgr::lineOverlapsWithPreviousContourLevel(
+    const cvf::Vec3d& lineCenter,
+    const RimContourMapProjection::ContourPolygons* previousLevel) const
+{
+    const int64_t jump = 50;
+    CVF_ASSERT(previousLevel);
+    double tolerance = 1.0e-3 * m_contourMapProjection->sampleSpacing();
+    for (const RimContourMapProjection::ContourPolygon& edgePolygon : *previousLevel)
+    {
+        std::pair<int64_t, double> closestIndex(0, std::numeric_limits<double>::infinity());
+        for (int64_t i = 0; i < (int64_t) edgePolygon.vertices.size(); i += jump)
+        {
+            const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
+            const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[(i + 1) % edgePolygon.vertices.size()];
+            double dist1 = cvf::GeometryTools::linePointSquareDist(edgeVertex1, edgeVertex2, lineCenter);
+            if (dist1 < tolerance)
+            {
+                return true;
+            }
+            if (dist1 < closestIndex.second)
+            {
+                closestIndex = std::make_pair(i, dist1);
+            }
+        }
+        for (int64_t i = std::max((int64_t)1, closestIndex.first - jump + 1); i < std::min((int64_t)edgePolygon.vertices.size(), closestIndex.first + jump); ++i)
+        {
+            const cvf::Vec3d& edgeVertex1 = edgePolygon.vertices[i];
+            const cvf::Vec3d& edgeVertex2 = edgePolygon.vertices[(i + 1) % edgePolygon.vertices.size()];
+            double            dist1       = cvf::GeometryTools::linePointSquareDist(edgeVertex1, edgeVertex2, lineCenter);
+            if (dist1 < tolerance)
+            {
+                return true;
+            }
+        }
+    }
+    return false;
+}

ApplicationCode/ModelVisualization/RivContourMapProjectionPartMgr.h

@@ -58,6 +58,9 @@ private:
     std::vector<std::vector<cvf::ref<cvf::Drawable>>>    createContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform, const std::vector<std::vector<cvf::BoundingBox>>& labelBBoxes) const;
     std::vector<cvf::ref<cvf::Drawable>>                 createContourLabels(const cvf::Camera* camera, const caf::DisplayCoordTransform* displayCoordTransform, std::vector<std::vector<cvf::BoundingBox>>* labelBBoxes) const;
     cvf::ref<cvf::DrawableGeo>                           createPickPointVisDrawable(const caf::DisplayCoordTransform* displayCoordTransform) const;
+    bool lineOverlapsWithPreviousContourLevel(const cvf::Vec3d& lineCenter,
+                                                                                              const RimContourMapProjection::ContourPolygons* previousLevel) const;
+
 private:
     caf::PdmPointer<RimContourMapProjection>  m_contourMapProjection;
     caf::PdmPointer<RimGridView>              m_parentContourMap;

ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp

@@ -710,6 +710,9 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
         }
     }
 
+    const double cellArea = m_sampleSpacing * m_sampleSpacing;
+    const double areaThreshold = 1.0e-5 * 0.5 * cellArea;
+
     std::vector<std::vector<std::vector<cvf::Vec3d>>> subtractPolygons;
     if (!m_contourPolygons.empty())
     {
@@ -722,12 +725,12 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
             }
         }
     }
-    std::vector<std::vector<cvf::Vec4d>> threadTriangles(omp_get_max_threads());
+    std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles(omp_get_max_threads());
 
 #pragma omp parallel
     {
         int myThread = omp_get_thread_num();
-        threadTriangles[myThread].reserve(faceList->size() / omp_get_num_threads());
+        threadTriangles[myThread].resize(std::max((size_t) 1, m_contourPolygons.size()));
 
         std::set<int64_t> excludedFaceIndices;
 
@@ -747,8 +750,8 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
 
             if (m_contourPolygons.empty())
             {
-                threadTriangles[myThread].insert(
+                threadTriangles[myThread][0].insert(
-                    threadTriangles[myThread].end(), triangleWithValues.begin(), triangleWithValues.end());
+                    threadTriangles[myThread][0].end(), triangleWithValues.begin(), triangleWithValues.end());
                 continue;
             }
 
@@ -778,7 +781,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
                     excludedFaceIndices.insert(i);
                     continue;
                 }
-
+                
                 std::vector<std::vector<cvf::Vec3d>> clippedPolygons;
 
                 if (!subtractPolygons[c].empty())
@@ -794,6 +797,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
                 {
                     clippedPolygons.swap(intersectPolygons);
                 }
+
                 {
                     std::vector<cvf::Vec4d> clippedTriangles;
                     for (std::vector<cvf::Vec3d>& clippedPolygon : clippedPolygons)
@@ -828,6 +832,10 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
                         }
                         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();
@@ -862,16 +870,37 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
                             }
                         }
                     }
-                    threadTriangles[myThread].insert(
+                    threadTriangles[myThread][c].insert(
-                        threadTriangles[myThread].end(), clippedTriangles.begin(), clippedTriangles.end());
+                        threadTriangles[myThread][c].end(), clippedTriangles.begin(), clippedTriangles.end());
                 }
             }
         }
     }
-    std::vector<cvf::Vec4d> finalTriangles;
+
-    for (size_t i = 0; i < threadTriangles.size(); ++i)
+    std::vector<std::vector<cvf::Vec4d>> trianglesPerLevel(std::max((size_t)1, contourLevels.size()));
+    for (size_t c = 0; c < trianglesPerLevel.size(); ++c)
     {
-        finalTriangles.insert(finalTriangles.end(), threadTriangles[i].begin(), threadTriangles[i].end());
+        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 = sumTriangleAreas(allTrianglesThisLevel);
+        if (triangleAreasThisLevel > 1.0e-3 * m_contourLevelCumulativeAreas[c])
+        {
+            trianglesPerLevel[c] = allTrianglesThisLevel;
+        }
+        else
+        {
+            m_contourPolygons[c].clear();
+        }
+    }
+
+    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;
@@ -929,26 +958,38 @@ void RimContourMapProjection::generateContourPolygons()
                     contourLevels.front() *= 0.5;
                 }
 
-                std::vector<caf::ContourLines::ClosedPolygons> closedContourLines = caf::ContourLines::create(
+                std::vector<caf::ContourLines::ListOfLineSegments> unorderedLineSegmentsPerLevel = caf::ContourLines::create(
-                    m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels, areaTreshold);
+                    m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels);
 
-                contourPolygons.resize(closedContourLines.size());
+                std::vector<ContourPolygons>(unorderedLineSegmentsPerLevel.size()).swap(contourPolygons);
-
+                
-                for (size_t i = 0; i < closedContourLines.size(); ++i)
+                for (size_t i = 0; i < unorderedLineSegmentsPerLevel.size(); ++i)
                 {
-                    for (size_t j = 0; j < closedContourLines[i].size(); ++j)
+                    std::vector<std::vector<cvf::Vec3d>> polygonsForThisLevel;
+                    RigCellGeometryTools::createPolygonFromLineSegments(unorderedLineSegmentsPerLevel[i], polygonsForThisLevel, 1.0e-8);
+                    for (size_t j = 0; j < polygonsForThisLevel.size(); ++j)
                     {
+                        double signedArea = cvf::GeometryTools::signedAreaPlanarPolygon(cvf::Vec3d::Z_AXIS, polygonsForThisLevel[j]);
                         ContourPolygon contourPolygon;
                         contourPolygon.value = contourLevels[i];
-                        contourPolygon.vertices.reserve(closedContourLines[i][j].size() / 2);
+                        if (signedArea < 0.0)
-
-                        for (size_t k = 0; k < closedContourLines[i][j].size(); k += 2)
                         {
-                            cvf::Vec3d contourPoint3d = cvf::Vec3d(closedContourLines[i][j][k], 0.0);
+                            contourPolygon.vertices.insert(contourPolygon.vertices.end(), polygonsForThisLevel[j].rbegin(), polygonsForThisLevel[j].rend());
-                            contourPolygon.vertices.push_back(contourPoint3d);
+                        }
-                            contourPolygon.bbox.add(contourPoint3d);
+                        else
+                        {
+                            contourPolygon.vertices = polygonsForThisLevel[j];
+                        }
+
+                        contourPolygon.area = cvf::GeometryTools::signedAreaPlanarPolygon(cvf::Vec3d::Z_AXIS, contourPolygon.vertices);
+                        if (contourPolygon.area > areaTreshold)
+                        {
+                            for (const cvf::Vec3d& vertex : contourPolygon.vertices)
+                            {
+                                contourPolygon.bbox.add(vertex);
+                            }
+                            contourPolygons[i].push_back(contourPolygon);
                         }
-                        contourPolygons[i].push_back(contourPolygon);
                     }
                 }
                 for (size_t i = 0; i < contourPolygons.size(); ++i)
@@ -959,6 +1000,13 @@ void RimContourMapProjection::generateContourPolygons()
                         smoothContourPolygons(&contourPolygons[i], clipBy, true);
                     }
                 }
+
+                m_contourLevelCumulativeAreas.resize(contourPolygons.size(), 0.0);
+                for (int64_t i = (int64_t) contourPolygons.size() - 1; i >= 0; --i)
+                {
+                    double levelOuterArea = sumPolygonArea(contourPolygons[i]);
+                    m_contourLevelCumulativeAreas[i] = levelOuterArea;
+                }
             }
         }
     }
@@ -1022,12 +1070,43 @@ void RimContourMapProjection::smoothContourPolygons(ContourPolygons*       conto
                 if (!intersections.empty())
                 {
                     polygon.vertices = intersections.front();
+                    polygon.area     = std::abs(cvf::GeometryTools::signedAreaPlanarPolygon(cvf::Vec3d::Z_AXIS, polygon.vertices));
                 }
             }
         }
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimContourMapProjection::sumPolygonArea(const ContourPolygons& contourPolygons)
+{
+    double sumArea = 0.0;
+    for (const ContourPolygon& polygon : contourPolygons)
+    {
+        sumArea += polygon.area;
+    }
+    return sumArea;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimContourMapProjection::sumTriangleAreas(const std::vector<cvf::Vec4d>& triangles)
+{
+    double sumArea = 0.0;
+    for (size_t i = 0; i < triangles.size(); i += 3)
+    {
+        cvf::Vec3d v1(triangles[i].x(), triangles[i].y(), triangles[i].z());
+        cvf::Vec3d v2(triangles[i + 1].x(), triangles[i + 1].y(), triangles[i + 1].z());
+        cvf::Vec3d v3(triangles[i + 2].x(), triangles[i + 2].y(), triangles[i + 2].z());
+        double area = 0.5 * ((v3 - v1) ^ (v2 - v1)).length();
+        sumArea += area;
+    }
+    return sumArea;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimContourMapProjection.h

@@ -49,6 +49,7 @@ public:
     {
         std::vector<cvf::Vec3d> vertices;
         double                  value;
+        double                  area;
         cvf::BoundingBox        bbox;
     };
 
@@ -146,6 +147,8 @@ protected:
     std::vector<cvf::Vec3d> generateVertices() const;
     void                    generateContourPolygons();
     void                    smoothContourPolygons(ContourPolygons* contourPolygons, const ContourPolygons* clipBy, bool favourExpansion);
+    static double           sumPolygonArea(const ContourPolygons& contourPolygons);
+    static double           sumTriangleAreas(const std::vector<cvf::Vec4d>& triangles);
 
     std::vector<CellIndexAndResult> cellOverlapVolumesAndResults(const cvf::Vec2d&          globalPos2d,
                                                                  const std::vector<double>& weightingResultValues) const;
@@ -208,6 +211,7 @@ protected:
     cvf::BoundingBox                      m_gridBoundingBox;
     double                                m_sampleSpacing;
     std::vector<ContourPolygons>          m_contourPolygons;
+    std::vector<double>                   m_contourLevelCumulativeAreas;
     std::vector<cvf::Vec4d>               m_trianglesWithVertexValues;
     int                                   m_currentResultTimestep;
 

ApplicationCode/ProjectDataModel/RimGeoMechContourMapProjection.cpp

@@ -141,7 +141,7 @@ void RimGeoMechContourMapProjection::ensureOnlyValidPorBarVisible(cvf::UByteArra
     const std::vector<float>& resultValues = resultCollection->resultValues(porBarAddr, 0, timeStep);
     for (int i = 0; i < static_cast<int>(visibility->size()); ++i)
     {
-        size_t resValueIdx = m_femPart->elementNodeResultIdx(i, 0);
+        size_t resValueIdx = m_femPart->elementNodeResultIdx((int) i, 0);
         double scalarValue = resultValues[resValueIdx];
         (*visibility)[i] &= scalarValue != std::numeric_limits<double>::infinity();
     }

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp

@@ -124,7 +124,9 @@ std::array<cvf::Vec3d, 8> RigCellGeometryTools::estimateHexOverlapWithBoundingBo
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons)
+void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>>& intersectionLineSegments,
+                                                         std::vector<std::vector<cvf::Vec3d>>&         polygons,
+                                                         double                                        tolerance)
 {
     bool startNewPolygon = true;
     while (!intersectionLineSegments.empty())
@@ -146,7 +148,6 @@ void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf
         //Search remaining list for next point...
 
         bool isFound = false;
-        float tolerance = 0.0001f;
 
         for (std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >::iterator lIt = intersectionLineSegments.begin(); lIt != intersectionLineSegments.end(); lIt++)
         {
@@ -191,9 +192,6 @@ void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf
             startNewPolygon = true;
         }
     }
-
-
-
 }
 
 //==================================================================================================

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h

@@ -37,7 +37,8 @@ public:
                                                                        cvf::BoundingBox*                overlapBoundingBox);
 
     static void createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>>& intersectionLineSegments,
-                                              std::vector<std::vector<cvf::Vec3d>>&         polygons);
+                                              std::vector<std::vector<cvf::Vec3d>>&         polygons,
+                                              double                                        tolerance = 1.0e-4);
 
     static void findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners,
                                  cvf::Vec3d&                      localXdirection,

ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -156,6 +156,26 @@ double GeometryTools::getAngle(const cvf::Vec3d& v1, const cvf::Vec3d& v2)
     return angle;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double GeometryTools::signedAreaPlanarPolygon(const cvf::Vec3d& planeNormal, const std::vector<cvf::Vec3d>& polygon)
+{
+    int Z = findClosestAxis(planeNormal);
+    int X = (Z + 1) % 3;
+    int Y = (Z + 2) % 3;
+
+    // Use Shoelace formula to calculate signed area.
+    // https://en.wikipedia.org/wiki/Shoelace_formula
+    double signedArea = 0.0;
+    for (size_t i = 0; i < polygon.size(); ++i)
+    {
+        signedArea += (polygon[(i + 1) % polygon.size()][X] - polygon[i][X]) *
+                      (polygon[(i + 1) % polygon.size()][Y] + polygon[i][Y]);
+    }
+    return signedArea;
+}
+
 /*
    Determine the intersection point of two line segments  
    From Paul Bourke, but modified to really handle coincident lines
@@ -178,7 +198,7 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect(
    numera = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
    numerb = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
 
-    double EPS = 1e-40;
+   double EPS = 1e-40;
 
    // Are the line coincident? 
    if (fabs(numera) < EPS && fabs(numerb) < EPS && fabs(denom) < EPS) 
@@ -199,7 +219,7 @@ GeometryTools::IntersectionStatus inPlaneLineIntersect(
 
        // Check if the p1 p2 line is a point
       
-       if (length12 < EPS )
+       if (length12 < EPS)
        {
            *x = x1;
            *y = y1;

ApplicationCode/ReservoirDataModel/cvfGeometryTools.h

@@ -58,6 +58,8 @@ public:
     static double     getAngle(const cvf::Vec3d& positiveNormalAxis, const cvf::Vec3d& v1, const cvf::Vec3d& v2);
     static double     getAngle(const cvf::Vec3d& v1, const cvf::Vec3d& v2);
 
+    static double signedAreaPlanarPolygon(const cvf::Vec3d& planeNormal, const std::vector<cvf::Vec3d>& polygon);
+
     static cvf::Vec3d polygonAreaNormal3D(const std::vector<cvf::Vec3d>& polygon);
 
     enum IntersectionStatus

Fwk/AppFwk/CommonCode/cafContourLines.cpp

@@ -214,17 +214,16 @@ void caf::ContourLines::create(const std::vector<double>& dataXY, const std::vec
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::vector<caf::ContourLines::ClosedPolygons> caf::ContourLines::create(const std::vector<double>& dataXY,
+std::vector<caf::ContourLines::ListOfLineSegments> caf::ContourLines::create(const std::vector<double>& dataXY,
-    const std::vector<double>& xPositions,
+                                                                             const std::vector<double>& xPositions,
-    const std::vector<double>& yPositions,
+                                                                             const std::vector<double>& yPositions,
-    const std::vector<double>& contourLevels,
+                                                                             const std::vector<double>& contourLevels)
-    double                     areaThreshold)
 {
     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());
+    std::vector<ListOfLineSegments> listOfSegmentsPerLevel(contourLevels.size());
 
     for (size_t i = 0; i < contourLevels.size(); ++i)
     {
@@ -232,84 +231,16 @@ std::vector<caf::ContourLines::ClosedPolygons> caf::ContourLines::create(const s
         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;
+            ListOfLineSegments 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]));
+                unorderedSegments.push_back(std::make_pair(cvf::Vec3d(contourLineSegments[i][j]), cvf::Vec3d(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. Use Shoelace formula to calculate signed area.
-                    // https://en.wikipedia.org/wiki/Shoelace_formula
-                    double signedArea = 0.0;
-                    for (size_t j = 0; j < closedPolygonDeque.size() - 1; ++j)
-                    {
-                        signedArea += (closedPolygonDeque[j + 1].x() - closedPolygonDeque[j].x()) *
-                                      (closedPolygonDeque[j + 1].y() + closedPolygonDeque[j].y());
-                    }
-                    if (std::abs(signedArea) > areaThreshold)
-                    {
-                        if (signedArea < 0.0)
-                        {
-                            closedPolygonsPerLevel[i].emplace_back(closedPolygonDeque.rbegin(), closedPolygonDeque.rend());
-                        }
-                        else
-                        {
-                            closedPolygonsPerLevel[i].emplace_back(closedPolygonDeque.begin(), closedPolygonDeque.end());
-                        }
-                    }
-                    closedPolygonDeque.clear();
-                }
             }
+            listOfSegmentsPerLevel[i] = unorderedSegments;
         }
     }
-    return closedPolygonsPerLevel;
+
+    return listOfSegmentsPerLevel;
 }
 
 //--------------------------------------------------------------------------------------------------

Fwk/AppFwk/CommonCode/cafContourLines.h

@@ -25,9 +25,11 @@
 
 #include "cvfBase.h"
 #include "cvfVector2.h"
+#include "cvfVector3.h"
 
 #include <deque>
 #include <limits>
+#include <list>
 #include <vector>
 
 namespace caf
@@ -35,14 +37,13 @@ namespace caf
 class ContourLines
 {
 public:
-    typedef std::vector<cvf::Vec2d> ClosedPolygon;
+    typedef std::pair<cvf::Vec3d, cvf::Vec3d> LineSegment;
-    typedef std::vector<ClosedPolygon> ClosedPolygons;
+    typedef std::list<LineSegment> ListOfLineSegments;
 
-    static std::vector<ClosedPolygons> create(const std::vector<double>&            dataXY,
+    static std::vector<ListOfLineSegments> create(const std::vector<double>&            dataXY,
-                                              const std::vector<double>&            xPositions,
+                                                  const std::vector<double>&            xPositions,
-                                              const std::vector<double>&            yPositions,
+                                                  const std::vector<double>&            yPositions,
-                                              const std::vector<double>&            contourLevels,
+                                                  const std::vector<double>&            contourLevels);
-                                              double                                areaTreshold = std::numeric_limits<double>::infinity());
     
 private:
     static void create(const std::vector<double>& dataXY,