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
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/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/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())
 
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()
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_Polygon.cpp b/src/core/geometry/codac2_Polygon.cpp
index a7f909549..52dec2431 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:
@@ -85,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()] });
           }
         }
         
@@ -164,12 +178,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]);
     }
 
@@ -193,6 +212,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;
@@ -215,18 +237,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 +263,30 @@ 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 +297,6 @@ namespace codac2
         transect = try_transect;
 
       i++;
-      assert(i < 10);
     } while(retry);
 
     // Now the number i of crossing can be computed.
@@ -294,25 +333,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/src/core/geometry/codac2_Segment.cpp b/src/core/geometry/codac2_Segment.cpp
index 5aa7826c0..781e244bf 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);
@@ -199,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/src/core/geometry/codac2_geometry.cpp b/src/core/geometry/codac2_geometry.cpp
index 66b8882cd..7381264a6 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;
@@ -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.))
@@ -69,7 +72,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 +108,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)
@@ -131,6 +135,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] });
@@ -143,7 +153,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();
         }
@@ -181,24 +191,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_Polygon.cpp b/tests/core/geometry/codac2_tests_Polygon.cpp
index dae41c7b3..9845b9c1a 100644
--- a/tests/core/geometry/codac2_tests_Polygon.cpp
+++ b/tests/core/geometry/codac2_tests_Polygon.cpp
@@ -29,6 +29,61 @@ 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}}));
+
+  { // 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));
+
+    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
+    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));
+  }
+
+  { // 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}}));
+  }
+
+  { // 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")
@@ -36,6 +91,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);
@@ -44,6 +100,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);
@@ -132,4 +189,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..9a6b4be90 100644
--- a/tests/core/geometry/codac2_tests_Polygon.py
+++ b/tests/core/geometry/codac2_tests_Polygon.py
@@ -27,11 +27,62 @@ 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))
+
+    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])
+    ]
+
+    edges2 = [
+      Segment([0,0],[1,0]), Segment([1,0],[0,2]), Segment([0,2],[0,0])
+    ]
+
+    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]]))
+
+    # 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()
     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)
@@ -39,6 +90,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)
@@ -116,5 +168,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
diff --git a/tests/core/geometry/codac2_tests_Segment.cpp b/tests/core/geometry/codac2_tests_Segment.cpp
index 926342d6e..437b95c1a 100644
--- a/tests/core/geometry/codac2_tests_Segment.cpp
+++ b/tests/core/geometry/codac2_tests_Segment.cpp
@@ -137,6 +137,47 @@ 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("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 3b9e53114..97852af4b 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.]))
@@ -129,6 +129,41 @@ 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())
+
+    # 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])
diff --git a/tests/core/geometry/codac2_tests_geometry.cpp b/tests/core/geometry/codac2_tests_geometry.cpp
index 232ca87ef..286632df3 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")
@@ -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")
@@ -150,8 +153,47 @@ 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++)
   //  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 2f5994d56..432031371 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):
     
@@ -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):
 
@@ -132,5 +135,39 @@ 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]))
+
+
+    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