Allow compute_intersections supports circular arcs

.gitignore

@@ -1,4 +1,5 @@
 bazel-*
+build/*
 .vscode
 checker/checker
 checker/logs/*

include/packingsolver/algorithms/common.hpp

@@ -2,6 +2,10 @@
 
 #define _USE_MATH_DEFINES
 
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
 #include "optimizationtools/utils/output.hpp"
 
 #include <cstdint>

include/packingsolver/irregular/shape.hpp

@@ -140,6 +140,9 @@ struct ShapeElement
     /** Length of the element. */
     LengthDbl length() const;
 
+    /** Check if a point is on the element. */
+    bool contains(const Point& point) const;
+
     ShapeElement rotate(Angle angle) const;
 
     ShapeElement axial_symmetry_identity_line() const;
@@ -219,6 +222,11 @@ struct Shape
             Angle angle = 0.0,
             bool mirror = false) const;
 
+    /** Check if the shape contains a given point. */
+    bool contains(
+            const Point& point,
+            bool strict = false) const;
+
     /* Check if the shape is connected and in anticlockwise direction. */
     bool check() const;
 

src/irregular/CMakeLists.txt

@@ -7,6 +7,7 @@ target_sources(PackingSolver_irregular PRIVATE
     algorithm_formatter.cpp
     optimize.cpp
     shape.cpp
+    shape_element_intersections.cpp
     shape_convex_hull.cpp
     shape_extract_borders.cpp
     shape_self_intersections_removal.cpp

src/irregular/shape.cpp

@@ -165,6 +165,68 @@ LengthDbl ShapeElement::length() const
     return -1;
 }
 
+bool ShapeElement::contains(const Point& point) const
+{
+    switch (type) {
+    case ShapeElementType::LineSegment: {
+        // Calculate the squared length of the line segment
+        LengthDbl line_length_squared = std::pow(end.x - start.x, 2) + std::pow(end.y - start.y, 2);
+
+        // If the line segment is actually a point
+        if (equal(line_length_squared, 0.0)) {
+            return equal(point.x, start.x) && equal(point.y, start.y);
+        }
+
+        // Calculate parameter t, representing the position of the point on the line segment (between 0 and 1 indicates on the segment)
+        LengthDbl t = ((point.x - start.x) * (end.x - start.x) + (point.y - start.y) * (end.y - start.y)) / line_length_squared;
+
+        // If t is outside the [0,1] range, the point is not on the line segment
+        if (strictly_lesser(t, 0.0) || strictly_greater(t, 1.0)) {
+            return false;
+        }
+
+        // Calculate the projection point
+        Point projection;
+        projection.x = start.x + t * (end.x - start.x);
+        projection.y = start.y + t * (end.y - start.y);
+
+        // Check if the distance to the projection point is small enough
+        LengthDbl distance_squared = std::pow(point.x - projection.x, 2) + std::pow(point.y - projection.y, 2);
+        return equal(distance_squared, 0.0);
+    } case ShapeElementType::CircularArc: {
+        // Check if point lies on circle
+        if (!equal(distance(point, this->center), distance(this->start, this->center))) {
+            return false;
+        }
+
+        // Calculate angles
+        Angle point_angle = angle_radian(point - this->center);
+        Angle start_angle = angle_radian(this->start - this->center);
+        Angle end_angle = angle_radian(this->end - this->center);
+
+        // Normalize angles
+        if (this->anticlockwise) {
+            if (end_angle <= start_angle) {
+                end_angle += 2 * M_PI;
+            }
+            if (point_angle < start_angle) {
+                point_angle += 2 * M_PI;
+            }
+            return point_angle >= start_angle && point_angle <= end_angle;
+        } else {
+            if (start_angle <= end_angle) {
+                start_angle += 2 * M_PI;
+            }
+            if (point_angle < end_angle) {
+                point_angle += 2 * M_PI;
+            }
+            return point_angle >= end_angle && point_angle <= start_angle;
+        }
+    }
+    }
+    return -1;
+}
+
 std::string ShapeElement::to_string() const
 {
     switch (type) {
@@ -537,6 +599,95 @@ std::pair<Point, Point> Shape::compute_min_max(
     return {{x_min, y_min}, {x_max, y_max}};
 }
 
+bool Shape::contains(
+        const Point& point,
+        bool strict) const
+{
+    if (this->elements.empty())
+        return false;
+
+    // First check if the point lies on any boundary.
+    for (const ShapeElement& element : this->elements)
+        if (element.contains(point))
+            return (strict)? false: true;
+
+    // Then use the ray-casting algorithm to check if the point is inside
+    int intersection_count = 0;
+    for (const ShapeElement& element: this->elements) {
+        if (element.type == ShapeElementType::LineSegment) {
+            // Handle the special case of horizontal line segments
+            if (equal(element.start.y, element.end.y)) {
+                // Horizontal line segment: if the point's y coordinate equals the segment's y coordinate,
+                // and the x coordinate is within the segment range, the point is on the segment
+                if (equal(point.y, element.start.y)
+                        && (!strictly_lesser(point.x, std::min(element.start.x, element.end.x)))
+                        && (!strictly_greater(point.x, std::max(element.start.x, element.end.x)))) {
+                    // Already checked in the previous section, no need to handle here
+                    continue;
+                }
+            }
+
+            // Standard ray-casting algorithm for line segment checking
+            // Cast a ray to the right from the point, count intersections with segments
+            bool cond1 = (strictly_greater(element.start.y, point.y) != strictly_greater(element.end.y, point.y));
+            bool cond2 = strictly_lesser(
+                    point.x,
+                    (element.end.x - element.start.x) * (point.y - element.start.y)
+                    / (element.end.y - element.start.y) + element.start.x);
+
+            if (cond1 && cond2) {
+                intersection_count++;
+            }
+        } else if (element.type == ShapeElementType::CircularArc) {
+            // Circular arc checking is more complex
+            LengthDbl dx = point.x - element.center.x;
+            LengthDbl dy = point.y - element.center.y;
+            LengthDbl distance = std::sqrt(dx * dx + dy * dy);
+
+            LengthDbl radius = std::sqrt(
+                std::pow(element.start.x - element.center.x, 2)
+                + std::pow(element.start.y - element.center.y, 2));
+
+            // If the point is inside the circle and to the left of the center, there may be intersections with a ray to the right
+            if (strictly_lesser(distance, radius) && strictly_lesser(point.x, element.center.x)) {
+                LengthDbl start_angle = angle_radian({element.start.x - element.center.x, element.start.y - element.center.y});
+                LengthDbl end_angle = angle_radian({element.end.x - element.center.x, element.end.y - element.center.y});
+
+                // Ensure angles are in the correct range
+                if (element.anticlockwise && end_angle <= start_angle) {
+                    end_angle += 2 * M_PI;
+                } else if (!element.anticlockwise && start_angle <= end_angle) {
+                    start_angle += 2 * M_PI;
+                }
+
+                // Calculate the point's line-of-sight angle (angle between the line from point to center and the horizontal)
+                LengthDbl point_angle = angle_radian({dx, dy});
+                if (strictly_lesser(point_angle, 0)) {
+                    point_angle += 2 * M_PI;  // Adjust angle to [0, 2π)
+                }
+
+                // Calculate the intersection angle of the ray to the right with the circle
+                LengthDbl ray_angle = 0; // Angle of ray to the right is 0
+
+                // Check if the ray intersects the arc
+                bool intersects_arc;
+                if (element.anticlockwise) {
+                    intersects_arc = (!strictly_lesser(ray_angle, start_angle) && !strictly_greater(ray_angle, end_angle));
+                } else {
+                    intersects_arc = (!strictly_greater(ray_angle, start_angle) && !strictly_lesser(ray_angle, end_angle));
+                }
+
+                if (intersects_arc) {
+                    intersection_count++;
+                }
+            }
+        }
+    }
+
+    // If the number of intersections is odd, the point is inside the shape
+    return (intersection_count % 2 == 1);
+}
+
 Shape& Shape::shift(
         LengthDbl x,
         LengthDbl y)