From 197c4e82bf4d4cf4cc98825e54cece18ec1acfd0 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 17:03:07 +0200
Subject: [PATCH 01/16] [doc] fix Doxygen API copy during manual build

---
 doc/CMakeLists.txt | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt
index 27eb92267..bce6a0a88 100644
--- a/doc/CMakeLists.txt
+++ b/doc/CMakeLists.txt
@@ -102,13 +102,15 @@
                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
                       COMMENT "Generating the manual website using Sphinx")
 
-    add_custom_command(
-      TARGET manual POST_BUILD
-      COMMAND ${CMAKE_COMMAND} -E copy_directory
-              ${CMAKE_BINARY_DIR}/doc/api/html
-              ${SPHINX_EXTRA_API}
-      COMMENT "Copying generated Doxygen HTML to Sphinx extra_html/api"
-    )
+    if(DOXYGEN_FOUND)
+      add_custom_command(
+        TARGET manual POST_BUILD
+        COMMAND ${CMAKE_COMMAND} -E copy_directory
+                ${CMAKE_CURRENT_BINARY_DIR}/api/html
+                ${SPHINX_BUILD}/extra_html/api
+        COMMENT "Copying generated Doxygen HTML to Sphinx extra_html/api"
+      )
+    endif()
     
     install(DIRECTORY ${SPHINX_BUILD}/
             DESTINATION share/codac_website

From 7f315a5b516e271dfe91bd57d3a24a328d496f3e Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 17:08:56 +0200
Subject: [PATCH 02/16] [doc] document tube_eval with scalar and vector tube
 examples

---
 .../functions/analytic/analytic_functions.rst |  4 +-
 doc/manual/manual/tubes/slicedtube.rst        | 38 +++++++++++++++++++
 doc/manual/manual/tubes/src.cpp               | 21 ++++++++++
 doc/manual/manual/tubes/src.py                | 19 ++++++++++
 4 files changed, 81 insertions(+), 1 deletion(-)

diff --git a/doc/manual/manual/functions/analytic/analytic_functions.rst b/doc/manual/manual/functions/analytic/analytic_functions.rst
index 60fb121fd..50555c8b5 100644
--- a/doc/manual/manual/functions/analytic/analytic_functions.rst
+++ b/doc/manual/manual/functions/analytic/analytic_functions.rst
@@ -25,7 +25,7 @@ This page provides an overview of the ``AnalyticFunction`` class, its key featur
   For defining a contractor based on an ``AnalyticFunction``, the ``CtcInverse`` class is available.
   It allows to solve constraints under the form :math:`\mathbf{f}(\mathbf{x})\in[\mathbf{y}]`.
 
-  See more.
+  :ref:`See more <sec-ctc-analytic-ctcinverse>`.
 
 
 Variables
@@ -226,6 +226,8 @@ Evaluations can be done by the ``.eval()`` method of ``AnalyticFunction``.
 
 For multivariate functions, it is possible to provide several inputs to the ``.eval()`` method, as in the above example with ``f2`` or ``f3``. In Codac C++, the arguments may be of different types (any mixed scalar, vector, or matrix types). In Python and Matlab however, the binding does not allow any combination of inputs, but several arguments of same structural types can be provided: a list of scalar values, or a list of vector values (``Vector``, ``IntervalVector``), or a list of matrix values (``Matrix``, ``IntervalMatrix``).
 
+Functions can also be evaluated on time-varying inputs represented by ``SlicedTube`` objects with ``f.tube_eval(x1, x2, ...)``. This produces a new tube by applying the interval evaluation to each temporal slice. For more information, see :ref:`sec-domains-tubes-slicedtube-function-eval`.
+
 
 Evaluation modes
 ^^^^^^^^^^^^^^^^
diff --git a/doc/manual/manual/tubes/slicedtube.rst b/doc/manual/manual/tubes/slicedtube.rst
index 87202c2ae..01c21cbac 100644
--- a/doc/manual/manual/tubes/slicedtube.rst
+++ b/doc/manual/manual/tubes/slicedtube.rst
@@ -274,6 +274,44 @@ For vector-valued codomains, pointwise arithmetic is also available. In particul
       :dedent: 4
 
 
+.. _sec-domains-tubes-slicedtube-function-eval:
+
+Evaluating analytic functions on tubes
+--------------------------------------
+
+An :class:`~codac.AnalyticFunction` can be evaluated on one or more sliced tubes using ``f.tube_eval(x1, x2, ...)``. The result is a new tube computed by applying the function to the codomains of each temporal slice:
+
+.. math::
+
+  [y]_k = [f]([x_1]_k, [x_2]_k, \ldots, [x_n]_k).
+
+The result encloses :math:`f(x_1(t),\ldots,x_n(t))` at every time :math:`t`.
+It shares the same ``TDomain`` as the input tubes, including its gates. All input tubes must share the same ``TDomain`` object, not merely have identical time bounds or discretizations. If the function is scalar-valued, ``tube_eval()`` returns a scalar tube.  If it is vector-valued, ``tube_eval()`` returns a vector tube. In Python, the number of input tubes supported by ``tube_eval()`` is currently limited. All input tubes must be of the same kind (scalar or vector), whereas C++ also supports mixed types.
+
+In the following example, a two-argument function maps heading and speed tubes to a 2D velocity tube :math:`(v_x,v_y)=(v\cos\theta,v\sin\theta)`. A second function computes its squared Euclidean norm :math:`v_x^2+v_y^2`, returning a scalar tube:
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [slicedtube-class-function-eval-beg]
+      :end-before: [slicedtube-class-function-eval-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [slicedtube-class-function-eval-beg]
+      :end-before: [slicedtube-class-function-eval-end]
+      :dedent: 4
+
+Unlike ``x(t)``, which evaluates a tube at a given time, ``f.tube_eval(x)`` evaluates an analytic function on an entire tube. It does not implicitly pass time to ``f``: if time is an argument, it must be supplied as an input tube.
+For functions of time alone, ``SlicedTube(tdomain, f)`` directly evaluates ``f(t)`` over the temporal slices (see `Creating sliced tubes`_).
+
+
 Integration and primitives
 --------------------------
 
diff --git a/doc/manual/manual/tubes/src.cpp b/doc/manual/manual/tubes/src.cpp
index 8e3a932f3..0eb2c62af 100644
--- a/doc/manual/manual/tubes/src.cpp
+++ b/doc/manual/manual/tubes/src.cpp
@@ -222,6 +222,27 @@ TEST_CASE("SlicedTube class - manual")
     // [slicedtube-class-7-end]
   }
 
+  {
+    // [slicedtube-class-function-eval-beg]
+    auto td = create_tdomain({0,5}, 0.1, true);
+
+    ScalarVar t;
+    SlicedTube heading(td, AnalyticFunction({t}, 0.2*sin(t)));
+    SlicedTube speed(td, Interval(1,2));
+
+    ScalarVar h, v;
+    AnalyticFunction f({h,v}, {v*cos(h), v*sin(h)});
+    auto velocity = f.tube_eval(heading,speed); // SlicedTube<IntervalVector>
+
+    VectorVar p(2);
+    AnalyticFunction g({p}, sqr(p[0])+sqr(p[1]));
+    auto velocity_norm_squared = g.tube_eval(velocity); // scalar tube: ||velocity||^2
+    // [slicedtube-class-function-eval-end]
+
+    CHECK(velocity.size() == 2);
+    CHECK(velocity_norm_squared.size() == 1);
+  }
+
   {
     auto td = create_tdomain({0,2});
     SlicedTube x(td, Interval());
diff --git a/doc/manual/manual/tubes/src.py b/doc/manual/manual/tubes/src.py
index 23d27c3e6..6d845dff8 100644
--- a/doc/manual/manual/tubes/src.py
+++ b/doc/manual/manual/tubes/src.py
@@ -187,6 +187,25 @@ def tests_SlicedTube_manual(test):
     y = A * b
     # [slicedtube-class-7-end]
 
+    # [slicedtube-class-function-eval-beg]
+    td = create_tdomain([0,5], 0.1, True)
+
+    t = ScalarVar()
+    heading = SlicedTube(td, AnalyticFunction([t], 0.2*sin(t)))
+    speed = SlicedTube(td, Interval(1,2))
+
+    h, v = ScalarVar(), ScalarVar()
+    f = AnalyticFunction([h,v], [v*cos(h), v*sin(h)])
+    velocity = f.tube_eval(heading,speed) # vector-valued tube
+
+    p = VectorVar(2)
+    g = AnalyticFunction([p], sqr(p[0])+sqr(p[1]))
+    velocity_norm_squared = g.tube_eval(velocity) # scalar tube: ||velocity||^2
+    # [slicedtube-class-function-eval-end]
+
+    test.assertTrue(velocity.size() == 2)
+    test.assertTrue(velocity_norm_squared.size() == 1)
+
     td = create_tdomain([0,2])
     x = SlicedTube(td, Interval())
 

From 9e251393d307f3a24a0ac8fb875cb2682329a350 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 17:11:30 +0200
Subject: [PATCH 03/16] [geom] Fix segment intersection with uncertain
 parallelism

---
 src/core/geometry/codac2_Segment.cpp         | 17 ++++++------
 tests/core/geometry/codac2_tests_Segment.cpp | 29 ++++++++++++++++++++
 tests/core/geometry/codac2_tests_Segment.py  | 25 +++++++++++++++++
 3 files changed, 63 insertions(+), 8 deletions(-)

diff --git a/src/core/geometry/codac2_Segment.cpp b/src/core/geometry/codac2_Segment.cpp
index 5aa7826c0..53e31d620 100644
--- a/src/core/geometry/codac2_Segment.cpp
+++ b/src/core/geometry/codac2_Segment.cpp
@@ -116,6 +116,9 @@ namespace codac2
   {
     IntervalVector b1 = e1.box(), b2 = e2.box();
 
+    if(e1.intersects(e2) == BoolInterval::FALSE)
+      return IntervalVector::empty(2);
+
     // Ill-conditioned cases
     // For tiny boxes, proj_intersection returns an overly pessimistic enclosure.
     {
@@ -169,15 +172,13 @@ namespace codac2
 
     if(c.contains(0.)) // e1 and e2 may be parallel
     {
-      // They are then either possibly colinear
-      //   => infinite intersection points
-      // Or not
-      //   => no intersection point => empty output
-
-      if((aligned(e1[0],e1[1],e2[0]) & BoolInterval::TRUE) == BoolInterval::TRUE)
-        return IntervalVector(2);
-      else
+      // If c is exactly zero, the supporting lines are certainly parallel.
+      // They are disjoint only when colinearity can be ruled out.
+      if(c == 0. && aligned(e1[0],e1[1],e2[0]) == BoolInterval::FALSE)
         return IntervalVector::empty(2);
+
+      // Otherwise, a non-parallel realization remains possible.
+      return IntervalVector(2);
     }
 
     Interval a = (x1*y2-y1*x2);
diff --git a/tests/core/geometry/codac2_tests_Segment.cpp b/tests/core/geometry/codac2_tests_Segment.cpp
index 926342d6e..ae13e99e6 100644
--- a/tests/core/geometry/codac2_tests_Segment.cpp
+++ b/tests/core/geometry/codac2_tests_Segment.cpp
@@ -137,6 +137,35 @@ TEST_CASE("Segment")
     CHECK(proj_intersection(Segment({2,0},{6,4}), Segment({6,5},{5,6})) == IntervalVector({6.5,4.5}));
   }
 
+  SECTION("intersection - uncertain parallelism")
+  {
+    Segment e1(
+      {0.,0.},
+      IntervalVector({1.,Interval(-1.,1.)})
+    );
+    Segment e2({0.,0.5},{1.,0.5});
+
+    CHECK(!proj_intersection(e1,e2).is_empty());
+    CHECK(e1.intersects(e2) == BoolInterval::UNKNOWN);
+
+    auto x = e1 & e2;
+    CHECK(!x.is_empty());
+    CHECK(x.is_superset(IntervalVector({0.5,0.5})));
+  }
+
+  SECTION("intersection - uncertain parallel but disjoint segments")
+  {
+    Segment e1(
+      {0.,0.},
+      IntervalVector({2.,Interval(2.-1e-10,2.+1e-10)})
+    );
+    Segment e2({0.,1.},{2.,3.});
+
+    CHECK(!proj_intersection(e1,e2).is_empty());
+    CHECK(e1.intersects(e2) == BoolInterval::FALSE);
+    CHECK((e1 & e2).is_empty());
+  }
+
   SECTION("intersection - near infinite cases")
   {
     Segment e1({-1,6},{-1,next_float(-oo)});
diff --git a/tests/core/geometry/codac2_tests_Segment.py b/tests/core/geometry/codac2_tests_Segment.py
index 3b9e53114..45c1f10b5 100644
--- a/tests/core/geometry/codac2_tests_Segment.py
+++ b/tests/core/geometry/codac2_tests_Segment.py
@@ -129,6 +129,31 @@ def test_intersection(self):
     self.assertTrue((Segment([2,0],[6,4]) & Segment([6,5],[5,6])) == IntervalVector.empty(2))
     self.assertTrue(proj_intersection(Segment([2,0],[6,4]), Segment([6,5],[5,6])) == IntervalVector([6.5,4.5]))
 
+    # Possible intersection under uncertain parallelism
+    e1 = Segment(
+      [0.,0.],
+      IntervalVector([1.,Interval(-1.,1.)])
+    )
+    e2 = Segment([0.,0.5],[1.,0.5])
+
+    self.assertFalse(proj_intersection(e1,e2).is_empty())
+    self.assertTrue(e1.intersects(e2) == BoolInterval.UNKNOWN)
+
+    x = e1 & e2
+    self.assertFalse(x.is_empty())
+    self.assertTrue(x.is_superset(IntervalVector([0.5,0.5])))
+
+    # Nearly parallel supporting lines, but certainly disjoint segments
+    e1 = Segment(
+      [0.,0.],
+      IntervalVector([2.,Interval(2.-1e-10,2.+1e-10)])
+    )
+    e2 = Segment([0.,1.],[2.,3.])
+
+    self.assertFalse(proj_intersection(e1,e2).is_empty())
+    self.assertTrue(e1.intersects(e2) == BoolInterval.FALSE)
+    self.assertTrue((e1 & e2).is_empty())
+
     # Near infinite cases
     e1 = Segment([-1,6],[-1,next_float(-oo)])
     e2 = Segment([-1,-1],[3,-6])

From 185aad0f5178c1d3e71924478ac23956c873674d Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 17:50:22 +0200
Subject: [PATCH 04/16] [cmake] Fix unnecessary Python binding rebuilds on
 CMake reconfiguration

---
 python/CMakeLists.txt | 88 +++++++++++++++++++++++--------------------
 1 file changed, 48 insertions(+), 40 deletions(-)

diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index e9272e5e4..609e5805f 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -16,31 +16,6 @@
   endif()
   # Adds pybind11::headers, pybind11::module, pybind11::embed
 
-  # The Python interpreter, which this directory, tests/ and examples/ run
-  # through PYTHON_EXECUTABLE: to generate the docstrings below, to build the
-  # wheel, and to run the Python tests and examples. pybind11 sets it where this
-  # directory sees it only with PYBIND11_FINDPYTHON=OFF, through
-  # FindPythonLibsNew.cmake. By default (COMPAT, or NEW) it uses FindPython,
-  # which sets Python_EXECUTABLE (or Python3_EXECUTABLE), and sets
-  # PYTHON_EXECUTABLE, in COMPAT, only within its own directory: here it stayed
-  # empty, execute_process() ran doxygen2docstring.py as a program, and the
-  # configuration stopped on "permission denied". It is initialized here with
-  # the interpreter pybind11 builds the modules for, unless it was given, and
-  # cached, so that tests/ and examples/, configured after this directory, get
-  # it too.
-  if(NOT PYTHON_EXECUTABLE)
-    if(Python_EXECUTABLE)
-      set(_codac_python_executable "${Python_EXECUTABLE}")
-    elseif(Python3_EXECUTABLE)
-      set(_codac_python_executable "${Python3_EXECUTABLE}")
-    else()
-      find_package(Python COMPONENTS Interpreter REQUIRED)
-      set(_codac_python_executable "${Python_EXECUTABLE}")
-    endif()
-    set(PYTHON_EXECUTABLE "${_codac_python_executable}" CACHE FILEPATH "Python interpreter" FORCE)
-  endif()
-  message(STATUS "Python interpreter: ${PYTHON_EXECUTABLE}")
-
   set(PYTHON_PACKAGE_NAME ${PROJECT_NAME})
   set(PYTHON_PACKAGE_DIR "${CMAKE_CURRENT_BINARY_DIR}/python_package")
 
@@ -83,18 +58,42 @@
 
     set(NO_DOCSTRING_UPDATE OFF CACHE BOOL "Disables the update of docstrings")
 
-    file(REMOVE_RECURSE ${CMAKE_CURRENT_BINARY_DIR}/docstring)
+    # Regenerate docstrings without changing unchanged files' timestamps.
+    set(_docstring_dir "${CMAKE_CURRENT_BINARY_DIR}/docstring")
+    set(_docstring_tmp "${CMAKE_CURRENT_BINARY_DIR}/docstring_tmp")
+    file(REMOVE_RECURSE "${_docstring_tmp}")
+    file(MAKE_DIRECTORY "${_docstring_dir}" "${_docstring_tmp}")
     message(STATUS "Generating docstrings from doxygen files...")
-    file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/docstring)
     execute_process(COMMAND ${PYTHON_EXECUTABLE}
       ${PROJECT_SOURCE_DIR}/scripts/doxygen/doxygen2docstring.py
       ${CMAKE_CURRENT_BINARY_DIR}/../doc/api/xml/
-      ${CMAKE_CURRENT_BINARY_DIR}/docstring
-      RESULT_VARIABLE DOCSTRING_RESULT)
-    if(NOT DOCSTRING_RESULT EQUAL 0)
-      message(FATAL_ERROR "doxygen2docstring.py failed (exit code ${DOCSTRING_RESULT}) while generating the Python docstring headers.")
+      "${_docstring_tmp}"
+      RESULT_VARIABLE _docstring_status)
+    if(NOT _docstring_status EQUAL 0)
+      message(FATAL_ERROR "Failed to generate Python docstrings (exit code ${_docstring_status})")
     endif()
 
+    file(GLOB _new_docstrings "${_docstring_tmp}/*")
+    foreach(_file IN LISTS _new_docstrings)
+      get_filename_component(_filename "${_file}" NAME)
+      execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different
+        "${_file}" "${_docstring_dir}/${_filename}"
+        RESULT_VARIABLE _copy_status)
+      if(NOT _copy_status EQUAL 0)
+        message(FATAL_ERROR "Failed to update Python docstring: ${_filename}")
+      endif()
+    endforeach()
+
+    # Remove headers that the generator no longer produces.
+    file(GLOB _old_docstrings "${_docstring_dir}/*")
+    foreach(_file IN LISTS _old_docstrings)
+      get_filename_component(_filename "${_file}" NAME)
+      if(NOT EXISTS "${_docstring_tmp}/${_filename}")
+        file(REMOVE "${_file}")
+      endif()
+    endforeach()
+    file(REMOVE_RECURSE "${_docstring_tmp}")
+
   endif()
 
 
@@ -124,14 +123,23 @@
       "${CMAKE_CURRENT_SOURCE_DIR}/../tests/codac2_tests_*.py"
       )
 
-    foreach(TEST_FILE ${TESTS_FILES})
-      string(REPLACE "codac2_tests_" "test_" TEST_FILE_UNITTEST ${TEST_FILE}) # unittest convention
-      get_filename_component(TEST_FILE_UNITTEST ${TEST_FILE_UNITTEST} NAME_WE)
-      configure_file(${TEST_FILE} ${PYTHON_PACKAGE_DIR}/${PYTHON_PACKAGE_NAME}/tests/${TEST_FILE_UNITTEST}.py COPYONLY)
-      #get_filename_component(TEST_NAME ${TEST_FILE} NAME_WE)
-      #set(TEST_NAME "python_${TEST_NAME}")
-      #add_test(NAME ${TEST_NAME} COMMAND ${PYTHON_EXECUTABLE} ${TEST_FILE})
-      #set_tests_properties(${TEST_NAME} PROPERTIES ENVIRONMENT "PYTHONPATH=${PYTHON_PACKAGE_DIR}")
+    # Copy Python tests at build time instead of reconfiguring CMake
+    # whenever an individual test file changes.
+    set(_python_test_outputs)
+    foreach(TEST_FILE IN LISTS TESTS_FILES)
+      string(REPLACE "codac2_tests_" "test_" TEST_FILE_UNITTEST "${TEST_FILE}") # unittest convention
+      get_filename_component(TEST_FILE_UNITTEST "${TEST_FILE_UNITTEST}" NAME_WE)
+      set(_test_output "${PYTHON_PACKAGE_DIR}/${PYTHON_PACKAGE_NAME}/tests/${TEST_FILE_UNITTEST}.py")
+      add_custom_command(
+        OUTPUT "${_test_output}"
+        COMMAND ${CMAKE_COMMAND} -E make_directory "${PYTHON_PACKAGE_DIR}/${PYTHON_PACKAGE_NAME}/tests"
+        COMMAND ${CMAKE_COMMAND} -E copy "${TEST_FILE}" "${_test_output}"
+        DEPENDS "${TEST_FILE}"
+        VERBATIM
+      )
+      list(APPEND _python_test_outputs "${_test_output}")
     endforeach()
+    add_custom_target(codac_python_test_files ALL DEPENDS ${_python_test_outputs})
+    add_dependencies(pip_package codac_python_test_files)
 
-  endif()
\ No newline at end of file
+  endif()

From 73f2403460d4fde08a234012c140a4c8180effeb Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:01:31 +0200
Subject: [PATCH 05/16] [geom] Handle vertex-aligned rays in polygon
 containment

---
 src/core/geometry/codac2_Polygon.cpp         | 32 +++++++++++++++-----
 tests/core/geometry/codac2_tests_Polygon.cpp |  8 +++++
 tests/core/geometry/codac2_tests_Polygon.py  |  7 +++++
 3 files changed, 39 insertions(+), 8 deletions(-)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index a7f909549..747cb0c84 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -215,18 +215,20 @@ namespace codac2
     // Odd number of crossing => point is inside
     // Even number of crossing => point is outside
 
+    const auto polygon_vertices = vertices();
     Index i = 0;
     // Some limit cases must be considered, for instance when the
-    // transect crosses exactly one vertex of the polygon, and/or 
+    // transect crosses exactly one vertex of the polygon, and/or
     // when the point to be tested is this very vertex. Below, we
-    // generate a convenient transect avoiding such configuration.
+    // first try the four axis-aligned transects, then oblique ones
+    // if necessary.
 
     do
     {
       retry = false;
 
       // Creating a ray candidate:
-      Segment try_transect = [i,&p]() -> Segment {
+      Segment try_transect = [i,&p,&box_hull,&polygon_vertices]() -> Segment {
         switch(i)
         {
           case 0:
@@ -239,14 +241,29 @@ namespace codac2
             return { {p[0],prev_float(oo)}, p };
           default:
           {
-            assert_release_unreachable("failed to test if the point is contained in polygon");
-            return Segment(IntervalVector(2),IntervalVector(2));
+            const Index k = i-4;
+            assert_release(k <= static_cast<Index>(polygon_vertices.size()) &&
+              "failed to find a vertex-free transect for polygon containment");
+            assert_release(
+              !box_hull.is_unbounded() &&
+              "failed to find a vertex-free transect for an unbounded polygon");
+
+            double span = std::max(1., box_hull.max_diam());
+            double x_out = prev_float(box_hull[0].lb()-span);
+            double dx = p[0].mid()-x_out;
+            double alpha =
+              -1. + 2.*static_cast<double>(k+1) /
+                    static_cast<double>(polygon_vertices.size()+2);
+            double y_out = p[1].mid()+alpha*dx;
+
+            return { {x_out,y_out}, p };
           }
         }
       }();
 
-      // The ray may pass through the vertices, we must double counting
-      for(const auto& pi : vertices())
+      // The ray may pass through a vertex, which would make crossing
+      // counting ambiguous.
+      for(const auto& pi : polygon_vertices)
         if(try_transect.contains(pi) != BoolInterval::FALSE)
         {
           retry = true;
@@ -257,7 +274,6 @@ namespace codac2
         transect = try_transect;
 
       i++;
-      assert(i < 10);
     } while(retry);
 
     // Now the number i of crossing can be computed.
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index dae41c7b3..296a41ced 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -132,4 +132,12 @@ TEST_CASE("Polygon contains - limit cases")
     CHECK(i == ConvexPolygon(IntervalVector({4.5,4.25})));
     CHECK(i.vertices().size() == 1);
   }
+
+  { // All four axis-aligned transects cross a vertex
+    Polygon p({{0,1},{1,0},{0,-1},{-1,0}});
+    CHECK(p.contains({0,0}) == BoolInterval::TRUE);
+
+    Polygon p_rev({{-1,0},{0,-1},{1,0},{0,1}});
+    CHECK(p_rev.contains({0,0}) == BoolInterval::TRUE);
+  }
 }
\ No newline at end of file
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index fd41c9b76..96279ae4e 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -116,5 +116,12 @@ def test_Polygon_limit_case(self):
     self.assertTrue(i == ConvexPolygon(IntervalVector([4.5,4.25])))
     self.assertTrue(len(i.vertices()) == 1)
 
+    # All four axis-aligned transects cross a vertex.
+    p = Polygon([[0,1],[1,0],[0,-1],[-1,0]])
+    self.assertTrue(p.contains([0,0]) == BoolInterval.TRUE)
+
+    p_rev = Polygon([[-1,0],[0,-1],[1,0],[0,1]])
+    self.assertTrue(p_rev.contains([0,0]) == BoolInterval.TRUE)
+
 if __name__ ==  '__main__':
   unittest.main()
\ No newline at end of file

From 6eb8ff19f6b7b85680957398be0e8ffa1da78d9a Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:09:59 +0200
Subject: [PATCH 06/16] [geom] Handle duplicate polygon vertices

---
 src/core/geometry/codac2_Polygon.cpp         | 16 ++++++++++++----
 tests/core/geometry/codac2_tests_Polygon.cpp | 10 ++++++++++
 tests/core/geometry/codac2_tests_Polygon.py  |  8 ++++++++
 3 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index 747cb0c84..e76037538 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -7,6 +7,7 @@
  *  \license    GNU Lesser General Public License (LGPL)
  */
 
+#include <algorithm>
 #include <iostream>
 #include "codac2_Polygon.h"
 #include "codac2_geometry.h"
@@ -35,10 +36,16 @@ namespace codac2
 
   Polygon::Polygon(const std::vector<IntervalVector>& v)
     : std::vector<Segment>(
-      [v]
+      [v = std::vector<IntervalVector>(v)]() mutable
       {
         vector<Segment> edges;
 
+        // Remove redundant consecutive vertices before simplifying
+        // colinear points.
+        v.erase(unique(v.begin(),v.end()),v.end());
+        if(v.size() > 1 && v.front() == v.back())
+          v.pop_back();
+
         switch(v.size())
         {
           case 0:
@@ -244,16 +251,17 @@ namespace codac2
             const Index k = i-4;
             assert_release(k <= static_cast<Index>(polygon_vertices.size()) &&
               "failed to find a vertex-free transect for polygon containment");
+
             assert_release(
               !box_hull.is_unbounded() &&
               "failed to find a vertex-free transect for an unbounded polygon");
 
-            double span = std::max(1., box_hull.max_diam());
+            double span = std::max(1.,box_hull.max_diam());
             double x_out = prev_float(box_hull[0].lb()-span);
             double dx = p[0].mid()-x_out;
             double alpha =
-              -1. + 2.*static_cast<double>(k+1) /
-                    static_cast<double>(polygon_vertices.size()+2);
+              -1. + 2.*static_cast<double>(k+1)/
+                static_cast<double>(polygon_vertices.size()+2);
             double y_out = p[1].mid()+alpha*dx;
 
             return { {x_out,y_out}, p };
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index 296a41ced..4f0104d18 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -29,6 +29,16 @@ TEST_CASE("Polygon base")
   CHECK(p3 == Polygon({{1,2},{5,1},{3,4}}));
   CHECK(p3 == Polygon({{5,1},{1,2},{3,4}}));
   CHECK(p3 == Polygon({{1,2},{3,4},{5,1}}));
+
+  { // Consecutive duplicate vertices
+    Polygon p({{0,0},{1,0},{1,0},{2,0},{2,1},{0,1}});
+    CHECK(p == Polygon({{0,0},{2,0},{2,1},{0,1}}));
+  }
+
+  { // Repeated duplicates and explicit closing vertex
+    Polygon p({{0,0},{2,0},{2,0},{2,0},{2,1},{0,1},{0,0}});
+    CHECK(p == Polygon({{0,0},{2,0},{2,1},{0,1}}));
+  }
 }
 
 TEST_CASE("Empty polygon")
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index 96279ae4e..c78624011 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -27,6 +27,14 @@ def test_polygon_base(self):
     self.assertTrue(p3 == Polygon([[5,1],[1,2],[3,4]]))
     self.assertTrue(p3 == Polygon([[1,2],[3,4],[5,1]]))
 
+    # Consecutive duplicate vertices.
+    p = Polygon([[0,0],[1,0],[1,0],[2,0],[2,1],[0,1]])
+    self.assertTrue(p == Polygon([[0,0],[2,0],[2,1],[0,1]]))
+
+    # Repeated duplicates and explicit closing vertex.
+    p = Polygon([[0,0],[2,0],[2,0],[2,0],[2,1],[0,1],[0,0]])
+    self.assertTrue(p == Polygon([[0,0],[2,0],[2,1],[0,1]]))
+
   def test_empty_polygon(self):
 
     p1 = Polygon.empty()

From 1ac9892aa8d534f8470d0f08f41cc5abbe6eb145 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:36:20 +0200
Subject: [PATCH 07/16] [geom] Fix orientation convention

---
 doc/manual/manual/geometry/src.cpp            |  8 ++++----
 doc/manual/manual/geometry/src.py             |  8 ++++----
 src/core/geometry/codac2_geometry.cpp         |  9 +++++----
 tests/core/geometry/codac2_tests_geometry.cpp | 18 +++++++++---------
 tests/core/geometry/codac2_tests_geometry.py  | 18 +++++++++---------
 5 files changed, 31 insertions(+), 30 deletions(-)

diff --git a/doc/manual/manual/geometry/src.cpp b/doc/manual/manual/geometry/src.cpp
index 88295740a..376bde3c0 100644
--- a/doc/manual/manual/geometry/src.cpp
+++ b/doc/manual/manual/geometry/src.cpp
@@ -27,15 +27,15 @@ TEST_CASE("geometry - manual")
 
     IntervalVector p1({0,0}), p2({1,1}), p3({-1,2});
     OrientationInterval c = orientation(p1,p2,p3);
-    // c == OrientationInterval::CLOCKWISE
+    // c == OrientationInterval::COUNTERCLOCKWISE
     OrientationInterval d = orientation(p3,p2,p1);
-    // d == OrientationInterval::COUNTERCLOCKWISE
+    // d == OrientationInterval::CLOCKWISE
     // [geometry-1-end]
     
     CHECK(a == OrientationInterval::COLINEAR);
     CHECK(b == OrientationInterval::UNKNOWN);
-    CHECK(c == OrientationInterval::CLOCKWISE);
-    CHECK(d == OrientationInterval::COUNTERCLOCKWISE);
+    CHECK(c == OrientationInterval::COUNTERCLOCKWISE);
+    CHECK(d == OrientationInterval::CLOCKWISE);
   }
 
   {
diff --git a/doc/manual/manual/geometry/src.py b/doc/manual/manual/geometry/src.py
index 609631e80..c69d91b6b 100644
--- a/doc/manual/manual/geometry/src.py
+++ b/doc/manual/manual/geometry/src.py
@@ -24,14 +24,14 @@ def tests_geometry_manual(test):
 
     p1, p2, p3 = IntervalVector([0,0]), IntervalVector([1,1]), IntervalVector([-1,2])
     c = orientation(p1,p2,p3)
-    # c == OrientationInterval.CLOCKWISE
+    # c == OrientationInterval.COUNTERCLOCKWISE
     d = orientation(p3,p2,p1)
-    # d == OrientationInterval.COUNTERCLOCKWISE
+    # d == OrientationInterval.CLOCKWISE
     # [geometry-1-end]
     test.assertTrue(a == OrientationInterval.COLINEAR)
     test.assertTrue(b == OrientationInterval.UNKNOWN)
-    test.assertTrue(c == OrientationInterval.CLOCKWISE)
-    test.assertTrue(d == OrientationInterval.COUNTERCLOCKWISE)
+    test.assertTrue(c == OrientationInterval.COUNTERCLOCKWISE)
+    test.assertTrue(d == OrientationInterval.CLOCKWISE)
 
     # [geometry-2-beg]
     a = aligned([0,0],[1,1],[10,10])
diff --git a/src/core/geometry/codac2_geometry.cpp b/src/core/geometry/codac2_geometry.cpp
index 66b8882cd..c299795cb 100644
--- a/src/core/geometry/codac2_geometry.cpp
+++ b/src/core/geometry/codac2_geometry.cpp
@@ -19,7 +19,7 @@ namespace codac2
   {
     assert_release(p1.size() == 2 && p2.size() == 2 && p3.size() == 2);
 
-    Interval o = (p2[1]-p1[1])*(p3[0]-p2[0])-(p2[0]-p1[0])*(p3[1]-p2[1]); 
+    Interval o = (p2[0]-p1[0])*(p3[1]-p1[1])-(p2[1]-p1[1])*(p3[0]-p1[0]); 
 
     if(o.is_empty())
       return OrientationInterval::EMPTY;
@@ -69,7 +69,8 @@ namespace codac2
       });
     pts.erase(unique(pts.begin(), pts.end()), pts.end());
 
-    if(pts.size() < 3) return pts;
+    if(pts.size() < 3)
+      return pts;
 
     // Implementation of a Graham scan method.
     // Based on some sources from OpenGenus Foundation.
@@ -104,7 +105,7 @@ namespace codac2
         {
           auto o = orientation(p1, p0, p2);
 
-          if(o == OrientationInterval::COUNTERCLOCKWISE)
+          if(o == OrientationInterval::CLOCKWISE)
             return true;
 
           else if((o == OrientationInterval::UNKNOWN || o == OrientationInterval::COLINEAR)
@@ -143,7 +144,7 @@ namespace codac2
         // points next-to-top, top, and pts[i] makes
         // a non-left turn
         while(l.size() > 1 &&
-          orientation(*std::prev(std::prev(l.end())), l.back(), pts[i]) == OrientationInterval::COUNTERCLOCKWISE)
+          orientation(*std::prev(std::prev(l.end())), l.back(), pts[i]) == OrientationInterval::CLOCKWISE)
         {
           l.pop_back();
         }
diff --git a/tests/core/geometry/codac2_tests_geometry.cpp b/tests/core/geometry/codac2_tests_geometry.cpp
index 232ca87ef..efc43f6c6 100644
--- a/tests/core/geometry/codac2_tests_geometry.cpp
+++ b/tests/core/geometry/codac2_tests_geometry.cpp
@@ -24,22 +24,22 @@ TEST_CASE("orientation")
   CHECK(orientation({0.,0.},{0.,0.},{0.,0.}) == OrientationInterval::COLINEAR);
 
   // Clockwise points
-  CHECK(orientation({0.,0.},{0.,-5.},{-2.,-5.}) == OrientationInterval::COUNTERCLOCKWISE);
-  CHECK(orientation({-999.,0.1},{0.,0.},{-999.,-0.1}) == OrientationInterval::COUNTERCLOCKWISE);
+  CHECK(orientation({0.,0.},{0.,-5.},{-2.,-5.}) == OrientationInterval::CLOCKWISE);
+  CHECK(orientation({-999.,0.1},{0.,0.},{-999.,-0.1}) == OrientationInterval::CLOCKWISE);
 
-  // Anticlockwise points
-  CHECK(orientation({0.,0.},{0.,-5.},{2.,-5.}) == OrientationInterval::CLOCKWISE);
-  CHECK(orientation({999.,0.1},{0.,0.},{999.,-0.1}) == OrientationInterval::CLOCKWISE);
+  // Counterclockwise points
+  CHECK(orientation({0.,0.},{0.,-5.},{2.,-5.}) == OrientationInterval::COUNTERCLOCKWISE);
+  CHECK(orientation({999.,0.1},{0.,0.},{999.,-0.1}) == OrientationInterval::COUNTERCLOCKWISE);
 
   // Other tests
-  CHECK(orientation({0,0},{1,1},{0,2}) == OrientationInterval::CLOCKWISE);
+  CHECK(orientation({0,0},{1,1},{0,2}) == OrientationInterval::COUNTERCLOCKWISE);
   CHECK(orientation({0,0},{1,1},{2,2}) == OrientationInterval::COLINEAR);
-  CHECK(orientation({0,0},{1,1},{2.2,2}) == OrientationInterval::COUNTERCLOCKWISE);
+  CHECK(orientation({0,0},{1,1},{2.2,2}) == OrientationInterval::CLOCKWISE);
   CHECK(orientation({0,0},{2,2},{2,2}) == OrientationInterval::COLINEAR);
   CHECK(orientation({0,0},{2,2},{4,4}) == OrientationInterval::COLINEAR);
   CHECK(orientation({0,0},{8,8},{4,4}) == OrientationInterval::COLINEAR);
-  CHECK(orientation({0,0},{10,1},{4,4}) == OrientationInterval::CLOCKWISE);
-  CHECK(orientation({0,0},{2,2},{10,1}) == OrientationInterval::COUNTERCLOCKWISE);
+  CHECK(orientation({0,0},{10,1},{4,4}) == OrientationInterval::COUNTERCLOCKWISE);
+  CHECK(orientation({0,0},{2,2},{10,1}) == OrientationInterval::CLOCKWISE);
 }
 
 TEST_CASE("Alignments")
diff --git a/tests/core/geometry/codac2_tests_geometry.py b/tests/core/geometry/codac2_tests_geometry.py
index 2f5994d56..90138608f 100644
--- a/tests/core/geometry/codac2_tests_geometry.py
+++ b/tests/core/geometry/codac2_tests_geometry.py
@@ -22,22 +22,22 @@ def test_orientation(self):
     self.assertTrue(orientation([0.,0.],[0.,0.],[0.,0.]) == OrientationInterval.COLINEAR)
 
     # Clockwise points
-    self.assertTrue(orientation([0.,0.],[0.,-5.],[-2.,-5.]) == OrientationInterval.COUNTERCLOCKWISE)
-    self.assertTrue(orientation([-999.,0.1],[0.,0.],[-999.,-0.1]) == OrientationInterval.COUNTERCLOCKWISE)
+    self.assertTrue(orientation([0.,0.],[0.,-5.],[-2.,-5.]) == OrientationInterval.CLOCKWISE)
+    self.assertTrue(orientation([-999.,0.1],[0.,0.],[-999.,-0.1]) == OrientationInterval.CLOCKWISE)
 
-    # Anticlockwise points
-    self.assertTrue(orientation([0.,0.],[0.,-5.],[2.,-5.]) == OrientationInterval.CLOCKWISE)
-    self.assertTrue(orientation([999.,0.1],[0.,0.],[999.,-0.1]) == OrientationInterval.CLOCKWISE)
+    # Counterclockwise points
+    self.assertTrue(orientation([0.,0.],[0.,-5.],[2.,-5.]) == OrientationInterval.COUNTERCLOCKWISE)
+    self.assertTrue(orientation([999.,0.1],[0.,0.],[999.,-0.1]) == OrientationInterval.COUNTERCLOCKWISE)
 
     # Other tests
-    self.assertTrue(orientation([0,0],[1,1],[0,2]) == OrientationInterval.CLOCKWISE)
+    self.assertTrue(orientation([0,0],[1,1],[0,2]) == OrientationInterval.COUNTERCLOCKWISE)
     self.assertTrue(orientation([0,0],[1,1],[2,2]) == OrientationInterval.COLINEAR)
-    self.assertTrue(orientation([0,0],[1,1],[2.2,2]) == OrientationInterval.COUNTERCLOCKWISE)
+    self.assertTrue(orientation([0,0],[1,1],[2.2,2]) == OrientationInterval.CLOCKWISE)
     self.assertTrue(orientation([0,0],[2,2],[2,2]) == OrientationInterval.COLINEAR)
     self.assertTrue(orientation([0,0],[2,2],[4,4]) == OrientationInterval.COLINEAR)
     self.assertTrue(orientation([0,0],[8,8],[4,4]) == OrientationInterval.COLINEAR)
-    self.assertTrue(orientation([0,0],[10,1],[4,4]) == OrientationInterval.CLOCKWISE)
-    self.assertTrue(orientation([0,0],[2,2],[10,1]) == OrientationInterval.COUNTERCLOCKWISE)
+    self.assertTrue(orientation([0,0],[10,1],[4,4]) == OrientationInterval.COUNTERCLOCKWISE)
+    self.assertTrue(orientation([0,0],[2,2],[10,1]) == OrientationInterval.CLOCKWISE)
 
   def test_alignments(self):
     

From 393d94cba2e82b2fb1c57813761c2acb4e5a4bf7 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:41:59 +0200
Subject: [PATCH 08/16] [geom] Handle two-point convex hull after collinear
 filtering

---
 src/core/geometry/codac2_geometry.cpp         | 6 ++++++
 tests/core/geometry/codac2_tests_geometry.cpp | 9 +++++++++
 tests/core/geometry/codac2_tests_geometry.py  | 9 +++++++++
 3 files changed, 24 insertions(+)

diff --git a/src/core/geometry/codac2_geometry.cpp b/src/core/geometry/codac2_geometry.cpp
index c299795cb..c70811819 100644
--- a/src/core/geometry/codac2_geometry.cpp
+++ b/src/core/geometry/codac2_geometry.cpp
@@ -132,6 +132,12 @@ namespace codac2
         m++; // Update size of modified array
       }
 
+      if(m < 3)
+      {
+        pts.resize(m);
+        return pts;
+      }
+
     // Create an empty stack and push first three points to it.
 
       list<IntervalVector> l({ pts[0], pts[1], pts[2] });
diff --git a/tests/core/geometry/codac2_tests_geometry.cpp b/tests/core/geometry/codac2_tests_geometry.cpp
index efc43f6c6..f654a23e1 100644
--- a/tests/core/geometry/codac2_tests_geometry.cpp
+++ b/tests/core/geometry/codac2_tests_geometry.cpp
@@ -150,6 +150,15 @@ TEST_CASE("convex_hull")
   CHECK(hull[6] == v_save[8]);
   CHECK(hull.size() == 7);
 
+  // Fully collinear points reduce to the two endpoints.
+  hull = convex_hull({
+    {0,0},{1,0},{2,0},{3,0}
+  });
+
+  CHECK(hull.size() == 2);
+  CHECK(hull[0] == IntervalVector({0,0}));
+  CHECK(hull[1] == IntervalVector({3,0}));
+
   //for(const auto& vi : v)
   //  DefaultFigure::draw_point(vi.mid(), Color::blue());
   //for(size_t i = 0 ; i < hull.size() ; i++)
diff --git a/tests/core/geometry/codac2_tests_geometry.py b/tests/core/geometry/codac2_tests_geometry.py
index 90138608f..1fcda460e 100644
--- a/tests/core/geometry/codac2_tests_geometry.py
+++ b/tests/core/geometry/codac2_tests_geometry.py
@@ -132,5 +132,14 @@ def test_convex_hull(self):
     self.assertTrue(hull[6] == v_save[8])
     self.assertTrue(len(hull) == 7)
 
+    # Fully collinear points reduce to the two endpoints.
+    hull = convex_hull([
+      [0,0],[1,0],[2,0],[3,0]
+    ])
+
+    self.assertTrue(len(hull) == 2)
+    self.assertTrue(hull[0] == IntervalVector([0,0]))
+    self.assertTrue(hull[1] == IntervalVector([3,0]))
+
 if __name__ ==  '__main__':
   unittest.main()
\ No newline at end of file

From 4557c6a58bb616a362f291545b8faf6bc10457b1 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:48:32 +0200
Subject: [PATCH 09/16] [geom] Fix merging of adjacent interval points

---
 .../src/core/geometry/codac2_py_geometry.cpp  | 11 ++++++-
 src/core/geometry/codac2_geometry.cpp         | 16 ++++------
 tests/core/geometry/codac2_tests_geometry.cpp | 32 ++++++++++++++++++-
 tests/core/geometry/codac2_tests_geometry.py  | 25 +++++++++++++++
 4 files changed, 72 insertions(+), 12 deletions(-)

diff --git a/python/src/core/geometry/codac2_py_geometry.cpp b/python/src/core/geometry/codac2_py_geometry.cpp
index 479ddd450..9e2264b26 100644
--- a/python/src/core/geometry/codac2_py_geometry.cpp
+++ b/python/src/core/geometry/codac2_py_geometry.cpp
@@ -70,7 +70,16 @@ void export_geometry(py::module& m)
     T_ROTATE_CONST_T_REF_CONST_INTERVAL_REF_CONST_INTERVALVECTOR_REF,
     "x"_a, "a"_a, "c"_a = Vector::zero(2));
 
-  m.def("merge_adjacent_points", &merge_adjacent_points,
+  m.def("merge_adjacent_points",
+    [](py::list l)
+    {
+      auto points = l.cast<list<IntervalVector>>();
+      merge_adjacent_points(points);
+
+      l.attr("clear")();
+      for(const auto& p : points)
+        l.append(p);
+    },
     VOID_MERGE_ADJACENT_POINTS_LIST_INTERVALVECTOR_REF,
     "l"_a);
 }
\ No newline at end of file
diff --git a/src/core/geometry/codac2_geometry.cpp b/src/core/geometry/codac2_geometry.cpp
index c70811819..7a467daa6 100644
--- a/src/core/geometry/codac2_geometry.cpp
+++ b/src/core/geometry/codac2_geometry.cpp
@@ -188,24 +188,20 @@ namespace codac2
     do
     {
       merged_something = false;
-      for(auto li = l.begin(); li != l.end();)
+
+      for(auto li = l.begin(); li != l.end() && !merged_something; ++li)
       {
-        for(auto lj = l.begin(); lj != l.end(); ++lj)
+        auto lj = std::next(li);
+        for(; lj != l.end(); ++lj)
         {
-          if(li != lj && li->intersects(*lj))
+          if(li->intersects(*lj))
           {
             *li |= *lj;
-            li = l.erase(li);
+            l.erase(lj);
             merged_something = true;
             break;
           }
-
-          else
-            ++li;
         }
-
-        if(merged_something)
-          break;
       }
     } while(merged_something);
   }
diff --git a/tests/core/geometry/codac2_tests_geometry.cpp b/tests/core/geometry/codac2_tests_geometry.cpp
index f654a23e1..970a856d4 100644
--- a/tests/core/geometry/codac2_tests_geometry.cpp
+++ b/tests/core/geometry/codac2_tests_geometry.cpp
@@ -163,4 +163,34 @@ TEST_CASE("convex_hull")
   //  DefaultFigure::draw_point(vi.mid(), Color::blue());
   //for(size_t i = 0 ; i < hull.size() ; i++)
   //  DefaultFigure::draw_line(hull[i].mid(), hull[(i+1)%hull.size()].mid(), Color::blue());
-}
\ No newline at end of file
+}
+
+TEST_CASE("geometry - merge adjacent points")
+{
+  list<IntervalVector> pts {
+    {{0,1},{0,1}},
+    {{0.5,2},{0.5,1.5}},
+    {{3,4},{3,4}}
+  };
+
+  merge_adjacent_points(pts);
+
+  REQUIRE(pts.size() == 2);
+  auto it = pts.begin();
+  CHECK(*it == IntervalVector({{0,2},{0,1.5}}));
+  ++it;
+  CHECK(*it == IntervalVector({{3,4},{3,4}}));
+
+  // Merging can create a new intersection, so the process must be iterative.
+  pts = {
+    {{0,1},{0,1}},
+    {{1,2},{0,1}},
+    {{2,3},{0,1}}
+  };
+
+  merge_adjacent_points(pts);
+
+  REQUIRE(pts.size() == 1);
+  CHECK(pts.front() == IntervalVector({{0,3},{0,1}}));
+}
+
diff --git a/tests/core/geometry/codac2_tests_geometry.py b/tests/core/geometry/codac2_tests_geometry.py
index 1fcda460e..29f4319d5 100644
--- a/tests/core/geometry/codac2_tests_geometry.py
+++ b/tests/core/geometry/codac2_tests_geometry.py
@@ -141,5 +141,30 @@ def test_convex_hull(self):
     self.assertTrue(hull[0] == IntervalVector([0,0]))
     self.assertTrue(hull[1] == IntervalVector([3,0]))
 
+
+    pts = [
+      IntervalVector([[0,1],[0,1]]),
+      IntervalVector([[0.5,2],[0.5,1.5]]),
+      IntervalVector([[3,4],[3,4]])
+    ]
+
+    merge_adjacent_points(pts)
+
+    self.assertTrue(len(pts) == 2)
+    self.assertTrue(pts[0] == IntervalVector([[0,2],[0,1.5]]))
+    self.assertTrue(pts[1] == IntervalVector([[3,4],[3,4]]))
+
+    # Merging can create a new intersection, so the process must be iterative.
+    pts = [
+      IntervalVector([[0,1],[0,1]]),
+      IntervalVector([[1,2],[0,1]]),
+      IntervalVector([[2,3],[0,1]])
+    ]
+
+    merge_adjacent_points(pts)
+
+    self.assertTrue(len(pts) == 1)
+    self.assertTrue(pts[0] == IntervalVector([[0,3],[0,1]]))
+
 if __name__ ==  '__main__':
   unittest.main()
\ No newline at end of file

From d799fd3231e08cf949cab7dfd1a61a30980d053a Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:52:23 +0200
Subject: [PATCH 10/16] [geom] Handle empty inputs in aligned

---
 src/core/geometry/codac2_geometry.cpp         | 5 ++++-
 tests/core/geometry/codac2_tests_geometry.cpp | 3 +++
 tests/core/geometry/codac2_tests_geometry.py  | 3 +++
 3 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/core/geometry/codac2_geometry.cpp b/src/core/geometry/codac2_geometry.cpp
index 7a467daa6..7381264a6 100644
--- a/src/core/geometry/codac2_geometry.cpp
+++ b/src/core/geometry/codac2_geometry.cpp
@@ -39,7 +39,10 @@ namespace codac2
 
     Interval cross_product = (p2[0]-p1[0])*(p3[1]-p1[1])-(p2[1]-p1[1])*(p3[0]-p1[0]);
 
-    if(cross_product == 0.)
+    if(cross_product.is_empty())
+      return BoolInterval::EMPTY;
+
+    else if(cross_product == 0.)
       return BoolInterval::TRUE;
 
     else if(!cross_product.contains(0.))
diff --git a/tests/core/geometry/codac2_tests_geometry.cpp b/tests/core/geometry/codac2_tests_geometry.cpp
index 970a856d4..286632df3 100644
--- a/tests/core/geometry/codac2_tests_geometry.cpp
+++ b/tests/core/geometry/codac2_tests_geometry.cpp
@@ -61,6 +61,9 @@ TEST_CASE("Alignments")
   CHECK(aligned({0,0},{1,1},{10,10}) == BoolInterval::TRUE);
   CHECK(aligned({{-oo,0},{-oo,0}},{1,1},{10,10}) == BoolInterval::UNKNOWN);
   CHECK(aligned(IntervalVector(2),{1,1},{10,10}) == BoolInterval::UNKNOWN);
+  CHECK(aligned(IntervalVector::empty(2),{1,1},{10,10}) == BoolInterval::EMPTY);
+  CHECK(aligned({1,1},IntervalVector::empty(2),{10,10}) == BoolInterval::EMPTY);
+  CHECK(aligned({1,1},{10,10},IntervalVector::empty(2)) == BoolInterval::EMPTY);
 }
 
 TEST_CASE("convex_hull")
diff --git a/tests/core/geometry/codac2_tests_geometry.py b/tests/core/geometry/codac2_tests_geometry.py
index 29f4319d5..432031371 100644
--- a/tests/core/geometry/codac2_tests_geometry.py
+++ b/tests/core/geometry/codac2_tests_geometry.py
@@ -59,6 +59,9 @@ def test_alignments(self):
     self.assertTrue(aligned([0,0],[1,1],[10,10]) == BoolInterval.TRUE)
     self.assertTrue(aligned([[-oo,0],[-oo,0]],[1,1],[10,10]) == BoolInterval.UNKNOWN)
     self.assertTrue(aligned(IntervalVector(2),[1,1],[10,10]) == BoolInterval.UNKNOWN)
+    self.assertTrue(aligned(IntervalVector.empty(2),[1,1],[10,10]) == BoolInterval.EMPTY)
+    self.assertTrue(aligned([1,1],IntervalVector.empty(2),[10,10]) == BoolInterval.EMPTY)
+    self.assertTrue(aligned([1,1],[10,10],IntervalVector.empty(2)) == BoolInterval.EMPTY)
 
   def test_convex_hull(self):
 

From 447981499c0e1b946fd4ac2551270f8e81374cdb Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 18:57:29 +0200
Subject: [PATCH 11/16] [geom] Fix polygon equality with multiple contours

---
 src/core/geometry/codac2_Polygon.cpp         | 31 +++++++++++---------
 tests/core/geometry/codac2_tests_Polygon.cpp | 26 ++++++++++++++++
 tests/core/geometry/codac2_tests_Polygon.py  | 24 +++++++++++++++
 3 files changed, 67 insertions(+), 14 deletions(-)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index e76037538..dd7ea6012 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -318,25 +318,28 @@ namespace codac2
 
   bool Polygon::operator==(const Polygon& p) const
   {
-    if(is_empty() || p.is_empty())
-      return is_empty() && p.is_empty();
-
-    size_t n = this->size();
-    if(p.size() != n)
+    if(p.size() != this->size())
       return false;
 
-    size_t i; // looking for same reference of first value
-    for(i = 0 ; i < n ; i++)
-      if((*this)[0] == p[i])
-        break;
+    vector<bool> matched(p.size(),false);
 
-    int way = 1;
-    if(n > 1)
-      way = ((*this)[1] == p[(i+1)%n]) ? 1 : -1;
+    for(const auto& e : *this)
+    {
+      bool found = false;
+
+      for(size_t i = 0 ; i < p.size() ; i++)
+      {
+        if(!matched[i] && e == p[i])
+        {
+          matched[i] = true;
+          found = true;
+          break;
+        }
+      }
 
-    for(size_t j = 0 ; j < n ; j++)
-      if((*this)[j] != p[(i+way*j+2*n)%n])
+      if(!found)
         return false;
+    }
 
     return true;
   }
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index 4f0104d18..8993c6852 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -30,6 +30,32 @@ TEST_CASE("Polygon base")
   CHECK(p3 == Polygon({{5,1},{1,2},{3,4}}));
   CHECK(p3 == Polygon({{1,2},{3,4},{5,1}}));
 
+  { // Multiple contours can have independent orientations
+    vector<Segment> edges1 {
+      {{0,0},{1,0}}, {{1,0},{0,1}}, {{0,1},{0,0}},
+      {{3,0},{4,0}}, {{4,0},{3,1}}, {{3,1},{3,0}}
+    };
+
+    vector<Segment> edges2 {
+      {{0,0},{0,1}}, {{0,1},{1,0}}, {{1,0},{0,0}},
+      {{3,0},{4,0}}, {{4,0},{3,1}}, {{3,1},{3,0}}
+    };
+
+    CHECK(Polygon(edges1) == Polygon(edges2));
+  }
+
+  { // Same number of edges, but a different polygon
+    vector<Segment> edges1 {
+      {{0,0},{1,0}}, {{1,0},{0,1}}, {{0,1},{0,0}}
+    };
+
+    vector<Segment> edges2 {
+      {{0,0},{1,0}}, {{1,0},{0,2}}, {{0,2},{0,0}}
+    };
+
+    CHECK_FALSE(Polygon(edges1) == Polygon(edges2));
+  }
+
   { // Consecutive duplicate vertices
     Polygon p({{0,0},{1,0},{1,0},{2,0},{2,1},{0,1}});
     CHECK(p == Polygon({{0,0},{2,0},{2,1},{0,1}}));
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index c78624011..effefa581 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -27,6 +27,30 @@ def test_polygon_base(self):
     self.assertTrue(p3 == Polygon([[5,1],[1,2],[3,4]]))
     self.assertTrue(p3 == Polygon([[1,2],[3,4],[5,1]]))
 
+    # Multiple contours can have independent orientations.
+    edges1 = [
+      Segment([0,0],[1,0]), Segment([1,0],[0,1]), Segment([0,1],[0,0]),
+      Segment([3,0],[4,0]), Segment([4,0],[3,1]), Segment([3,1],[3,0])
+    ]
+
+    edges2 = [
+      Segment([0,0],[0,1]), Segment([0,1],[1,0]), Segment([1,0],[0,0]),
+      Segment([3,0],[4,0]), Segment([4,0],[3,1]), Segment([3,1],[3,0])
+    ]
+
+    self.assertTrue(Polygon(edges1) == Polygon(edges2))
+
+    # Same number of edges, but a different polygon.
+    edges1 = [
+      Segment([0,0],[1,0]), Segment([1,0],[0,1]), Segment([0,1],[0,0])
+    ]
+
+    edges2 = [
+      Segment([0,0],[1,0]), Segment([1,0],[0,2]), Segment([0,2],[0,0])
+    ]
+
+    self.assertFalse(Polygon(edges1) == Polygon(edges2))
+
     # Consecutive duplicate vertices.
     p = Polygon([[0,0],[1,0],[1,0],[2,0],[2,1],[0,1]])
     self.assertTrue(p == Polygon([[0,0],[2,0],[2,1],[0,1]]))

From 683c93528191209bef6dd371b6d2a969ab7ea193 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 19:01:02 +0200
Subject: [PATCH 12/16] [geom] Fix duplicate Segment test name

---
 tests/core/geometry/codac2_tests_Segment.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/core/geometry/codac2_tests_Segment.py b/tests/core/geometry/codac2_tests_Segment.py
index 45c1f10b5..b1a027002 100644
--- a/tests/core/geometry/codac2_tests_Segment.py
+++ b/tests/core/geometry/codac2_tests_Segment.py
@@ -21,7 +21,7 @@ def hull(x):
 
 class TestSegment(unittest.TestCase):
 
-  def test_contains(self):
+  def test_contains_horizontal_segment(self):
 
     p1 = Vector([0,-10])
     p2 = Vector([-10,-10])
@@ -43,7 +43,7 @@ def test_intersects(self):
     self.assertTrue(Segment(IntervalVector([hull(-oo),2]),Vector([0,2])).intersects(Segment(Vector([0,0]),Vector([0,1]))) == BoolInterval.FALSE)
     self.assertTrue(Segment(IntervalVector([hull(-oo),2]),Vector([0,2])).intersects(Segment(Vector([0,0]),Vector([0,2]))) == BoolInterval.TRUE)
 
-  def test_contains(self):
+  def test_contains_oblique_segment(self):
 
     p1,p2,p3 = Vector([0.,0.]), Vector([1.5,3.]), Vector([0.,2.])
     e1 = Segment(Vector([-10.,-10.]), Vector([10.,10.]))

From a5d295eae3df1fa523e59595d57c2ad883193f85 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 19:10:11 +0200
Subject: [PATCH 13/16] [geom] Fix vertices of multi-contour polygons

---
 src/core/geometry/codac2_Polygon.cpp         | 9 +++++++--
 tests/core/geometry/codac2_tests_Polygon.cpp | 9 +++++++++
 tests/core/geometry/codac2_tests_Polygon.py  | 9 +++++++++
 3 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index dd7ea6012..9f61fb44c 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -171,12 +171,17 @@ namespace codac2
 
     if(!is_empty())
     {
-      size_t n = this->size();
+      size_t n = this->size(), ring_begin = 0;
 
       for(size_t i = 0 ; i < n ; i++)
+      {
         v.push_back((*this)[i][0]);
 
-      if((*this)[0][0] != (*this)[n-1][1])
+        if(i+1 < n && !(*this)[i][1].intersects((*this)[i+1][0]))
+          ring_begin = i+1;
+      }
+
+      if((*this)[ring_begin][0] != (*this)[n-1][1])
         v.push_back((*this)[n-1][1]);
     }
 
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index 8993c6852..41dd82e74 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -42,6 +42,15 @@ TEST_CASE("Polygon base")
     };
 
     CHECK(Polygon(edges1) == Polygon(edges2));
+
+    auto vertices = Polygon(edges1).vertices();
+    REQUIRE(vertices.size() == 6);
+    CHECK(vertices[0] == IntervalVector({0,0}));
+    CHECK(vertices[1] == IntervalVector({1,0}));
+    CHECK(vertices[2] == IntervalVector({0,1}));
+    CHECK(vertices[3] == IntervalVector({3,0}));
+    CHECK(vertices[4] == IntervalVector({4,0}));
+    CHECK(vertices[5] == IntervalVector({3,1}));
   }
 
   { // Same number of edges, but a different polygon
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index effefa581..908db1b42 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -40,6 +40,15 @@ def test_polygon_base(self):
 
     self.assertTrue(Polygon(edges1) == Polygon(edges2))
 
+    vertices = Polygon(edges1).vertices()
+    self.assertTrue(len(vertices) == 6)
+    self.assertTrue(vertices[0] == IntervalVector([0,0]))
+    self.assertTrue(vertices[1] == IntervalVector([1,0]))
+    self.assertTrue(vertices[2] == IntervalVector([0,1]))
+    self.assertTrue(vertices[3] == IntervalVector([3,0]))
+    self.assertTrue(vertices[4] == IntervalVector([4,0]))
+    self.assertTrue(vertices[5] == IntervalVector([3,1]))
+
     # Same number of edges, but a different polygon.
     edges1 = [
       Segment([0,0],[1,0]), Segment([1,0],[0,1]), Segment([0,1],[0,0])

From b39ae7742edeaed2d18f5e5aaa1e910c2f19a4dd Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 19:17:57 +0200
Subject: [PATCH 14/16] [geom] Propagate empty points in polygon containment

---
 src/core/geometry/codac2_Polygon.cpp         | 3 +++
 tests/core/geometry/codac2_tests_Polygon.cpp | 2 ++
 tests/core/geometry/codac2_tests_Polygon.py  | 2 ++
 3 files changed, 7 insertions(+)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index 9f61fb44c..e76bf701e 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -205,6 +205,9 @@ namespace codac2
   {
     assert_release(p.size() == 2);
 
+    if(p.is_empty())
+      return BoolInterval::EMPTY;
+
     IntervalVector box_hull = box();
     if(!p.intersects(box_hull))
       return BoolInterval::FALSE;
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index 41dd82e74..4df2def5e 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -81,6 +81,7 @@ TEST_CASE("Empty polygon")
   Polygon p1 = Polygon::empty();
   CHECK(p1.contains(IntervalVector({1,1})) == BoolInterval::FALSE);
   CHECK(p1.contains(IntervalVector(2)) == BoolInterval::FALSE);
+  CHECK(p1.contains(IntervalVector::empty(2)) == BoolInterval::EMPTY);
   CHECK(p1.is_empty());
   CHECK(p1.edges().size() == 0);
   CHECK(p1.vertices().size() == 0);
@@ -89,6 +90,7 @@ TEST_CASE("Empty polygon")
 TEST_CASE("Polygon contains")
 {
   Polygon p1({{3,-1},{3,4},{5,6},{-1,1}});
+  CHECK(p1.contains(IntervalVector::empty(2)) == BoolInterval::EMPTY);
   CHECK(p1.contains({3.1,3}) == BoolInterval::FALSE);
   CHECK(p1.contains({2.9,3}) == BoolInterval::TRUE);
   CHECK(p1.contains({3,3}) == BoolInterval::TRUE);
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index 908db1b42..9e6522327 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -73,6 +73,7 @@ def test_empty_polygon(self):
     p1 = Polygon.empty()
     self.assertTrue(p1.contains(IntervalVector([1,1])) == BoolInterval.FALSE)
     self.assertTrue(p1.contains(IntervalVector(2)) == BoolInterval.FALSE)
+    self.assertTrue(p1.contains(IntervalVector.empty(2)) == BoolInterval.EMPTY)
     self.assertTrue(p1.is_empty())
     self.assertTrue(len(p1.edges()) == 0)
     self.assertTrue(len(p1.vertices()) == 0)
@@ -80,6 +81,7 @@ def test_empty_polygon(self):
   def test_Polygon(self):
 
     p1 = Polygon([[3,-1],[3,4],[5,6],[-1,1]])
+    self.assertTrue(p1.contains(IntervalVector.empty(2)) == BoolInterval.EMPTY)
     self.assertTrue(p1.contains([3.1,3]) == BoolInterval.FALSE)
     self.assertTrue(p1.contains([2.9,3]) == BoolInterval.TRUE)
     self.assertTrue(p1.contains([3,3]) == BoolInterval.TRUE)

From ec3a5253eee9ca2be4c2ca5c940cbf92e04c33b4 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 19:22:06 +0200
Subject: [PATCH 15/16] [geom] Handle empty segments in colinear

---
 src/core/geometry/codac2_Segment.cpp         |  5 ++++-
 tests/core/geometry/codac2_tests_Segment.cpp | 12 ++++++++++++
 tests/core/geometry/codac2_tests_Segment.py  | 10 ++++++++++
 3 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/src/core/geometry/codac2_Segment.cpp b/src/core/geometry/codac2_Segment.cpp
index 53e31d620..781e244bf 100644
--- a/src/core/geometry/codac2_Segment.cpp
+++ b/src/core/geometry/codac2_Segment.cpp
@@ -200,7 +200,10 @@ namespace codac2
     Interval c = ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
     auto a = aligned(e1[0],e1[1],e2[0]);
 
-    if(c == 0 && a == BoolInterval::TRUE)
+    if(c.is_empty() || a == BoolInterval::EMPTY)
+      return BoolInterval::EMPTY;
+
+    else if(c == 0 && a == BoolInterval::TRUE)
       return BoolInterval::TRUE;
 
     else if(c.contains(0.) // e1 and e2 may be parallel
diff --git a/tests/core/geometry/codac2_tests_Segment.cpp b/tests/core/geometry/codac2_tests_Segment.cpp
index ae13e99e6..437b95c1a 100644
--- a/tests/core/geometry/codac2_tests_Segment.cpp
+++ b/tests/core/geometry/codac2_tests_Segment.cpp
@@ -166,6 +166,18 @@ TEST_CASE("Segment")
     CHECK((e1 & e2).is_empty());
   }
 
+  SECTION("colinear - empty segments")
+  {
+    Segment empty(
+      IntervalVector::empty(2),
+      IntervalVector::empty(2)
+    );
+    Segment e({0,0},{1,1});
+
+    CHECK(colinear(empty,e) == BoolInterval::EMPTY);
+    CHECK(colinear(e,empty) == BoolInterval::EMPTY);
+  }
+
   SECTION("intersection - near infinite cases")
   {
     Segment e1({-1,6},{-1,next_float(-oo)});
diff --git a/tests/core/geometry/codac2_tests_Segment.py b/tests/core/geometry/codac2_tests_Segment.py
index b1a027002..97852af4b 100644
--- a/tests/core/geometry/codac2_tests_Segment.py
+++ b/tests/core/geometry/codac2_tests_Segment.py
@@ -154,6 +154,16 @@ def test_intersection(self):
     self.assertTrue(e1.intersects(e2) == BoolInterval.FALSE)
     self.assertTrue((e1 & e2).is_empty())
 
+    # Empty segments propagate EMPTY in colinearity tests
+    empty = Segment(
+      IntervalVector.empty(2),
+      IntervalVector.empty(2)
+    )
+    e = Segment([0,0],[1,1])
+
+    self.assertTrue(colinear(empty,e) == BoolInterval.EMPTY)
+    self.assertTrue(colinear(e,empty) == BoolInterval.EMPTY)
+
     # Near infinite cases
     e1 = Segment([-1,6],[-1,next_float(-oo)])
     e2 = Segment([-1,-1],[3,-6])

From 6b01e64d69f0ce1dfbac04e87346ac1b8594d553 Mon Sep 17 00:00:00 2001
From: Simon Rohou <simon.rohou@gmail.com>
Date: Wed, 14 Oct 2026 19:26:21 +0200
Subject: [PATCH 16/16] [geom] Handle degenerate polygons after vertex
 simplification

---
 src/core/geometry/codac2_Polygon.cpp         | 11 +++++++++--
 tests/core/geometry/codac2_tests_Polygon.cpp | 10 ++++++++++
 tests/core/geometry/codac2_tests_Polygon.py  |  9 +++++++++
 3 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/src/core/geometry/codac2_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index e76bf701e..52dec2431 100644
--- a/src/core/geometry/codac2_Polygon.cpp
+++ b/src/core/geometry/codac2_Polygon.cpp
@@ -92,8 +92,15 @@ namespace codac2
               }
             }
 
-            for(size_t i = 0 ; i < v_noaligned.size() ; i++)
-              edges.push_back({ v_noaligned[i], v_noaligned[(i+1)%v_noaligned.size()] });
+            if(v_noaligned.size() == 1)
+              edges = { Segment(v_noaligned[0],v_noaligned[0]) };
+
+            else if(v_noaligned.size() == 2)
+              edges = { Segment(v_noaligned[0],v_noaligned[1]) };
+
+            else
+              for(size_t i = 0 ; i < v_noaligned.size() ; i++)
+                edges.push_back({ v_noaligned[i], v_noaligned[(i+1)%v_noaligned.size()] });
           }
         }
         
diff --git a/tests/core/geometry/codac2_tests_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index 4df2def5e..9845b9c1a 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -65,6 +65,16 @@ TEST_CASE("Polygon base")
     CHECK_FALSE(Polygon(edges1) == Polygon(edges2));
   }
 
+  { // Fully collinear vertices reduce to one segment
+    Polygon p({{0,0},{1,0},{2,0},{3,0}});
+    Polygon expected({{0,0},{3,0}});
+
+    CHECK(p == expected);
+    REQUIRE(p.edges().size() == 1);
+    CHECK(p.vertices().size() == 2);
+    CHECK(p.contains({1.5,0}) == BoolInterval::TRUE);
+  }
+
   { // Consecutive duplicate vertices
     Polygon p({{0,0},{1,0},{1,0},{2,0},{2,1},{0,1}});
     CHECK(p == Polygon({{0,0},{2,0},{2,1},{0,1}}));
diff --git a/tests/core/geometry/codac2_tests_Polygon.py b/tests/core/geometry/codac2_tests_Polygon.py
index 9e6522327..9a6b4be90 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -60,6 +60,15 @@ def test_polygon_base(self):
 
     self.assertFalse(Polygon(edges1) == Polygon(edges2))
 
+    # Fully collinear vertices reduce to one segment.
+    p = Polygon([[0,0],[1,0],[2,0],[3,0]])
+    expected = Polygon([[0,0],[3,0]])
+
+    self.assertTrue(p == expected)
+    self.assertTrue(len(p.edges()) == 1)
+    self.assertTrue(len(p.vertices()) == 2)
+    self.assertTrue(p.contains([1.5,0]) == BoolInterval.TRUE)
+
     # Consecutive duplicate vertices.
     p = Polygon([[0,0],[1,0],[1,0],[2,0],[2,1],[0,1]])
     self.assertTrue(p == Polygon([[0,0],[2,0],[2,1],[0,1]]))