diff --git a/apps/attributes/gradient_magnitude.hpp b/apps/attributes/gradient_magnitude.hpp
index 35baa969b6796e54a3515d3d690eb37fd4dd6243..3b5953e995c7ff78acb9a1702e2f9ac72c9dc6ca 100644
--- a/apps/attributes/gradient_magnitude.hpp
+++ b/apps/attributes/gradient_magnitude.hpp
@@ -3,8 +3,8 @@
 #include <apps/tos/croutines.hpp>
 
 #include <mln/accu/accumulators/mean.hpp>
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/image/image2d.hpp>
-#include <mln/core/math_ops.hpp>
 #include <mln/morpho/component_tree/component_tree.hpp>
 
 
diff --git a/apps/attributes/meaningfullness.hpp b/apps/attributes/meaningfullness.hpp
index 41d47976afe8be1656052698f31bdc898da19f26..5e0793a958aca5ff332f2c6f73face394135855e 100644
--- a/apps/attributes/meaningfullness.hpp
+++ b/apps/attributes/meaningfullness.hpp
@@ -4,11 +4,11 @@
 
 #include <mln/core/colors.hpp>
 #include <mln/core/extension/fill.hpp>
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/image/image2d.hpp>
-#include <mln/core/math_ops.hpp>
 #include <mln/core/neighb2d.hpp>
 #include <mln/core/trace.hpp>
-#include <mln/core/vec/vec_math_ops.hpp>
+#include <mln/core/vec/vec_functional_ops.hpp>
 #include <mln/core/wrt_offset.hpp>
 
 
diff --git a/doc/source/core.rst b/doc/source/core.rst
index 319d69010938809f0fbf771a4943b877b1486f92..7cbd3a086c96006fc450c85bbfc4fb0607386571 100644
--- a/doc/source/core.rst
+++ b/doc/source/core.rst
@@ -12,6 +12,7 @@ Core Module
    core/ranges
    core/neighborhood
    core/stl
+   core/functional
 
 
 
@@ -63,18 +64,18 @@ See :doc:`core/images` for a description of the image concepts and image basics.
     :widths: auto
     :class: full
 
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
+    | :cpp:func:`view::transform(ima, f) <mln::view::transform>`    | Views the image with a function applied to the values.                                     |
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
+    | :cpp:func:`view::filter(ima, pred) <mln::view::filter>`       | Views the image restricted to pixels whose values pass a predicate.                        |
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
+    | :cpp:func:`view::mask(ima, mask) <mln::view::mask>`           | Views the image restricted to pixels in a binary mask.                                     |
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
+    | :cpp:func:`view::clip(ima, roi) <mln::view::clip>`            | Views the image restricted to a sub-region.                                                |
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
+    | :cpp:func:`view::zip(ima1, ima2, ..., imaN) <mln::view::zip>` | Views a list of images as a single image whose values are tuples of the each image values. |
+    +---------------------------------------------------------------+--------------------------------------------------------------------------------------------+
 
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
-    | :cpp:func:`view::transform(ima, f) <mln::view::transform>` | Views the image with a function applied to the values.                                     |
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
-    | ``view::filter(ima, pred)``                                | Views the image restricted to pixels whose values pass a predicate.                        |
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
-    | ``view::mask(ima, mask)``                                  | Views the image restricted to pixels in a binary mask.                                     |
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
-    | ``view::clip(ima, roi)``                                   | Views the image restricted to a sub-region.                                                |
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
-    | ``view::zip(ima1, ima2, ..., imaN)``                       | Views a list of images as a single image whose values are tuples of the each image values. |
-    +------------------------------------------------------------+--------------------------------------------------------------------------------------------+
 
 .. topic:: Common image operators (views)
 
@@ -82,22 +83,33 @@ See :doc:`core/images` for a description of the image concepts and image basics.
     :widths: auto
     :class: full
 
-    +--------------+-----------------------------------------------------------+
-    | Arithmetical | `a + b`, `a - b`, `a * b`, `a / b`, `a % b`               |
-    +--------------+-----------------------------------------------------------+
-    | Logical      | `a && b`, `a || b`, `! a`                                 |
-    +--------------+-----------------------------------------------------------+
-    | Comparison   | `a < b`, `a <= b`, `a == b`, `a != b`, `a >= b`, `a > b`  |
-    +--------------+-----------------------------------------------------------+
-    | Conditional  | ``view::ifelse(a, b, c)``                                 |
-    +--------------+-----------------------------------------------------------+
+    +--------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Arithmetical | :cpp:func:`a + b <mln::view::ops::operator+>`, :cpp:func:`a - b <mln::view::ops::operator->`, :cpp:func:`a * b <mln::view::ops::operator*>`, :cpp:func:`a / b <mln::view::ops::operator/>`, :cpp:func:`a % b <mln::view::ops::operator%>`                                                                                                                    |
+    +--------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Logical      | :cpp:func:`a && b <mln::view::ops::operator&&>`, :cpp:func:`a || b <mln::view::ops::operator||>`, :cpp:func:`\!a <mln::view::ops::operator!>`                                                                                                                                                                                                                |
+    +--------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Comparison   | :cpp:func:`a \< b <mln::view::ops::operator<>`, :cpp:func:`a \<= b <mln::view::ops::operator<=>`, :cpp:func:`a == b <mln::view::ops::operator==>`, :cpp:func:`a != b <mln::view::ops::operator!=>`, :cpp:func:`a >= b <mln::view::ops::operator>=>`, :cpp:func:`a > b <mln::view::ops::operator>>`, :cpp:func:`equalFP(a, b, eps) <mln::view::ops::equalFP>` |
+    +--------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Conditional  | :cpp:func:`view::ifelse(a, b, c) <mln::view::ifelse>`                                                                                                                                                                                                                                                                                                        |
+    +--------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 
 
 .. topic:: Mathematical image operator (views)
 
-  +-------+-------+
-  | FIXME | FIXME |
-  +-------+-------+
+  .. table::
+    :widths: auto
+    :class: full
+
+    +----------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Component-wise | :cpp:func:`abs(a) <mln::view::maths::abs>`, :cpp:func:`sqr(a) <mln::view::maths::sqr>`, :cpp:func:`pow(a, p) <mln::view::maths::pow>`, :cpp:func:`cbrt(a) <mln::view::maths::cbrt>`, :cpp:func:`sqrt(a) <mln::view::maths::sqrt>`                                                                                                                                     |
+    +----------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Reduction      | :cpp:func:`sum(a) <mln::view::maths::sum>`, :cpp:func:`prod(a) <mln::view::maths::prod>`, :cpp:func:`min(a) <mln::view::maths::min>`, :cpp:func:`max(a) <mln::view::maths::max>`, :cpp:func:`dot(a, b) <mln::view::maths::dot>`, :cpp:func:`cross(a, b) <mln::view::maths::cross>`                                                                                    |
+    +----------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Norms          | :cpp:func:`l0norm(a) <mln::view::maths::l0norm>`, :cpp:func:`l1norm(a) <mln::view::maths::l1norm>`, :cpp:func:`l2norm(a) <mln::view::maths::l2norm>`, :cpp:func:`l2norm_sqr(a) <mln::view::maths::l2norm_sqr>`, :cpp:func:`linfnorm(a) <mln::view::maths::linfnorm>`, :cpp:func:`lpnorm\<p>(a) <template <unsigned V> mln::view::maths::lpnorm<V>>`                   |
+    +----------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+    | Distances      | :cpp:func:`l0dist(a, b) <mln::view::maths::l0dist>`, :cpp:func:`l1dist(a, b) <mln::view::maths::l1dist>`, :cpp:func:`l2dist(a, b) <mln::view::maths::l2dist>`, :cpp:func:`l2dist_sqr(a, b) <mln::view::maths::l2dist_sqr>`, :cpp:func:`linfdist(a, b) <mln::view::maths::linfdist>`, :cpp:func:`lpdist\<p>(a, b) <template <unsigned V> mln::view::maths::lpdist<V>>` |
+    +----------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+
 
 .. topic:: Other image operators (views)
 
@@ -105,15 +117,15 @@ See :doc:`core/images` for a description of the image concepts and image basics.
     :widths: auto
     :class: full
     
-    +---------------------------+---------------------------------------------------+
-    | ``view::cast<V>(ima)``    | Views the image with the values casted into ``V`` |
-    +---------------------------+---------------------------------------------------+
-    | ``view::channel(ima, k)`` | Views the k-th channel of the image               |
-    +---------------------------+---------------------------------------------------+
-    | ``view::red(ima)``        |                                                   |
-    | ``view::green(ima)``      |                                                   |
-    | ``view::blue(ima)``       | Views the red, green or blue channel of the image |
-    +---------------------------+---------------------------------------------------+
+    +----------------------------------------------------------------------------+---------------------------------------------------+
+    | :cpp:func:`view::cast\<T>(ima) <template <typename T> mln::view::cast<T>>` | Views the image with the values casted into ``T`` |
+    +----------------------------------------------------------------------------+---------------------------------------------------+
+    | :cpp:func:`view::channel(ima, k) <mln::view::channel>`                     | Views the k-th channel of the image               |
+    +----------------------------------------------------------------------------+---------------------------------------------------+
+    | :cpp:func:`view::red(ima) <mln::view::red>`                                |                                                   |
+    | :cpp:func:`view::green(ima) <mln::view::green>`                            | Views the red, green or blue channel of the image |
+    | :cpp:func:`view::blue(ima) <mln::view::blue>`                              |                                                   |
+    +----------------------------------------------------------------------------+---------------------------------------------------+
 
 
 .. raw:: html
diff --git a/doc/source/core/functional.rst b/doc/source/core/functional.rst
new file mode 100644
index 0000000000000000000000000000000000000000..061b64de45f0757f43dbd07576534f76fd23f85b
--- /dev/null
+++ b/doc/source/core/functional.rst
@@ -0,0 +1,438 @@
+Functional maths operators
+##########################
+
+
+Include :file:`<mln/core/functional_ops.hpp>`
+
+.. cpp:namespace:: mln::functional
+
+Component-wise
+==============
+
+
+1. .. cpp:function:: auto abs(Scalar s)
+2. .. cpp:function:: auto abs(Vector v)
+
+    1. Returns the absolute value of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the absolute value of each component of ``v``: :math:`r_i = \lvert v_i \rvert`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = -2;
+        mln::signed_rgb8 v = {1, -2, 3};
+        auto r_s = abs(s); // 2
+        auto r_v = abs(v); // {1, 2, 3}
+
+
+1. .. cpp:function:: auto pow(Scalar s, Scalar p)
+2. .. cpp:function:: auto pow(Vector v, Scalar p)
+
+    1. Returns ``s`` exponent ``p``: :math:`r = s^p`.
+    2. Returns each component of ``v`` exponent ``p``: :math:`r_i = v_i^p`.
+
+    :param s: Input scalar
+    :param v: Input vector
+    :param p: Input exponent
+
+    .. code::
+
+        mln::int8 s = -2;
+        mln::signed_rgb8 v = {1, -2, 3};
+        auto r_s = pow(s, 3); // -8
+        auto r_v = pow(v, 3); // {1, -8, 27}
+
+
+1. .. cpp:function:: auto sqr(Scalar s)
+2. .. cpp:function:: auto sqr(Vector v)
+
+    1. Returns the square value of ``s``: :math:`r = s^2`.
+    2. Returns the square value of each component of ``v``: :math:`r_i = v_i^2`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = -2;
+        mln::signed_rgb8 v = {1, -2, 3};
+        auto r_s = sqr(s); // 4
+        auto r_v = sqr(v); // {1, 4, 9}
+
+
+1. .. cpp:function:: auto cbrt(Scalar s)
+2. .. cpp:function:: auto cbrt(Vector v)
+
+    1. Returns the cubic root of ``s``: :math:`r = \sqrt[3]{s}`.
+    2. Returns the cubic root of each component of ``v``: :math:`r_i = \sqrt[3]{v_i}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 729;
+        mln::signed_rgb8 v = {-0.125, 8, 729};
+        auto r_s = cbrt(s); // 9
+        auto r_v = cbrt(v); // {-0.5, 2, 9}
+
+
+1. .. cpp:function:: auto sqrt(Scalar s)
+2. .. cpp:function:: auto sqrt(Vector v)
+
+    1. Returns the square root of ``s``: :math:`r = \sqrt{s}`.
+    2. Returns the square root of each component of ``v``: :math:`r_i = \sqrt{v_i}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 81;
+        mln::signed_rgb8 v = {-0.25, 4, 9};
+        auto r_s = sqrt(s); // 9
+        auto r_v = sqrt(v); // {-0.5, 2, 3}
+
+
+Reduction
+=========
+
+
+1. .. cpp:function:: auto sum(Scalar s)
+2. .. cpp:function:: auto sum(Vector v)
+
+    1. Returns ``s``: :math:`r = s`.
+    2. Returns the sum of all components of ``v``: :math:`r = \sum{v_i}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = sum(s); // 9
+        auto r_v = sum(v); // 8
+
+
+1. .. cpp:function:: auto prod(Scalar s)
+2. .. cpp:function:: auto prod(Vector v)
+
+    1. Returns ``s``: :math:`r = s`.
+    2. Returns the product of all components of ``v``: :math:`r = \prod{v_i}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = prod(s); // 9
+        auto r_v = prod(v); // -18
+
+
+1. .. cpp:function:: auto min(Scalar s)
+2. .. cpp:function:: auto min(Vector v)
+
+    1. Returns ``s``: :math:`r = s`.
+    2. Returns the minimum of all components of ``v``: :math:`r = \min(v_1, v_2, \dots, v_n)`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = min(s); // 9
+        auto r_v = min(v); // -2
+
+
+1. .. cpp:function:: auto max(Scalar s)
+2. .. cpp:function:: auto max(Vector v)
+
+    1. Returns ``s``: :math:`r = s`.
+    2. Returns the maximum of all components of ``v``: :math:`r = \max(v_1, v_2, \dots, v_n)`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = max(s); // 9
+        auto r_v = max(v); // 9
+
+
+1. .. cpp:function:: auto dot(Vector v1, Vector v2)
+
+    2. Returns the scalar product of ``v1`` by ``v2``: :math:`r = v_1 \cdot v_2`.
+
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::signed_rgb8 v1= {9, 1, -2};
+        mln::signed_rgb8 v2 = {1, -2, 9};
+        auto v = dot(v1, v2); // -11
+
+
+1. .. cpp:function:: auto cross(Vector v1, Vector v2)
+
+    2. Returns the cross product of ``v1`` by ``v2``: :math:`r = v_1 \wedge v_2`.
+
+
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::signed_rgb8 v1= {9, 1, -2};
+        mln::signed_rgb8 v2 = {1, -2, 9};
+        auto v = cross(v1, v2); // {5, -83, -19}
+
+
+    .. warning:: The cross product is only defined for vectors of dimension 3 and 7.
+
+
+Norms
+=====
+
+
+1. .. cpp:function:: auto l0norm(Scalar s)
+2. .. cpp:function:: auto l0norm(Vector v)
+
+    1. Returns the absolute value norm of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the absolute value (min) norm of ``v``: :math:`r = \lVert v \rVert_0 = \min(\lvert v_1 \rvert, \lvert v_2 \rvert, \dots, \lvert v_n \rvert)`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = l0norm(s); // 9
+        auto r_v = l0norm(v); // 1
+
+
+1. .. cpp:function:: auto l1norm(Scalar s)
+2. .. cpp:function:: auto l1norm(Vector v)
+
+    1. Returns the absolute value norm of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the absolute value (sum) norm of ``v``: :math:`r = \lVert v \rVert_1 = \sum{\lvert v_i \rvert}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = l1norm(s); // 9
+        auto r_v = l1norm(v); // 12
+
+
+1. .. cpp:function:: auto l2norm(Scalar s)
+2. .. cpp:function:: auto l2norm(Vector v)
+
+    1. Returns the euclidean norm of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the euclidean norm of ``v``: :math:`r = \lVert v \rVert_2 = \sqrt{\sum{v_i^2}}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {2, -4, 4};
+        auto r_s = l2norm(s); // 9
+        auto r_v = l2norm(v); // 6
+
+
+1. .. cpp:function:: auto l2norm_sqr(Scalar s)
+2. .. cpp:function:: auto l2norm_sqr(Vector v)
+
+    1. Returns the squared euclidean norm of ``s``: :math:`r = \lvert s \rvert^2`.
+    2. Returns the squared euclidean norm of ``v``: :math:`r = \lVert v \rVert_2^2 = \sum{v_i^2}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {2, -4, 4};
+        auto r_s = l2norm_sqr(s); // 81
+        auto r_v = l2norm_sqr(v); // 36
+
+
+1. .. cpp:function:: auto linfnorm(Scalar s)
+2. .. cpp:function:: auto linfnorm(Vector v)
+
+    1. Returns the absolute value norm of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the absolute value (max) norm of ``v``: :math:`r = \lVert v \rVert_{inf} = \max(\lvert v_1 \rvert, \lvert v_2 \rvert, \dots, \lvert v_n \rvert)`.
+
+    :param s: Input scalar
+    :param v: Input vector
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {1, -2, 9};
+        auto r_s = linfnorm(s); // 9
+        auto r_v = linfnorm(v); // 9
+
+
+1. .. cpp:function:: template <unsigned p> auto lpnorm<p>(Scalar s)
+2. .. cpp:function:: template <unsigned p> auto lpnorm<p>(Vector v)
+
+    1. Returns the p-euclidean norm of ``s``: :math:`r = \lvert s \rvert`.
+    2. Returns the p-euclidean norm of ``v``: :math:`r = \lVert v \rVert_p = \sqrt[p]{\sum{\lvert v_i \rvert^p}}`.
+
+    :param s: Input scalar
+    :param v: Input vector
+    :tparam p: Norm exponent
+
+    .. code::
+
+        mln::int8 s = 9;
+        mln::signed_rgb8 v = {2, -4, 9};
+        auto r_s = lpnorm<4>(s); // 9
+        auto r_v = lpnorm<4>(v); // 9.09186...
+
+
+Distances
+=========
+
+
+1. .. cpp:function:: auto l0dist(Scalar s1, Scalar s2)
+2. .. cpp:function:: auto l0dist(Vector v1, Vector v2)
+
+    1. Returns the distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert`.
+    2. Returns the distance normalized (min of absolute values) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_0`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = l0dist(s1, s2); // 11
+        auto d_v = l0dist(v1, v2); // 2
+
+
+1. .. cpp:function:: auto l1dist(Scalar s1, Scalar s2)
+2. .. cpp:function:: auto l1dist(Vector v1, Vector v2)
+
+    1. Returns the distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert`.
+    2. Returns the distance normalized (sum of absolute values) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_1`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = l1dist(s1, s2); // 11
+        auto d_v = l1dist(v1, v2); // 10
+
+
+1. .. cpp:function:: auto l2dist(Scalar s1, Scalar s2)
+2. .. cpp:function:: auto l2dist(Vector v1, Vector v2)
+
+    1. Returns the distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert`.
+    2. Returns the distance normalized (euclidean) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_2`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = l2dist(s1, s2); // 11
+        auto d_v = l2dist(v1, v2); // 6
+
+
+1. .. cpp:function:: auto l2dist_sqr(Scalar s1, Scalar s2)
+2. .. cpp:function:: auto l2dist_sqr(Vector v1, Vector v2)
+
+    1. Returns the squared distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert^2`.
+    2. Returns the squared distance normalized (euclidean) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_2^2`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = l2dist_sqr(s1, s2); // 121
+        auto d_v = l2dist_sqr(v1, v2); // 36
+
+
+1. .. cpp:function:: auto linfdist(Scalar s1, Scalar s2)
+2. .. cpp:function:: auto linfdist(Vector v1, Vector v2)
+
+    1. Returns the distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert`.
+    2. Returns the distance normalized (max of absolute values) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_{inf}`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+
+    .. code::
+
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = linfdist(s1, s2); // 11
+        auto d_v = linfdist(v1, v2); // 4
+
+
+1. .. cpp:function:: template <unsigned p> auto lpdist<p>(Scalar s1, Scalar s2)
+2. .. cpp:function:: template <unsigned p> auto lpdist<p>(Vector v1, Vector v2)
+
+    1. Returns the distance between ``s1`` and ``s2``: :math:`d = \lvert s_1 - s_2 \rvert`.
+    2. Returns the distance normalized (p-euclidean) between ``v1`` and ``v2``: :math:`d = \lVert v_1 - v_2 \rVert_p`.
+
+    :param s1: Input scalar
+    :param s2: Input scalar
+    :param v1: Input vector
+    :param v2: Input vector
+    :tparam p: Norm Exponent
+
+    .. code::
+    
+        mln::int8 s1 = 9;
+        mln::int8 s2 = -2;
+        mln::signed_rgb8 v1 = {1, 2, 2};
+        mln::signed_rgb8 v2 = {3, 6, 6};
+        auto d_s = lpdist<4>(s1, s2); // 11
+        auto d_v = lpdist<4>(v1, v2); // 4.79356
diff --git a/doc/source/core/images.rst b/doc/source/core/images.rst
index a0b355b865bd340c872a0bb58094c227ea01dbe1..7a4a30fe40ef64896b25a2dae29e7b39a2756403 100644
--- a/doc/source/core/images.rst
+++ b/doc/source/core/images.rst
@@ -8,8 +8,17 @@ Forsee to add pdim (point) and vdim (value) for dimension constant in iamge conc
 
 .. toctree::
    :hidden:
-
+   
+   view/cast
+   view/clip
+   view/channel
+   view/filter
+   view/mask
+   view/maths
+   view/operators
+   view/rgb
    view/transform
+   view/zip
 
 
 .. contents::
diff --git a/doc/source/core/view/cast.rst b/doc/source/core/view/cast.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c139c9221dc10c11d8b95065c4b0704ee7f12c66
--- /dev/null
+++ b/doc/source/core/view/cast.rst
@@ -0,0 +1,53 @@
+Cast (View)
+###########
+
+Include :file:`<mln/core/image/view/cast.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: template <typename T> auto cast<T>(Image ima)
+
+    1. Makes a view from `ima` where for each pixel value evals to :cpp:expr:`out(p) = static_cast<T>(ima(p))`
+    
+    :param ima: Input range
+    :tparam T: Type to cast to
+
+    .. code::
+    
+        mln::image2d<double> ima = ...;
+        auto g1 = mln::view::cast<int>(ima);
+
++----------------+---------------+-------+----------------+
+|                |               | *ima* | *cast<T>(ima)* |
++================+===============+=======+================+
+|                | Forward       | X     | X              |
++                +---------------+-------+----------------+
+| **Category**   | Bidirectional | X     | X              |
++                +---------------+-------+----------------+
+|                | Raw           | X     |                |
++----------------+---------------+-------+----------------+
+|                | Writable      | X     |                |
++                +---------------+-------+----------------+
+| **Properties** | Accessible    | X     | X              |
++                +---------------+-------+----------------+
+|                | Indexable     | X     | X              |
++----------------+---------------+-------+----------------+
+
+
+Examples
+--------
+   
+* Double image to int:
+
+.. code::
+
+    mln::image2d<double> input = { {0.1, 1.1, 2.1}, {3.5, 4.5, 5.5} }; 
+    auto g = mln::view::cast<int>(input);
+    mln::io::imprint(g);
+
+Outputs:
+
+.. code-block:: text
+
+    0 1 2
+    3 4 5
diff --git a/doc/source/core/view/channel.rst b/doc/source/core/view/channel.rst
new file mode 100644
index 0000000000000000000000000000000000000000..99d1f56a13ca5b5b0f5b14046cd1a62b2c373ef9
--- /dev/null
+++ b/doc/source/core/view/channel.rst
@@ -0,0 +1,59 @@
+Channel (View)
+##############
+
+Include :file:`<mln/core/image/view/rgb.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto channel(Image ima, unsigned k)
+
+    2. Makes a view that project the k-th channel of a multi-channelled colored image, when k is only known at runtime.
+    
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        auto ima_red = mln::view::channel(ima, 0);
+
+
++----------------+---------------+-------+-------------------+
+|                |               | *ima* | *channel(ima, k)* |
++================+===============+=======+===================+
+|                | Forward       | X     | X                 |
++                +---------------+-------+-------------------+
+| **Category**   | Bidirectional | X     | X                 |
++                +---------------+-------+-------------------+
+|                | Raw           | X     |                   |
++----------------+---------------+-------+-------------------+
+|                | Writable      | X     | X                 |
++                +---------------+-------+-------------------+
+| **Properties** | Accessible    | X     | X                 |
++                +---------------+-------+-------------------+
+|                | Indexable     | X     | X                 |
++----------------+---------------+-------+-------------------+
+
+
+Examples
+--------
+   
+* Set the red channel to max:
+
+.. code::
+
+    mln::image2d<mln::rgb8> ima = {
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} },
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} },
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} }
+      };
+    auto ima_red = mln::view::channel(ima, 0);
+    mln::fill(ima_red, 255);
+    mln::io::imprint(ima);
+
+Outputs:
+
+.. code-block:: text
+
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
diff --git a/doc/source/core/view/clip.rst b/doc/source/core/view/clip.rst
new file mode 100644
index 0000000000000000000000000000000000000000..ba3fbc6d896cf413d1a9841fb381622051ee2e84
--- /dev/null
+++ b/doc/source/core/view/clip.rst
@@ -0,0 +1,76 @@
+Clip (View)
+###########
+
+Include :file:`<mln/core/image/view/clip.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto clip(Image ima, Domain dom)
+
+    1. Makes a view that clip the original image to a sub-domain. A pixel is kept if and only if :cpp:expr:`dom.has(p.point())` is true.
+    
+    :param ima: Input range
+    :param dom: Input range of point
+
+    .. code::
+    
+        mln::image2d<int> ima = ...;
+        mln::image2d<bool> mask = ...;
+        auto g = mln::view::mask(ima, mask);
+
+
++----------------+---------------+-------+------------------+
+|                |               | *ima* | *clip(ima, dom)* |
++================+===============+=======+==================+
+|                | Forward       | X     | X                |
++                +---------------+-------+------------------+
+| **Category**   | Bidirectional | X     | X                |
++                +---------------+-------+------------------+
+|                | Raw           | X     |                  |
++----------------+---------------+-------+------------------+
+|                | Writable      | X     | X                |
++                +---------------+-------+------------------+
+| **Properties** | Accessible    | X     | X                |
++                +---------------+-------+------------------+
+|                | Indexable     | X     | X                |
++----------------+---------------+-------+------------------+
+
+
+Examples
+--------
+   
+* Set a sub-domain to 42:
+
+.. code::
+
+    mln::image2d<int> ima = { {0, 1, 2}, {3, 4, 5}, {6, 7, 8} };
+    box2d d = {{0, 0}, {1, 1}};
+    auto g = mln::view::clip(ima, d);
+    mln::fill(g, 42);
+    mln::io::imprint(ima);
+
+Outputs:
+
+.. code-block:: text
+
+    42 42 2
+    42 42 5
+    6  7  8
+
+* Set all odd values to 42:
+
+.. code::
+
+    mln::image2d<int> ima = { {0, 1, 2}, {3, 4, 5}, {6, 7, 8} };
+    auto odd_d = mln::where(ima % 2 == 1);
+    auto g = mln::view::clip(ima, odd_d);
+    mln::fill(g, 42);
+    mln::io::imprint(ima);
+
+Outputs:
+
+.. code-block:: text
+
+    0  42 2
+    42 4  42
+    6  42 8
diff --git a/doc/source/core/view/filter.rst b/doc/source/core/view/filter.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6aec691ce0ca685d8912d91d4e8cfe71690eee49
--- /dev/null
+++ b/doc/source/core/view/filter.rst
@@ -0,0 +1,54 @@
+Filter (View)
+#############
+
+Include :file:`<mln/core/image/view/filter.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto filter(Image ima, Predicate pred)
+
+    1. Makes a view where each pixel is kept if and only if the predicate :cpp:expr:`pred(ima(p))` is true.
+    
+    :param ima: Input range
+    :param pred: Predicate to apply on pixel values
+
+    .. code::
+    
+        mln::image2d<int> ima = ...;
+        auto g1 = mln::view::filter(ima, [](int v) { return v > 125; });
+
+
++----------------+---------------+-------+---------------------+
+|                |               | *ima* | *filter(ima, pred)* |
++================+===============+=======+=====================+
+|                | Forward       | X     | X                   |
++                +---------------+-------+---------------------+
+| **Category**   | Bidirectional | X     | X                   |
++                +---------------+-------+---------------------+
+|                | Raw           | X     |                     |
++----------------+---------------+-------+---------------------+
+|                | Writable      | X     |                     |
++                +---------------+-------+---------------------+
+| **Properties** | Accessible    | X     | X                   |
++                +---------------+-------+---------------------+
+|                | Indexable     | X     | X                   |
++----------------+---------------+-------+---------------------+
+
+
+Examples
+--------
+   
+* Keep all add values:
+
+.. code::
+
+    mln::image2d<int> input = { {0, 1, 2}, {3, 4, 5} };
+    auto g = mln::view::transform(input, [](int v) { return v % 2 != 0; });
+    mln::io::imprint(g);
+
+Outputs:
+
+.. code-block:: text
+
+    0 1
+    3 5
diff --git a/doc/source/core/view/mask.rst b/doc/source/core/view/mask.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8ea7a9ca0ac4f4df1166253aca47957f7cedfb63
--- /dev/null
+++ b/doc/source/core/view/mask.rst
@@ -0,0 +1,57 @@
+Mask (View)
+###########
+
+Include :file:`<mln/core/image/view/mask.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto mask(Image ima, BoolImage m)
+
+    1. Makes a view where each pixel is kept if and only if the masked boolean :cpp:expr:`m(p)` is true.
+    
+    :param ima: Input range
+    :param mask: Input range
+
+    .. code::
+    
+        mln::image2d<int> ima = ...;
+        mln::image2d<bool> m = ...;
+        auto g = mln::view::mask(ima, m);
+
+
++----------------+---------------+-------+-------------------+
+|                |               | *ima* | *mask(ima, mask)* |
++================+===============+=======+===================+
+|                | Forward       | X     | X                 |
++                +---------------+-------+-------------------+
+| **Category**   | Bidirectional | X     | X                 |
++                +---------------+-------+-------------------+
+|                | Raw           | X     |                   |
++----------------+---------------+-------+-------------------+
+|                | Writable      | X     | X                 |
++                +---------------+-------+-------------------+
+| **Properties** | Accessible    | X     | X                 |
++                +---------------+-------+-------------------+
+|                | Indexable     | X     | X                 |
++----------------+---------------+-------+-------------------+
+
+
+Examples
+--------
+   
+* Set all odd values to 42:
+
+.. code::
+
+    mln::image2d<int> ima = { {0, 1, 2}, {3, 4, 5} };
+    mln::image2d<bool> m = ima % 2 == 1;
+    auto g = mln::view::mask(ima, m);
+    mln::fill(g, 42);
+    mln::io::imprint(ima);
+
+Outputs:
+
+.. code-block:: text
+
+    0  42 2
+    42 4  42
diff --git a/doc/source/core/view/maths.rst b/doc/source/core/view/maths.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1b40cc6c037a920851ee82df352bfb36ca4a9147
--- /dev/null
+++ b/doc/source/core/view/maths.rst
@@ -0,0 +1,363 @@
+Maths (view)
+############
+
+Include :file:`<mln/core/image/view/maths.hpp>`
+
+.. cpp:namespace:: mln::view::maths
+
+Component-wise
+==============
+
+
+#. .. cpp:function:: auto abs(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = abs(ima(p))`.
+       Refer to :cpp:func:`functional::abs <mln::functional::abs>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::abs(ima);
+
+
+#. .. cpp:function:: auto pow(Image ima, Scalar s)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = pow(ima(p), s)`.
+       Refer to :cpp:func:`functional::pow <mln::functional::pow>` documentation for details.
+
+    :param ima: Input range
+    :param s: Exponent
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::pow(ima, 3);
+
+
+#. .. cpp:function:: auto sqr(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = sqr(ima(p))`.
+       Refer to :cpp:func:`functional::sqr <mln::functional::sqr>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::sqr(ima);
+
+
+#. .. cpp:function:: auto cbrt(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = cbrt(ima(p))`. 
+       Refer to :cpp:func:`functional::cbrt <mln::functional::cbrt>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::cbrt(ima);
+
+
+#. .. cpp:function:: auto sqrt(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = sqrt(ima(p))`. 
+       Refer to :cpp:func:`functional::sqrt <mln::functional::sqrt>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::sqrt(ima);
+
+
+Reduction
+=========
+
+
+#. .. cpp:function:: auto sum(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = sum(ima(p))`. 
+       Refer to :cpp:func:`functional::sum <mln::functional::sum>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::sum(ima);
+
+
+#. .. cpp:function:: auto prod(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = prod(ima(p))`. 
+       Refer to :cpp:func:`functional::prod <mln::functional::prod>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::prod(ima);
+
+
+#. .. cpp:function:: auto min(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = min(ima(p))`. 
+       Refer to :cpp:func:`functional::min <mln::functional::min>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::min(ima);
+
+
+#. .. cpp:function:: auto max(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = max(ima(p))`. 
+       Refer to :cpp:func:`functional::max <mln::functional::max>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::max(ima);
+
+
+#. .. cpp:function:: auto dot(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = dot(ima1(p), ima2(p))`. 
+       Refer to :cpp:func:`functional::dot <mln::functional::dot>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::dot(ima1, ima2);
+
+
+#. .. cpp:function:: auto cross(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = cross(ima1(p), ima2(p))`. 
+       Refer to :cpp:func:`functional::cross <mln::functional::cross>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::rgb8> g = mln::view::maths::cross(ima1, ima2);
+
+
+Norms
+=====
+
+
+#. .. cpp:function:: auto l0norm(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l0norm(ima(p))`. 
+       Refer to :cpp:func:`functional::l0norm <mln::functional::l0norm>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l0norm(ima);
+
+
+#. .. cpp:function:: auto l1norm(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l1norm(ima(p))`. 
+       Refer to :cpp:func:`functional::l1norm <mln::functional::l1norm>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l1norm(ima);
+
+
+#. .. cpp:function:: auto l2norm(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l2norm(ima(p))`. 
+       Refer to :cpp:func:`functional::l2norm <mln::functional::l2norm>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l2norm(ima);
+
+
+#. .. cpp:function:: auto l2norm_sqr(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l2norm_sqr(ima(p))`. 
+       Refer to :cpp:func:`functional::l2norm_sqr <mln::functional::l2norm_sqr>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l2norm_sqr(ima);
+
+
+#. .. cpp:function:: auto linfnorm(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = linfnorm(ima(p))`. 
+       Refer to :cpp:func:`functional::linfnorm <mln::functional::linfnorm>` documentation for details.
+
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::linfnorm(ima);
+
+
+#. .. cpp:function:: template <unsigned V> auto lpnorm<V>(Image ima)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = lpnorm<V>(ima(p))`. 
+       Refer to :cpp:func:`functional::lpnorm <template <unsigned p> mln::functional::lpnorm<p>>` documentation for details.
+
+    :param ima: Input range
+    :tparam V: Norm exponent
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::lpnorm<4>(ima);
+
+
+
+
+Distances
+=========
+
+
+#. .. cpp:function:: auto l0dist(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l0dist(ima1, ima2)`. 
+       Refer to :cpp:func:`functional::l0dist <mln::functional::l0dist>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l0dist(ima1, ima2);
+
+
+#. .. cpp:function:: auto l1dist(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l1dist(ima1, ima2)`. 
+       Refer to :cpp:func:`functional::l1dist <mln::functional::l1dist>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l1dist(ima1, ima2);
+
+
+#. .. cpp:function:: auto l2dist(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l2dist(ima1, ima2)`. 
+       Refer to :cpp:func:`functional::l2dist <mln::functional::l2dist>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l2dist(ima1, ima2);
+
+
+#. .. cpp:function:: auto l2dist_sqr(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = l2dist_sqr(ima1, ima2)`. 
+       Refer to :cpp:func:`functional::l2dist_sqr <mln::functional::l2dist_sqr>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::l2dist_sqr(ima1, ima2);
+
+
+
+#. .. cpp:function:: auto linfdist(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = linfdist(ima1, ima2)`. 
+       Refer to :cpp:func:`functional::linfdist <mln::functional::linfdist>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::linfdist(ima1, ima2);
+
+
+#. .. cpp:function:: template <unsigned V> auto lpdist<V>(Image ima1, Image ima2)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = lpdist<V>(ima1(p), ima2(p))`. 
+       Refer to :cpp:func:`functional::lpdist <template <unsigned p> mln::functional::lpdist<p>>` documentation for details.
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :tparam V: Norm Exponent
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        mln::image2d<mln::uint8> g = mln::view::maths::lpdist<4>(ima1, ima2);
+
+
++----------------+---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+|                |               | *ima, ima1, ima2* | *abs, pow, sqr, cbrt, sqrt, sum, prod, min, max, dot, cross, l0norm, l1norm, l2norm, l2norm_sqr, linfnorm, lpnorm, l0dist, l1dist, l2dist, l2dist_sqr, linfdist, lpdist* |
++================+===============+===================+==========================================================================================================================================================================+
+|                | Forward       | X                 | X                                                                                                                                                                        |
++                +---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+| **Category**   | Bidirectional | X                 | X                                                                                                                                                                        |
++                +---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+|                | Raw           | X                 |                                                                                                                                                                          |
++----------------+---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+|                | Writable      | X                 |                                                                                                                                                                          |
++                +---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+| **Properties** | Accessible    | X                 | X                                                                                                                                                                        |
++                +---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+|                | Indexable     | X                 | X                                                                                                                                                                        |
++----------------+---------------+-------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
+
+
+Examples
+--------
+
+TODO
diff --git a/doc/source/core/view/operators.rst b/doc/source/core/view/operators.rst
new file mode 100644
index 0000000000000000000000000000000000000000..4a95aea3c9d98af413818732ced85a35e04816da
--- /dev/null
+++ b/doc/source/core/view/operators.rst
@@ -0,0 +1,429 @@
+Operators (View)
+################
+
+Include :file:`<mln/core/image/view/operators.hpp>`
+
+.. cpp:namespace:: mln::view::ops
+
+
+Arithmetical
+============
+
+
+1. .. cpp:function:: auto operator+(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator+(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator+(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) + ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) + s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s + ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 + ima2;
+        double s = 5.5;
+        auto g1 = ima1 + s;
+        auto g1 = s + ima1;
+
+
+1. .. cpp:function:: auto operator-(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator-(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator-(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) - ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) - s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s - ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 - ima2;
+        double s = 5.5;
+        auto g1 = ima1 - s;
+        auto g1 = s - ima1;
+
+
+1. .. cpp:function:: auto operator*(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator*(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator*(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) * ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) * s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s * ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 * ima2;
+        double s = 5.5;
+        auto g1 = ima1 * s;
+        auto g1 = s * ima1;
+
+
+1. .. cpp:function:: auto operator/(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator/(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator/(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) / ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) / s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s / ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. warning:: Mind the division by zero.
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 / ima2;
+        double s = 5.5;
+        auto g1 = ima1 / s;
+        auto g1 = s / ima1;
+
+
+1. .. cpp:function:: auto operator%(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator%(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator%(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) % ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) % s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s % ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 % ima2;
+        double s = 5.5;
+        auto g1 = ima1 % s;
+        auto g1 = s % ima1;
+
+
+Logical
+=======
+
+
+1. .. cpp:function:: auto operator&&(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator&&(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator&&(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) && ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) && s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s && ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 && ima2;
+        double s = 5.5;
+        auto g1 = ima1 && s;
+        auto g1 = s && ima1;
+
+
+1. .. cpp:function:: auto operator||(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator||(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator||(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) || ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) || s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s || ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 || ima2;
+        double s = 5.5;
+        auto g1 = ima1 || s;
+        auto g1 = s || ima1;
+
+
+1. .. cpp:function:: auto operator!(Image ima)
+
+    1. Makes a view where for each pixel value evals to `out(p) = !ima(p)`
+
+
+    :param ima: Input range
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        auto g = !ima;
+
+
+Comparison
+==========
+
+
+1. .. cpp:function:: auto operator<(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator<(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator<(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) < ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) < s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s < ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 < ima2;
+        double s = 5.5;
+        auto g1 = ima1 < s;
+        auto g1 = s < ima1;
+
+
+1. .. cpp:function:: auto operator<=(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator<=(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator<=(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) <= ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) <= s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s <= ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 <= ima2;
+        double s = 5.5;
+        auto g1 = ima1 <= s;
+        auto g1 = s <= ima1;
+
+
+1. .. cpp:function:: auto operator==(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator==(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator==(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) == ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) == s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s == ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 == ima2;
+        double s = 5.5;
+        auto g1 = ima1 == s;
+        auto g1 = s == ima1;
+
+
+1. .. cpp:function:: auto equalFP(Image ima1, Image ima2, double eps)
+2. .. cpp:function:: auto equalFP(Image ima1, Scalar s, double eps)
+3. .. cpp:function:: auto equalFP(Scalar s, Image ima1, double eps)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = abs(ima1(p) - ima2(p)) < eps`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = abs(ima1(p) - s) < eps`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = abs(s - ima1(p)) < eps`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+    :param eps: Epsilon precision
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = equalFP(ima1, ima2, 10e-6);
+        double s = 5.5;
+        auto g1 = equalFP(ima1, s, 10e-6);
+        auto g1 = equalFP(s, ima1, 10e-6);
+
+
+1. .. cpp:function:: auto operator!=(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator!=(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator!=(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) != ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) != s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s != ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 != ima2;
+        double s = 5.5;
+        auto g1 = ima1 != s;
+        auto g1 = s != ima1;
+
+
+#. .. cpp:function:: auto operator>(Image ima1, Image ima2)
+#. .. cpp:function:: auto operator>(Image ima1, Scalar s)
+#. .. cpp:function:: auto operator>(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) > ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) > s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s > ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 > ima2;
+        double s = 5.5;
+        auto g1 = ima1 > s;
+        auto g1 = s > ima1;
+
+
+1. .. cpp:function:: auto operator>=(Image ima1, Image ima2)
+2. .. cpp:function:: auto operator>=(Image ima1, Scalar s)
+3. .. cpp:function:: auto operator>=(Scalar s, Image ima1)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) >= ima2(p)`
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima1(p) >= s`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = s >= ima1(p)`.
+
+
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima1 = ...;
+        mln::image2d<mln::rgb8> ima2 = ...;
+        auto g1 = ima1 >= ima2;
+        double s = 5.5;
+        auto g1 = ima1 >= s;
+        auto g1 = s >= ima1;
+
+
+.. cpp:namespace:: mln::view
+
+Conditional
+===========
+
+
+1. .. cpp:function:: auto ifelse(Image ima, Image ima1, Image ima2)
+2. .. cpp:function:: auto ifelse(Image ima, Image ima1, Scalar s_else)
+3. .. cpp:function:: auto ifelse(Image ima, Scalar s_if, Image ima2)
+4. .. cpp:function:: auto ifelse(Image ima, Scalar s_if, Scalar s_else)
+
+    1. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima(p) ? ima1(p) : ima2(p)`.
+    2. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima(p) ? ima1(p) : s_else`.
+    3. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima(p) ? s_if : ima2(p)`.
+    4. Makes a view where for each pixel value evals to :cpp:expr:`out(p) = ima(p) ? s_if : s_else`.
+
+
+    :param ima: Input range
+    :param ima1: Input range
+    :param ima2: Input range
+    :param s: Scalar value
+
+
+    .. code::
+    
+        mln::image2d<int> ima = ...;
+        mln::image2d<int> ima1 = ...;
+        mln::image2d<int> ima2 = ...;
+        auto g1 = view::ifelse(ima == 0, ima1, ima2);
+        auto g2 = view::ifelse(ima < 125, 125, ima);
+        auto g3 = view::ifelse(ima > 125, ima, 125);
+        auto g4 = view::ifelse(ima < 125, 0, 255);
+
+
++----------------+---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+|                |               | *ima, ima1, ima2* | *+, -, \*, /, %, <, <=, ==, !=, >, >=* | *ifelse(ima, ima1, ima2)* | *ifelse(ima, ima1, s_else) / ifelse(ima, s_if, ima2) / ifelse(ima, s_if, s_else)*   |
++================+===============+===================+========================================+===========================+=====================================================================================+
+|                | Forward       | X                 | X                                      | X                         | X                                                                                   |
++                +---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+| **Category**   | Bidirectional | X                 | X                                      | X                         | X                                                                                   |
++                +---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+|                | Raw           | X                 |                                        |                           |                                                                                     |
++----------------+---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+|                | Writable      | X                 |                                        | X                         |                                                                                     |
++                +---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+| **Properties** | Accessible    | X                 | X                                      | X                         | X                                                                                   |
++                +---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+|                | Indexable     | X                 | X                                      | X                         | X                                                                                   |
++----------------+---------------+-------------------+----------------------------------------+---------------------------+-------------------------------------------------------------------------------------+
+
+
+Examples
+--------
+
+TODO
diff --git a/doc/source/core/view/rgb.rst b/doc/source/core/view/rgb.rst
new file mode 100644
index 0000000000000000000000000000000000000000..26c6b025624a5843203c6a80a4c395962d2a58b2
--- /dev/null
+++ b/doc/source/core/view/rgb.rst
@@ -0,0 +1,66 @@
+RGB (View)
+##########
+
+Include :file:`<mln/core/image/view/rgb.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto red(Image ima)
+#. .. cpp:function:: auto green(Image ima)
+#. .. cpp:function:: auto blue(Image ima)
+
+    1. Makes a view that project the red channel of an rgb colored image.
+    2. Makes a view that project the green channel of an rgb colored image.
+    3. Makes a view that project the blue channel of an rgb colored image.
+    
+    :param ima: Input range
+
+    .. code::
+    
+        mln::image2d<mln::rgb8> ima = ...;
+        auto ima_r = mln::view::red(ima);
+        auto ima_g = mln::view::green(ima);
+        auto ima_b = mln::view::blue(ima);
+
+
++----------------+---------------+-------+-----------------------------------+
+|                |               | *ima* | *red(ima), green(ima), blue(ima)* |
++================+===============+=======+===================================+
+|                | Forward       | X     | X                                 |
++                +---------------+-------+-----------------------------------+
+| **Category**   | Bidirectional | X     | X                                 |
++                +---------------+-------+-----------------------------------+
+|                | Raw           | X     |                                   |
++----------------+---------------+-------+-----------------------------------+
+|                | Writable      | X     | X                                 |
++                +---------------+-------+-----------------------------------+
+| **Properties** | Accessible    | X     | X                                 |
++                +---------------+-------+-----------------------------------+
+|                | Indexable     | X     | X                                 |
++----------------+---------------+-------+-----------------------------------+
+
+
+Examples
+--------
+   
+* Set the red channel to max:
+
+.. code::
+
+    mln::image2d<mln::rgb8> ima = {
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} },
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} },
+      { {0, 0, 255}, {0, 0, 255}, {0, 0, 255} }
+      };
+    auto ima_r = mln::view::red(ima);
+    mln::fill(ima_r, 255);
+    mln::io::imprint(ima);
+
+Outputs:
+
+.. code-block:: text
+
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
+    (255, 0, 255) (255, 0, 255) (255, 0, 255)
+
diff --git a/doc/source/core/view/transform.rst b/doc/source/core/view/transform.rst
index 008b0507cec13c862c55745755df8dc388805fcb..173e86a2cc8264560b002327c7de99006047ec8b 100644
--- a/doc/source/core/view/transform.rst
+++ b/doc/source/core/view/transform.rst
@@ -8,18 +8,17 @@ Include :file:`<mln/core/image/view/transform.hpp>`
 #. .. cpp:function:: auto transform(Image ima, UnaryFunction f)
 #. .. cpp:function:: auto transform(Image ima, Image ima2, BinaryFunction f)
 
-    1. Makes a view from `ima` where for each pixel value evals to `out(p) = f(ima(p))`
-    2. Makes a view from `ima` where for each pixel value evals to `out(p) = f(ima(p), ima2(p))`. For this overload, the domain of both image must be equal.
+    1. Makes a view from `ima` where for each pixel value evals to :cpp:expr:`out(p) = f(ima(p))`
+    2. Makes a view from `ima` where for each pixel value evals to :cpp:expr:`out(p) = f(ima(p), ima2(p))`. For this overload, the domain of both image must be equal.
     
-
     :param ima: Input range
     :param ima2: Input range
     :param f: Function to apply on pixel values
 
 
-
 .. warning:: The function f must not create dangling references on the input.
 
+
     .. code::
     
         mln::image2d<mln::rgb8> ima = ...;
@@ -27,23 +26,21 @@ Include :file:`<mln/core/image/view/transform.hpp>`
         auto g2 = mln::view::transform(g1, [](const rgb8& v) -> const uint8& { return v[0]; }); // KO: create a dangling reference to g1(p)
 
 
-
-+----------------+---------------+------------+---------------------+---------------------------+
-|                |               | *ima,ima2* | *transform(ima, f)* | *transform(ima, ima2, f)* |
-+================+===============+============+=====================+===========================+
-|                | Forward       | X          | X                   | X                         |
-+                +---------------+------------+---------------------+---------------------------+
-| **Category**   | Bidirectional | X          | X                   | X                         |
-+                +---------------+------------+---------------------+---------------------------+
-|                | Raw           | X          |                     |                           |
-+----------------+---------------+------------+---------------------+---------------------------+
-|                | Writable      | X          | X*                  | X*                        |
-+                +---------------+------------+---------------------+---------------------------+
-| **Properties** | Accessible    | X          | X                   | X                         |
-+                +---------------+------------+---------------------+---------------------------+
-|                | Indexable     | X          | X                   | X                         |
-+----------------+---------------+------------+---------------------+---------------------------+
-
++----------------+---------------+-------------+---------------------+---------------------------+
+|                |               | *ima, ima2* | *transform(ima, f)* | *transform(ima, ima2, f)* |
++================+===============+=============+=====================+===========================+
+|                | Forward       | X           | X                   | X                         |
++                +---------------+-------------+---------------------+---------------------------+
+| **Category**   | Bidirectional | X           | X                   | X                         |
++                +---------------+-------------+---------------------+---------------------------+
+|                | Raw           | X           |                     |                           |
++----------------+---------------+-------------+---------------------+---------------------------+
+|                | Writable      | X           | X*                  | X*                        |
++                +---------------+-------------+---------------------+---------------------------+
+| **Properties** | Accessible    | X           | X                   | X                         |
++                +---------------+-------------+---------------------+---------------------------+
+|                | Indexable     | X           | X                   | X                         |
++----------------+---------------+-------------+---------------------+---------------------------+
 
 (*): only if `f` is a projection.
 
@@ -55,7 +52,7 @@ Examples
 
 .. code::
 
-    mln::image2d<uint8_t> input = { {0,1,2}, {3,4,5} }; 
+    mln::image2d<uint8_t> input = { {0, 1, 2}, {3, 4, 5} };
     auto g = mln::view::transform(input, [](mln::uint8 v) { return v * 2; });
     mln::io::imprint(g);
 
@@ -70,9 +67,10 @@ Outputs:
 
 .. code::
 
-    mln::image2d<mln::rgb8> input = {{{0,0,0}, {1,1,1}, {2,2,2}}, {{3,3,3}, {4,4,4}, {5,5,5}}}; 
+    mln::image2d<mln::rgb8> input = {{{0, 0, 0}, {1, 1, 1}, {2, 2, 2}}, {{3, 3, 3}, {4, 4, 4}, {5, 5, 5}}};
     auto g = mln::view::transform(input, [](mln::rgb8& x) -> uint8_t& { return x[0]; });
     mln::fill(g, 42);
+    mln::io::imprint(g);
 
 Outputs:
 
@@ -86,7 +84,7 @@ Outputs:
 
 .. code-block:: cpp
 
-    mln::image2d<uint8_t> input1 = ...; 
-    mln::image2d<uint8_t> input2 = ...; 
+    mln::image2d<uint8_t> input1 = ...;
+    mln::image2d<uint8_t> input2 = ...;
     auto f = [](uint8_t x, uint8_t y) { return alpha * x + (1 - alpha) * y; };
     auto g = mln::view::transform(input1, input2, f);
diff --git a/doc/source/core/view/zip.rst b/doc/source/core/view/zip.rst
new file mode 100644
index 0000000000000000000000000000000000000000..393f4516ceb38fed51bf56d122b56a58a4649837
--- /dev/null
+++ b/doc/source/core/view/zip.rst
@@ -0,0 +1,67 @@
+Zip (View)
+##########
+
+Include :file:`<mln/core/image/view/zip.hpp>`
+
+.. cpp:namespace:: mln::view
+
+#. .. cpp:function:: auto zip(Image... ima)
+
+    1. Makes a view from `ima` where for each pixel value evals to :cpp:expr:`out(p) = std::make_tuple(ima(p)...)`
+    
+    :param ima: Input ranges
+
+    .. code::
+    
+        mln::image2d<int> ima1 = ...;
+        mln::image2d<double> ima2 = ...;
+        mln::image2d<mln::rgb8> ima3 = ...;
+        auto g = mln::view::zip(ima1, ima2, ima3);
+
+
++----------------+---------------+----------+---------------+
+|                |               | *ima...* | *zip(ima...)* |
++================+===============+==========+===============+
+|                | Forward       | X        | X             |
++                +---------------+----------+---------------+
+| **Category**   | Bidirectional | X        | X             |
++                +---------------+----------+---------------+
+|                | Raw           | X        |               |
++----------------+---------------+----------+---------------+
+|                | Writable      | X        | X             |
++                +---------------+----------+---------------+
+| **Properties** | Accessible    | X        | X             |
++                +---------------+----------+---------------+
+|                | Indexable     | X        | X             |
++----------------+---------------+----------+---------------+
+
+
+Examples
+--------
+   
+* Traversing both input and output image at the same time:
+
+.. code::
+
+    mln::image2d<uint8_t> input = { {0, 1, 2}, {3, 4, 5} };
+    mln::image2d<uint8_t> output{input.domain()};
+    mln::fill(output, 0);
+    auto zipped_ima = mln::view::zip(input, output);
+    for (auto&& [v_in, v_out] : zipped_ima.values())
+      v_out = v_in * v_in;
+    
+    auto vals = zipped_ima.values();
+    for (auto r : mln::ranges::rows(vals))
+    {
+      for (auto&& [_, v_out] : r)
+        std::cout << v_out << " ";
+      std::cout << std::endl;
+    }
+
+
+Outputs:
+
+.. code-block:: text
+
+    0 1 4 
+    9 16 25 
diff --git a/fixtures/ImageCompare/include/fixtures/ImageCompare/image_compare.hpp b/fixtures/ImageCompare/include/fixtures/ImageCompare/image_compare.hpp
index 9053edfea3119d1a86f13f02396ae2f1bff21163..398cc26bb0a8b260e117be3e12e10429d5dac032 100644
--- a/fixtures/ImageCompare/include/fixtures/ImageCompare/image_compare.hpp
+++ b/fixtures/ImageCompare/include/fixtures/ImageCompare/image_compare.hpp
@@ -1,8 +1,14 @@
 #pragma once
 
 #include <mln/core/image/image.hpp>
+#include <mln/core/rangev3/rows.hpp>
+#include <mln/core/rangev3/view/zip.hpp>
 #include <mln/io/imprint.hpp>
 
+#include <range/v3/begin_end.hpp>
+
+#include <string>
+
 #include <gtest/gtest.h>
 
 
@@ -17,6 +23,12 @@ namespace fixtures::ImageCompare
   template <class I, class J>
   bool compare(const mln::Image<I>& f, const mln::Image<J>& g);
 
+  namespace experimental
+  {
+    template <class ImageLhs, class ImageRhs>
+    bool compare(ImageLhs f, ImageRhs g);
+  }
+
   /******************************************/
   /****          Implementation          ****/
   /******************************************/
@@ -51,8 +63,63 @@ namespace fixtures::ImageCompare
     }
   } // namespace impl
 
+
+  namespace experimental
+  {
+    template <class ImageLhs, class ImageRhs>
+    bool compare(ImageLhs f, ImageRhs g)
+    {
+      static_assert(mln::is_a<ImageLhs, mln::experimental::Image>());
+      static_assert(mln::is_a<ImageRhs, mln::experimental::Image>());
+
+      auto f_dom = f.domain();
+      auto g_dom = f.domain();
+
+      auto f_dom_b = ::ranges::begin(f_dom);
+      auto g_dom_b = ::ranges::begin(g_dom);
+      auto f_dom_e = ::ranges::end(f_dom);
+      auto g_dom_e = ::ranges::end(g_dom);
+      for (; (f_dom_b != f_dom_e) || (g_dom_b != g_dom_e); ++f_dom_b, ++g_dom_b)
+        if (*f_dom_b != *g_dom_b)
+          return false;
+
+      if ((f_dom_b != f_dom_e) || (g_dom_b != g_dom_e))
+        return false;
+
+      auto zipped_vals = mln::ranges::view::zip(f.new_values(), g.new_values());
+      for (auto&& r : mln::ranges::rows(zipped_vals))
+        for (auto&& [f_v, g_v] : r)
+          if (f_v != g_v)
+            return false;
+
+      return true;
+    }
+
+    namespace impl
+    {
+      template <class ImageLhs, class ImageRhs>
+      std::string err_compare_msg(ImageLhs f, ImageRhs g)
+      {
+        static_assert(mln::is_a<ImageLhs, mln::experimental::Image>());
+        static_assert(mln::is_a<ImageRhs, mln::experimental::Image>());
+
+        std::stringstream msg;
+        msg << "The following images differs:\n";
+        mln::io::experimental::imprint(f, msg);
+        msg << " and\n:";
+        mln::io::experimental::imprint(g, msg);
+        return msg.str();
+      }
+    } // namespace impl
+  }   // namespace experimental
+
+
 } // namespace fixtures::ImageCompare
 
 /// \brief GTest macro for image equality test
 #define ASSERT_IMAGES_EQ(f, g)                                                                                         \
   ASSERT_TRUE(::fixtures::ImageCompare::compare(f, g)) << ::fixtures::ImageCompare::impl::err_compare_msg(f, g)
+
+#define ASSERT_IMAGES_EQ_EXP(f, g)                                                                                     \
+  ASSERT_TRUE(::fixtures::ImageCompare::experimental::compare(f, g))                                                   \
+      << ::fixtures::ImageCompare::experimental::impl::err_compare_msg(f, g)
diff --git a/pylene/include/mln/accu/accumulators/moment_of_inertia.hpp b/pylene/include/mln/accu/accumulators/moment_of_inertia.hpp
index 6543f54ad8aa11733ec7e5d0e13b98d747b9bd13..f193783383368a8002e0e4dae226dafb853d16db 100644
--- a/pylene/include/mln/accu/accumulators/moment_of_inertia.hpp
+++ b/pylene/include/mln/accu/accumulators/moment_of_inertia.hpp
@@ -4,7 +4,7 @@
 #include <mln/accu/accumulators/mean.hpp>
 #include <mln/accu/accumulators/sum.hpp>
 #include <mln/accu/composite_accumulator.hpp>
-#include <mln/core/math_ops.hpp>
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/value/value_traits.hpp>
 
 #include <boost/type_traits/promote.hpp>
diff --git a/pylene/include/mln/accu/accumulators/variance.hpp b/pylene/include/mln/accu/accumulators/variance.hpp
index c6b0e29c4217acd6a0b3dbd55b8af56a30382ccf..554e3e576d5be76bed84420d906e5fe315178663 100644
--- a/pylene/include/mln/accu/accumulators/variance.hpp
+++ b/pylene/include/mln/accu/accumulators/variance.hpp
@@ -4,7 +4,7 @@
 #include <mln/accu/accumulators/mean.hpp>
 #include <mln/accu/accumulators/sum.hpp>
 #include <mln/accu/composite_accumulator.hpp>
-#include <mln/core/math_ops.hpp>
+#include <mln/core/functional_ops.hpp>
 
 #include <boost/type_traits/promote.hpp>
 
diff --git a/pylene/include/mln/core/colors.hpp b/pylene/include/mln/core/colors.hpp
index 8dcf0d0ac0f09a610f6715d5dae420cb5f43365e..367912578c87367bc99edd4772c76fb3d61f3270 100644
--- a/pylene/include/mln/core/colors.hpp
+++ b/pylene/include/mln/core/colors.hpp
@@ -44,14 +44,18 @@ namespace mln
   template <typename T>
   using rgb = internal::vec_base<T, 3, rgb_tag>;
 
-  typedef rgb<uint8>  rgb8;
-  typedef rgb<uint16> rgb16;
+  using rgb8         = rgb<uint8>;
+  using signed_rgb8  = rgb<int8>;
+  using rgb16        = rgb<uint16>;
+  using signed_rgb16 = rgb<int16>;
 
   template <typename T>
   using bgr = internal::vec_base<T, 3, bgr_tag>;
 
-  typedef bgr<uint8>  bgr8;
-  typedef bgr<uint16> bgr16;
+  using bgr8         = bgr<uint8>;
+  using signed_bgr8  = bgr<int8>;
+  using bgr16        = bgr<uint16>;
+  using signed_bgr16 = bgr<int16>;
 
   template <class I>
   internal::transformed_image<I, dyn_getter, false> channel(Image<I>&& f, int k)
diff --git a/pylene/include/mln/core/functional_ops.hpp b/pylene/include/mln/core/functional_ops.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0871530e97ac42cba3d2e911483c5b68567343fd
--- /dev/null
+++ b/pylene/include/mln/core/functional_ops.hpp
@@ -0,0 +1,127 @@
+#pragma once
+
+#include <mln/core/private/maths_ops.hpp>
+#include <mln/core/vec/vec_maths_ops.hpp>
+
+#include <cmath>
+
+
+#include <mln/core/private/def_functional_ops_gen_code_macros.hpp>
+
+namespace mln
+{
+  namespace functional
+  {
+    template <typename T = void>
+    struct sqrt_t;
+    template <typename T = void>
+    struct cbrt_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct pow_t;
+    template <typename T = void>
+    struct sqr_t;
+    template <typename T = void>
+    struct abs_t;
+    template <typename T = void>
+    struct sum_t;
+    template <typename T = void>
+    struct prod_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct dot_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct cross_t;
+    template <typename T = void>
+    struct minimum_t;
+    template <typename T = void>
+    struct maximum_t;
+
+    template <typename T = void>
+    struct l0norm_t;
+    template <typename T = void>
+    struct l1norm_t;
+    template <typename T = void>
+    struct l2norm_t;
+    template <typename T = void>
+    struct l2norm_sqr_t;
+    template <typename T = void>
+    struct linfnorm_t;
+    template <unsigned p, typename T = void>
+    struct lpnorm_t;
+
+    template <typename T1 = void, typename T2 = T1>
+    struct l0dist_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct l1dist_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct l2dist_t;
+    template <unsigned p, typename T1 = void, typename T2 = T1>
+    struct lpdist_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct l2dist_sqr_t;
+    template <typename T1 = void, typename T2 = T1>
+    struct linfdist_t;
+  }; // namespace functional
+
+
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(sqrt);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(cbrt);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(pow);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(sqr);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(abs);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(sum);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(prod);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(dot);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(cross);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(maximum);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(minimum);
+
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(l0norm);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(l1norm);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(l2norm);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(l2norm_sqr);
+  MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(linfnorm);
+  namespace functional
+  {
+    template <unsigned p, typename T>
+    struct lpnorm_t
+    {
+      typedef decltype(lpnorm<p>(std::declval<T>())) result_type;
+      auto operator()(const T& x) const -> decltype(lpnorm<p>(x)) { return lpnorm<p>(x); }
+    };
+    template <unsigned p>
+    struct lpnorm_t<p, void>
+    {
+      template <typename T>
+      auto operator()(const T& x) const -> decltype(lpnorm<p>(x))
+      {
+        return lpnorm<p>(x);
+      }
+    };
+  } // namespace functional
+
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(l0dist);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(l1dist);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(l2dist);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(l2dist_sqr);
+  MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(linfdist);
+  namespace functional
+  {
+    template <unsigned p, typename T1, typename T2>
+    struct lpdist_t
+    {
+      typedef decltype(lpdist<p>(std::declval<T1>(), std::declval<T2>())) result_type;
+      auto operator()(const T1& x, const T2& y) const -> decltype(lpdist<p>(x, y)) { return lpdist<p>(x, y); }
+    };
+    template <unsigned p>
+    struct lpdist_t<p, void, void>
+    {
+      template <typename T1, typename T2>
+      auto operator()(const T1& x, const T2& y) const -> decltype(lpdist<p>(x, y))
+      {
+        return lpdist<p>(x, y);
+      }
+    };
+  } // namespace functional
+} // namespace mln
+
+#include <mln/core/private/undef_functional_ops_gen_code_macros.hpp>
diff --git a/pylene/include/mln/core/image/image.hpp b/pylene/include/mln/core/image/image.hpp
index 6f732690f01c42ed6dfd7a6dc3ba69b5de0c94f3..ba5532b6a13ad67cbbf5eae9a4313fb7a6ef068e 100644
--- a/pylene/include/mln/core/image/image.hpp
+++ b/pylene/include/mln/core/image/image.hpp
@@ -22,7 +22,7 @@
 #include <mln/core/image/private/image_traits.hpp>
 
 #include <mln/core/image/image_expr.hpp>
-#include <mln/core/image/image_math_ops.hpp>
+#include <mln/core/image/image_maths_ops.hpp>
 #include <mln/core/image/image_ops.hpp>
 
 #include <mln/core/image/image_routines.hpp>
diff --git a/pylene/include/mln/core/image/image_math_ops.hpp b/pylene/include/mln/core/image/image_maths_ops.hpp
similarity index 96%
rename from pylene/include/mln/core/image/image_math_ops.hpp
rename to pylene/include/mln/core/image/image_maths_ops.hpp
index 27f2015ab28f20d985a64a1d1a17c5f3f6f1e63d..12a9f643c7440c17f2685c49603964b65ebfaead 100644
--- a/pylene/include/mln/core/image/image_math_ops.hpp
+++ b/pylene/include/mln/core/image/image_maths_ops.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/image/image_expr.hpp>
 #include <mln/core/image/image_ops.hpp>
-#include <mln/core/math_ops.hpp>
 
 namespace mln
 {
diff --git a/pylene/include/mln/core/image/private/def_gen_code_macros.hpp b/pylene/include/mln/core/image/private/def_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..80f9ca1eb74b943693d70e2c84ee1b20cca2bfe8
--- /dev/null
+++ b/pylene/include/mln/core/image/private/def_gen_code_macros.hpp
@@ -0,0 +1,48 @@
+#include <mln/core/image/image.hpp>
+#include <mln/core/image/view/transform.hpp>
+#include <mln/core/image/view/zip.hpp>
+
+#include <type_traits>
+
+
+#define MLN_PRIVATE_DEFINE_UNARY_OPERATOR(op, f)                                                                       \
+  template <class I, class = std::enable_if_t<::mln::is_a<I, ::mln::experimental::Image>::value>>                      \
+  auto op(const I& ima)                                                                                                \
+  {                                                                                                                    \
+    return ::mln::view::transform(static_cast<const I&>(ima), f);                                                      \
+  }
+
+#define MLN_PRIVATE_DEFINE_BINARY_OPERATOR(op, f)                                                                      \
+  namespace impl                                                                                                       \
+  {                                                                                                                    \
+    template <class I1, class I2>                                                                                      \
+    auto op(const ::mln::experimental::Image<I1>& ima1, const ::mln::experimental::Image<I2>& ima2)                    \
+    {                                                                                                                  \
+      return ::mln::view::transform(static_cast<const I1&>(ima1), static_cast<const I2&>(ima2), f);                    \
+    }                                                                                                                  \
+                                                                                                                       \
+    template <class I, class Scalar,                                                                                   \
+              class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>                       \
+    auto op(const ::mln::experimental::Image<I>& ima1, Scalar s)                                                       \
+    {                                                                                                                  \
+      auto g = [f_ = f, s_ = s](auto&& arg) { return f_(arg, s_); };                                                   \
+      return ::mln::view::transform(static_cast<const I&>(ima1), g);                                                   \
+    }                                                                                                                  \
+                                                                                                                       \
+    template <class Scalar, class I,                                                                                   \
+              class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>                       \
+    auto op(Scalar s, const ::mln::experimental::Image<I>& ima2)                                                       \
+    {                                                                                                                  \
+      auto g = [f_ = f, s_ = s](auto&& arg) { return f_(s_, arg); };                                                   \
+      return ::mln::view::transform(static_cast<const I&>(ima2), g);                                                   \
+    }                                                                                                                  \
+  } /* namespace impl */                                                                                               \
+                                                                                                                       \
+  /* This overload is there to be a best match wrt old API impl */                                                     \
+  template <class A, class B,                                                                                          \
+            class = std::enable_if_t<(::mln::is_a<A, ::mln::experimental::Image>::value ||                             \
+                                      ::mln::is_a<B, ::mln::experimental::Image>::value)>>                             \
+  auto op(const A& lhs, const B& rhs)                                                                                  \
+  {                                                                                                                    \
+    return impl::op(lhs, rhs);                                                                                         \
+  }
diff --git a/pylene/include/mln/core/image/private/undef_gen_code_macros.hpp b/pylene/include/mln/core/image/private/undef_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3156bf117ff7d648c7e8312aae66380806184e03
--- /dev/null
+++ b/pylene/include/mln/core/image/private/undef_gen_code_macros.hpp
@@ -0,0 +1,3 @@
+
+#undef MLN_PRIVATE_DEFINE_UNARY_OPERATOR
+#undef MLN_PRIVATE_DEFINE_BINARY_OPERATOR
diff --git a/pylene/include/mln/core/image/view/cast.hpp b/pylene/include/mln/core/image/view/cast.hpp
index 9dd883b8b5e8d3f06dec6d1f73a3435edbc99f00..546d3ee113f345beb9cd3feb966e1f3ff7b27121 100644
--- a/pylene/include/mln/core/image/view/cast.hpp
+++ b/pylene/include/mln/core/image/view/cast.hpp
@@ -24,7 +24,7 @@ namespace mln::view
   template <class V, class InputImage>
   auto cast(InputImage ima)
   {
-    static_assert(is_a<InputImage, experimental::Image>());
+    static_assert(is_a<InputImage, mln::experimental::Image>());
 
     return view::transform(std::move(ima), detail::cast_to<V>{});
   }
diff --git a/pylene/include/mln/core/image/view/channel.hpp b/pylene/include/mln/core/image/view/channel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..47d418b1fe64df2445e2633747bf040eac757916
--- /dev/null
+++ b/pylene/include/mln/core/image/view/channel.hpp
@@ -0,0 +1,26 @@
+#pragma once
+
+#include <mln/core/image/image.hpp>
+#include <mln/core/image/view/transform.hpp>
+#include <mln/core/ops.hpp>
+#include <mln/core/vec/vec_ops.hpp>
+
+
+namespace mln::view
+{
+  template <typename InputImage>
+  auto channel(InputImage ima, int k)
+  {
+    static_assert(is_a<InputImage, mln::experimental::Image>());
+
+    return mln::view::transform(std::move(ima), dyn_getter{k});
+  }
+
+  template <unsigned k, typename InputImage>
+  auto channel(InputImage ima)
+  {
+    static_assert(is_a<InputImage, mln::experimental::Image>());
+
+    return mln::view::transform(std::move(ima), getter<k>{});
+  }
+} // namespace mln::view
diff --git a/pylene/include/mln/core/image/view/maths.hpp b/pylene/include/mln/core/image/view/maths.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..84dfb9bc04395973c0e2a271248e48b4b6425238
--- /dev/null
+++ b/pylene/include/mln/core/image/view/maths.hpp
@@ -0,0 +1,78 @@
+#pragma once
+
+#include <mln/core/functional_ops.hpp>
+
+#include <mln/core/image/private/def_gen_code_macros.hpp>
+
+namespace mln::view
+{
+
+  namespace maths
+  {
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(abs, functional::abs_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(sqr, functional::sqr_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(pow, functional::pow_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(sqrt, functional::sqrt_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(cbrt, functional::cbrt_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(sum, functional::sum_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(prod, functional::prod_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(dot, functional::dot_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(cross, functional::cross_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(min, functional::minimum_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(max, functional::maximum_t<>());
+
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(l0norm, functional::l0norm_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(l1norm, functional::l1norm_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(l2norm, functional::l2norm_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(l2norm_sqr, functional::l2norm_sqr_t<>());
+    MLN_PRIVATE_DEFINE_UNARY_OPERATOR(linfnorm, functional::linfnorm_t<>());
+    template <unsigned p, class I, class = std::enable_if_t<::mln::is_a<I, ::mln::experimental::Image>::value>>
+    auto lpnorm(const I& ima)
+    {
+      return ::mln::view::transform(static_cast<const I&>(ima), functional::lpnorm_t<p>{});
+    }
+
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(l0dist, functional::l0dist_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(l1dist, functional::l1dist_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(l2dist, functional::l2dist_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(l2dist_sqr, functional::l2dist_sqr_t<>());
+    MLN_PRIVATE_DEFINE_BINARY_OPERATOR(linfdist, functional::linfdist_t<>());
+    namespace impl
+    {
+      template <unsigned p, class I1, class I2>
+      auto lpdist(const ::mln::experimental::Image<I1>& ima1, const ::mln::experimental::Image<I2>& ima2)
+      {
+        return ::mln::view::transform(static_cast<const I1&>(ima1), static_cast<const I2&>(ima2),
+                                      functional::lpdist_t<p>{});
+      }
+
+      template <unsigned p, class I, class Scalar,
+                class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>
+      auto lpdist(const ::mln::experimental::Image<I>& ima1, Scalar s)
+      {
+        auto g = [f_ = functional::lpdist_t<p>{}, s_ = s](auto&& arg) { return f_(arg, s_); };
+        return ::mln::view::transform(static_cast<const I&>(ima1), g);
+      }
+
+      template <unsigned p, class Scalar, class I,
+                class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>
+      auto lpdist(Scalar s, const ::mln::experimental::Image<I>& ima2)
+      {
+        auto g = [f_ = functional::lpdist_t<p>{}, s_ = s](auto&& arg) { return f_(s_, arg); };
+        return ::mln::view::transform(static_cast<const I&>(ima2), g);
+      }
+    } /* namespace impl */
+
+    /* This overload is there to be a best match wrt old API impl */
+    template <unsigned p, class A, class B,
+              class = std::enable_if_t<(::mln::is_a<A, ::mln::experimental::Image>::value ||
+                                        ::mln::is_a<B, ::mln::experimental::Image>::value)>>
+    auto lpdist(const A& lhs, const B& rhs)
+    {
+      return impl::lpdist<p>(lhs, rhs);
+    }
+
+  } // namespace maths
+} // namespace mln::view
+
+#include <mln/core/image/private/undef_gen_code_macros.hpp>
diff --git a/pylene/include/mln/core/image/view/operators.hpp b/pylene/include/mln/core/image/view/operators.hpp
index 59618a1eb592a0f67d2434a3efa6abd2cb7831ed..a4d6519ae046fb8d138f4feebee12fd56b586772 100644
--- a/pylene/include/mln/core/image/view/operators.hpp
+++ b/pylene/include/mln/core/image/view/operators.hpp
@@ -1,62 +1,24 @@
 #pragma once
 
+#include <mln/core/image/image.hpp>
 #include <mln/core/image/view/transform.hpp>
 #include <mln/core/image/view/zip.hpp>
 #include <mln/core/rangev3/foreach.hpp>
 
 #include <functional>
+#include <type_traits>
+#include <utility>
 
+#include <mln/core/image/private/def_gen_code_macros.hpp>
 
 namespace mln::view
 {
   // Prevent ADL on those operators so that A+B uses old API
-  // and using mln::experimental::ops A+B uses new API
+  // and using mln::view::ops A+B uses new API
   // to be reverted when the new API is validated
   namespace ops
   {
 
-#define MLN_PRIVATE_DEFINE_UNARY_OPERATOR(op, f)                                                                       \
-  template <class I, class = std::enable_if_t<mln::is_a<I, mln::experimental::Image>::value>>                          \
-  auto op(const I& ima)                                                                                                \
-  {                                                                                                                    \
-    return view::transform(static_cast<const I&>(ima), f);                                                             \
-  }
-
-
-#define MLN_PRIVATE_DEFINE_BINARY_OPERATOR(op, f)                                                                      \
-  namespace impl                                                                                                       \
-  {                                                                                                                    \
-    template <class I1, class I2>                                                                                      \
-    auto op(const mln::experimental::Image<I1>& ima1, const mln::experimental::Image<I2>& ima2)                        \
-    {                                                                                                                  \
-      return view::transform(static_cast<const I1&>(ima1), static_cast<const I2&>(ima2), f);                           \
-    }                                                                                                                  \
-                                                                                                                       \
-    template <class I, class Scalar, class = std::enable_if_t<!is_a<Scalar, mln::experimental::Image>::value>>         \
-    auto op(const mln::experimental::Image<I>& ima1, Scalar s)                                                         \
-    {                                                                                                                  \
-      auto g = [f_ = f, s_ = s](auto&& arg) { return f_(arg, s_); };                                                   \
-      return view::transform(static_cast<const I&>(ima1), g);                                                          \
-    }                                                                                                                  \
-                                                                                                                       \
-    template <class Scalar, class I, class = std::enable_if_t<!is_a<Scalar, mln::experimental::Image>::value>>         \
-    auto op(Scalar s, const mln::experimental::Image<I>& ima2)                                                         \
-    {                                                                                                                  \
-      auto g = [f_ = f, s_ = s](auto&& arg) { return f_(s_, arg); };                                                   \
-      return view::transform(static_cast<const I&>(ima2), g);                                                          \
-    }                                                                                                                  \
-  }                                                                                                                    \
-                                                                                                                       \
-  /* This overload is there to be a best match wrt old API impl */                                                     \
-  template <class A, class B,                                                                                          \
-            class = std::enable_if_t<(is_a<A, mln::experimental::Image>::value ||                                      \
-                                      is_a<B, mln::experimental::Image>::value)>>                                      \
-  auto op(const A& lhs, const B& rhs)                                                                                  \
-  {                                                                                                                    \
-    return mln::view::ops::impl::op(lhs, rhs);                                                                         \
-  }
-
-
     namespace details
     {
       template <class Base>
@@ -88,6 +50,64 @@ namespace mln::view
     MLN_PRIVATE_DEFINE_BINARY_OPERATOR(operator*, details::to_common_type<std::multiplies<>>());
     MLN_PRIVATE_DEFINE_BINARY_OPERATOR(operator/, details::to_common_type<std::divides<>>());
     MLN_PRIVATE_DEFINE_BINARY_OPERATOR(operator%, details::to_common_type<std::modulus<>>());
+
+    // FIXME: integrate a transient mechanism to do this comparison : implement operators on vec_base
+    namespace impl
+    {
+      struct equalFP_t
+      {
+        equalFP_t(double eps)
+          : eps_(eps)
+        {
+        }
+
+        bool operator()(double l, double r) const { return std::abs(l - r) < eps_; }
+
+        template <unsigned dim, class tag>
+        bool operator()(const mln::internal::vec_base<double, dim, tag>& l,
+                        const mln::internal::vec_base<double, dim, tag>& r) const
+        {
+          for (unsigned i = 0; i < dim; ++i)
+            if (std::abs(l[i] - r[i]) >= eps_)
+              return false;
+          return true;
+        }
+
+      private:
+        double eps_;
+      };
+
+      template <class I1, class I2>
+      auto equalFP(const ::mln::experimental::Image<I1>& ima1, const ::mln::experimental::Image<I2>& ima2, double eps)
+      {
+        return ::mln::view::transform(static_cast<const I1&>(ima1), static_cast<const I2&>(ima2), equalFP_t{eps});
+      }
+
+      template <class I, class Scalar,
+                class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>
+      auto equalFP(const ::mln::experimental::Image<I>& ima1, Scalar s, double eps)
+      {
+        auto g = [f_ = equalFP_t{eps}, s_ = s](auto&& arg) { return f_(arg, s_); };
+        return ::mln::view::transform(static_cast<const I&>(ima1), g);
+      }
+
+      template <class Scalar, class I,
+                class = std::enable_if_t<!::mln::is_a<Scalar, ::mln::experimental::Image>::value>>
+      auto equalFP(Scalar s, const ::mln::experimental::Image<I>& ima2, double eps)
+      {
+        auto g = [f_ = equalFP_t{eps}, s_ = s](auto&& arg) { return f_(s_, arg); };
+        return ::mln::view::transform(static_cast<const I&>(ima2), g);
+      }
+    } /* namespace impl */
+
+    /* This overload is there to be a best match wrt old API impl */
+    template <class A, class B,
+              class = std::enable_if_t<(::mln::is_a<A, ::mln::experimental::Image>::value ||
+                                        ::mln::is_a<B, ::mln::experimental::Image>::value)>>
+    auto equalFP(const A& lhs, const B& rhs, double eps)
+    {
+      return impl::equalFP(lhs, rhs, eps);
+    }
   } // namespace ops
 
 
@@ -98,9 +118,9 @@ namespace mln::view
 
 
     template <class ICond, class ITrue, class IFalse>
-    struct ifelse_fn<
-        ICond, ITrue, IFalse,
-        std::enable_if_t<is_a<ITrue, mln::experimental::Image>::value && is_a<IFalse, mln::experimental::Image>::value>>
+    struct ifelse_fn<ICond, ITrue, IFalse,
+                     std::enable_if_t<(is_a<ITrue, mln::experimental::Image>::value &&
+                                       is_a<IFalse, mln::experimental::Image>::value)>>
     {
       auto operator()(const ICond& cond, ITrue iftrue, IFalse iffalse) const
       {
@@ -115,8 +135,8 @@ namespace mln::view
 
     template <class ICond, class ITrue, class IFalse>
     struct ifelse_fn<ICond, ITrue, IFalse,
-                     std::enable_if_t<!is_a<ITrue, mln::experimental::Image>::value &&
-                                      is_a<IFalse, mln::experimental::Image>::value>>
+                     std::enable_if_t<(!is_a<ITrue, mln::experimental::Image>::value &&
+                                       is_a<IFalse, mln::experimental::Image>::value)>>
     {
       auto operator()(const ICond& cond, ITrue vtrue, IFalse iffalse) const
       {
@@ -130,8 +150,8 @@ namespace mln::view
 
     template <class ICond, class ITrue, class IFalse>
     struct ifelse_fn<ICond, ITrue, IFalse,
-                     std::enable_if_t<is_a<ITrue, mln::experimental::Image>::value &&
-                                      !is_a<IFalse, mln::experimental::Image>::value>>
+                     std::enable_if_t<(is_a<ITrue, mln::experimental::Image>::value &&
+                                       !is_a<IFalse, mln::experimental::Image>::value)>>
     {
       auto operator()(const ICond& cond, ITrue iftrue, IFalse vfalse) const
       {
@@ -146,8 +166,8 @@ namespace mln::view
 
     template <class ICond, class ITrue, class IFalse>
     struct ifelse_fn<ICond, ITrue, IFalse,
-                     std::enable_if_t<!is_a<ITrue, mln::experimental::Image>::value &&
-                                      !is_a<IFalse, mln::experimental::Image>::value>>
+                     std::enable_if_t<(!is_a<ITrue, mln::experimental::Image>::value &&
+                                       !is_a<IFalse, mln::experimental::Image>::value)>>
     {
       auto operator()(const ICond& cond, ITrue vtrue, IFalse vfalse) const
       {
@@ -164,7 +184,6 @@ namespace mln::view
                                                       std::move(iffalse));
   }
 
-#undef MLN_PRIVATE_DEFINE_UNARY_OPERATOR
-#undef MLN_PRIVATE_DEFINE_BINARY_OPERATOR
+#include <mln/core/image/private/undef_gen_code_macros.hpp>
 
 } // namespace mln::view
diff --git a/pylene/include/mln/core/image/view/rgb.hpp b/pylene/include/mln/core/image/view/rgb.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b432d9168cf848438b0610013a229ad41c02ea5
--- /dev/null
+++ b/pylene/include/mln/core/image/view/rgb.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <mln/core/image/image.hpp>
+#include <mln/core/image/view/channel.hpp>
+
+
+namespace mln::view
+{
+  template <typename InputImage>
+  auto red(InputImage ima)
+  {
+    static_assert(is_a<InputImage, mln::experimental::Image>());
+
+    return mln::view::channel<0>(std::move(ima));
+  }
+
+  template <typename InputImage>
+  auto green(InputImage ima)
+  {
+    static_assert(is_a<InputImage, mln::experimental::Image>());
+
+    return mln::view::channel<1>(std::move(ima));
+  }
+
+  template <typename InputImage>
+  auto blue(InputImage ima)
+  {
+    static_assert(is_a<InputImage, mln::experimental::Image>());
+
+    return mln::view::channel<2>(std::move(ima));
+  }
+} // namespace mln::view
diff --git a/pylene/include/mln/core/image/views.hpp b/pylene/include/mln/core/image/views.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..44918f794f7d3856ee77b79f5015bac56d056cd2
--- /dev/null
+++ b/pylene/include/mln/core/image/views.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+#include <mln/core/image/view/adaptor.hpp>
+#include <mln/core/image/view/cast.hpp>
+#include <mln/core/image/view/channel.hpp>
+#include <mln/core/image/view/clip.hpp>
+#include <mln/core/image/view/filter.hpp>
+#include <mln/core/image/view/mask.hpp>
+#include <mln/core/image/view/maths.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/rgb.hpp>
+#include <mln/core/image/view/transform.hpp>
+#include <mln/core/image/view/zip.hpp>
diff --git a/pylene/include/mln/core/math_ops.hpp b/pylene/include/mln/core/math_ops.hpp
deleted file mode 100644
index e3e096cd4ddd55db14885e3140500d1b89962a5c..0000000000000000000000000000000000000000
--- a/pylene/include/mln/core/math_ops.hpp
+++ /dev/null
@@ -1,475 +0,0 @@
-#pragma once
-
-/**
- * \file
- * \brief Define fundamental mathematical and statistical operators
- *
- */
-
-#include <cmath>
-#include <cstdlib>
-#include <utility>
-
-
-namespace mln
-{
-  /* Element wise operators */
-  int                sqr(int x);
-  long               sqr(long n);
-  long long          sqr(long long n);
-  unsigned int       sqr(unsigned int x);
-  unsigned long      sqr(unsigned long n);
-  unsigned long long sqr(unsigned long long n);
-  float              sqr(float x);
-  double             sqr(double x);
-  long double        sqr(long double x);
-
-  int                abs(int x);
-  long               abs(long n);
-  long long          abs(long long n);
-  unsigned int       abs(unsigned int x);
-  unsigned long      abs(unsigned long n);
-  unsigned long long abs(unsigned long long n);
-  float              abs(float x);
-  double             abs(double x);
-  long double        abs(long double x);
-
-  using std::cbrt;
-  using std::pow;
-  using std::sqrt;
-
-  /* Reduction operators */
-  int                sum(int x);
-  long               sum(long n);
-  long long          sum(long long n);
-  unsigned int       sum(unsigned int x);
-  unsigned long      sum(unsigned long n);
-  unsigned long long sum(unsigned long long n);
-  float              sum(float x);
-  double             sum(double x);
-  long double        sum(long double x);
-
-  int                prod(int x);
-  long               prod(long n);
-  long long          prod(long long n);
-  unsigned int       prod(unsigned int x);
-  unsigned long      prod(unsigned long n);
-  unsigned long long prod(unsigned long long n);
-  float              prod(float x);
-  double             prod(double x);
-  long double        prod(long double x);
-
-  int                minimum(int x);
-  long               minimum(long n);
-  long long          minimum(long long n);
-  unsigned int       minimum(unsigned int x);
-  unsigned long      minimum(unsigned long n);
-  unsigned long long minimum(unsigned long long n);
-  float              minimum(float x);
-  double             minimum(double x);
-  long double        minimum(long double x);
-
-  int                maximum(int x);
-  long               maximum(long n);
-  long long          maximum(long long n);
-  unsigned int       maximum(unsigned int x);
-  unsigned long      maximum(unsigned long n);
-  unsigned long long maximum(unsigned long long n);
-  float              maximum(float x);
-  double             maximum(double x);
-  long double        maximum(long double x);
-
-  int                l0norm(int x);
-  long               l0norm(long n);
-  long long          l0norm(long long n);
-  unsigned int       l0norm(unsigned int x);
-  unsigned long      l0norm(unsigned long n);
-  unsigned long long l0norm(unsigned long long n);
-  float              l0norm(float x);
-  double             l0norm(double x);
-  long double        l0norm(long double x);
-
-  int                l1norm(int x);
-  long               l1norm(long n);
-  long long          l1norm(long long n);
-  unsigned int       l1norm(unsigned int x);
-  unsigned long      l1norm(unsigned long n);
-  unsigned long long l1norm(unsigned long long n);
-  float              l1norm(float x);
-  double             l1norm(double x);
-  long double        l1norm(long double x);
-
-  int                l2norm(int x);
-  long               l2norm(long n);
-  long long          l2norm(long long n);
-  unsigned int       l2norm(unsigned int x);
-  unsigned long      l2norm(unsigned long n);
-  unsigned long long l2norm(unsigned long long n);
-  float              l2norm(float x);
-  double             l2norm(double x);
-  long double        l2norm(long double x);
-
-  int                l2dist(int x, int y);
-  long               l2dist(long x, long y);
-  long long          l2dist(long long x, long long y);
-  unsigned int       l2dist(unsigned int x, unsigned int y);
-  unsigned long      l2dist(unsigned long x, unsigned long y);
-  unsigned long long l2dist(unsigned long long x, unsigned long long y);
-  float              l2dist(float x, float y);
-  double             l2dist(double x, double y);
-  long double        l2dist(long double x, long double y);
-
-  int                l1dist(int x, int y);
-  long               l1dist(long x, long y);
-  long long          l1dist(long long x, long long y);
-  unsigned int       l1dist(unsigned int x, unsigned int y);
-  unsigned long      l1dist(unsigned long x, unsigned long y);
-  unsigned long long l1dist(unsigned long long x, unsigned long long y);
-  float              l1dist(float x, float y);
-  double             l1dist(double x, double y);
-  long double        l1dist(long double x, long double y);
-
-  int                l2norm_sqr(int x);
-  long               l2norm_sqr(long n);
-  long long          l2norm_sqr(long long n);
-  unsigned int       l2norm_sqr(unsigned int x);
-  unsigned long      l2norm_sqr(unsigned long n);
-  unsigned long long l2norm_sqr(unsigned long long n);
-  float              l2norm_sqr(float x);
-  double             l2norm_sqr(double x);
-  long double        l2norm_sqr(long double x);
-
-  int                l2dist_sqr(int x, int y);
-  long               l2dist_sqr(long x, long y);
-  long long          l2dist_sqr(long long x, long long y);
-  unsigned int       l2dist_sqr(unsigned int x, unsigned int y);
-  unsigned long      l2dist_sqr(unsigned long x, unsigned long y);
-  unsigned long long l2dist_sqr(unsigned long long x, unsigned long long y);
-  float              l2dist_sqr(float x, float y);
-  double             l2dist_sqr(double x, double y);
-  long double        l2dist_sqr(long double x, long double y);
-
-  int                linfnorm(int x);
-  long               linfnorm(long n);
-  long long          linfnorm(long long n);
-  unsigned int       linfnorm(unsigned int x);
-  unsigned long      linfnorm(unsigned long n);
-  unsigned long long linfnorm(unsigned long long n);
-  float              linfnorm(float x);
-  double             linfnorm(double x);
-  long double        linfnorm(long double x);
-
-  template <unsigned p>
-  int lpnorm(int x);
-  template <unsigned p>
-  long lpnorm(long x);
-  template <unsigned p>
-  long long lpnorm(long long x);
-  template <unsigned p>
-  unsigned int lpnorm(unsigned int x);
-  template <unsigned p>
-  unsigned long lpnorm(unsigned long n);
-  template <unsigned p>
-  unsigned long long lpnorm(unsigned long long n);
-  template <unsigned p>
-  float lpnorm(float x);
-  template <unsigned p>
-  double lpnorm(double x);
-  template <unsigned p>
-  long double lpnorm(long double x);
-
-  namespace functional
-  {
-    template <typename T = void>
-    struct sqrt_t;
-    template <typename T = void>
-    struct cbrt_t;
-    template <typename T1 = void, typename T2 = T1>
-    struct pow_t;
-    template <typename T = void>
-    struct sqr_t;
-    template <typename T = void>
-    struct abs_t;
-
-    template <typename T = void>
-    struct sum_t;
-    template <typename T = void>
-    struct prod_t;
-    template <typename T = void>
-    struct minimum_t;
-    template <typename T = void>
-    struct maximum_t;
-    template <typename T = void>
-    struct l0norm_t;
-    template <typename T = void>
-    struct l1norm_t;
-    template <typename T = void>
-    struct l2norm_t;
-    template <typename T = void>
-    struct l2norm_sqr_t;
-    template <typename T = void>
-    struct linfnorm_t;
-    template <unsigned p, typename T = void>
-    struct lpnorm_t;
-
-    template <typename T1 = void, typename T2 = T1>
-    struct l0dist_t;
-    template <typename T1 = void, typename T2 = T1>
-    struct l1dist_t;
-    template <typename T1 = void, typename T2 = T1>
-    struct l2dist_t;
-    template <typename T1 = void, typename T2 = T1>
-    struct linfdist_t;
-  }; // namespace functional
-
-    /********************************/
-    /*** Implementation            **/
-    /********************************/
-
-#define MLN_GEN_CODE(FUN, TYPE, OP)                                                                                    \
-  inline TYPE FUN(TYPE x) { return OP; }
-
-#define MLN_GEN_BINARY_CODE(FUN, TYPE, OP)                                                                             \
-  inline TYPE FUN(TYPE x, TYPE y) { return OP; }
-
-#define MLN_GEN_CODE_2(TEMPLATE, FUN, TYPE, OP)                                                                        \
-  TEMPLATE inline TYPE FUN(TYPE x) { return OP; }
-
-#define MLN_FUNCTIONAL_GEN_UNARY_CODE(FUN)                                                                             \
-  namespace functional                                                                                                 \
-  {                                                                                                                    \
-    using mln::FUN;                                                                                                    \
-    template <typename T>                                                                                              \
-    struct FUN##_t                                                                                                     \
-    {                                                                                                                  \
-      typedef decltype(FUN(std::declval<T>())) result_type;                                                            \
-      auto                                     operator()(const T& x) const -> decltype(FUN(x)) { return FUN(x); }     \
-    };                                                                                                                 \
-    template <>                                                                                                        \
-    struct FUN##_t<void>                                                                                               \
-    {                                                                                                                  \
-      template <typename T>                                                                                            \
-      auto operator()(const T& x) const -> decltype(FUN(x))                                                            \
-      {                                                                                                                \
-        return FUN(x);                                                                                                 \
-      }                                                                                                                \
-    };                                                                                                                 \
-  }
-
-#define MLN_FUNCTIONAL_GEN_BINARY_CODE(FUN)                                                                            \
-  namespace functional                                                                                                 \
-  {                                                                                                                    \
-    using mln::FUN;                                                                                                    \
-    template <typename T1, typename T2>                                                                                \
-    struct FUN##_t                                                                                                     \
-    {                                                                                                                  \
-      typedef decltype(FUN(std::declval<T1>(), std::declval<T2>())) result_type;                                       \
-      auto operator()(const T1& x, const T2& y) const -> decltype(FUN(x, y)) { return FUN(x, y); }                     \
-    };                                                                                                                 \
-    template <>                                                                                                        \
-    struct FUN##_t<void, void>                                                                                         \
-    {                                                                                                                  \
-      template <typename T1, typename T2>                                                                              \
-      auto operator()(const T1& x, const T2& y) const -> decltype(FUN(x, y))                                           \
-      {                                                                                                                \
-        return FUN(x, y);                                                                                              \
-      }                                                                                                                \
-    };                                                                                                                 \
-  }
-
-#define MLN_FUNCTIONAL_GEN_CODE_2(TTYPE, TNAME, FUN)                                                                   \
-  namespace functional                                                                                                 \
-  {                                                                                                                    \
-    template <TTYPE TNAME, typename T>                                                                                 \
-    struct FUN##_t                                                                                                     \
-    {                                                                                                                  \
-      typedef decltype(FUN<TTYPE>(std::declval<T>())) result_type;                                                     \
-      auto operator()(const T& x) const -> decltype(FUN<TTYPE>(x)) { return FUN<TTYPE>(x); }                           \
-    };                                                                                                                 \
-    template <TTYPE TNAME>                                                                                             \
-    struct FUN##_t<TNAME, void>                                                                                        \
-    {                                                                                                                  \
-      template <typename T>                                                                                            \
-      auto operator()(const T& x) const -> decltype(FUN<TTYPE>(x))                                                     \
-      {                                                                                                                \
-        return FUN<TTYPE>(x);                                                                                          \
-      }                                                                                                                \
-    };                                                                                                                 \
-  }
-
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(sqrt);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(cbrt);
-  MLN_FUNCTIONAL_GEN_BINARY_CODE(pow);
-
-  MLN_GEN_CODE(sqr, int, x* x);
-  MLN_GEN_CODE(sqr, long, x* x);
-  MLN_GEN_CODE(sqr, long long, x* x);
-  MLN_GEN_CODE(sqr, unsigned int, x* x);
-  MLN_GEN_CODE(sqr, unsigned long, x* x);
-  MLN_GEN_CODE(sqr, unsigned long long, x* x);
-  MLN_GEN_CODE(sqr, float, x* x);
-  MLN_GEN_CODE(sqr, double, x* x);
-  MLN_GEN_CODE(sqr, long double, x* x);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(sqr);
-
-  MLN_GEN_CODE(abs, int, std::abs(x));
-  MLN_GEN_CODE(abs, long, std::abs(x));
-  MLN_GEN_CODE(abs, long long, std::abs(x));
-  MLN_GEN_CODE(abs, unsigned int, x);
-  MLN_GEN_CODE(abs, unsigned long, x);
-  MLN_GEN_CODE(abs, unsigned long long, x);
-  MLN_GEN_CODE(abs, float, std::abs(x));
-  MLN_GEN_CODE(abs, double, std::abs(x));
-  MLN_GEN_CODE(abs, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(abs);
-
-  MLN_GEN_CODE(l0norm, int, std::abs(x));
-  MLN_GEN_CODE(l0norm, long, std::abs(x));
-  MLN_GEN_CODE(l0norm, long long, std::abs(x));
-  MLN_GEN_CODE(l0norm, unsigned int, x);
-  MLN_GEN_CODE(l0norm, unsigned long, x);
-  MLN_GEN_CODE(l0norm, unsigned long long, x);
-  MLN_GEN_CODE(l0norm, float, std::abs(x));
-  MLN_GEN_CODE(l0norm, double, std::abs(x));
-  MLN_GEN_CODE(l0norm, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(l0norm);
-
-  MLN_GEN_CODE(l1norm, int, std::abs(x));
-  MLN_GEN_CODE(l1norm, long, std::abs(x));
-  MLN_GEN_CODE(l1norm, long long, std::abs(x));
-  MLN_GEN_CODE(l1norm, unsigned int, x);
-  MLN_GEN_CODE(l1norm, unsigned long, x);
-  MLN_GEN_CODE(l1norm, unsigned long long, x);
-  MLN_GEN_CODE(l1norm, float, std::abs(x));
-  MLN_GEN_CODE(l1norm, double, std::abs(x));
-  MLN_GEN_CODE(l1norm, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(l1norm);
-
-  MLN_GEN_CODE(l2norm, int, std::abs(x));
-  MLN_GEN_CODE(l2norm, long, std::abs(x));
-  MLN_GEN_CODE(l2norm, long long, std::abs(x));
-  MLN_GEN_CODE(l2norm, unsigned int, x);
-  MLN_GEN_CODE(l2norm, unsigned long, x);
-  MLN_GEN_CODE(l2norm, unsigned long long, x);
-  MLN_GEN_CODE(l2norm, float, std::abs(x));
-  MLN_GEN_CODE(l2norm, double, std::abs(x));
-  MLN_GEN_CODE(l2norm, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(l2norm);
-
-  MLN_GEN_BINARY_CODE(l2dist, int, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, long long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, unsigned int, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, unsigned long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, unsigned long long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, float, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, double, abs(x - y));
-  MLN_GEN_BINARY_CODE(l2dist, long double, abs(x - y));
-  MLN_FUNCTIONAL_GEN_BINARY_CODE(l2dist);
-
-  MLN_GEN_BINARY_CODE(l1dist, int, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, long long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, unsigned int, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, unsigned long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, unsigned long long, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, float, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, double, abs(x - y));
-  MLN_GEN_BINARY_CODE(l1dist, long double, abs(x - y));
-  MLN_FUNCTIONAL_GEN_BINARY_CODE(l1dist);
-
-  MLN_GEN_CODE(l2norm_sqr, int, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, long, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, long long, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, unsigned int, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, unsigned long, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, unsigned long long, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, float, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, double, sqr(x));
-  MLN_GEN_CODE(l2norm_sqr, long double, sqr(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(l2norm_sqr);
-
-  MLN_GEN_BINARY_CODE(l2dist_sqr, int, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, long, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, long long, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, unsigned int, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, unsigned long, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, unsigned long long, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, float, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, double, sqr(x - y));
-  MLN_GEN_BINARY_CODE(l2dist_sqr, long double, sqr(x - y));
-  MLN_FUNCTIONAL_GEN_BINARY_CODE(l2dist_sqr);
-
-  MLN_GEN_CODE(linfnorm, int, std::abs(x));
-  MLN_GEN_CODE(linfnorm, long, std::abs(x));
-  MLN_GEN_CODE(linfnorm, long long, std::abs(x));
-  MLN_GEN_CODE(linfnorm, unsigned int, x);
-  MLN_GEN_CODE(linfnorm, unsigned long, x);
-  MLN_GEN_CODE(linfnorm, unsigned long long, x);
-  MLN_GEN_CODE(linfnorm, float, std::abs(x));
-  MLN_GEN_CODE(linfnorm, double, std::abs(x));
-  MLN_GEN_CODE(linfnorm, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(linfnorm);
-
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, int, std::abs(x));
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, long, std::abs(x));
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, long long, std::abs(x));
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, unsigned int, x);
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, unsigned long, x);
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, unsigned long long, x);
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, float, std::abs(x));
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, double, std::abs(x));
-  MLN_GEN_CODE_2(template <unsigned p>, lpnorm, long double, std::abs(x));
-  MLN_FUNCTIONAL_GEN_CODE_2(unsigned, p, lpnorm);
-
-  MLN_GEN_CODE(maximum, int, x);
-  MLN_GEN_CODE(maximum, long, x);
-  MLN_GEN_CODE(maximum, long long, x);
-  MLN_GEN_CODE(maximum, unsigned int, x);
-  MLN_GEN_CODE(maximum, unsigned long, x);
-  MLN_GEN_CODE(maximum, unsigned long long, x);
-  MLN_GEN_CODE(maximum, float, x);
-  MLN_GEN_CODE(maximum, double, x);
-  MLN_GEN_CODE(maximum, long double, x);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(maximum);
-
-  MLN_GEN_CODE(minimum, int, x);
-  MLN_GEN_CODE(minimum, long, x);
-  MLN_GEN_CODE(minimum, long long, x);
-  MLN_GEN_CODE(minimum, unsigned int, x);
-  MLN_GEN_CODE(minimum, unsigned long, x);
-  MLN_GEN_CODE(minimum, unsigned long long, x);
-  MLN_GEN_CODE(minimum, float, x);
-  MLN_GEN_CODE(minimum, double, x);
-  MLN_GEN_CODE(minimum, long double, x);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(minimum);
-
-  MLN_GEN_CODE(sum, int, x);
-  MLN_GEN_CODE(sum, long, x);
-  MLN_GEN_CODE(sum, long long, x);
-  MLN_GEN_CODE(sum, unsigned int, x);
-  MLN_GEN_CODE(sum, unsigned long, x);
-  MLN_GEN_CODE(sum, unsigned long long, x);
-  MLN_GEN_CODE(sum, float, x);
-  MLN_GEN_CODE(sum, double, x);
-  MLN_GEN_CODE(sum, long double, x);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(sum);
-
-  MLN_GEN_CODE(prod, int, x);
-  MLN_GEN_CODE(prod, long, x);
-  MLN_GEN_CODE(prod, long long, x);
-  MLN_GEN_CODE(prod, unsigned int, x);
-  MLN_GEN_CODE(prod, unsigned long, x);
-  MLN_GEN_CODE(prod, unsigned long long, x);
-  MLN_GEN_CODE(prod, float, x);
-  MLN_GEN_CODE(prod, double, x);
-  MLN_GEN_CODE(prod, long double, x);
-  MLN_FUNCTIONAL_GEN_UNARY_CODE(prod);
-
-#undef MLN_GEN_CODE
-#undef MLN_GEN_BINARY_CODE
-#undef MLN_GEN_CODE_2
-#undef MLN_FUNCTIONAL_GEN_UNARY_CODE
-#undef MLN_FUNCTIONAL_GEN_BINARY_CODE
-#undef MLN_FUNCTIONAL_GEN_CODE_2
-} // namespace mln
diff --git a/pylene/include/mln/core/ops.hpp b/pylene/include/mln/core/ops.hpp
index ba621223f87fc2c8e7054a2e1eb86df2f26a5aab..49246d400439d45cf553d25bd082af314e13dda4 100644
--- a/pylene/include/mln/core/ops.hpp
+++ b/pylene/include/mln/core/ops.hpp
@@ -2,6 +2,8 @@
 
 #include <algorithm>
 #include <functional>
+#include <tuple>
+#include <utility>
 
 
 /**
diff --git a/pylene/include/mln/core/private/def_functional_ops_gen_code_macros.hpp b/pylene/include/mln/core/private/def_functional_ops_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c283c62ebe0bbfacbf2f2440beb4d69c5a25d61
--- /dev/null
+++ b/pylene/include/mln/core/private/def_functional_ops_gen_code_macros.hpp
@@ -0,0 +1,43 @@
+
+#define MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE(FUN)                                                                     \
+  namespace functional                                                                                                 \
+  {                                                                                                                    \
+    using mln::FUN;                                                                                                    \
+    template <typename T>                                                                                              \
+    struct FUN##_t                                                                                                     \
+    {                                                                                                                  \
+      typedef decltype(FUN(std::declval<T>())) result_type;                                                            \
+                                                                                                                       \
+      auto operator()(const T& x) const -> decltype(FUN(x)) { return FUN(x); }                                         \
+    };                                                                                                                 \
+    template <>                                                                                                        \
+    struct FUN##_t<void>                                                                                               \
+    {                                                                                                                  \
+      template <typename T>                                                                                            \
+      auto operator()(const T& x) const -> decltype(FUN(x))                                                            \
+      {                                                                                                                \
+        return FUN(x);                                                                                                 \
+      }                                                                                                                \
+    };                                                                                                                 \
+  }
+
+#define MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE(FUN)                                                                    \
+  namespace functional                                                                                                 \
+  {                                                                                                                    \
+    using mln::FUN;                                                                                                    \
+    template <typename T1, typename T2>                                                                                \
+    struct FUN##_t                                                                                                     \
+    {                                                                                                                  \
+      typedef decltype(FUN(std::declval<T1>(), std::declval<T2>())) result_type;                                       \
+      auto operator()(const T1& x, const T2& y) const -> decltype(FUN(x, y)) { return FUN(x, y); }                     \
+    };                                                                                                                 \
+    template <>                                                                                                        \
+    struct FUN##_t<void, void>                                                                                         \
+    {                                                                                                                  \
+      template <typename T1, typename T2>                                                                              \
+      auto operator()(const T1& x, const T2& y) const -> decltype(FUN(x, y))                                           \
+      {                                                                                                                \
+        return FUN(x, y);                                                                                              \
+      }                                                                                                                \
+    };                                                                                                                 \
+  }
diff --git a/pylene/include/mln/core/private/def_maths_ops_gen_code_macros.hpp b/pylene/include/mln/core/private/def_maths_ops_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..338911a159febf2b140e8dd33b2320c495e84f45
--- /dev/null
+++ b/pylene/include/mln/core/private/def_maths_ops_gen_code_macros.hpp
@@ -0,0 +1,12 @@
+
+#define MLN_PRIVATE_GEN_UNARY_CODE(FUN, TYPE, OP)                                                                      \
+  inline TYPE FUN(TYPE x) { return OP; }
+
+#define MLN_PRIVATE_GEN_BINARY_CODE(FUN, TYPE, OP)                                                                     \
+  inline TYPE FUN(TYPE x, TYPE y) { return OP; }
+
+#define MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(TEMPLATE, FUN, TYPE, OP)                                                  \
+  TEMPLATE inline TYPE FUN(TYPE x) { return OP; }
+
+#define MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(TEMPLATE, FUN, TYPE, OP)                                                 \
+  TEMPLATE inline TYPE FUN(TYPE x, TYPE y) { return OP; }
diff --git a/pylene/include/mln/core/private/maths_ops.hpp b/pylene/include/mln/core/private/maths_ops.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..263e73659e70ede781e753f33adf52b4a9bdf54e
--- /dev/null
+++ b/pylene/include/mln/core/private/maths_ops.hpp
@@ -0,0 +1,388 @@
+#pragma once
+
+/**
+ * \file
+ * \brief Define fundamental mathematical and statistical operators
+ *
+ */
+
+#include <cmath>
+#include <cstdlib>
+#include <utility>
+
+
+#include <mln/core/private/def_maths_ops_gen_code_macros.hpp>
+
+namespace mln
+{
+  using std::cbrt;
+  using std::pow;
+  using std::sqrt;
+
+  /* Element wise operators */
+  int                sqr(int x);
+  long               sqr(long n);
+  long long          sqr(long long n);
+  unsigned int       sqr(unsigned int x);
+  unsigned long      sqr(unsigned long n);
+  unsigned long long sqr(unsigned long long n);
+  float              sqr(float x);
+  double             sqr(double x);
+  long double        sqr(long double x);
+
+  int                abs(int x);
+  long               abs(long n);
+  long long          abs(long long n);
+  unsigned int       abs(unsigned int x);
+  unsigned long      abs(unsigned long n);
+  unsigned long long abs(unsigned long long n);
+  float              abs(float x);
+  double             abs(double x);
+  long double        abs(long double x);
+
+  /* Reduction operators */
+  int                sum(int x);
+  long               sum(long n);
+  long long          sum(long long n);
+  unsigned int       sum(unsigned int x);
+  unsigned long      sum(unsigned long n);
+  unsigned long long sum(unsigned long long n);
+  float              sum(float x);
+  double             sum(double x);
+  long double        sum(long double x);
+
+  int                prod(int x);
+  long               prod(long n);
+  long long          prod(long long n);
+  unsigned int       prod(unsigned int x);
+  unsigned long      prod(unsigned long n);
+  unsigned long long prod(unsigned long long n);
+  float              prod(float x);
+  double             prod(double x);
+  long double        prod(long double x);
+
+  int                minimum(int x);
+  long               minimum(long n);
+  long long          minimum(long long n);
+  unsigned int       minimum(unsigned int x);
+  unsigned long      minimum(unsigned long n);
+  unsigned long long minimum(unsigned long long n);
+  float              minimum(float x);
+  double             minimum(double x);
+  long double        minimum(long double x);
+
+  int                maximum(int x);
+  long               maximum(long n);
+  long long          maximum(long long n);
+  unsigned int       maximum(unsigned int x);
+  unsigned long      maximum(unsigned long n);
+  unsigned long long maximum(unsigned long long n);
+  float              maximum(float x);
+  double             maximum(double x);
+  long double        maximum(long double x);
+
+  int                l0norm(int x);
+  long               l0norm(long n);
+  long long          l0norm(long long n);
+  unsigned int       l0norm(unsigned int x);
+  unsigned long      l0norm(unsigned long n);
+  unsigned long long l0norm(unsigned long long n);
+  float              l0norm(float x);
+  double             l0norm(double x);
+  long double        l0norm(long double x);
+
+  int                l1norm(int x);
+  long               l1norm(long n);
+  long long          l1norm(long long n);
+  unsigned int       l1norm(unsigned int x);
+  unsigned long      l1norm(unsigned long n);
+  unsigned long long l1norm(unsigned long long n);
+  float              l1norm(float x);
+  double             l1norm(double x);
+  long double        l1norm(long double x);
+
+  int                l2norm(int x);
+  long               l2norm(long n);
+  long long          l2norm(long long n);
+  unsigned int       l2norm(unsigned int x);
+  unsigned long      l2norm(unsigned long n);
+  unsigned long long l2norm(unsigned long long n);
+  float              l2norm(float x);
+  double             l2norm(double x);
+  long double        l2norm(long double x);
+
+  int                l2dist(int x, int y);
+  long               l2dist(long x, long y);
+  long long          l2dist(long long x, long long y);
+  unsigned int       l2dist(unsigned int x, unsigned int y);
+  unsigned long      l2dist(unsigned long x, unsigned long y);
+  unsigned long long l2dist(unsigned long long x, unsigned long long y);
+  float              l2dist(float x, float y);
+  double             l2dist(double x, double y);
+  long double        l2dist(long double x, long double y);
+
+  int                l1dist(int x, int y);
+  long               l1dist(long x, long y);
+  long long          l1dist(long long x, long long y);
+  unsigned int       l1dist(unsigned int x, unsigned int y);
+  unsigned long      l1dist(unsigned long x, unsigned long y);
+  unsigned long long l1dist(unsigned long long x, unsigned long long y);
+  float              l1dist(float x, float y);
+  double             l1dist(double x, double y);
+  long double        l1dist(long double x, long double y);
+
+  int                l2norm_sqr(int x);
+  long               l2norm_sqr(long n);
+  long long          l2norm_sqr(long long n);
+  unsigned int       l2norm_sqr(unsigned int x);
+  unsigned long      l2norm_sqr(unsigned long n);
+  unsigned long long l2norm_sqr(unsigned long long n);
+  float              l2norm_sqr(float x);
+  double             l2norm_sqr(double x);
+  long double        l2norm_sqr(long double x);
+
+  int                l2dist_sqr(int x, int y);
+  long               l2dist_sqr(long x, long y);
+  long long          l2dist_sqr(long long x, long long y);
+  unsigned int       l2dist_sqr(unsigned int x, unsigned int y);
+  unsigned long      l2dist_sqr(unsigned long x, unsigned long y);
+  unsigned long long l2dist_sqr(unsigned long long x, unsigned long long y);
+  float              l2dist_sqr(float x, float y);
+  double             l2dist_sqr(double x, double y);
+  long double        l2dist_sqr(long double x, long double y);
+
+  int                linfnorm(int x);
+  long               linfnorm(long n);
+  long long          linfnorm(long long n);
+  unsigned int       linfnorm(unsigned int x);
+  unsigned long      linfnorm(unsigned long n);
+  unsigned long long linfnorm(unsigned long long n);
+  float              linfnorm(float x);
+  double             linfnorm(double x);
+  long double        linfnorm(long double x);
+
+  template <unsigned p>
+  int lpnorm(int x);
+  template <unsigned p>
+  long lpnorm(long x);
+  template <unsigned p>
+  long long lpnorm(long long x);
+  template <unsigned p>
+  unsigned int lpnorm(unsigned int x);
+  template <unsigned p>
+  unsigned long lpnorm(unsigned long n);
+  template <unsigned p>
+  unsigned long long lpnorm(unsigned long long n);
+  template <unsigned p>
+  float lpnorm(float x);
+  template <unsigned p>
+  double lpnorm(double x);
+  template <unsigned p>
+  long double lpnorm(long double x);
+
+  template <unsigned p>
+  int lpdist(int x, int y);
+  template <unsigned p>
+  long lpdist(long x, long y);
+  template <unsigned p>
+  long long lpdist(long long x, long long y);
+  template <unsigned p>
+  unsigned int lpdist(unsigned int x, unsigned int y);
+  template <unsigned p>
+  unsigned long lpdist(unsigned long x, unsigned long y);
+  template <unsigned p>
+  unsigned long long lpdist(unsigned long long x, unsigned long long y);
+  template <unsigned p>
+  float lpdist(float x, float y);
+  template <unsigned p>
+  double lpdist(double x, double y);
+  template <unsigned p>
+  long double lpdist(long double x, long double y);
+
+
+  /********************************/
+  /*** Implementation            **/
+  /********************************/
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, int, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, long, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, long long, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, unsigned int, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, unsigned long, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, unsigned long long, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, float, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, double, x* x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sqr, long double, x* x);
+
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, int, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, long long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, float, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, double, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(abs, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, int, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, long long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, float, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, double, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l0norm, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, int, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, long long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, float, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, double, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l1norm, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, int, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, long long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, float, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, double, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, int, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, long, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, long long, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, unsigned int, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, unsigned long, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, unsigned long long, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, float, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, double, abs(x - y));
+  MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE(template <unsigned p>, lpdist, long double, abs(x - y));
+
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, unsigned int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, unsigned long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, unsigned long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, float, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, double, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist, long double, abs(x - y));
+
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, unsigned int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, unsigned long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, unsigned long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, float, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, double, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l1dist, long double, abs(x - y));
+
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, unsigned int, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, unsigned long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, unsigned long long, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, float, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, double, abs(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l0dist, long double, abs(x - y));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, int, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, long, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, long long, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, unsigned int, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, unsigned long, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, unsigned long long, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, float, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, double, sqr(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(l2norm_sqr, long double, sqr(x));
+
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, int, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, long long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, unsigned int, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, unsigned long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, unsigned long long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, float, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, double, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, long double, sqr(x - y));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, int, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, long long, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, float, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, double, std::abs(x));
+  MLN_PRIVATE_GEN_UNARY_CODE(linfnorm, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, int, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, long long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, unsigned int, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, unsigned long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, unsigned long long, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, float, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, double, sqr(x - y));
+  MLN_PRIVATE_GEN_BINARY_CODE(linfdist, long double, sqr(x - y));
+
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, int, std::abs(x));
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, long, std::abs(x));
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, long long, std::abs(x));
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, unsigned int, x);
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, unsigned long, x);
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, unsigned long long, x);
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, float, std::abs(x));
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, double, std::abs(x));
+  MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE(template <unsigned p>, lpnorm, long double, std::abs(x));
+
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, float, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, double, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(maximum, long double, x);
+
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, float, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, double, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(minimum, long double, x);
+
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, float, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, double, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(sum, long double, x);
+
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, unsigned int, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, unsigned long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, unsigned long long, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, float, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, double, x);
+  MLN_PRIVATE_GEN_UNARY_CODE(prod, long double, x);
+
+} // namespace mln
+
+#include <mln/core/private/undef_maths_ops_gen_code_macros.hpp>
diff --git a/pylene/include/mln/core/private/undef_functional_ops_gen_code_macros.hpp b/pylene/include/mln/core/private/undef_functional_ops_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff5b83875fa373f8d5da2ff7f48a18d798f38cec
--- /dev/null
+++ b/pylene/include/mln/core/private/undef_functional_ops_gen_code_macros.hpp
@@ -0,0 +1,3 @@
+
+#undef MLN_PRIVATE_FUNCTIONAL_GEN_UNARY_CODE
+#undef MLN_PRIVATE_FUNCTIONAL_GEN_BINARY_CODE
diff --git a/pylene/include/mln/core/private/undef_maths_ops_gen_code_macros.hpp b/pylene/include/mln/core/private/undef_maths_ops_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..76927dc59982867696af845916de9de21482cd15
--- /dev/null
+++ b/pylene/include/mln/core/private/undef_maths_ops_gen_code_macros.hpp
@@ -0,0 +1,5 @@
+
+#undef MLN_PRIVATE_GEN_UNARY_CODE
+#undef MLN_PRIVATE_GEN_BINARY_CODE
+#undef MLN_PRIVATE_GEN_TEMPLATED_UNARY_CODE
+#undef MLN_PRIVATE_GEN_TEMPLATED_BINARY_CODE
diff --git a/pylene/include/mln/core/rangev3/views.hpp b/pylene/include/mln/core/rangev3/views.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f751e4fdeaeb6ef5e41dccd93c0508d757fed4b
--- /dev/null
+++ b/pylene/include/mln/core/rangev3/views.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+#include <mln/core/rangev3/view/filter.hpp>
+#include <mln/core/rangev3/view/remove_if.hpp>
+#include <mln/core/rangev3/view/reverse.hpp>
+#include <mln/core/rangev3/view/transform.hpp>
+#include <mln/core/rangev3/view/zip.hpp>
+#include <mln/core/rangev3/view/zip_with.hpp>
diff --git a/pylene/include/mln/core/vec.hpp b/pylene/include/mln/core/vec.hpp
index 2477814ed5be229a0883b4e5ea57ff6db3a84aa9..a7e779e8cc70c9e5b513bf8aaf62adf1933a93f8 100644
--- a/pylene/include/mln/core/vec.hpp
+++ b/pylene/include/mln/core/vec.hpp
@@ -3,7 +3,7 @@
 #include <mln/core/value/value_traits.hpp>
 #include <mln/core/vec/compare.hpp>
 #include <mln/core/vec/vec_io.hpp>
-#include <mln/core/vec/vec_math_ops.hpp>
+#include <mln/core/vec/vec_maths_ops.hpp>
 #include <mln/core/vec/vec_ops.hpp>
 #include <mln/core/vec_base.hpp>
 
diff --git a/pylene/include/mln/core/vec/private/def_gen_code_macros.hpp b/pylene/include/mln/core/vec/private/def_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6bc95ed42139103291cac6b1e3bbc470ffbf096d
--- /dev/null
+++ b/pylene/include/mln/core/vec/private/def_gen_code_macros.hpp
@@ -0,0 +1,26 @@
+
+#define MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, fun) internal::vec_base<decltype(fun(std::declval<T>())), dim, tag>
+
+#define MLN_VEC_PRIVATE_PROMOTE(T1, T2) decltype(std::declval<T1>() + std::declval<T2>())
+
+#define MLN_VEC_PRIVATE_PROMOTE_FUN(T, fun) decltype(fun(std::declval<T>()))
+
+#define MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, f, g) decltype(f(g(std::declval<T>())))
+
+#define MLN_VEC_PRIVATE_GEN_UNARY_CODE(FUN)                                                                            \
+  template <class T, unsigned dim, class tag>                                                                          \
+  inline MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, FUN) FUN(const internal::vec_base<T, dim, tag>& x)               \
+  {                                                                                                                    \
+    MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, FUN) res;                                                             \
+    for (unsigned i = 0; i < dim; ++i)                                                                                 \
+      res[i] = FUN(x[i]);                                                                                              \
+    return res;                                                                                                        \
+  }
+
+#define MLN_VEC_PRIVATE_GEN_BINARY_CODE(FUN, FUNBASE, EXPR)                                                            \
+  template <class T, unsigned dim, class tag>                                                                          \
+  inline auto FUN(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)                  \
+      ->decltype(FUNBASE(EXPR))                                                                                        \
+  {                                                                                                                    \
+    return FUNBASE(EXPR);                                                                                              \
+  }
diff --git a/pylene/include/mln/core/vec/private/undef_gen_code_macros.hpp b/pylene/include/mln/core/vec/private/undef_gen_code_macros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ddf412c3ac9a1be012e1a7efb9f23044df1202b0
--- /dev/null
+++ b/pylene/include/mln/core/vec/private/undef_gen_code_macros.hpp
@@ -0,0 +1,7 @@
+
+#undef MLN_VEC_PRIVATE_PROMOTE_FUN
+#undef MLN_VEC_PRIVATE_PROMOTE
+#undef MLN_VEC_PRIVATE_PROMOTE_FUN
+#undef MLN_VEC_PRIVATE_PROMOTE_FUNCOMP
+#undef MLN_VEC_PRIVATE_GEN_UNARY_CODE
+#undef MLN_VEC_PRIVATE_GEN_BINARY_CODE
diff --git a/pylene/include/mln/core/vec/vec_math_ops.hpp b/pylene/include/mln/core/vec/vec_math_ops.hpp
deleted file mode 100644
index f0dfb5eca53cc09a27148b49778606fe73dbbc7f..0000000000000000000000000000000000000000
--- a/pylene/include/mln/core/vec/vec_math_ops.hpp
+++ /dev/null
@@ -1,276 +0,0 @@
-#pragma once
-
-#include <mln/core/math_ops.hpp>
-#include <mln/core/vec.hpp>
-
-
-/**
- * \file Define vectorial overloads for fundamental mathematical and
- * statistical operators
- *
- */
-
-#define MLN_VEC_PROMOTE_FUN(T, dim, tag, fun) internal::vec_base<decltype(fun(std::declval<T>())), dim, tag>
-
-#define MLN_PROMOTE(T1, T2) decltype(std::declval<T1>() + std::declval<T2>())
-
-#define MLN_PROMOTE_FUN(T, fun) decltype(fun(std::declval<T>()))
-
-#define MLN_PROMOTE_FUNCOMP(T, f, g) decltype(f(g(std::declval<T>())))
-
-namespace mln
-{
-
-  /* Element wise operators */
-  template <class U, class V, unsigned dim, class tag>
-  internal::vec_base<decltype(pow(std::declval<U>(), std::declval<V>())), dim, tag>
-      pow(const internal::vec_base<U, dim, tag>& x, V exp);
-
-  template <class T, unsigned dim, class tag>
-  MLN_VEC_PROMOTE_FUN(T, dim, tag, sqr)
-  sqr(const internal::vec_base<T, dim, tag>& x);
-
-  template <class T, unsigned dim, class tag>
-  MLN_VEC_PROMOTE_FUN(T, dim, tag, abs)
-  abs(const internal::vec_base<T, dim, tag>& x);
-
-  template <class T, unsigned dim, class tag>
-  MLN_VEC_PROMOTE_FUN(T, dim, tag, sqrt)
-  sqrt(const internal::vec_base<T, dim, tag>& x);
-
-  template <class T, unsigned dim, class tag>
-  MLN_VEC_PROMOTE_FUN(T, dim, tag, cqrt)
-  cbrt(const internal::vec_base<T, dim, tag>& x);
-
-  /* algebraic operators */
-  template <class U, class V, unsigned dim, class tag>
-  MLN_PROMOTE(U, V)
-  dot(const internal::vec_base<U, dim, tag>& x, const internal::vec_base<V, dim, tag>& y);
-
-  template <class U, class V, class tag>
-  internal::vec_base<MLN_PROMOTE(U, V), 3, tag> cross(const internal::vec_base<U, 3, tag>& x,
-                                                      const internal::vec_base<V, 3, tag>& y);
-
-  /* Reduction operators */
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE(T, T)
-  sum(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE(T, T)
-  prod(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  T maximum(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  T minimum(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUN(T, abs)
-  l0norm(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUN(T, abs)
-  linfnorm(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUN(T, abs)
-  l1norm(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUNCOMP(T, sqrt, sqr)
-  l2norm(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUN(T, sqr) l2norm_sqr(const internal::vec_base<T, dim, tag>& x);
-
-  template <unsigned p, typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUNCOMP(T, pow, abs)
-  lpnorm(const internal::vec_base<T, dim, tag>& x);
-
-  template <typename T, unsigned dim, class tag>
-  inline auto l1dist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
-      -> decltype(l1norm(x - y));
-
-  template <typename T, unsigned dim, class tag>
-  inline auto l2dist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
-      -> decltype(l2norm(x - y));
-
-  template <typename T, unsigned dim, class tag>
-  inline auto l2dist_sqr(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
-      -> decltype(l2norm_sqr(x - y));
-
-  /*****************************/
-  /**   Implementation       ***/
-  /*****************************/
-
-#define MLN_GEN_CODE(FUN)                                                                                              \
-  template <class T, unsigned dim, class tag>                                                                          \
-  inline MLN_VEC_PROMOTE_FUN(T, dim, tag, FUN) FUN(const internal::vec_base<T, dim, tag>& x)                           \
-  {                                                                                                                    \
-    MLN_VEC_PROMOTE_FUN(T, dim, tag, FUN) res;                                                                         \
-    for (unsigned i = 0; i < dim; ++i)                                                                                 \
-      res[i] = FUN(x[i]);                                                                                              \
-    return res;                                                                                                        \
-  }
-
-#define MLN_GEN_BINARY_CODE(FUN, FUNBASE, EXPR)                                                                        \
-  template <class T, unsigned dim, class tag>                                                                          \
-  inline auto FUN(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)                  \
-      ->decltype(FUNBASE(EXPR))                                                                                        \
-  {                                                                                                                    \
-    return FUNBASE(EXPR);                                                                                              \
-  }
-
-  MLN_GEN_CODE(sqr);
-  MLN_GEN_CODE(sqrt);
-  MLN_GEN_CODE(abs);
-  MLN_GEN_CODE(cbrt);
-
-  MLN_GEN_BINARY_CODE(l1dist, l1norm, x - y);
-  MLN_GEN_BINARY_CODE(l2dist, l2norm, x - y);
-  MLN_GEN_BINARY_CODE(l2dist_sqr, l2norm_sqr, x - y);
-
-  template <class U, class V, unsigned dim, class tag>
-  internal::vec_base<decltype(pow(std::declval<U>(), std::declval<V>())), dim, tag>
-      pow(const internal::vec_base<U, dim, tag>& x, V exp)
-  {
-    typedef decltype(pow(std::declval<U>(), std::declval<V>())) R;
-    internal::vec_base<R, dim, tag>                             res;
-    for (unsigned i = 0; i < dim; ++i)
-      res[i] = pow(x[i], exp);
-    return res;
-  }
-
-  template <class U, class V, unsigned dim, class tag>
-  MLN_PROMOTE(U, V)
-  dot(const internal::vec_base<U, dim, tag>& x, const internal::vec_base<V, dim, tag>& y)
-  {
-    MLN_PROMOTE(U, V) res = 0;
-    for (int i = 0; i < dim; ++i)
-      res += x[i] * y[i];
-    return res;
-  }
-
-  template <class U, class V, class tag>
-  internal::vec_base<MLN_PROMOTE(U, V), 3, tag> cross(const internal::vec_base<U, 3, tag>& u,
-                                                      const internal::vec_base<V, 3, tag>& v)
-  {
-    return {u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]};
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE(T, T) sum(const internal::vec_base<T, dim, tag>& x)
-  {
-    MLN_PROMOTE(T, T) res = x[0];
-    for (unsigned i = 1; i < dim; ++i)
-      res += x[i];
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE(T, T) prod(const internal::vec_base<T, dim, tag>& x)
-  {
-    MLN_PROMOTE(T, T) res = x[0];
-    for (unsigned i = 1; i < dim; ++i)
-      res *= x[i];
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline T maximum(const internal::vec_base<T, dim, tag>& x)
-  {
-    T res = x[0];
-    for (unsigned i = 1; i < dim; ++i)
-      if (res < x[i])
-        res = x[i];
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline T minimum(const internal::vec_base<T, dim, tag>& x)
-  {
-    T res = x[0];
-    for (unsigned i = 1; i < dim; ++i)
-      if (x[i] < res)
-        res = x[i];
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUN(T, abs) l0norm(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUN(T, abs) U;
-    U res = abs(x[0]);
-    for (unsigned i = 1; i < dim; ++i)
-    {
-      U v = abs(x[i]);
-      if (v < res)
-        res = v;
-    }
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUN(T, abs) linfnorm(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUN(T, abs) U;
-    U res = abs(x[0]);
-    for (unsigned i = 1; i < dim; ++i)
-    {
-      U v = abs(x[i]);
-      if (res < v)
-        res = v;
-    }
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUN(T, abs) l1norm(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUN(T, abs) U;
-    U res = abs(x[0]);
-    for (unsigned i = 1; i < dim; ++i)
-      res += abs(x[i]);
-    return res;
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUNCOMP(T, sqrt, sqr) l2norm(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUNCOMP(T, sqrt, sqr) U;
-    U res = sqr(x[0]);
-    for (unsigned i = 1; i < dim; ++i)
-      res += sqr(x[i]);
-    return sqrt(res);
-  }
-
-  template <typename T, unsigned dim, class tag>
-  inline MLN_PROMOTE_FUN(T, sqr) l2norm_sqr(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUNCOMP(T, sqrt, sqr) U;
-    U res = sqr(x[0]);
-    for (unsigned i = 1; i < dim; ++i)
-      res += sqr(x[i]);
-    return res;
-  }
-
-  template <unsigned p, typename T, unsigned dim, class tag>
-  MLN_PROMOTE_FUNCOMP(T, pow, abs)
-  lpnorm(const internal::vec_base<T, dim, tag>& x)
-  {
-    typedef MLN_PROMOTE_FUNCOMP(T, pow, abs) U;
-    U res = pow(abs(x[0]), p);
-    for (unsigned i = 1; i < dim; ++i)
-      res += pow(abs(x[i]), p);
-    return pow(res, 1 / p);
-  }
-} // namespace mln
-
-#undef MLN_GEN_CODE
-#undef MLN_GEN_BINARY_CODE
-#undef MLN_VEC_PROMOTE_FUN
-#undef MLN_VEC_PROMOTE
-#undef MLN_PROMOTE_FUNCOMP
diff --git a/pylene/include/mln/core/vec/vec_maths_ops.hpp b/pylene/include/mln/core/vec/vec_maths_ops.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..60119ac997125c03bd8ddd3528a84d7e42a5fc3e
--- /dev/null
+++ b/pylene/include/mln/core/vec/vec_maths_ops.hpp
@@ -0,0 +1,295 @@
+#pragma once
+
+#include <mln/core/private/maths_ops.hpp>
+#include <mln/core/vec/vec_ops.hpp>
+
+/**
+ * \file Define vectorial overloads for fundamental mathematical and
+ * statistical operators
+ *
+ */
+
+#include <mln/core/vec/private/def_gen_code_macros.hpp>
+
+namespace mln
+{
+  /* Element wise operators */
+  template <class U, class V, unsigned dim, class tag>
+  internal::vec_base<decltype(pow(std::declval<U>(), std::declval<V>())), dim, tag>
+      pow(const internal::vec_base<U, dim, tag>& x, V exp);
+
+  template <class T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, sqr)
+  sqr(const internal::vec_base<T, dim, tag>& x);
+
+  template <class T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, abs)
+  abs(const internal::vec_base<T, dim, tag>& x);
+
+  template <class T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, sqrt)
+  sqrt(const internal::vec_base<T, dim, tag>& x);
+
+  template <class T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_VEC_PROMOTE_FUN(T, dim, tag, cqrt)
+  cbrt(const internal::vec_base<T, dim, tag>& x);
+
+  /* algebraic operators */
+  template <class U, class V, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE(U, V)
+  dot(const internal::vec_base<U, dim, tag>& x, const internal::vec_base<V, dim, tag>& y);
+
+  template <class U, class V, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), 3, tag> cross(const internal::vec_base<U, 3, tag>& x,
+                                                                  const internal::vec_base<V, 3, tag>& y);
+  template <class U, class V, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), 7, tag> cross(const internal::vec_base<U, 7, tag>& x,
+                                                                  const internal::vec_base<V, 7, tag>& y);
+  template <class U, class V, unsigned dim, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), dim, tag> cross(const internal::vec_base<U, dim, tag>& x,
+                                                                    const internal::vec_base<V, dim, tag>& y);
+
+  /* Reduction operators */
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE(T, T)
+  sum(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE(T, T)
+  prod(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  T maximum(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  T minimum(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs)
+  l0norm(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs)
+  linfnorm(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs)
+  l1norm(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, sqrt, sqr)
+  l2norm(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUN(T, sqr) l2norm_sqr(const internal::vec_base<T, dim, tag>& x);
+
+  template <unsigned p, typename T, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, pow, abs)
+  lpnorm(const internal::vec_base<T, dim, tag>& x);
+
+  template <typename T, unsigned dim, class tag>
+  inline auto l0dist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(l0norm(x - y));
+
+  template <typename T, unsigned dim, class tag>
+  inline auto l1dist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(l1norm(x - y));
+
+  template <typename T, unsigned dim, class tag>
+  inline auto l2dist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(l2norm(x - y));
+
+  template <unsigned p, typename T, unsigned dim, class tag>
+  inline auto lpdist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(lpnorm<p>(x - y));
+
+  template <typename T, unsigned dim, class tag>
+  inline auto l2dist_sqr(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(l2norm_sqr(x - y));
+
+  template <typename T, unsigned dim, class tag>
+  inline auto linfdist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(linfnorm(x - y));
+
+  /*****************************/
+  /**   Implementation       ***/
+  /*****************************/
+
+  MLN_VEC_PRIVATE_GEN_UNARY_CODE(sqr);
+  MLN_VEC_PRIVATE_GEN_UNARY_CODE(sqrt);
+  MLN_VEC_PRIVATE_GEN_UNARY_CODE(abs);
+  MLN_VEC_PRIVATE_GEN_UNARY_CODE(cbrt);
+
+  MLN_VEC_PRIVATE_GEN_BINARY_CODE(l0dist, l0norm, x - y);
+  MLN_VEC_PRIVATE_GEN_BINARY_CODE(l1dist, l1norm, x - y);
+  MLN_VEC_PRIVATE_GEN_BINARY_CODE(l2dist, l2norm, x - y);
+  MLN_VEC_PRIVATE_GEN_BINARY_CODE(l2dist_sqr, l2norm_sqr, x - y);
+  MLN_VEC_PRIVATE_GEN_BINARY_CODE(linfdist, linfnorm, x - y);
+
+  template <unsigned p, class T, unsigned dim, class tag>
+  inline auto lpdist(const internal::vec_base<T, dim, tag>& x, const internal::vec_base<T, dim, tag>& y)
+      -> decltype(lpnorm<p>(x - y))
+  {
+    return lpnorm<p>(x - y);
+  }
+
+
+  template <class U, class V, unsigned dim, class tag>
+  internal::vec_base<decltype(pow(std::declval<U>(), std::declval<V>())), dim, tag>
+      pow(const internal::vec_base<U, dim, tag>& x, V exp)
+  {
+    typedef decltype(pow(std::declval<U>(), std::declval<V>())) R;
+
+    internal::vec_base<R, dim, tag> res;
+    for (unsigned i = 0; i < dim; ++i)
+      res[i] = pow(x[i], exp);
+    return res;
+  }
+
+  template <class U, class V, unsigned dim, class tag>
+  MLN_VEC_PRIVATE_PROMOTE(U, V)
+  dot(const internal::vec_base<U, dim, tag>& x, const internal::vec_base<V, dim, tag>& y)
+  {
+    MLN_VEC_PRIVATE_PROMOTE(U, V) res = 0;
+    for (unsigned i = 0; i < dim; ++i)
+      res += x[i] * y[i];
+    return res;
+  }
+
+  template <class U, class V, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), 3, tag> cross(const internal::vec_base<U, 3, tag>& u,
+                                                                  const internal::vec_base<V, 3, tag>& v)
+  {
+    return {u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]};
+  }
+
+  template <class U, class V, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), 7, tag> cross(const internal::vec_base<U, 7, tag>& u,
+                                                                  const internal::vec_base<V, 7, tag>& v)
+  {
+    // cf. https://en.wikipedia.org/wiki/Seven-dimensional_cross_product#Coordinate_expressions
+    return {u[1] * v[3] - u[3] * v[1] + u[2] * v[6] - u[6] * v[2] + u[4] * v[5] - u[5] * v[4],
+            u[2] * v[4] - u[4] * v[2] + u[3] * v[0] - u[0] * v[3] + u[5] * v[6] - u[6] * v[5],
+            u[3] * v[5] - u[5] * v[3] + u[4] * v[1] - u[1] * v[4] + u[6] * v[0] - u[0] * v[6],
+            u[4] * v[6] - u[6] * v[4] + u[5] * v[2] - u[2] * v[5] + u[0] * v[1] - u[1] * v[0],
+            u[5] * v[0] - u[0] * v[5] + u[6] * v[3] - u[3] * v[6] + u[1] * v[2] - u[2] * v[1],
+            u[6] * v[1] - u[1] * v[6] + u[0] * v[4] - u[4] * v[0] + u[2] * v[3] - u[3] * v[2],
+            u[0] * v[2] - u[2] * v[0] + u[1] * v[5] - u[5] * v[1] + u[3] * v[4] - u[4] * v[3]};
+  }
+
+  template <class U, class V, unsigned dim, class tag>
+  internal::vec_base<MLN_VEC_PRIVATE_PROMOTE(U, V), dim, tag>
+      cross([[maybe_unused]] const internal::vec_base<U, dim, tag>& u,
+            [[maybe_unused]] const internal::vec_base<V, dim, tag>& v)
+  {
+    static_assert(dim != 3 && dim != 7, "cross product is only defined for vector of size 3 or 7");
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE(T, T) sum(const internal::vec_base<T, dim, tag>& x)
+  {
+    MLN_VEC_PRIVATE_PROMOTE(T, T) res = x[0];
+    for (unsigned i = 1; i < dim; ++i)
+      res += x[i];
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE(T, T) prod(const internal::vec_base<T, dim, tag>& x)
+  {
+    MLN_VEC_PRIVATE_PROMOTE(T, T) res = x[0];
+    for (unsigned i = 1; i < dim; ++i)
+      res *= x[i];
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline T maximum(const internal::vec_base<T, dim, tag>& x)
+  {
+    T res = x[0];
+    for (unsigned i = 1; i < dim; ++i)
+      if (res < x[i])
+        res = x[i];
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline T minimum(const internal::vec_base<T, dim, tag>& x)
+  {
+    T res = x[0];
+    for (unsigned i = 1; i < dim; ++i)
+      if (x[i] < res)
+        res = x[i];
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) l0norm(const internal::vec_base<T, dim, tag>& x)
+  {
+    typedef MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) U;
+    U res = abs(x[0]);
+    for (unsigned i = 1; i < dim; ++i)
+    {
+      U v = abs(x[i]);
+      if (v < res)
+        res = v;
+    }
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) linfnorm(const internal::vec_base<T, dim, tag>& x)
+  {
+    typedef MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) U;
+    U res = abs(x[0]);
+    for (unsigned i = 1; i < dim; ++i)
+    {
+      U v = abs(x[i]);
+      if (res < v)
+        res = v;
+    }
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) l1norm(const internal::vec_base<T, dim, tag>& x)
+  {
+    typedef MLN_VEC_PRIVATE_PROMOTE_FUN(T, abs) U;
+    U res = abs(x[0]);
+    for (unsigned i = 1; i < dim; ++i)
+      res += abs(x[i]);
+    return res;
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, sqrt, sqr) l2norm(const internal::vec_base<T, dim, tag>& x)
+  {
+    typedef MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, sqrt, sqr) U;
+    U res = sqr(x[0]);
+    for (unsigned i = 1; i < dim; ++i)
+      res += sqr(x[i]);
+    return sqrt(res);
+  }
+
+  template <typename T, unsigned dim, class tag>
+  inline MLN_VEC_PRIVATE_PROMOTE_FUN(T, sqr) l2norm_sqr(const internal::vec_base<T, dim, tag>& x)
+  {
+    typedef MLN_VEC_PRIVATE_PROMOTE_FUNCOMP(T, sqrt, sqr) U;
+    U res = sqr(x[0]);
+    for (unsigned i = 1; i < dim; ++i)
+      res += sqr(x[i]);
+    return res;
+  }
+
+  template <unsigned p, typename T, unsigned dim, class tag>
+  auto lpnorm(const internal::vec_base<T, dim, tag>& x)
+  {
+    double res = pow(abs(x[0]), p);
+    for (unsigned i = 1; i < dim; ++i)
+      res += pow(abs(x[i]), p);
+    return pow(res, 1. / static_cast<double>(p));
+  }
+} // namespace mln
+
+#include <mln/core/vec/private/undef_gen_code_macros.hpp>
diff --git a/pylene/include/mln/io/imprint.hpp b/pylene/include/mln/io/imprint.hpp
index 403d854bb9d238820f7b8b50840d4e5666c64308..a0ac4871d76c834c5482da3f9f33ad4d0e7c3500 100644
--- a/pylene/include/mln/io/imprint.hpp
+++ b/pylene/include/mln/io/imprint.hpp
@@ -2,6 +2,7 @@
 
 #include <mln/core/domain/box.hpp>
 #include <mln/core/grays.hpp>
+#include <mln/core/rangev3/rows.hpp>
 #include <mln/core/value/value_traits.hpp>
 #include <mln/io/format.hpp>
 
@@ -23,6 +24,16 @@ namespace mln
     template <typename I>
     void imprint_with_border(const Image<I>& ima, std::ostream& os = std::cout);
 
+
+    namespace experimental
+    {
+      template <typename InputImage>
+      void imprint(InputImage ima, std::ostream& os = std::cout);
+
+      template <typename InputImage>
+      void imprint_with_border(InputImage ima, std::ostream& os = std::cout);
+    } // namespace experimental
+
     /******************************************/
     /****          Implementation          ****/
     /******************************************/
@@ -36,7 +47,14 @@ namespace mln
         for (auto p : domain)
         {
           os << '{' << p << ",";
-          format(os, ima(p)) << "}," << std::endl;
+          if constexpr (mln::is_a<Image, mln::experimental::Image>())
+          {
+            format(os, (*const_cast<Image*>(&ima))(p)) << "}," << std::endl;
+          }
+          else
+          {
+            format(os, ima(p)) << "}," << std::endl;
+          }
         }
       }
 
@@ -45,10 +63,24 @@ namespace mln
           imprint(const Image& ima, Domain domain, std::ostream& os)
       {
         (void)domain;
-        mln_foreach (const auto& pix, ima.pixels())
+
+        if constexpr (mln::is_a<Image, mln::experimental::Image>())
         {
-          os << '{' << pix.point() << ",";
-          format(os, pix.val()) << "}," << std::endl;
+          auto pixels = const_cast<Image*>(&ima)->new_pixels();
+          for (auto&& r : mln::ranges::rows(pixels))
+            for (auto&& pix : r)
+            {
+              os << '{' << pix.point() << ",";
+              format(os, pix.val()) << "}," << std::endl;
+            }
+        }
+        else
+        {
+          mln_foreach (const auto& pix, ima.pixels())
+          {
+            os << '{' << pix.point() << ",";
+            format(os, pix.val()) << "}," << std::endl;
+          }
         }
       }
 
@@ -63,8 +95,21 @@ namespace mln
 
         int               wtext = 0;
         frmt_max_width<V> frmter;
-        mln_foreach (V v, ima.values())
-          wtext = std::max(wtext, frmter(v));
+
+        if constexpr (mln::is_a<Image, mln::experimental::Image>())
+        {
+          auto vals = const_cast<Image*>(&ima)->new_values();
+          for (auto&& r : mln::ranges::rows(vals))
+            for (auto&& v : r)
+            {
+              wtext = std::max(wtext, frmter(v));
+            }
+        }
+        else
+        {
+          mln_foreach (V v, ima.values())
+            wtext = std::max(wtext, frmter(v));
+        }
 
         os << domain << "(" << typeid(V).name() << ")" << std::endl;
         os.width(4);
@@ -75,7 +120,14 @@ namespace mln
           for (p[1] = domain.pmin[1]; p[1] < domain.pmax[1]; ++p[1])
           {
             os << std::setw(wtext);
-            format(os, ima(p)) << " ";
+            if constexpr (mln::is_a<Image, mln::experimental::Image>())
+            {
+              format(os, (*const_cast<Image*>(&ima))(p)) << " ";
+            }
+            else
+            {
+              format(os, ima(p)) << " ";
+            }
           }
           os << std::endl;
         }
@@ -93,8 +145,21 @@ namespace mln
 
         int               wtext = 0;
         frmt_max_width<V> frmter;
-        mln_foreach (V v, ima.values())
-          wtext = std::max(wtext, frmter(v));
+
+        if constexpr (mln::is_a<Image, mln::experimental::Image>())
+        {
+          auto vals = const_cast<Image*>(&ima)->new_values();
+          for (auto&& r : mln::ranges::rows(vals))
+            for (auto&& v : r)
+            {
+              wtext = std::max(wtext, frmter(v));
+            }
+        }
+        else
+        {
+          mln_foreach (V v, ima.values())
+            wtext = std::max(wtext, frmter(v));
+        }
 
         os << domain << "(" << typeid(V).name() << ")" << std::endl;
         os.width(4);
@@ -107,7 +172,14 @@ namespace mln
             for (p[2] = domain.pmin[2]; p[2] < domain.pmax[2]; ++p[2])
             {
               os << std::setw(wtext);
-              format(os, ima(p)) << " ";
+              if constexpr (mln::is_a<Image, mln::experimental::Image>())
+              {
+                format(os, (*const_cast<Image*>(&ima))(p)) << " ";
+              }
+              else
+              {
+                format(os, ima(p)) << " ";
+              }
             }
             os << std::endl;
           }
@@ -127,8 +199,21 @@ namespace mln
 
         int               wtext = 0;
         frmt_max_width<V> frmter;
-        mln_foreach (V v, ima.values())
-          wtext = std::max(wtext, frmter(v));
+
+        if constexpr (mln::is_a<Image, mln::experimental::Image>())
+        {
+          auto vals = ima.new_values();
+          for (auto&& r : mln::ranges::rows(vals))
+            for (auto&& v : r)
+            {
+              wtext = std::max(wtext, frmter(v));
+            }
+        }
+        else
+        {
+          mln_foreach (V v, ima.values())
+            wtext = std::max(wtext, frmter(v));
+        }
 
         os << domain << "(" << typeid(V).name() << ")" << std::endl;
         os.width(4);
@@ -140,7 +225,14 @@ namespace mln
           for (p[1] = domain.pmin[1] - border; p[1] < domain.pmax[1] + border; ++p[1])
           {
             os << std::setw(wtext);
-            format(os, ima.at(p)) << " ";
+            if constexpr (mln::is_a<Image, mln::experimental::Image>())
+            {
+              format(os, const_cast<Image*>(&ima)->at(p)) << " ";
+            }
+            else
+            {
+              format(os, ima.at(p)) << " ";
+            }
           }
           os << std::endl;
         }
@@ -159,5 +251,24 @@ namespace mln
     {
       internal::imprint_with_border(exact(ima), exact(ima).domain(), os);
     }
-  } // namespace io
+
+    namespace experimental
+    {
+      template <typename InputImage>
+      void imprint(InputImage ima, std::ostream& os)
+      {
+        static_assert(mln::is_a<InputImage, mln::experimental::Image>());
+
+        internal::imprint(ima, ima.domain(), os);
+      }
+
+      template <typename InputImage>
+      void imprint_with_border(InputImage ima, std::ostream& os)
+      {
+        static_assert(mln::is_a<InputImage, mln::experimental::Image>());
+
+        internal::imprint_with_border(ima, ima.domain(), os);
+      }
+    } // namespace experimental
+  }   // namespace io
 } // namespace mln
diff --git a/pylene/include/mln/io/saver.hpp b/pylene/include/mln/io/saver.hpp
index a00469ebf0f38742d6569a05b20e8cb67008427a..5b553748dc4a90d051ef6d96ed8b03f6c1c7a7c6 100644
--- a/pylene/include/mln/io/saver.hpp
+++ b/pylene/include/mln/io/saver.hpp
@@ -69,7 +69,7 @@ namespace mln
     void Saver<I>::save_experimental(I ima, PluginWriter* plugin, bool permissive) const
     {
       static_assert(is_a<I, mln::experimental::Image>());
-      
+
       if (not plugin->can_write(typeid(mln_value(I))))
       {
         std::string msg = "The plugin does not support writing " + internal::demangle(typeid(mln_value(I)).name());
@@ -89,7 +89,7 @@ namespace mln
 
       // FIXME: Oh my! So dirty!
       // To be fixed with new io facility based on type-erased image for python
-      if constexpr (is_a<I, experimental::Image>{})
+      if constexpr (mln::is_a<I, mln::experimental::Image>{})
       {
         for (auto v : const_cast<I*>(&ima)->new_values())
         {
diff --git a/pylene/include/mln/morpho/structural/gradient.hpp b/pylene/include/mln/morpho/structural/gradient.hpp
index 9b9313d3b5d6a0ec5d64c2eea911e17dd0a6d075..1de42be5f3ac1316eb8c2dee6ae3b14f1406312d 100644
--- a/pylene/include/mln/morpho/structural/gradient.hpp
+++ b/pylene/include/mln/morpho/structural/gradient.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
 #include <mln/core/algorithm/transform.hpp>
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/image/image.hpp>
-#include <mln/core/math_ops.hpp>
 #include <mln/core/trace.hpp>
 #include <mln/morpho/se/se.hpp>
 #include <mln/morpho/structural/dilate.hpp>
diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt
index 5cb7e60e49ef14950df30cf39478cc18c2b7e916..90b97882b04708edaf1c8c4b5cd09da68fa6dfbb 100644
--- a/tests/core/CMakeLists.txt
+++ b/tests/core/CMakeLists.txt
@@ -18,10 +18,18 @@ add_core_test(${test_prefix}range_zip           range/zip.cpp)
 
 # test Views
 add_core_test(${test_prefix}view_adaptor     image/view/adaptor.cpp)
-add_core_test(${test_prefix}view_transform   image/view/transform.cpp)
-add_core_test(${test_prefix}view_filter      image/view/filter.cpp)
+add_core_test(${test_prefix}view_cast        image/view/cast.cpp)
+add_core_test(${test_prefix}view_channel     image/view/channel.cpp)
 add_core_test(${test_prefix}view_clip        image/view/clip.cpp)
+add_core_test(${test_prefix}view_filter      image/view/filter.cpp)
 add_core_test(${test_prefix}view_mask        image/view/mask.cpp)
+add_core_test(${test_prefix}view_maths_comps image/view/maths_components.cpp)
+add_core_test(${test_prefix}view_maths_reds  image/view/maths_reductions.cpp)
+add_core_test(${test_prefix}view_maths_norms image/view/maths_norms.cpp)
+add_core_test(${test_prefix}view_maths_dists image/view/maths_distances.cpp)
+add_core_test(${test_prefix}view_ops         image/view/operators.cpp)
+add_core_test(${test_prefix}view_rgb         image/view/rgb.cpp)
+add_core_test(${test_prefix}view_transform   image/view/transform.cpp)
 add_core_test(${test_prefix}view_zip         image/view/zip.cpp)
 
 
@@ -33,7 +41,6 @@ add_core_test(${test_prefix}transformed_image       image/morphers/transformed_i
 add_core_test(${test_prefix}filtered_image          image/morphers/filtered_image.cpp)
 add_core_test(${test_prefix}zip_image               image/zip_image.cpp)
 add_core_test(${test_prefix}sub_image               image/sub_image.cpp)
-add_core_test(${test_prefix}image_ops               image/image_ops.cpp)
 add_core_test(${test_prefix}image_routines          image/image_routines.cpp)
 add_core_test(${test_prefix}image_ndimage           image/ndimage.cpp)
 add_core_test(${test_prefix}image_private_ndpixel   image/private/ndimage_pixel.cpp)
diff --git a/tests/core/image/view/cast.cpp b/tests/core/image/view/cast.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da3ed71892d1fea68e56534763cde7466bea94c3
--- /dev/null
+++ b/tests/core/image/view/cast.cpp
@@ -0,0 +1,30 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/cast.hpp>
+#include <mln/core/image/view/operators.hpp>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, cast)
+{
+  using namespace mln;
+  using namespace view::ops;
+
+  image2d<double> ima = {{0.1, 1.2, 2.3, 3.4, 4.5}, //
+                         {5.6, 6.7, 4.8, 8.7, 9.6}, //
+                         {10.5, 11.4, 12.3, 13.2, 14.1}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  auto casted_ima = view::cast<int>(ima);
+  ASSERT_TRUE(all_of(casted_ima == ref));
+
+  image2d<double> ima2 = {{0.1, 0.2, 0.3, 0.4, 0.5}, //
+                          {0.6, 0.7, 0.8, 0.7, 0.6}, //
+                          {0.5, 0.4, 0.3, 0.2, 0.1}};
+
+  ASSERT_TRUE(all_of(ima == casted_ima + ima2));
+}
diff --git a/tests/core/image/view/channel.cpp b/tests/core/image/view/channel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..38a78c0f3d143824555f01b86fd71dfec568ab78
--- /dev/null
+++ b/tests/core/image/view/channel.cpp
@@ -0,0 +1,59 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/algorithm/fill.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/channel.hpp>
+#include <mln/core/image/view/operators.hpp>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, channel_static)
+{
+  using namespace mln::view::ops;
+
+  mln::image2d<mln::rgb8> ima = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima_ch1 = mln::view::channel<1>(ima);
+
+  mln::image2d<unsigned> ref = {{2, 5, 8}, //
+                                {2, 5, 8}, //
+                                {2, 5, 8}};
+
+  ASSERT_TRUE(mln::all_of(ref == ima_ch1));
+
+  mln::fill(ima_ch1, 42);
+
+  mln::image2d<mln::rgb8> ref2 = {{{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}, //
+                                  {{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}, //
+                                  {{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}};
+
+  ASSERT_TRUE(mln::all_of(ref2 == ima));
+}
+
+TEST(View, channel_dynamic)
+{
+  using namespace mln::view::ops;
+
+  mln::image2d<mln::rgb8> ima = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima_ch2 = mln::view::channel(ima, 2);
+
+  mln::image2d<unsigned> ref = {{3, 6, 9}, //
+                                {3, 6, 9}, //
+                                {3, 6, 9}};
+
+  ASSERT_TRUE(mln::all_of(ref == ima_ch2));
+
+  mln::fill(ima_ch2, 42);
+
+  mln::image2d<mln::rgb8> ref2 = {{{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}, //
+                                  {{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}, //
+                                  {{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}};
+
+  ASSERT_TRUE(mln::all_of(ref2 == ima));
+}
diff --git a/tests/core/image/view/clip.cpp b/tests/core/image/view/clip.cpp
index a42d82f5bb2e3cc72dae4a92a771e1cdbc073a91..199377dc4296c890a1b58808d4a2326005b797e9 100644
--- a/tests/core/image/view/clip.cpp
+++ b/tests/core/image/view/clip.cpp
@@ -7,8 +7,6 @@
 #include <mln/core/image/view/clip.hpp>
 #include <mln/core/image/view/operators.hpp>
 
-#include <fixtures/ImageCompare/image_compare.hpp>
-
 #include <gtest/gtest.h>
 
 
@@ -180,6 +178,27 @@ TEST(Core, Clip_wherex2_joints)
   ASSERT_TRUE(all_of(ima == ref));
 }
 
+TEST(Core, Clip_wherex2_disjoints)
+{
+  using namespace mln;
+  using namespace mln::view::ops;
+
+  image2d<int> ima(5, 5);
+
+  image2d<int> ref = {
+      {0, 1, 2, 3, 4},      //
+      {5, 6, 7, 8, 9},      //
+      {10, 42, 42, 42, 42}, //
+      {42, 42, 42, 42, 42}, //
+      {20, 21, 22, 23, 24}  //
+  };
+
+  mln::iota(ima, 0);
+  auto clipped_ima = view::clip(ima, experimental::where(ima > 10));
+  mln::fill(view::clip(clipped_ima, experimental::where(clipped_ima < 20)), 42);
+  ASSERT_TRUE(all_of(ima == ref));
+}
+
 TEST(Core, Clip_where_and)
 {
   using namespace mln;
diff --git a/tests/core/image/view/filter.cpp b/tests/core/image/view/filter.cpp
index 9b84c23f8f80256f1b93688265ec36570871ec68..4b5f5bee8b2ef4dfb456f88a899601886f28bca2 100644
--- a/tests/core/image/view/filter.cpp
+++ b/tests/core/image/view/filter.cpp
@@ -1,18 +1,13 @@
-#include <mln/core/image/image2d.hpp>
-#include <mln/core/image/view/filter.hpp>
-#include <mln/core/image/view/operators.hpp>
-
 #include <mln/core/algorithm/all_of.hpp>
 #include <mln/core/algorithm/fill.hpp>
 #include <mln/core/algorithm/iota.hpp>
-
+#include <mln/core/concept/new/concepts.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/filter.hpp>
+#include <mln/core/image/view/operators.hpp>
 #include <mln/core/rangev3/foreach.hpp>
 #include <mln/core/rangev3/view/zip.hpp>
 
-#include <mln/core/concept/new/concepts.hpp>
-
-#include <fixtures/ImageCompare/image_compare.hpp>
-
 #include <gtest/gtest.h>
 
 
diff --git a/tests/core/image/view/mask.cpp b/tests/core/image/view/mask.cpp
index 612060e7414b072b2630b83d733c9363d1a7898d..676a9adb90bf1230a66ec1ce457bb709a6ba2048 100644
--- a/tests/core/image/view/mask.cpp
+++ b/tests/core/image/view/mask.cpp
@@ -1,14 +1,11 @@
-#include <mln/core/image/view/mask.hpp>
-
 #include <mln/core/algorithm/all_of.hpp>
 #include <mln/core/algorithm/fill.hpp>
 #include <mln/core/concept/new/archetype/image.hpp>
 #include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/mask.hpp>
 #include <mln/core/image/view/operators.hpp>
 #include <mln/core/rangev3/foreach.hpp>
 
-#include <fixtures/ImageCompare/image_compare.hpp>
-
 #include <gtest/gtest.h>
 
 
diff --git a/tests/core/image/view/maths_components.cpp b/tests/core/image/view/maths_components.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..916999189a7d23db893558d6f7bf5f69f62b1e99
--- /dev/null
+++ b/tests/core/image/view/maths_components.cpp
@@ -0,0 +1,143 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/maths.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/transform.hpp>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, maths_abs)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == abs(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ref2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  ASSERT_TRUE(all_of(ref2 == abs(ima2)));
+}
+
+TEST(View, maths_sqr)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  image2d<int> ref = {{0, 1, 4, 9, 16},     //
+                      {25, 36, 16, 64, 81}, //
+                      {100, 121, 144, 169, 196}};
+
+  ASSERT_TRUE(all_of(ref == sqr(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ref2 = {{{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}, //
+                             {{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}, //
+                             {{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}};
+
+  ASSERT_TRUE(all_of(ref2 == sqr(ima2)));
+}
+
+TEST(View, maths_pow)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  auto ima = view::transform(ref, [](auto x) { return x * x * x; });
+
+  ASSERT_TRUE(all_of(pow(ref, 3) == ima));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb16> ref2 = {{{1, 8, 27}, {64, 125, 216}, {343, 512, 729}}, //
+                              {{1, 8, 27}, {64, 125, 216}, {343, 512, 729}}, //
+                              {{1, 8, 27}, {64, 125, 216}, {343, 512, 729}}};
+
+  ASSERT_TRUE(all_of(ref2 == pow(ima2, 3)));
+}
+
+TEST(View, maths_sqrt)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, 1, 4, 9, 16},     //
+                      {25, 36, 16, 64, 81}, //
+                      {100, 121, 144, 169, 196}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == sqrt(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}, //
+                             {{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}, //
+                             {{1, 4, 9}, {16, 25, 36}, {49, 64, 81}}};
+
+  image2d<mln::rgb8> ref2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  ASSERT_TRUE(all_of(ref2 == sqrt(ima2)));
+}
+
+TEST(View, maths_cbrt)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<double> ref = {{0, 1, 2, 3, 4}, //
+                         {5, 6, 4, 8, 9}, //
+                         {10, 11, 12, 13, 14}};
+
+  auto ima = pow(ref, 3);
+
+  ASSERT_TRUE(all_of(equalFP(ref, cbrt(ima), 10e-6)));
+
+
+  image2d<mln::rgb<double>> ref2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                    {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                    {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima2 = pow(ref2, 3);
+
+  ASSERT_TRUE(all_of(equalFP(ref2, cbrt(ima2), 10e-6)));
+}
diff --git a/tests/core/image/view/maths_distances.cpp b/tests/core/image/view/maths_distances.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a3278e0830c8c69061d3ce26fc76a1f52336d59
--- /dev/null
+++ b/tests/core/image/view/maths_distances.cpp
@@ -0,0 +1,227 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/maths.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/transform.hpp>
+
+#include <cmath>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, maths_l0dist)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 2, 2, 12, 12},    //
+                      {30, 30, 20, 56, 90}, //
+                      {90, 132, 132, 182, 182}};
+
+  ASSERT_TRUE(all_of(ref == l0dist(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}, //
+                             {{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}, //
+                             {{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}};
+
+  image2d<mln::rgb8> ima4 = {{{3, 6, 6}, {2, 4, 4}, {1, 2, 2}}, //
+                             {{3, 6, 6}, {2, 4, 4}, {1, 2, 2}}, //
+                             {{3, 6, 6}, {2, 4, 4}, {1, 2, 2}}};
+
+  image2d<int> ref2 = {{2, 0, 2}, //
+                       {2, 0, 2}, //
+                       {2, 0, 2}};
+
+  ASSERT_TRUE(all_of(ref2 == l0dist(ima3, ima4)));
+}
+
+TEST(View, maths_l1dist)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 2, 2, 12, 12},    //
+                      {30, 30, 20, 56, 90}, //
+                      {90, 132, 132, 182, 182}};
+
+  ASSERT_TRUE(all_of(ref == l1dist(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima4 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<int> ref2 = {{18, 0, 18}, //
+                       {18, 0, 18}, //
+                       {18, 0, 18}};
+
+  ASSERT_TRUE(all_of(ref2 == l1dist(ima3, ima4)));
+}
+
+TEST(View, maths_l2dist)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<double> ref = {{0, 2, 2, 12, 12},    //
+                         {30, 30, 20, 56, 90}, //
+                         {90, 132, 132, 182, 182}};
+
+  ASSERT_TRUE(all_of(ref == l2dist(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima4 = {{{5, 6, 5}, {4, 5, 6}, {5, 4, 5}}, //
+                             {{5, 6, 5}, {4, 5, 6}, {5, 4, 5}}, //
+                             {{5, 6, 5}, {4, 5, 6}, {5, 4, 5}}};
+
+  image2d<double> ref2 = {{6, 0, 6}, //
+                          {6, 0, 6}, //
+                          {6, 0, 6}};
+
+  ASSERT_TRUE(all_of(ref2 == l2dist(ima3, ima4)));
+}
+
+TEST(View, maths_lpdist)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 2, 2, 12, 12},    //
+                      {30, 30, 20, 56, 90}, //
+                      {90, 132, 132, 182, 182}};
+
+  ASSERT_TRUE(all_of(ref == lpdist<4>(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima4 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<double> ref2 = {{std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.), 0, std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.)}, //
+                          {std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.), 0, std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.)}, //
+                          {std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.), 0, std::pow(6 * 6 * 6 * 6 * 3, 1. / 4.)}};
+
+  ASSERT_TRUE(all_of(ref2 == lpdist<4>(ima3, ima4)));
+}
+
+TEST(View, maths_l2dist_sqr)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 4, 4, 144, 144},         //
+                      {900, 900, 400, 3136, 8100}, //
+                      {8100, 17424, 17424, 33124, 33124}};
+
+  ASSERT_TRUE(all_of(ref == l2dist_sqr(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima4 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<int> ref2 = {{108, 0, 108}, //
+                       {108, 0, 108}, //
+                       {108, 0, 108}};
+
+  ASSERT_TRUE(all_of(ref2 == l2dist_sqr(ima3, ima4)));
+}
+
+TEST(View, maths_linfdist)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima1 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ima2 = {{0, -1, 2, -3, 4},  //
+                       {-5, 6, -4, 8, -9}, //
+                       {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 4, 4, 144, 144},         //
+                      {900, 900, 400, 3136, 8100}, //
+                      {8100, 17424, 17424, 33124, 33124}};
+
+  ASSERT_TRUE(all_of(ref == linfdist(ima1, sqr(ima2))));
+
+
+  image2d<mln::rgb8> ima3 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima4 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<int> ref2 = {{6, 0, 6}, //
+                       {6, 0, 6}, //
+                       {6, 0, 6}};
+
+  ASSERT_TRUE(all_of(ref2 == linfdist(ima3, ima4)));
+}
diff --git a/tests/core/image/view/maths_norms.cpp b/tests/core/image/view/maths_norms.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7273f1a3cfa43055d626ae84fd0c7e1ef3f21fc
--- /dev/null
+++ b/tests/core/image/view/maths_norms.cpp
@@ -0,0 +1,179 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/maths.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/transform.hpp>
+
+#include <cmath>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, maths_l0norm)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == l0norm(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{1, 4, 7}, //
+                       {1, 4, 7}, //
+                       {1, 4, 7}};
+
+  ASSERT_TRUE(all_of(ref2 == l0norm(ima2)));
+}
+
+TEST(View, maths_l1norm)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == l1norm(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{6, 15, 24}, //
+                       {6, 15, 24}, //
+                       {6, 15, 24}};
+
+  ASSERT_TRUE(all_of(ref2 == l1norm(ima2)));
+}
+
+TEST(View, maths_l2norm)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == l2norm(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}, //
+                             {{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}, //
+                             {{1, 2, 2}, {2, 4, 4}, {3, 6, 6}}};
+
+  image2d<int> ref2 = {{3, 6, 9}, //
+                       {3, 6, 9}, //
+                       {3, 6, 9}};
+
+  ASSERT_TRUE(all_of(ref2 == l2norm(ima2)));
+}
+
+TEST(View, maths_lpnorm)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == lpnorm<4>(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<double> ref2 = {{std::pow(98, 1. / 4.), std::pow(2177, 1. / 4.), std::pow(13058, 1. / 4.)}, //
+                          {std::pow(98, 1. / 4.), std::pow(2177, 1. / 4.), std::pow(13058, 1. / 4.)}, //
+                          {std::pow(98, 1. / 4.), std::pow(2177, 1. / 4.), std::pow(13058, 1. / 4.)}};
+
+  ASSERT_TRUE(all_of(ref2 == lpnorm<4>(ima2)));
+}
+
+TEST(View, maths_l2norm_sqr)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 4, 9, 16},     //
+                      {25, 36, 16, 64, 81}, //
+                      {100, 121, 144, 169, 196}};
+
+  ASSERT_TRUE(all_of(ref == l2norm_sqr(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{14, 77, 194}, //
+                       {14, 77, 194}, //
+                       {14, 77, 194}};
+
+  ASSERT_TRUE(all_of(ref2 == l2norm_sqr(ima2)));
+}
+
+TEST(View, maths_linfnorm)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, 1, 2, 3, 4}, //
+                      {5, 6, 4, 8, 9}, //
+                      {10, 11, 12, 13, 14}};
+
+  ASSERT_TRUE(all_of(ref == linfnorm(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{3, 6, 9}, //
+                       {3, 6, 9}, //
+                       {3, 6, 9}};
+
+  ASSERT_TRUE(all_of(ref2 == linfnorm(ima2)));
+}
diff --git a/tests/core/image/view/maths_reductions.cpp b/tests/core/image/view/maths_reductions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..27f325ebf629087265362322a62c197218c817fe
--- /dev/null
+++ b/tests/core/image/view/maths_reductions.cpp
@@ -0,0 +1,163 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/maths.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/transform.hpp>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, maths_sum)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  ASSERT_TRUE(all_of(ref == sum(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{6, 15, 24}, //
+                       {6, 15, 24}, //
+                       {6, 15, 24}};
+
+  ASSERT_TRUE(all_of(ref2 == sum(ima2)));
+}
+
+TEST(View, maths_prod)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  ASSERT_TRUE(all_of(ref == prod(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{6, 120, 504}, //
+                       {6, 120, 504}, //
+                       {6, 120, 504}};
+
+  ASSERT_TRUE(all_of(ref2 == prod(ima2)));
+}
+
+TEST(View, maths_cross)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<mln::rgb8> ima1 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima2 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<mln::signed_rgb8> ref = {{{-6, 12, -6}, {0, 0, 0}, {6, -12, 6}}, //
+                                   {{-6, 12, -6}, {0, 0, 0}, {6, -12, 6}}, //
+                                   {{-6, 12, -6}, {0, 0, 0}, {6, -12, 6}}};
+
+  ASSERT_TRUE(all_of(ref == cross(ima1, ima2)));
+}
+
+TEST(View, maths_dot)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<mln::rgb8> ima1 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<mln::rgb8> ima2 = {{{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}, //
+                             {{7, 8, 9}, {4, 5, 6}, {1, 2, 3}}};
+
+  image2d<int> ref = {{50, 77, 50}, //
+                      {50, 77, 50}, //
+                      {50, 77, 50}};
+
+  ASSERT_TRUE(all_of(ref == dot(ima1, ima2)));
+}
+
+TEST(View, maths_min)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  ASSERT_TRUE(all_of(ref == min(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{1, 4, 7}, //
+                       {1, 4, 7}, //
+                       {1, 4, 7}};
+
+  ASSERT_TRUE(all_of(ref2 == min(ima2)));
+}
+
+TEST(View, maths_max)
+{
+  using namespace mln;
+  using namespace view::ops;
+  using namespace view::maths;
+
+  image2d<int> ima = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  image2d<int> ref = {{0, -1, 2, -3, 4},  //
+                      {-5, 6, -4, 8, -9}, //
+                      {10, -11, 12, -13, 14}};
+
+  ASSERT_TRUE(all_of(ref == max(ima)));
+
+
+  image2d<mln::rgb8> ima2 = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  image2d<int> ref2 = {{3, 6, 9}, //
+                       {3, 6, 9}, //
+                       {3, 6, 9}};
+
+  ASSERT_TRUE(all_of(ref2 == max(ima2)));
+}
diff --git a/tests/core/image/image_ops.cpp b/tests/core/image/view/operators.cpp
similarity index 94%
rename from tests/core/image/image_ops.cpp
rename to tests/core/image/view/operators.cpp
index b9c5c1ca3cbbf000c67213b68951b176dd783297..3717df3fe0b5ef7a073c57a57b32ba279ca82f47 100644
--- a/tests/core/image/image_ops.cpp
+++ b/tests/core/image/view/operators.cpp
@@ -7,7 +7,6 @@
 #include <mln/core/image/image_ops.hpp>
 #include <mln/core/image/view/operators.hpp>
 #include <mln/core/image/view/transform.hpp>
-
 #include <mln/io/imprint.hpp>
 
 #include <boost/tuple/tuple_io.hpp>
@@ -119,6 +118,18 @@ TEST(Core, BinaryOperator_SameTypes)
   ASSERT_TRUE(mln::all_of(g3 == ref));
 }
 
+TEST(Core, BinaryOperator_EqualFP)
+{
+  using namespace mln;
+  using namespace mln::view::ops;
+
+  image2d<double> ima = {{1.015, 2.015, 3.015}, {4.015, 5.015, 6.015}};
+  image2d<double> ref = {{1.016, 2.016, 3.016}, {4.016, 5.016, 6.016}};
+
+  ASSERT_TRUE(mln::all_of(equalFP(ima, ref, 10e-3)));
+  ASSERT_TRUE(mln::all_of(not equalFP(ima, ref, 10e-5)));
+}
+
 TEST(Core, BinaryOperators_MixedTypes)
 {
   using namespace mln;
diff --git a/tests/core/image/view/rgb.cpp b/tests/core/image/view/rgb.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c16f54fad677b377e0c1c7040dc7efdbf1a3cf54
--- /dev/null
+++ b/tests/core/image/view/rgb.cpp
@@ -0,0 +1,84 @@
+#include <mln/core/algorithm/all_of.hpp>
+#include <mln/core/algorithm/fill.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/image/view/rgb.hpp>
+
+#include <gtest/gtest.h>
+
+
+TEST(View, rgb_red)
+{
+  using namespace mln::view::ops;
+
+  mln::image2d<mln::rgb8> ima = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima_red = mln::view::red(ima);
+
+  mln::image2d<unsigned> ref = {{1, 4, 7}, //
+                                {1, 4, 7}, //
+                                {1, 4, 7}};
+
+  ASSERT_TRUE(mln::all_of(ref == ima_red));
+
+  mln::fill(ima_red, 42);
+
+  mln::image2d<mln::rgb8> ref2 = {{{42, 2, 3}, {42, 5, 6}, {42, 8, 9}}, //
+                                  {{42, 2, 3}, {42, 5, 6}, {42, 8, 9}}, //
+                                  {{42, 2, 3}, {42, 5, 6}, {42, 8, 9}}};
+
+  ASSERT_TRUE(mln::all_of(ref2 == ima));
+}
+
+TEST(View, rgb_green)
+{
+  using namespace mln::view::ops;
+
+  mln::image2d<mln::rgb8> ima = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima_green = mln::view::green(ima);
+
+  mln::image2d<unsigned> ref = {{2, 5, 8}, //
+                                {2, 5, 8}, //
+                                {2, 5, 8}};
+
+  ASSERT_TRUE(mln::all_of(ref == ima_green));
+
+  mln::fill(ima_green, 42);
+
+  mln::image2d<mln::rgb8> ref2 = {{{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}, //
+                                  {{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}, //
+                                  {{1, 42, 3}, {4, 42, 6}, {7, 42, 9}}};
+
+  ASSERT_TRUE(mln::all_of(ref2 == ima));
+}
+
+TEST(View, rgb_blue)
+{
+  using namespace mln::view::ops;
+
+  mln::image2d<mln::rgb8> ima = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //
+                                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}};
+
+  auto ima_blue = mln::view::blue(ima);
+
+  mln::image2d<unsigned> ref = {{3, 6, 9}, //
+                                {3, 6, 9}, //
+                                {3, 6, 9}};
+
+  ASSERT_TRUE(mln::all_of(ref == ima_blue));
+
+  mln::fill(ima_blue, 42);
+
+  mln::image2d<mln::rgb8> ref2 = {{{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}, //
+                                  {{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}, //
+                                  {{1, 2, 42}, {4, 5, 42}, {7, 8, 42}}};
+
+  ASSERT_TRUE(mln::all_of(ref2 == ima));
+}
diff --git a/tests/core/image/view/transform.cpp b/tests/core/image/view/transform.cpp
index 5e66c3fb357ab4eabcb333d26b24b6397d84b6d9..7ae1b021276b5a775e86a1d27a037fc4ec68b44f 100644
--- a/tests/core/image/view/transform.cpp
+++ b/tests/core/image/view/transform.cpp
@@ -12,8 +12,6 @@
 #include <mln/core/rangev3/foreach.hpp>
 #include <mln/core/rangev3/view/zip.hpp>
 
-#include <fixtures/ImageCompare/image_compare.hpp>
-
 #include <vector>
 
 #include <gtest/gtest.h>
diff --git a/tests/core/vec/vec_math_ops.cpp b/tests/core/vec/vec_math_ops.cpp
index feee37f58899bf792a6a45f157c2473067de4f2f..d92050fad46d5190561af9f1ac33bbd4132f5f4c 100644
--- a/tests/core/vec/vec_math_ops.cpp
+++ b/tests/core/vec/vec_math_ops.cpp
@@ -1,5 +1,5 @@
+#include <mln/core/functional_ops.hpp>
 #include <mln/core/vec.hpp>
-#include <mln/core/vec/vec_math_ops.hpp>
 
 #include <gtest/gtest.h>