diff --git a/CMakeLists.txt b/CMakeLists.txt
index 270fdfceb833cb4770462d3e4457e6efb97eafdd..b4523db3404406294ff5894a82b579bbc5686dca 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,7 +23,6 @@ project(intersection)
 add_executable(
 	tests
 	tests/unit_tests/simple_triangle_tests.cpp
-	tests/unit_tests/contourize_tests.cpp
 	tests/unit_tests/point_in_triangle_tests.cpp
 	tests/unit_tests/triangulation_tests.cpp
 	tests/unit_tests/ear_removal_tests.cpp
diff --git a/include/contourize.h b/include/contourize.h
index 2660423519ae5f8062622b2b18056674359c0eb6..cb605f7fc51edc95bbcf6d8de759e0892e6bc471 100644
--- a/include/contourize.h
+++ b/include/contourize.h
@@ -5,6 +5,5 @@ struct Triangle;
 struct Point;
 
 
-std::vector<Point> contourize(const Triangle &t1, const Triangle &t2);
-
-std::vector<Point> sortPoints(std::vector<Point> &points);
\ No newline at end of file
+// creates a contour from a list of points 
+std::vector<Point> contourize(const std::vector<Point> &points);
\ No newline at end of file
diff --git a/include/edge.h b/include/edge.h
index df43ce0af78da01cf16ccd5757e20c487ad15b3f..144fe3f5b96b796b208d64df39e63d738e1b247a 100644
--- a/include/edge.h
+++ b/include/edge.h
@@ -5,4 +5,5 @@ struct Edge {
 	Point p1, p2;
 
     bool positiveSide(const Point &p) const;
+
 };
\ No newline at end of file
diff --git a/include/intersections.h b/include/intersections.h
index 029dd406123a4a5515bde8ef23bfedaebd54b95f..83326ed4fe6cd25a72c00d5da0b275cd4a63f9c6 100644
--- a/include/intersections.h
+++ b/include/intersections.h
@@ -11,6 +11,8 @@ struct TrigTrigInterResults {
     std::vector<Point> results;
 };
 
+std::optional<Point> intersectionWithinEdge(const Edge &e1, const Edge &e2);
 std::optional<Point> intersection(const Edge &e1, const Edge &e2);
 
-std::vector<Point> intersections(const Triangle &t1, const Triangle &t2);
\ No newline at end of file
+std::vector<Point> intersections(const Triangle &t1, const Triangle &t2);
+std::optional<Point> intersectionWithinEdgeDirection(const Edge &e1, const Edge &e2);
\ No newline at end of file
diff --git a/src/union.cpp b/src/union.cpp
index 9e17540353678311fcd6380f9516aaa309e4a353..0061093ca5170c1501e1264022fb4d4b7be418c3 100644
--- a/src/union.cpp
+++ b/src/union.cpp
@@ -1,39 +1,70 @@
 #include <vector>
-#include <optional>
-#include <union.h>
-#include <contourize.h>
+#include <triangle.h>
+#include <triangle_edges.h>
 #include <intersections.h>
+#include <contourize.h>
 #include <triangulation.h>
 
-std::vector<Triangle> unionize2(const Triangle &t1, const Triangle &t2) {
-	// if neighbours, do nothing
-	if (t1.neighbours(t2)) {
-		return std::vector{t1, t2};
-	}
-
-	// look for points in triangles
-	// if 0 points in triangles
-	// 	check if triangles intersect
-	// if 1 point in triangle
-	// 	
-}
+std::vector<Point> getPointsOnSide(const Edge &e, const std::vector<Point> intr, const Triangle &bottom)
+{
+    std::vector<Point> result;
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        if (e.positiveSide(bottom.points[i]))
+        {
+            result.push_back(bottom.points[i]);
+        }
+    }
 
-std::vector<Triangle> unionize(const Triangle &t1, const Triangle &t2) {
+    for (const auto &p : intr)
+    {
+        if (e.positiveSide(p))
+        {
+            result.push_back(p);
+        }
+    }
+    return result;
+}
 
-	// if neighbours, do nothing
-	if (t1.neighbours(t2)) {
-		return std::vector{t1, t2};
-	}
+std::vector<Triangle> unionizeTopAndBottom(const Triangle &top, const Triangle &bottom)
+{
+    if (intersections(top, bottom).empty()) {
+        return {};
+    }
+    std::vector<Triangle> result;
+    TriangleEdges topEdges = TriangleEdges(top);
+    TriangleEdges botEdges = TriangleEdges(bottom);
 
-	// at most 3? if infinite -> consider no intersections
-	std::vector<Point> newIntersections = intersections(t1, t2);
+    std::vector<Point> intrPoints;
 
-	if (newIntersections.empty()) {
-		// either completely envelops each other, no intersections	
-		return std::vector{t1, t2};
-	}
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        for (int j = 0; j < NB_TRIANGLE_SIDES; j++)
+        {
+            auto cand = intersectionWithinEdgeDirection(topEdges.edges[i], botEdges.edges[i]);
+            if (cand.has_value())
+            {
+                intrPoints.push_back(cand.value());
+            }
+        }
+    }
 
-	std::vector<Point> contour = contourize(t1, t2);
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        std::vector<Point> currentShape = getPointsOnSide(topEdges.edges[i], intrPoints, bottom);
+        std::vector<Point> orderedPoints = contourize(currentShape);
+        std::vector<Triangle> newTriangles = triangulate(orderedPoints);
+        result.insert(result.begin(), newTriangles.begin(), newTriangles.end());
+    }
 
-	return triangulate(contour);
+    return result;
 }
+
+std::vector<Triangle> unionize(const Triangle &t1, const Triangle &t2)
+{
+    if (t1.depth < t2.depth)
+    {
+        return unionizeTopAndBottom(t1, t2);
+    }
+    return unionizeTopAndBottom(t2, t1);
+}
\ No newline at end of file
diff --git a/src/utilities/complete_triangle.cpp b/src/utilities/complete_triangle.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ae71ebac8098051a7c788122d6a3b50a62e2481b
--- /dev/null
+++ b/src/utilities/complete_triangle.cpp
@@ -0,0 +1,6 @@
+#include <vector>
+#include <triangle.h>
+
+std::vector<Triangle> completeTriangles() {
+    // 
+}
\ No newline at end of file
diff --git a/src/utilities/contourize.cpp b/src/utilities/contourize.cpp
index 62528522be1a3a92b5e6dbbb81953f3c8b80f3bd..078d0bf2b8c514df711bb3c3b828baa1428fc8d4 100644
--- a/src/utilities/contourize.cpp
+++ b/src/utilities/contourize.cpp
@@ -1,141 +1,33 @@
-#include <vector>
-#include <triangle.h>
+#include <contourize.h>
+#include <point.h>
 #include <orientation.h>
-#include <triangle_edges.h>
-#include <stack_vector.h>
-#include <intersections.h>
+#include <list>
 
-Triangle clockwiseToCounterclockwise(const Triangle &t)
-{
-    return Triangle(t.points[0], t.points[2], t.points[1], t.depth);
-}
-
-
-
-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 &t1, const Triangle &t2)
-{
-    int leftmostInd = 0;
-    bool t1Leftmost = true;
-    for (int i = 0; i< NB_TRIANGLE_SIDES; i++) {
-        if (t1.points[i].x < (t1Leftmost ? t1 : t2).points[leftmostInd].x) {
-            leftmostInd = i;
-            t1Leftmost = true;
-        }
-        if (t2.points[i].x < (t1Leftmost ? t1 : t2).points[leftmostInd].x) {
-            leftmostInd = i;
-            t1Leftmost = false;
-        }
-    }
-
-   
-    std::vector<Point> result;
-
-    const Triangle &leftmostTriangle =t1Leftmost ? t1 : t2;
-    const Triangle &otherTriangle = t1Leftmost ? t2 : t1;  
-
-    const TriangleEdges lftEdges(leftmostTriangle);
-    const TriangleEdges otherEdges(otherTriangle);
-
-    bool usingLeftTriangle = true;
-
-    StaticVector<IntersectionAndEdge> intersectionVector;
-    int prev = leftmostInd;
-    do
-    {
-        const Triangle &currentTriangle = usingLeftTriangle ? leftmostTriangle : otherTriangle;
-        const TriangleEdges &other = !usingLeftTriangle ? lftEdges : otherEdges;
-        result.push_back(currentTriangle.points[prev]);
-
-        int nextPointInTriangle = (prev + 1) % NB_TRIANGLE_SIDES;
-        const Edge &e1 = Edge{currentTriangle.points[prev], currentTriangle.points[nextPointInTriangle]};
-
-        intersectionVector.clear();
-        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});
-            }
-        }
-
-        if (intersectionVector.getSize() == 0)
-        {
-            prev = nextPointInTriangle;
+bool isCounterClockwiseForAll(const Point &previous, const Point &candidate, const std::vector<Point> & points) {
+    for (const auto &point : points) {
+        if (previous == point || candidate == point) {
             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);
-            }
-
-            const Edge &e2 = other.edges[edgeIndex];
-            const std::pair<int, int> edgePointIndexes = TriangleEdges::otherPoint(edgeIndex);
-
-            if (e1.positiveSide(e2.p1))
-            {
-                prev = edgePointIndexes.first;
-            }
-            else
-            {
-                prev = edgePointIndexes.second;
-            }
+        if (orientation(previous, point, candidate) != Counterclockwise) {
+            return false;
         }
-
-    } while (!(prev == leftmostInd && usingLeftTriangle));
-    return result;
+    }
+    return true;
 }
 
-/*
-
-std::vector<Point> contourize(const Triangle &t1, const Triangle &t2, const std::vector<Point> &newIntersections) {
+std::vector<Point> contourize(const std::vector<Point> &points) {
+    std::list<Point> candidates(points.begin(), points.end());
     std::vector<Point> result;
-    // Go through each event point.
-    // If an even point is in a triangle, it is not part of the final Triangle
-    for (auto &p : t1.points) {
-        if (!t2.pointInTriangle(p)) {
-            result.push_back(p);
-        }
-    }
 
-    for (auto &p : t2.points) {
-        if (!t2.pointInTriangle(p)) {
-            result.push_back(p);
+    Point previous = points[0];
+    while (!candidates.empty()) {
+        if (isCounterClockwiseForAll(previous, candidates.front(), points)) {
+            previous = candidates.front();
+            candidates.pop_front();
+        } else {
+            candidates.push_back(candidates.front());
+            candidates.pop_front();
         }
     }
-
-    result.insert(result.end(), newIntersections.begin(), newIntersections.end());
-
     return result;
-}
-
-*/
\ No newline at end of file
+}
\ No newline at end of file
diff --git a/src/utilities/intersections.cpp b/src/utilities/intersections.cpp
index bf6c299d23370f956e37a46e1465250fb62595d8..8c22dd79f833aa9a02e1866f25161c279086a4db 100644
--- a/src/utilities/intersections.cpp
+++ b/src/utilities/intersections.cpp
@@ -1,9 +1,9 @@
 #include "intersections.h"
 #include <optional>
-#include <edge.h>
 #include <triangle_edges.h>
 #include <triangle.h>
 
+#include <edge.h>
 
 
 std::optional<float> getB(std::optional<float> slope, Point p) {
@@ -32,12 +32,39 @@ bool withinEdge(Edge e, Point p) {
 }
 
 void intersection(Edge e1, Edge e2, std::vector<Point> &results) {
-    auto point = intersection(e1, e2);
+    auto point = intersectionWithinEdge(e1, e2);
     if (point.has_value()) {
         results.push_back(point.value());
     }
 }
 
+
+std::optional<Point> intersectionWithinEdge(const Edge &e1, const Edge &e2) {
+    std::optional<Point> candPoint = intersection(e1, e2);
+    if (candPoint.has_value() && 
+    withinEdge(e1, candPoint.value()) &&
+     withinEdge(e2, candPoint.value())) {
+        return candPoint;
+    }
+    return {};
+}
+
+// returns intersection with e1 being extended infinitly in its direction
+std::optional<Point> intersectionWithinEdgeDirection(const Edge &e1, const Edge &e2) {
+    auto candPoint = intersection(e1, e2);
+
+    if (candPoint.has_value()) {
+        auto s1 = getSlope(Edge{e1.p1, candPoint.value()});
+        auto s2 = getSlope(e1);
+
+        // check if both slopes have the same sign
+        if (s1.value() * s2.value() >= 0) {
+            return candPoint;
+        }
+    }
+    return {};
+}
+
 std::optional<Point> intersection(const Edge &e1, const Edge &e2) {
     auto slope1 = getSlope(e1);
     auto slope2 = getSlope(e2);
@@ -61,10 +88,7 @@ std::optional<Point> intersection(const Edge &e1, const Edge &e2) {
     float candX = (b2.value() - b1.value()) / (slope1.value() - slope2.value());
     auto candPoint = Point{ candX, slope1.value() * candX + b1.value()};
 
-    if (withinEdge(e1, candPoint) && withinEdge(e2, candPoint)) {
         return candPoint;
-    }
-    return {};
 }
 void intersections(const Edge &e1, const TriangleEdges &te, std::vector<Point> &results) {
     for (int i = 0; i < NB_TRIANGLE_SIDES; i++) {
diff --git a/tests/unit_tests/intersections_tests.cpp b/tests/unit_tests/intersections_tests.cpp
index b0f12893cf63239c4d4f208aa1d2d56b4e6fa9df..6d888d28f64c68a3c9fa1770a0c212e86066c673 100644
--- a/tests/unit_tests/intersections_tests.cpp
+++ b/tests/unit_tests/intersections_tests.cpp
@@ -6,7 +6,7 @@ TEST (IntersectionTests, EdgeToEdgeTest) {
     Edge e1 {{5,0}, {2,6}};
     Edge e2 {{1,1}, {4,2}};
 
-    auto p = intersection(e1, e2);
+    auto p = intersectionWithinEdge(e1, e2);
     Point actual_point = p.value();
     Point expected_point = Point{4,2};
     EXPECT_TRUE(actual_point == expected_point);