include/triangle.h

 #pragma once
 #include "point.h"
+#include "print_triangle.h"
+#include <vector>
+#include <memory>
 
-// points specified counterclockwise
-struct Triangle {
-	Point points[3];
-	int depth;
-    bool neighbours(const Triangle &other) const;
-    Triangle(Point p1, Point p2, Point p3, int depth);
+// note : Triangle equality is only specified for tests
+
+// points specified in counterclockwise or collinear (line / point) order
+struct Triangle
+{
+    Point points[3];
+    int mainTriangleId;
+    int fragmentId;
+    std::shared_ptr<std::vector<int>> neighbours;
+
+
+    Triangle(const Point &p1, const Point &p2, const Point &p3, int mainTriangleId = 0, const std::vector<int> neighbours = {});
+    Triangle(const Point &p1, const Point &p2, const Point &p3, int mainTriangleId, std::shared_ptr<std::vector<int>> neighbours);
     bool pointInTriangle(const Point &p) const;
     Point nextPoint(int pointIndex) const;
+    bool operator==(const Triangle &other) const;
 };

include/triangle_fragment_id_assigner.h

+#pragma once
+
+class TriangleFragmentIdAssigner {
+   public:
+        static TriangleFragmentIdAssigner& getInstance()
+        {
+            static TriangleFragmentIdAssigner instance; 
+            return instance;
+        }
+        int generateUniqueId() {
+            currentInd++;
+            return currentInd;
+        };
+    private:
+        int currentInd = 0;
+        TriangleFragmentIdAssigner() {}       
+
+    public:
+        TriangleFragmentIdAssigner(TriangleFragmentIdAssigner const&) = delete;
+        void operator=(TriangleFragmentIdAssigner const&) = delete;
+
+};
\ No newline at end of file

include/triangulation.h

-#pragma once
-#include <vector>
-#include <optional>
-#include <triangle.h>
-
-class PointList;
-struct Point;
-
-std::optional<Triangle> removeEar(int &index, PointList &pointList, const std::vector<Point> &allPoints);
-
-std::vector<Triangle> triangulate(std::vector<Point> points);
-

src/data_structures/box.cpp

+#include <box.h>
+#include <triangle.h>
+#include <constants.h>
+
+Box::Box(Point topLeft, Point bottomRight) : 
+topLeft{topLeft},
+ bottomRight{bottomRight},
+  points{bottomRight, {topLeft.x, bottomRight.y}, topLeft, {bottomRight.x, topLeft.y}},
+  edges{{points}} {}
+
+
+bool Box::intersects(const Triangle &t) const {
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++) {
+        float x = t.points[i].x;
+        float y = t.points[i].y;
+        if (intersects({x,y})) {
+            return true;
+        }
+    }
+
+    for (int i = 0; i < BOX_NB_POINTS; i++) {
+        if (t.pointInTriangle(points[i])) {
+            return true;
+        }
+    }
+    return false;
+}
+
+bool Box::intersects(const Point &p) const {
+    return p.x <= topLeft.x && p.y <= topLeft.y && bottomRight.x <= p.x && bottomRight.y <= p.y; 
+}
+
+
+Box Box::firstQuadrant() const {
+    return Box({topLeft.x / 2, topLeft.y}, {bottomRight.x, topLeft.y / 2});
+}
+
+Box Box::secondQuadrant() const {
+    return Box(topLeft, {topLeft.x / 2, topLeft.y / 2});
+}
+
+Box Box::thirdQuadrant() const {
+    return Box({topLeft.x, topLeft.y / 2}, {topLeft.x/2,bottomRight.y});
+}
+
+Box Box::fourthQuadrant() const {
+    return Box({topLeft.x / 2, topLeft.y / 2}, bottomRight);
+}
\ No newline at end of file

src/data_structures/box_edges.cpp

+#include <box_edges.h>
+#include <box.h>
+
+BoxEdges::BoxEdges(const Point p [size]) {
+    for (int i = 0; i < size; i++) {
+        e[i] = {p[i], p[(i + 1) % size]};
+    }
+}
\ No newline at end of file

src/data_structures/quad_tree.cpp

+#include <quad_tree.h>
+#include <union.h>
+#include <iterator>
+#include <set>
+
+QuadTree::QuadTree(Box b) : b{b}, level{0} {}
+
+QuadTree::QuadTree(Box b, int level = 0) : b{b}, level{level} {}
+
+void QuadTree::collectUniqueTriangleFragments(const Triangle &t, std::set<int> &seen, std::vector<Triangle> &result) const {
+    if (seen.count(t.fragmentId) > 0) {
+        return;
+    }
+
+    result.push_back(t);
+    seen.insert(t.fragmentId);
+}
+
+std::vector<Triangle> QuadTree::visibleSurface() const {
+    std::set<int> seen;
+    return visibleSurface(seen);
+}
+
+std::vector<Triangle> QuadTree::visibleSurface(std::set<int> & seen) const
+{
+    std::vector<Triangle> result;
+
+    for (const Triangle &t : triangles) {
+        collectUniqueTriangleFragments(t, seen, result);
+    }
+    for (const QuadTree &c : children)
+    {
+        std::vector<Triangle> childResult = c.visibleSurface(seen);
+        result.insert(result.end(), childResult.begin(), childResult.end());
+    }
+
+    return result;
+}
+
+void QuadTree::split()
+{
+    const int nextLevel = level + 1;
+    children = {
+        QuadTree(b.firstQuadrant(), nextLevel),
+        QuadTree(b.secondQuadrant(), nextLevel),
+        QuadTree(b.thirdQuadrant(), nextLevel),
+        QuadTree(b.fourthQuadrant(), nextLevel)};
+}
+
+void QuadTree::addTriangle(const Triangle &triangle)
+{
+    if (!b.intersects(triangle))
+    {
+        return;
+    }
+
+    if (triangles.empty())
+    {
+        triangles.push_back(triangle);
+        return;
+    }
+
+    if (!children.empty())
+    {
+        for (QuadTree &child : children)
+        {
+            child.addTriangle(triangle);
+        }
+        return;
+    }
+
+    std::vector<Triangle> currentTriangleFragments{triangle};
+
+    // triangles that are not the current triangle and that have already been unioned
+    std::vector<Triangle> otherTriangleFragments;
+    for (const auto &otherTriangle : triangles)
+    {
+        for (const auto &fragment : currentTriangleFragments)
+        {
+            auto newTriangles = unionize(triangle, otherTriangle);
+            std::vector<Triangle> bottoms;
+            if (newTriangles.size() > 1)
+            {
+                bottoms = {newTriangles.begin(), newTriangles.begin() + newTriangles.size() - 2};
+            }
+            const Triangle &top = newTriangles.back();
+            otherTriangleFragments.push_back(top);
+            if (top.mainTriangleId == triangle.mainTriangleId)
+            {
+                otherTriangleFragments.insert(otherTriangleFragments.begin(), bottoms.begin(), bottoms.end());
+            }
+            else
+            {
+                currentTriangleFragments.insert(currentTriangleFragments.begin(), bottoms.begin(), bottoms.end());
+            }
+        }
+    }
+
+    // split
+    if (level >= QUADTREE_MAX_DEPTH) {
+        triangles.insert(triangles.end(), currentTriangleFragments.begin(), currentTriangleFragments.end());
+        triangles.insert(triangles.end(), otherTriangleFragments.begin(), otherTriangleFragments.end());
+        return;
+    }
+    split();
+    for (QuadTree &q : children)
+    {
+        for (const Triangle &t : otherTriangleFragments)
+        {
+            q.addNonIntersectingTriangle(t);
+        }
+        for (const Triangle &t : currentTriangleFragments)
+        {
+            q.addNonIntersectingTriangle(t);
+        }
+    }
+
+    triangles.clear();
+}
+
+int QuadTree::pointIntersection(const Point &p) const
+{
+    if (b.intersects(p))
+    {
+        for (const Triangle &triangle : triangles)
+        {
+            if (triangle.pointInTriangle(p))
+            {
+                return triangle.mainTriangleId;
+            }
+        }
+
+        for (const QuadTree &child : children)
+        {
+            int id = child.pointIntersection(p);
+            if (id != POINT_NOT_IN_QUADTREE)
+            {
+                return id;
+            }
+        }
+    }
+
+    return POINT_NOT_IN_QUADTREE;
+}
+
+void QuadTree::addNonIntersectingTriangle(const Triangle &t)
+{
+    if (b.intersects(t))
+    {
+        triangles.push_back(t);
+    }
+}
\ No newline at end of file

src/debug_utilities/mostly_equal.cpp

+#include <stdlib.h>
+#include <point.h>
+#include <triangle.h>
+
+#define DEBUG_EPSILON 0.0001
+
+bool mostlyEqual(float a, float b)
+{
+    return abs(a - b) < DEBUG_EPSILON;
+}
+
+bool mostlyEqual(const Point &a, const Point &b)
+{
+    return mostlyEqual(a.x, b.x) && mostlyEqual(a.y, b.y) && mostlyEqual(a.z, b.z);
+}
+
+// mostly equal
+bool Triangle::operator==(const Triangle &other) const {
+	return mostlyEqual(points[0], other.points[0]) &&
+		mostlyEqual(points[1], other.points[1]) &&
+		mostlyEqual(points[2], other.points[2]) &&
+		mainTriangleId == other.mainTriangleId;
+}
\ No newline at end of file

src/debug_utilities/print_triangles.cpp

+#include <print_triangle.h>
+#include <triangle.h>
+#include <point.h>
+
+std::ostream &operator<<(std::ostream &os, const Triangle &t) {
+	return os << "Triangle" << "(" << t.points[0] << ", " << t.points[1] << ", " << t.points[2] << ", " << t.mainTriangleId << ")";
+}
+/*
+std::ostream &operator<<(std::ostream &os, const std::vector<Triangle> &triangles) {
+	for (const auto &t : triangles) {
+
+	}
+	return os;
+}*/

src/debug_utilities/shift_triangle.cpp

+#include <shift_triangle.h>
+#include <triangle.h>
+#include <constants.h>
+
+void shiftZ(Triangle &t, float z) {
+    for (int i =0; i < NB_TRIANGLE_SIDES; i++) {
+        t.points[i].z += z;
+    }
+}
\ No newline at end of file

src/shapes/edge.cpp

@@ -4,4 +4,19 @@ bool Edge::positiveSide(const Point &p) const {
     float h = p1.y - p2.y;
     float g = p1.x - p2.x; 
     return -h * (p.x - p1.x) + g * (p.y - p1.y) >= 0;
+}
+
+// Output: produces a list of edges in counterclockwise order
+// Requirements: points are provided in counterclockwise order
+std::vector<Edge> makeEdges(const std::vector<Point> &points) {
+    std::vector<Edge> result;
+    if (points.size() < 2) {
+        return result;
+    }
+
+    for (int i = 0; i < points.size() - 1; i++) {
+        result.push_back(Edge{points[i], points[i + 1]});
+    }
+    result.push_back(Edge{points[points.size() - 1], points[0]});
+    return result;
 }
\ No newline at end of file

src/shapes/point.cpp

 #include "point.h"
 
+
+
 bool Point::operator==(const Point &other) const {
-    return x == other.x && y == other.y;
-}
\ No newline at end of file
+    return (x == other.x) && (y == other.y) && (z == other.z);
+}
+
+std::ostream &operator<<(std::ostream &os, Point const &p) { 
+    return os << "{" << p.x << "," << p.y << "," << p.z << "}";
+}
+
+Point Point::operator-() const {
+    return {-x, -y,-z};
+}
+
+
+Point Point::operator+(const Point &other) const {
+    return {x + other.x, y + other.y, z + other.z};
+}
+
+
+Point Point::operator-(const Point &other) const {
+    return *this + (- other);
+}

src/shapes/triangle.cpp

 #include "triangle.h"
 #include "triangle_edges.h"
 #include "edge.h"
+#include <triangle_fragment_id_assigner.h>
 
-bool Triangle::neighbours(const Triangle &other) const {
-    return false;
-}
+Triangle::Triangle(const Point &p1, const Point &p2, const Point &p3, int mainTriangleId, std::shared_ptr<std::vector<int>> neighbours = {}) : 
+points{p1, p2, p3}, mainTriangleId{mainTriangleId}, fragmentId{TriangleFragmentIdAssigner::getInstance().generateUniqueId()}, neighbours{neighbours} {}
 
-Triangle::Triangle(Point p1, Point p2, Point p3, int depth) : points{p1,p2,p3}, depth{depth} {}
+Triangle::Triangle(const Point &p1, const Point &p2, const Point &p3, int mainTriangleId, const std::vector<int> neighbours) : 
+points{p1, p2, p3}, mainTriangleId{mainTriangleId}, fragmentId{TriangleFragmentIdAssigner::getInstance().generateUniqueId()}, neighbours{make_shared<std::vector<int>>(neighbours)} {}
 
-bool Triangle::pointInTriangle(const Point &p) const {
+bool Triangle::pointInTriangle(const Point &p) const
+{
     // all tests must be positive
     auto edges = TriangleEdges(*this);
 
-    return !edges.e1.positiveSide(p) &&
-    !edges.e2.positiveSide(p) &&
-    !edges.e3.positiveSide(p);
-}
\ No newline at end of file
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        if (edges.edges[i].positiveSide(p))
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+int nextPoint(int pointIndex)
+{
+    return (pointIndex + 1) % 3;
+}

src/union.cpp

-#include <vector>
-#include <optional>
 #include <union.h>
 #include <contourize.h>
+#include <triangle.h>
+#include <constants.h>
+#include <list>
+#include <triangle_edges.h>
 #include <intersections.h>
-#include <triangulation.h>
+#include <split_triangle.h>
+#include <convex_triangulation.h>
+#include <orientation.h>
 
-std::vector<Triangle> unionize(const Triangle &t1, const Triangle &t2) {
+std::vector<Triangle> unionizeTopAndBottom(const Triangle &top, const Triangle &bot)
+{
+    // degenerate line / point case
+    if (orientation(top) == Collinear) {
+        return {bot, top};
+    }
 
-	// if neighbours, do nothing
-	if (t1.neighbours(t2)) {
-		return std::vector{t1, t2};
-	}
+    std::vector<Triangle> result;
 
-	// at most 3? if infinite -> consider no intersections
-	std::vector<Point> newIntersections = intersections(t1, t2);
+    TriangleEdges topEdges = TriangleEdges(top);
 
-	if (newIntersections.empty()) {
-		// either completely envelops each other, no intersections	
-		return std::vector{t1, t2};
-	}
+    // keep track of relevant triangles
+    std::vector<Point> relv = {bot.points[0], bot.points[1], bot.points[2]};
 
-	std::vector<Point> contour = contourize(t1, t2);
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        const Edge &e = topEdges.edges[i];
 
-	return triangulate(contour);
+        // index 0: On the "outside" of the edge 
+        // index 1: On the "inside" of the edge
+        auto shapes = splitShape(relv, e);
+        // split triangle if exists
+        // currently relevant triangles
+        // add these to result
+        if (!shapes[0].empty())
+        {
+            std::vector<Triangle> relvTriangles = convexTriangulation(shapes[0], bot.mainTriangleId, bot.neighbours);
+            result.insert(result.end(), relvTriangles.begin(), relvTriangles.end());
+        }
+        // future relevant triangles
+        relv = shapes[1];
+    }
+
+    result.push_back(top);
+    return result;
+}
+
+std::vector<Triangle> unionize(const Triangle &t1, const Triangle &t2)
+{
+    if (t1.mainTriangleId == t2.mainTriangleId) {
+        return {t1, t2};
+    }
+
+	if (std::find(t1.neighbours->begin(), t1.neighbours->end(), t2.mainTriangleId) != t1.neighbours->end()) {
+		return {t1, t2};
+	}
+    if (t1.points[0].z < t2.points[0].z)
+    {
+        return unionizeTopAndBottom(t1, t2);
+    }
+    return unionizeTopAndBottom(t2, t1);
 }

src/utilities/contourize.cpp

-#include <vector>
-#include <triangle.h>
+#include <contourize.h>
+#include <point.h>
 #include <orientation.h>
+#include <list>
 
-Triangle clockwiseToCounterclockwise(const Triangle &t) {
-    return Triangle(t.points[0], t.points[2], t.points[1], t.depth);
+const char* ContourizeException::what() const throw() {
+    return "Infinite Loop detected";
 }
 
-/*
 bool isCounterClockwiseForAll(const Point &previous, const Point &candidate, const std::vector<Point> & points) {
     for (const auto &point : points) {
         if (previous == point || candidate == point) {
             continue;
         }
-        if (orientation(previous, point, candidate) != Counterclockwise) {
+	auto o = orientation(previous, candidate, point);
+        if (o == Collinear) {
+            continue;
+        }
+
+        if (o != Counterclockwise) {
             return false;
         }
     }
     return true;
 }
-*/
-
-/*
-Point nextPoint(const Point &previous, std::vector<Point> &candidates, const std::vector<Point> &points) {
-    Point mostCounterclockwise = candidates[0];
-    for (int i = 1; i < candidates.size(); i++) {
-        if (isCounterClockwiseForAll(previous, candidates[i], points)) {
-            Point p = candidates[i];
-            candidates.erase(i + candidates.begin());
-            return p;
-        }
-    }
-    
-    return Point{-69,-69};
-} 
-*/
 
-std::vector<Point> contourize(const Triangle &t1, const Triangle &t2, const std::vector<Point> &intersections) {
-    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;
-
-    points.insert(points.end(), t1.points, t1.points + NUM_TRIANGLE_POINTS);
-    points.insert(points.end(), t2.points, t2.points + NUM_TRIANGLE_POINTS);
-
-    int minX = 0;
-    for (int i = 0; i < points.size(); i++) {
-        if (points[minX].x > points[i].x) {
-            minX = i;
-        }
+std::vector<Point> contourize(const std::vector<Point> &points) {
+    if (points.size() < 3) {
+        return points;
     }
 
-    // add intersection points
-    points.insert(points.end(), intersections.begin(), intersections.end());
 
-    int previous = minX;
-    int candidate;
 
-    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;
-            }
-        }
-        previous = candidate;
-    } while (previous != minX);
-
-    return result;
-}
+    std::vector<Point> result;
 
+    Point previous = points[0];
+    result.push_back(previous);
+    std::list<Point> candidates(points.begin() + 1, points.end());
+    const int size = candidates.size();
 
-/*
+    // infinite loop detection
+    int seen = 0;
 
-std::vector<Point> contourize(const Triangle &t1, const Triangle &t2, const std::vector<Point> &newIntersections) {
-    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);
+    while (candidates.size() > 1) {
+        // detect infinite loop
+        if (seen >= candidates.size()) {
+            throw ContourizeException();
         }
-    }
-
-    for (auto &p : t2.points) {
-        if (!t2.pointInTriangle(p)) {
-            result.push_back(p);
+  
+
+        if (isCounterClockwiseForAll(previous, candidates.front(), points)) {
+            previous = candidates.front();
+            candidates.pop_front();
+	        result.push_back(previous);
+            seen = 0;
+        } else {
+            candidates.push_back(candidates.front());
+            candidates.pop_front();
+            seen++;
         }
     }
 
-    result.insert(result.end(), newIntersections.begin(), newIntersections.end());
-
+    result.push_back(candidates.front());
     return result;
 }
-
-*/
\ No newline at end of file

src/utilities/convex_triangulation.cpp

+#include <convex_triangulation.h>
+#include <orientation.h>
+
+// filters out degenerate triangles
+void appendResults(std::vector<Triangle> &results, const Triangle &t) {
+	if (orientation(t) == Collinear) {
+		return;
+	}
+	results.push_back(t);
+}
+
+std::vector<Triangle> convexTriangulation(const std::vector<Point> &points, int triangleId, std::shared_ptr<std::vector<int>> neighbours) {
+	const int size = points.size();
+	if (size < 3) {
+		return std::vector<Triangle>{};
+	}
+
+	const Triangle t1 = Triangle(points[0], points[1], points[2], triangleId, neighbours);
+	std::vector<Triangle> results;
+
+	appendResults(results, t1);
+
+	
+	Point first = points[0];
+	for (int i = size - 1; i >= 3; i--) {
+		const Triangle t = Triangle(first, points[2], points[i], triangleId, neighbours);
+		appendResults(results, t);
+		
+		first = points[i];
+	}
+	return results;
+
+}
+

src/utilities/interpolate_z.cpp

+#include <interpolate_z.h>
+#include <limits>
+#include <point.h>
+#include <edge.h>
+#include <cmath>
+
+float dist2D(const Point &p) {
+	return sqrt(p.x * p.x + p.y * p.y);
+}
+
+float interpolateZ(const Edge &e, const Point &p) {
+	const float p1Dist = dist2D(e.p1);
+	const float p2Dist = dist2D(e.p2);
+    const float pDist = dist2D(p);
+
+    if (p1Dist - p2Dist == 0) {
+        return e.p1.z;
+    }
+
+    const float m = (e.p1.z - e.p2.z) / (p1Dist - p2Dist);
+    const float b = e.p1.z - p1Dist * m;
+
+
+	return m * pDist + b;
+}

src/utilities/intersections.cpp

 #include "intersections.h"
 #include <optional>
-#include <edge.h>
 #include <triangle_edges.h>
+#include <triangle.h>
+#include <constants.h>
+#include <edge.h>
+#include <interpolate_z.h>
 
-
-
-std::optional<float> getB(std::optional<float> slope, Point p) {
-    if (!slope.has_value()) {
+std::optional<float> getB(std::optional<float> slope, Point p)
+{
+    if (!slope.has_value())
+    {
         return {};
     }
 
     return p.y - slope.value() * p.x;
 }
 
-std::optional<float> getSlope(Edge e) {
-    if (e.p2.x - e.p1.x == 0) {
+std::optional<float> getSlope(Edge e)
+{
+    if (e.p2.x - e.p1.x == 0)
+    {
         return {};
     }
-    return (e.p2.y - e.p1.y) / (e.p2.x - e.p1.x) ;
+    return (e.p2.y - e.p1.y) / (e.p2.x - e.p1.x);
+}
+
+bool withinEdge(Edge e, Point p)
+{
+    float minX = std::min(e.p1.x, e.p2.x) - EPSILON;
+    float maxX = std::max(e.p1.x, e.p2.x) + EPSILON;
+
+    float minY = std::min(e.p1.y, e.p2.y) - EPSILON;
+    float maxY = std::max(e.p1.y, e.p2.y) + EPSILON;
+
+    return minX <= p.x && p.x <= maxX && minY <= p.y && p.y <= maxY;
+}
+
+void intersection(Edge e1, Edge e2, std::vector<Point> &results)
+{
+    auto point = intersectionWithinEdge(e1, e2);
+    if (point.has_value())
+    {
+        results.push_back(point.value());
+    }
 }
 
-bool withinEdge(Edge e, Point p) {
-    float minX = std::min(e.p1.x, e.p2.x);
-    float maxX = std::max(e.p1.x, e.p2.x);
+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 {};
+}
 
-    float minY = std::min(e.p1.y, e.p2.y);
-    float maxY = std::max(e.p1.y, e.p2.y);
+// intersection with first edge being considered infinite in length
+std::optional<Point> intersectionFirstSide(const Edge &e1, const Edge &e2) {
+    auto candPoint = intersection(e1, e2);
 
-    return minX <= p.x && p.x <= maxX && minY <= p.y && p.y <= maxY; 
+    if (candPoint.has_value() && withinEdge(e2, candPoint.value())) {
+        float z = interpolateZ(e2, candPoint.value());
+        candPoint.value().z = z;
+        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() && withinEdge(e2, candPoint.value()))
+    {
+
+        // check if in direction e1 is pointing at
+        // case where e1 is pointing towards positive x axis
+        if (e1.p1.x < e1.p2.x && candPoint.value().x >= e1.p1.x)
+        {
+            return candPoint;
+        }
+
+        // case where e1 is pointing towards negative x axis
+        if (e1.p1.x > e1.p2.x && candPoint.value().x <= e1.p1.x)
+        {
+            return candPoint;
+        }
+
+        // edge case where p1 and p2 form a vertical line
+        // case where e1 is pointing towards positive y axis
+        if (e1.p1.y < e1.p2.y && candPoint.value().y >= e1.p1.y)
+        {
+            return candPoint;
+        }
+
+        // case where e1 is pointing towards negative axis
+        if (e1.p1.y > e1.p2.y && candPoint.value().y <= e1.p1.y)
+        {
+            return candPoint;
+        }
+    }
+    return {};
+}
 
-void intersection(Edge e1, Edge e2, std::vector<Point> &results) {
+std::optional<Point> intersection(const Edge &e1, const Edge &e2)
+{
     auto slope1 = getSlope(e1);
     auto slope2 = getSlope(e2);
 
     auto b1 = getB(slope1, e1.p1);
     auto b2 = getB(slope2, e2.p1);
 
-
     // ignore overlapping case
-    if (!slope1.has_value() && !slope2.has_value()) {
-        return;
+    if (!slope1.has_value() && !slope2.has_value())
+    {
+        return {};
     }
 
-    if (!slope1.has_value()) {
-        results.push_back(Point{e1.p1.x, slope2.value() * e1.p1.x + b2.value()});
-        return;
+    if (!slope1.has_value())
+    {
+        return Point{e1.p1.x, slope2.value() * e1.p1.x + b2.value()};
     }
 
-    if (!slope2.has_value()) {
-        results.push_back(Point{e2.p1.x, slope1.value() * e2.p1.x + b1.value()});
-        return;
+    if (!slope2.has_value())
+    {
+        return Point{e2.p1.x, slope1.value() * e2.p1.x + b1.value()};
     }
     float candX = (b2.value() - b1.value()) / (slope1.value() - slope2.value());
-    auto candPoint = Point{ candX, slope1.value() * candX + b1.value()};
+    auto candPoint = Point{candX, slope1.value() * candX + b1.value()};
 
-    if (withinEdge(e1, candPoint) && withinEdge(e2, candPoint)) {
-        results.push_back(candPoint);
+    // ignore case where intersection is at an endpoint
+    /*
+    if (candPoint == e1.p1 || candPoint == e1.p2 || candPoint == e2.p1 || candPoint == e2.p2) {
+        return {};
+    }*/
+    if (candPoint.x > MAX_POINT || candPoint.y > MAX_POINT || candPoint.x < MIN_POINT || candPoint.y < MIN_POINT)
+    {
+        return {};
     }
+    return candPoint;
 }
-void intersections(Edge e1, TriangleEdges te, std::vector<Point> &results) {
-   intersection(e1, te.e1, results);
-   intersection(e1, te.e2, results);
-   intersection(e1, te.e3, results); 
+void intersections(const Edge &e1, const TriangleEdges &te, std::vector<Point> &results)
+{
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        intersection(e1, te.edges[i], results);
+    }
 }
-std::vector<Point> intersections(Triangle t1, Triangle t2) {
+std::vector<Point> intersections(const Triangle &t1, const Triangle &t2)
+{
     TriangleEdges t1Edges = TriangleEdges(t1);
 
     TriangleEdges t2Edges = TriangleEdges(t2);
     std::vector<Point> results;
 
-    intersections(t1Edges.e1, t2Edges, results);
-    intersections(t1Edges.e2, t2Edges, results);
-    intersections(t1Edges.e3, t2Edges, results);
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        intersections(t1Edges.edges[i], t2Edges, results);
+    }
 
     return results;
 }
 
+std::vector<Point> intersections(const std::vector<Point> &points, const Edge &line)
+{
+    std::vector<Point> result;
+    std::vector<Edge> es = makeEdges(points);
+    for (int i = 0; i < es.size(); i++)
+    {
+        auto cand = intersectionFirstSide(line, es[i]);
+        if (cand.has_value())
+        {
+            result.push_back(cand.value());
+        }
+    }
+    return result;
+}

src/utilities/orientation.cpp

 #include "orientation.h"
 #include <point.h>
+#include <triangle.h>
+
+float edgeValue(const Point &p1, const Point &p2) {
+    return (p2.x - p1.x) * (p2.y + p1.y);
+}
 
 Orientation orientation(const Point &p1, const Point &p2, const Point &p3) {
-    int val = (p2.y - p1.y) * (p3.x - p2.x) - (p2.x - p1.x) * (p3.y - p2.y);
+    float val = edgeValue(p1, p2) + edgeValue(p2, p3) + edgeValue(p3, p1);
     if (val == 0) {
         return Collinear;
     }
 
     return (val > 0) ? Clockwise : Counterclockwise; 
+}
+
+
+Orientation orientation(const Triangle &t) {
+    return orientation(t.points[0], t.points[1], t.points[2]);
 }
\ No newline at end of file

src/utilities/pointList.cpp

@@ -4,14 +4,16 @@ bool PointList::empty() const {
     return size == 0;
 }
 
-int PointList::prev(const PointNode &p) const {
-    const int index = (p.prevIndex + size) % size;
-    return index;
+int PointList::getSize() const {
+    return size;
 }
 
-int PointList::next(const PointNode &p) const {
-    const int index = p.nextIndex % size;
-    return index;
+int PointNode::next() const {
+    return nextIndex;
+}
+
+int PointNode::prev() const {
+    return prevIndex;
 }
 
 void PointList::remove(int i) {
@@ -20,10 +22,8 @@ void PointList::remove(int i) {
         return;
     }
     const PointNode &p = points[i];
-    int prevInd = prev(p);
-    int nextInd = next(p);
-    points[prevInd].nextIndex = nextInd;
-    points[nextInd].prevIndex = prevInd;
+    points[p.prev()].nextIndex = p.next();
+    points[p.next()].prevIndex = p.prev();
     size--;
 }
 
@@ -31,4 +31,6 @@ PointList::PointList(const std::vector<Point> &points) : size {points.size()} {
     for (int i = 0; i < points.size(); i++) {
         this->points[i] = PointNode{i - 1, i + 1, points[i]};
     }
+    this->points[0].prevIndex = size - 1;
+    this->points[size - 1].nextIndex = 0; 
 }
\ No newline at end of file

src/utilities/split_triangle.cpp

+#include <split_triangle.h>
+#include <constants.h>
+#include <contourize.h>
+#include <triangle.h>
+#include <edge.h>
+#include <intersections.h>
+
+std::vector<std::vector<Point>> splitShape(const std::vector<Point> &points, const Edge &line)
+{
+    // get intersections
+    auto intr = intersections(points, line);
+
+    std::vector<Point> pos;
+    std::vector<Point> neg;
+    for (int i = 0; i < points.size(); i++)
+    {
+        if (line.positiveSide(points[i]))
+        {
+            pos.push_back(points[i]);
+        }
+        else
+        {
+            neg.push_back(points[i]);
+        }
+    }
+
+    if (!(pos.empty() || neg.empty()))
+    {
+        
+            pos.insert(pos.end(), intr.begin(), intr.end());
+            neg.insert(neg.end(), intr.begin(), intr.end());
+    }    // Build Positive side
+
+    pos = contourize(pos);
+
+    // Build Negative side
+    neg = contourize(neg);
+    return std::vector<std::vector<Point>>{pos, neg};
+}
\ No newline at end of file