diff --git a/src/utilities/contourize.cpp b/src/utilities/contourize.cpp
index 3cee1fc744b8f71efe0ee60eddecf0e099cc9827..192877555e911f0b5544bacf22cbf87bd97bd0d9 100644
--- a/src/utilities/contourize.cpp
+++ b/src/utilities/contourize.cpp
@@ -1,8 +1,12 @@
 #include <vector>
 #include <triangle.h>
 #include <orientation.h>
+#include <triangle_edges.h>
+#include <stack_vector.h>
+#include <intersections.h>
 
-Triangle clockwiseToCounterclockwise(const Triangle &t) {
+Triangle clockwiseToCounterclockwise(const Triangle &t)
+{
     return Triangle(t.points[0], t.points[2], t.points[1], t.depth);
 }
 
@@ -30,53 +34,96 @@ Point nextPoint(const Point &previous, std::vector<Point> &candidates, const std
             return p;
         }
     }
-    
+
     return Point{-69,-69};
-} 
+}
 */
 
-std::vector<Point> contourize(const Triangle &t1, const Triangle &t2, const std::vector<Point> &intersections) {
+float distanceSquared(const Point &p1, const Point &p2)
+{
+    return (p2.x - p1.x) * (p2.x - p1.x) + (p2.y - p1.y) * (p2.y - p1.y);
+}
+
+bool closestPoint(const Point &refPoint, const Point &p1, const Point &p2)
+{
+    float d1 = distanceSquared(refPoint, p1);
+    float d2 = distanceSquared(refPoint, p2);
+    return d1 < d2;
+}
+
+struct IntersectionAndEdge
+{
+    Point p;
+    int edge;
+};
+
+std::vector<Point> contourize(const Triangle &leftmostTriangle, const Triangle &otherTriangle, int leftmostInd)
+{
     std::vector<Point> result;
-    const int NUM_TRIANGLE_POINTS = 3;
-    // using convex hull algorithm
-    // start at leftmost point
 
-    // get points
-    // initially ignoring intersection points because they cannot be the leftmost
-    std::vector<Point> points;
+    const TriangleEdges lftEdges(leftmostTriangle);
+    const TriangleEdges otherEdges(otherTriangle);
+    StaticVector<IntersectionAndEdge> intersectionVector;
+
+    bool usingLeftTriangle = true;
+
+    int prev = leftmostInd;
+    do
+    {
+        const Triangle &currentTriangle = usingLeftTriangle ? leftmostTriangle : otherTriangle;
+        const TriangleEdges &other = !usingLeftTriangle ? lftEdges : otherEdges;
+        result.push_back(currentTriangle.points[prev]);
 
-    points.insert(points.end(), t1.points, t1.points + NUM_TRIANGLE_POINTS);
-    points.insert(points.end(), t2.points, t2.points + NUM_TRIANGLE_POINTS);
+        int nextPointInTriangle = (prev + 1) % NB_TRIANGLE_SIDES;
+        const Edge &e1 = Edge{currentTriangle.points[prev], currentTriangle.points[nextPointInTriangle]};
 
-    int minX = 0;
-    for (int i = 0; i < points.size(); i++) {
-        if (points[minX].x > points[i].x) {
-            minX = i;
+        for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+        {
+            auto intrs = intersection(e1, other.edges[i]);
+            if (intrs.has_value())
+            {
+                intersectionVector.push_back(IntersectionAndEdge{intrs.value(), i});
+            }
         }
-    }
 
-    // add intersection points
-    points.insert(points.end(), intersections.begin(), intersections.end());
+        if (intersectionVector.getSize() == 0)
+        {
+            prev = nextPointInTriangle;
+            continue;
+        }
+        else
+        {
+            // intersections found. Must switch to next triangle
+            usingLeftTriangle = !usingLeftTriangle;
+            int edgeIndex = intersectionVector.items[0].edge;
+            if (intersectionVector.getSize() == 2 && !closestPoint(currentTriangle.points[prev], intersectionVector.items->p, intersectionVector.items[1].p))
+            {
+                // closest edge
+                edgeIndex = intersectionVector.items[1].edge;
+                result.push_back(intersectionVector.items[1].p);
+            }
+            else
+            {
+                result.push_back(intersectionVector.items[0].p);
+            }
 
-    int previous = minX;
-    int candidate;
+            const Edge &e2 = other.edges[edgeIndex];
+            const std::pair<int, int> edgePointIndexes = TriangleEdges::otherPoint(edgeIndex);
 
-    do {
-        
-        result.push_back(points[previous]);
-        candidate = (previous + 1) % points.size();
-        for (int i = 0; i < points.size(); i++) {
-            if (orientation(points[previous], points[candidate], points[i]) == Counterclockwise) {
-                candidate = i;
+            if (!e1.positiveSide(e2.p1))
+            {
+                prev = edgePointIndexes.first;
+            }
+            else
+            {
+                prev = edgePointIndexes.second;
             }
         }
-        previous = candidate;
-    } while (previous != minX);
 
+    } while (!(prev == leftmostInd));
     return result;
 }
 
-
 /*
 
 std::vector<Point> contourize(const Triangle &t1, const Triangle &t2, const std::vector<Point> &newIntersections) {
diff --git a/tests/unit_tests/contourize_tests.cpp b/tests/unit_tests/contourize_tests.cpp
index 98b279b3c8fa45db6b8309f538c1aa45bf33387c..00e1d5e1685b541631fe4599f41c4f8b63c1e162 100644
--- a/tests/unit_tests/contourize_tests.cpp
+++ b/tests/unit_tests/contourize_tests.cpp
@@ -21,4 +21,20 @@ TEST(ContourizeTest, TwoTriangles) {
     Point{2,3}};
     EXPECT_EQ(contour, expectedContour);
 
+}
+
+TEST (ContourizeTests, TriangleSingleDipTest) {
+    auto t1 = Triangle({0,0}, {5,0}, {2,6}, 0);
+    auto t2 = Triangle({1,1}, {7,3}, {8,6}, 0);
+
+    auto newIntersections = intersections(t1, t2);
+
+    auto contour = contourize(t1, t2, newIntersections);
+
+    std::vector<Point> expectedContour{
+        {0,0}, {5,0}, {5,2}, {7,3}, {8,6}, {4,3}, {2,6}
+    };
+
+    EXPECT_EQ(contour, expectedContour);
+    
 }
\ No newline at end of file