quadtree implementation

CMakeLists.txt

@@ -23,11 +23,20 @@ project(intersection)
 add_executable(
 	tests
 	tests/unit_tests/split_triangle_tests.cpp
-	src/debug_utilities/mostly_equal.cpp
 	tests/unit_tests/simple_triangle_tests.cpp
 	tests/unit_tests/point_in_triangle_tests.cpp
-	src/utilities/intersections.cpp
 	tests/union_tests.cpp
+	tests/unit_tests/orientation_tests.cpp
+	tests/unit_tests/contourize_tests.cpp
+	tests/unit_tests/intersection_tests/intersections_tests.cpp
+	tests/edge_union_cases.cpp
+	tests/unit_tests/orientation_tests.cpp
+	tests/union_neighbours_tests.cpp
+	tests/quad_tree_tests.cpp
+	tests/unit_tests/interpolate_z_tests.cpp
+	tests/unit_tests/convex_triangulation_tests.cpp
+	src/debug_utilities/mostly_equal.cpp
+	src/utilities/intersections.cpp
 	src/shapes/triangle.cpp
 	src/shapes/triangle_edges.cpp
 	src/utilities/contourize.cpp
@@ -36,20 +45,14 @@ add_executable(
 	src/shapes/point.cpp
 	src/union.cpp
 	src/utilities/pointList.cpp
-	tests/unit_tests/orientation_tests.cpp
-	tests/unit_tests/intersections_tests.cpp
 	src/data_structures/box.cpp
-	tests/edge_union_cases.cpp
-	tests/unit_tests/convex_triangulation_tests.cpp
 	src/utilities/convex_triangulation.cpp
-	tests/unit_tests/contourize_tests.cpp
 	src/debug_utilities/print_triangles.cpp
 	src/utilities/split_triangle.cpp
-	tests/unit_tests/orientation_tests.cpp
-	tests/union_neighbours_tests.cpp
 	src/debug_utilities/shift_triangle.cpp
 	src/utilities/interpolate_z.cpp
-	tests/unit_tests/interpolate_z_tests.cpp
+	src/data_structures/quad_tree.cpp
+	src/data_structures/box_edges.cpp
 )
 
 include_directories(include)

include/box.h

 #pragma once
+#include <point.h>
+#include "box_edges.h"
 
 struct Triangle;
-struct Point;
+#define BOX_NB_POINTS 4
 
+// left -> +x
+// up -> +y
 struct Box {
-    float minX, maxX, minY, maxY;
-    bool isIn(const Triangle &t) const;
-    bool isIn(const Point &p) const;
+    Point topLeft;
+    Point bottomRight;
+    Point points [BOX_NB_POINTS];
+    BoxEdges edges;
+    Box(Point topLeft, Point bottomRight);
+    bool intersects(const Triangle &t) const;
+    bool intersects(const Point &p) const;
+
+    Box firstQuadrant() const;
+    Box secondQuadrant() const;
+    Box thirdQuadrant() const;
+    Box fourthQuadrant() const;
 };
\ No newline at end of file

include/box_edges.h

+#pragma once
+#include <edge.h>
+struct Box;
+
+struct BoxEdges {
+    static const int size = 4;
+    Edge e [4];
+    BoxEdges(const Point p[size]);
+};
\ No newline at end of file

include/convex_triangulation.h

@@ -5,4 +5,4 @@
 #include <vector>
 
 // Requires a convex shape
-std::vector<Triangle> convexTriangulation(const std::vector<Point> &points, int triangleId);
+std::vector<Triangle> convexTriangulation(const std::vector<Point> &points, int triangleId, std::shared_ptr<std::vector<int>> neighbours);

include/quad_tree.h

@@ -3,22 +3,30 @@
 #include <memory>
 #include <vector>
 #include "box.h"
+#include <set>
 
 const int POINT_NOT_IN_QUADTREE = -1;
+const int QUADTREE_MAX_DEPTH = 5;
 
 #define QUADTREE_NODE_MAX_SHAPE 4
 
-class QuadTree {
+class QuadTree
+{
     Box b;
-    std::vector<Triangle> triangles; 
+    int level;
+    std::vector<Triangle> triangles;
 
     std::vector<QuadTree> children;
+    QuadTree(Box b, int level);
+    void collectUniqueTriangleFragments(const Triangle &t, std::set<int> &seen, std::vector<Triangle> &result) const;
+    void split();
+    void addNonIntersectingTriangle(const Triangle &t);
 
-    public:
-        QuadTree(Box b);
-        void addTriangle(Triangle triangle);
-        std::vector<Triangle> visibleSurface() const;
-        // returns triangle id
-        int pointIntersection(Point p) const;
-
+    std::vector<Triangle> visibleSurface(std::set<int> &seen) const;
+public:
+    QuadTree(Box b);
+    void addTriangle(const Triangle &triangle);
+    std::vector<Triangle> visibleSurface() const;
+    // returns triangle id
+    int pointIntersection(const Point &p) const;
 };
\ No newline at end of file

include/triangle.h

@@ -2,16 +2,21 @@
 #include "point.h"
 #include "print_triangle.h"
 #include <vector>
+#include <memory>
+
+// note : Triangle equality is only specified for tests
 
 // points specified in counterclockwise or collinear (line / point) order
 struct Triangle
 {
     Point points[3];
-    int id;
-    std::vector<int> neighbours;
+    int mainTriangleId;
+    int fragmentId;
+    std::shared_ptr<std::vector<int>> neighbours;
 
 
-    Triangle(const Point &p1, const Point &p2, const Point &p3, int id = 0, const 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

src/data_structures/box.cpp

@@ -2,20 +2,47 @@
 #include <triangle.h>
 #include <constants.h>
 
-bool Box::isIn(const Triangle &t) const {
+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 (x > minX && x < maxX && y > minY && y < maxY) {
+        if (intersects({x,y})) {
             return true;
         }
     }
-    return false;
-}
 
-bool Box::isIn(const Point &p) const {
-    if (p.x > minX && p.x < maxX && p.y > minY && p.y < maxY) {
-        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

@@ -19,5 +19,5 @@ 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]) &&
-		id == other.id;
+		mainTriangleId == other.mainTriangleId;
 }
\ No newline at end of file

src/debug_utilities/print_triangles.cpp

@@ -3,7 +3,7 @@
 #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.id << ")";
+	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) {

src/shapes/triangle.cpp

 #include "triangle.h"
 #include "triangle_edges.h"
 #include "edge.h"
+#include <triangle_fragment_id_assigner.h>
 
-Triangle::Triangle(const Point &p1, const Point &p2, const Point &p3, int id, std::vector<int> neighbours) : points{p1,p2,p3}, id{id}, neighbours{neighbours} {}
+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} {}
 
-bool Triangle::pointInTriangle(const Point &p) const {
+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
+{
     // all tests must be positive
     auto edges = TriangleEdges(*this);
 
-    for (int i = 0; i < NB_TRIANGLE_SIDES; i++) {
-        if (edges.edges[i].positiveSide(p)) {
+    for (int i = 0; i < NB_TRIANGLE_SIDES; i++)
+    {
+        if (edges.edges[i].positiveSide(p))
+        {
             return false;
         }
     }
     return true;
 }
 
-
-int nextPoint(int pointIndex)  {
+int nextPoint(int pointIndex)
+{
     return (pointIndex + 1) % 3;
 }

src/union.cpp

@@ -35,7 +35,7 @@ std::vector<Triangle> unionizeTopAndBottom(const Triangle &top, const Triangle &
         // add these to result
         if (!shapes[0].empty())
         {
-            std::vector<Triangle> relvTriangles = convexTriangulation(shapes[0], bot.id);
+            std::vector<Triangle> relvTriangles = convexTriangulation(shapes[0], bot.mainTriangleId, bot.neighbours);
             result.insert(result.end(), relvTriangles.begin(), relvTriangles.end());
         }
         // future relevant triangles
@@ -48,7 +48,11 @@ std::vector<Triangle> unionizeTopAndBottom(const Triangle &top, const Triangle &
 
 std::vector<Triangle> unionize(const Triangle &t1, const Triangle &t2)
 {
-	if (std::find(t1.neighbours.begin(), t1.neighbours.end(), t2.id) != t1.neighbours.end()) {
+    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 (intersections(t1, t2).empty())

src/utilities/convex_triangulation.cpp

@@ -9,13 +9,13 @@ void appendResults(std::vector<Triangle> &results, const Triangle &t) {
 	results.push_back(t);
 }
 
-std::vector<Triangle> convexTriangulation(const std::vector<Point> &points, int triangleId) {
+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);
+	const Triangle t1 = Triangle(points[0], points[1], points[2], triangleId, neighbours);
 	std::vector<Triangle> results;
 
 	appendResults(results, t1);
@@ -23,7 +23,7 @@ std::vector<Triangle> convexTriangulation(const std::vector<Point> &points, int
 	
 	Point first = points[0];
 	for (int i = size - 1; i >= 3; i--) {
-		const Triangle t = Triangle(first, points[2], points[i], triangleId);
+		const Triangle t = Triangle(first, points[2], points[i], triangleId, neighbours);
 		appendResults(results, t);
 		
 		first = points[i];

tests/quad_tree_tests.cpp

+#include <gtest/gtest.h>
+#include <quad_tree.h>
+
+TEST(QuadTreeTests, SimplePointQueryTest) {
+    QuadTree q({{10,10}, {0,0}});
+
+    Triangle t({0,0}, {5,0}, {2,3}, 42);
+    q.addTriangle(t);
+
+
+    Point p{1,1};
+    int id = q.pointIntersection(p);
+
+    EXPECT_EQ(id, t.mainTriangleId);
+}
+
+TEST (QuadTreeTests, NotInQueryPointTest) {
+
+    QuadTree q({{10,10}, {0,0}});
+    Triangle t({0,0}, {5,0}, {2,3});
+    Point p{9,9};
+    int id = q.pointIntersection(p);
+    EXPECT_EQ(id, POINT_NOT_IN_QUADTREE);
+}
+
+TEST (QuadTreeTests, OverlapTest) {
+
+    QuadTree q({{10,10}, {0,0}});
+    Triangle t1({0, 0, 0}, {5, 0, 0}, {2, 3, 0}, 1);
+    Triangle t2({3, 1, 1}, {6, 1, 1}, {4, 3, 1}, 2);
+
+    q.addTriangle(t1);
+    q.addTriangle(t2);
+
+    Point p{1,1};
+    int id = q.pointIntersection(p);
+    EXPECT_EQ(id, t1.mainTriangleId);
+}
+
+TEST (QuadTreeTests, NonOverlappingTest) {
+
+    QuadTree q({{10,10}, {0,0}});
+    Triangle t1({0, 0, 0}, {5, 0, 0}, {2, 3, 0}, 1);
+    Triangle t2{{6, 0, 1}, {10, 0, 1}, {7, 5, 1}, 2};
+
+    q.addTriangle(t1);
+    q.addTriangle(t2);
+}
+
+TEST (QuadTreeTests, VisibleSurfaceTestNonOverlapping) {
+
+    QuadTree q({{10,10}, {0,0}});
+    Triangle t1({0, 0, 0}, {5, 0, 0}, {2, 3, 0}, 1);
+    Triangle t2{{6, 0, 1}, {10, 0, 1}, {7, 5, 1}, 2};
+
+    q.addTriangle(t1);
+    q.addTriangle(t2);
+
+    auto surface = q.visibleSurface();
+
+    std::vector<Triangle> expected_surface {{{6,0,1}, {10, 0, 1}, {7, 5, 1}, 2}, {{0,0,0}, {5,0,0},{2,3,0}, 1}};
+
+    EXPECT_EQ(surface, expected_surface);
+}
+
+TEST (QuadTreeTests, NoVisibleSurfaceTest) {
+    QuadTree q({{1,1}, {0,0}});
+    auto surface = q.visibleSurface();
+    EXPECT_EQ(surface, std::vector<Triangle>{});
+}
+
+TEST (QuadTreeTests, FiveTrianglesTest) {
+
+    QuadTree q({{1,1}, {0,0}});
+
+    Triangle t1{{0,0,0}, {0,0,0}, {0,0,0}, 1};
+    Triangle t2{{0,0,0}, {0,0,0}, {0,0,0}, 1};
+    Triangle t3{{0,0,0}, {0,0,0}, {0,0,0}, 1};
+    Triangle t4{{0,0,0}, {0,0,0}, {0,0,0}, 1};
+}
+
+TEST (QuadTreeTests, DegenerateTrianglesTest) {
+    QuadTree q({{1,1}, {0,0}});
+    Triangle t1{{0,0,0}, {0,0,0}, {0,0,0}, 1};
+    Triangle t2{{0.2, 0.1, 1}, {0.2, 0.1, 1}, {0.3, 0.2, 1}, 2};
+    Triangle t3{{1,0,0}, {1,0,0}, {0,0,0}, 3};
+
+    q.addTriangle(t3);
+    q.addTriangle(t2);
+    q.addTriangle(t1);
+
+    auto surface = q.visibleSurface();
+
+    std::vector<Triangle> expected_surface {t1, t3, t2};
+
+    EXPECT_EQ(surface, expected_surface);
+}
\ No newline at end of file

tests/unit_tests/convex_triangulation_tests.cpp

@@ -14,7 +14,7 @@ TEST(ConvexTriangulationTests, Simple) {
 		p1, p2, p3, p4, p5
 	};	
 
-	auto results = convexTriangulation(points, 2);
+	auto results = convexTriangulation(points, 2, std::make_shared<std::vector<int>>());
 
 	std::vector<Triangle> expected_result {
 		Triangle(p1, p2, p3, 2),
@@ -33,7 +33,7 @@ TEST(ConvexTriangulationTests, TrivialTriangle) {
 		{1.5, 2}
 	};
 
-	auto results = convexTriangulation(points, 1);
+	auto results = convexTriangulation(points, 1, std::make_shared<std::vector<int>>());
 	Triangle expected_triangle = Triangle({0,0}, {2,0}, {1.5, 2}, 1);
 
 	EXPECT_EQ(results.size(), 1);