[sep] SepImage

doc/manual/index.rst

@@ -283,6 +283,7 @@ User manual
       * SepNot
       * SepAction
    * Analytic separators
+      * :ref:`sec-ctc-analytic-sepimage`
       * SepInverse
       * SepTransform
    * Geometrical separators

doc/manual/manual/contractors/analytic/index.rst

@@ -4,5 +4,6 @@ Analytic contractors
 .. toctree::
 
    Directed operators <http://codac.io>
+   sepimage.rst
    ctcinverse.rst
-   CtcInverseNotIn <http://codac.io>
+   CtcInverseNotIn <http://codac.io>

doc/manual/manual/contractors/analytic/sepimage.rst

@@ -0,0 +1,176 @@
+.. _sec-ctc-analytic-sepimage:
+
+The SepImage separator
+======================
+
+  Main author: `Maël Godard <https://godardma.github.io>`_
+
+Definition
+----------
+
+Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. 
+The ``SepImage`` separator allows one to handle constraints of the form :math:`\mathbf{y}\in\mathbf{f}(\mathbb{X})` 
+by separating input boxes :math:`[\mathbf{y}]\in\mathbb{IR}^p`.
+
+Construction and basic usage
+----------------------------
+
+The ``SepImage`` relies on a boundary approach to compute the image set. The boundary of the initial set needs to be 
+covered by a :ref:`gnomonic atlas <subsec-functions-peibos-gnomonic-atlas>`. To compute the image of the boundary,
+the initial box :math:`\left[-1,1\right]^m` can be subdivided up to a size :math:`\epsilon` and the image of each resulting box 
+is computed. Note that if one chooses to take :math:`\epsilon=2`, only one computation will be done per chart of the atlas. 
+
+Once this image of the boundary has been computed, the ``SepImage`` needs to a way to characterize if a point is inside or 
+outside of the image set. To do so, it looks for an antecedent of the point in the initial set. An additionnal contractor on 
+the initial set is then required.
+
+The typical workflow is:
+
+1. Define analytic variables (scalar, vector, matrix) associated with the domain of the function.
+2. Build an :class:`~codac2.AnalyticFunction`.
+3. Define a gnomonic atlas on the boundary of the initial set
+4. Define a contractor on this same initial set
+5. Instantiate ``SepImage`` with the atlas, the function, the resolution :math:`\epsilon` and the contractor.
+6. Contract an input box :math:`[\mathbf{y}]` or pave the separator in the image space.
+
+Example
+-------
+
+Consider that we want to construct a separator on the image of the unit disk by the function
+
+.. math:: 
+  \mathbf{f}(\mathbf{x})=
+  \left(
+  \begin{array}{c}
+  3(x_{1}+1)\\
+  x_{2}+ 0.5\sin (3 x_{1})
+  \end{array}
+  \right)
+
+This function can be constructed in codac as follows
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [sepimage-1-beg]
+      :end-before: [sepimage-1-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [sepimage-1-beg]
+      :end-before: [sepimage-1-end]
+      :dedent: 4
+
+We then need a contractor for the unit disk, and a gnomonic atlas for its boundary (the unit circle).
+For the atlas, the image of :math:`\left[-1,1\right]` by the function
+
+.. math:: 
+  \psi_{0} =
+  \left(
+  \begin{array}{c}
+  \cos\left(x\cdot\frac{\pi}{2}\right)\\
+  \sin\left(x\cdot\frac{\pi}{2}\right)
+  \end{array}
+  \right)
+
+is half of the unit cicle. The other half can be obtained with a rotation of :math:`\pi` rad. Such atlas 
+is constucted in codac as  
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [sepimage-2-beg]
+      :end-before: [sepimage-2-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [sepimage-2-beg]
+      :end-before: [sepimage-2-end]
+      :dedent: 4
+
+In this example, the constraint on the initial set can be seen as a distance constraint. It can be treated as an inversion as follows :
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [sepimage-3-beg]
+      :end-before: [sepimage-3-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [sepimage-3-beg]
+      :end-before: [sepimage-3-end]
+      :dedent: 4
+
+The separator can then be constructed and used, for example with a paver to get both an inner and an outer approximation of the image set.
+The resulting paving is showed in the next figure.
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [sepimage-4-beg]
+      :end-before: [sepimage-4-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [sepimage-4-beg]
+      :end-before: [sepimage-4-end]
+      :dedent: 4
+
+.. figure:: ./sep_image.png
+  :width: 400px
+  :align: center
+
+With a lower resolution, the performances of the separator improve as shown in the following figure.
+
+.. tabs::
+
+  .. group-tab:: Python
+
+    .. literalinclude:: src.py
+      :language: py
+      :start-after: [sepimage-5-beg]
+      :end-before: [sepimage-5-end]
+      :dedent: 4
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src.cpp
+      :language: c++
+      :start-after: [sepimage-5-beg]
+      :end-before: [sepimage-5-end]
+      :dedent: 4
+
+.. figure:: ./sep_image_fine.png
+  :width: 400px
+  :align: center
+
+A more complex example is available on the public github repository. It treats the topic of the explored area, 
+which is a classical problem in robotics.
+
+* `Python version <https://github.com/codac-team/codac/blob/codac2/examples/04_explored_area/main_peibos.py>`_
+* `C++ version <https://github.com/codac-team/codac/blob/codac2/examples/04_explored_area/main_peibos.cpp>`_

doc/manual/manual/contractors/analytic/src.cpp

@@ -12,6 +12,7 @@
 #include <codac2_CtcFixpoint.h>
 #include <codac2_SepInverse.h>
 #include <codac2_CtcInverseNotIn.h>
+#include <codac2_SepImage.h>
 #include <codac2_Approx.h>
 #include <codac2_Figure2D.h>
 
@@ -86,4 +87,40 @@ TEST_CASE("CtcInverse - manual")
     CHECK(c.fnc().input_size() == 2);
     CHECK(c.fnc().output_size() == 1);
   }
+}
+
+TEST_CASE("SepImage - manual")
+{
+  {
+    // [sepimage-1-beg]
+    VectorVar y (2);
+    AnalyticFunction f ({y},{3.*(y[0]+1),y[1]+0.5*sin(3.*y[0])});
+    // [sepimage-1-end]
+
+    // [sepimage-2-beg]
+    // {psi0,Sigma} is a gnomonic atlas of the unit circle
+    VectorVar X(1);
+    AnalyticFunction psi0 ({X},{cos(X[0]*PI/2.),sin(X[0]*PI/2.)});
+
+    OctaSym id ({1, 2});
+    OctaSym s ({-1, -2});
+
+    vector<OctaSym> Sigma ({id,s});
+    // [sepimage-2-end]
+
+    // [sepimage-3-beg]
+    AnalyticFunction h ({y},sqrt(sqr(y[0])+sqr(y[1])));
+    CtcInverse ctc_in (h,Interval(0,1));    
+    // [sepimage-3-end]
+
+    // [sepimage-4-beg]
+    SepImage sep1 (f,psi0,Sigma,0.125,ctc_in);
+    DefaultFigure::pave({{-0.5,6.5},{-1.5,1.5}},sep1,0.05);
+    // [sepimage-4-end]
+
+    // [sepimage-5-beg]
+    SepImage sep2 (f,psi0,Sigma,0.0625,ctc_in);
+    DefaultFigure::pave({{-0.5,6.5},{-1.5,1.5}},sep2,0.05);
+    // [sepimage-5-end]
+  }
 }

doc/manual/manual/contractors/analytic/src.py

@@ -106,5 +106,38 @@ def tests_CtcInverse_manual(test):
     assert c.fnc().output_size() == 1
     # [ctcinv-7-end]
 
+def tests_SepImage_manual(test):
+
+    # [sepimage-1-beg]
+    y = VectorVar(2)
+    f = AnalyticFunction([y], [3.*(y[0]+1),y[1]+0.5*sin(3.*y[0])])
+    # [sepimage-1-end]
+
+    # [sepimage-2-beg]
+    # {psi0,Sigma} is a gnomonic atlas of the box [-1,1]^2
+    X = VectorVar(1)
+    psi0 = AnalyticFunction([X],[cos(X[0]*PI/2.),sin(X[0]*PI/2.)])
+
+    id = OctaSym([1,2])
+    s = OctaSym([-1,-2])
+
+    Sigma = [id,s]
+    # [sepimage-2-end]
+
+    # [sepimage-3-beg]
+    h = AnalyticFunction([y],sqrt(sqr(y[0])+sqr(y[1])))
+    ctc_in = CtcInverse(h,Interval(0,1))
+    # [sepimage-3-end]
+
+    # [sepimage-4-beg]
+    sep = SepImage(f,psi0,Sigma,0.125,ctc_in)
+    DefaultFigure.pave([[-0.5,6.5],[-1.5,1.5]],sep,0.05)
+    # [sepimage-4-end]
+
+    # [sepimage-5-beg]
+    sep = SepImage(f,psi0,Sigma,0.0625,ctc_in)
+    DefaultFigure.pave([[-0.5,6.5],[-1.5,1.5]],sep,0.05)
+    # [sepimage-5-end]
+
 if __name__ ==  '__main__':
   unittest.main()

doc/manual/manual/contractors/index.rst

@@ -5,6 +5,7 @@ Contractors, separators
 .. toctree::
 
   CtcInter <set/ctcinter>
+  SepImage <analytic/sepimage>
   CtcInverse <analytic/ctcinverse>
   CtcLohner <dynamic/ctclohner>
   CtcDist <geometric/ctcdist>
@@ -77,6 +78,9 @@ Overview of contractors and separators
    * - ``CtcAction``
      - ``SepAction``
 
+   * - —
+     - :ref:`SepImage <sec-ctc-analytic-sepimage>`
+
    * - :ref:`CtcInverse <sec-ctc-analytic-ctcinverse>`
      - ``SepInverse``
 

doc/manual/manual/extensions/capd/peibos_capd.rst

@@ -15,7 +15,7 @@ Considering a dynamical system :math:`\dot{\mathbf{x}}=\gamma(\mathbf{x})`, the
 Gnomonic atlas
 --------------
 
-To handle the boundary of the initial set :math:`\mathbb{X}_0`, the PEIBOS tool relies on a gnomonic atlas. See :ref:`subsec-functions-peibos-gnomonic-atals`.
+To handle the boundary of the initial set :math:`\mathbb{X}_0`, the PEIBOS tool relies on a gnomonic atlas. See :ref:`subsec-functions-peibos-gnomonic-atlas`.
 
 Use
 ---

doc/manual/manual/functions/parallelepiped/parallelepiped_eval.rst

@@ -9,7 +9,7 @@ Parallelepiped inclusion function
 Use case
 --------
 
-Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to\mathbb{R}^m`. In the case where :math:`0<n<m`,
+Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to\mathbb{R}^m`. In the case where :math:`0<n\leq m`,
 a parallelepiped inclusion function is available in the library.
 
 A parallelepiped inclusion function of :math:`\mathbf{f}` is noted :math:`\langle\mathbf{f}\rangle`.

doc/manual/manual/functions/peibos/peibos.rst

@@ -12,7 +12,7 @@ The PEIBOS tool provides a way to compute the Parallelepipedic Enclosure of the
 Let us consider an initial set :math:`\mathbb{X}_0 \subset \mathbb{R}^n` with its boundary :math:`\partial \mathbb{X}_0`. 
 Considering a function :math:`\mathbf{f}:\mathbb{R}^n \to \mathbb{R}^p`, :math:`n \leq p`, the PEIBOS tool allows to compute the set :math:`\mathbf{Y}=\left\{ \mathbf{f}(\mathbf{x}) \mid \mathbf{x} \in \partial \mathbb{X}_0 \right\}`.
 
-.. _subsec-functions-peibos-gnomonic-atals:
+.. _subsec-functions-peibos-gnomonic-atlas:
 
 Gnomonic atlas
 --------------

examples/04_explored_area/CMakeLists.txt

@@ -31,4 +31,9 @@
   add_executable(${PROJECT_NAME} main.cpp)
   target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS})
   target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
-  target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES})
+  target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES})
+
+  add_executable(main_peibos main_peibos.cpp)
+  target_compile_options(main_peibos PUBLIC ${CODAC_CXX_FLAGS})
+  target_include_directories(main_peibos SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
+  target_link_libraries(main_peibos PUBLIC ${CODAC_LIBRARIES})

examples/04_explored_area/main_peibos.cpp

@@ -0,0 +1,45 @@
+#include <codac>
+
+using namespace std;
+using namespace codac2;
+
+int main()
+{
+  set_nb_threads(max_threads());
+
+  // {psi0,Sigma} is a gnomonic atlas of the box [-1,1]^2
+  VectorVar X(1);
+  AnalyticFunction psi0 ({X},{X[0],1});
+
+  OctaSym id ({1, 2});
+  OctaSym s ({-2, 1});
+
+  vector<OctaSym> Sigma ({id,s,s*s,s.invert()});
+
+  // In this box, we consider that the x-axis is the width of the linear sensor
+  // and the y-axis is the time
+  VectorVar y (2);
+  auto L = 0.1*y[0];
+  auto t = 1.2*y[1];
+
+  // We construct the trajectory of the robot
+  auto traj  = vec(pow(t,3)-t, 1-sqr(t));
+  // We need its derivative to compute the orthogonal to the trajectory (for the sensor)
+  auto dtraj =  vec(3*sqr(t)-1, -2*t);
+  auto dtraj_norm = sqrt(sqr(dtraj[0])+sqr(dtraj[1]));
+
+  // The image of the box [-1,1]^2 by f is the swept area
+  AnalyticFunction f ({y},{traj[0]-L*dtraj[1]/dtraj_norm, traj[1]+L*dtraj[0]/dtraj_norm});
+
+  // For the SepImage, we need a contractor on the initial set (here a simple box)
+  IntervalVector X0 = IntervalVector::constant(2,{-1,1});
+  CtcWrapper ctc_in (X0);
+
+  // Separator on the area seen by a robot
+  SepImage sep (f,psi0,Sigma,0.01,ctc_in);
+
+  // Visualizing the separator
+  IntervalVector Y0 ({{-0.8,0.8},{-0.7,1.3}});
+
+  DefaultFigure::pave(Y0,sep,0.01);
+}