Commit dbb3e872 authored by Edwin Carlinet's avatar Edwin Carlinet

Add ruminations for doc.

parent 7fe0be81
......@@ -13,7 +13,8 @@ Contents:
Introduction <intro>
tutorial
Réference <reference>
Reference <reference>
Rumination <ruminations>
Indices and tables
......@@ -24,3 +25,4 @@ Indices and tables
* :ref:`search`
Ruminations
###########
.. toctree::
ruminations/neighborhood
ruminations/iterators
Cleaning iterators
##################
Histoire
========
Les traitements de base en TI s'expriment généralement sous deux formes:
1. Pour chaque pixels, effectuer une action. Par exemple, un changement de
contraste `g` s'écrira::
void change_contrast(int* buffer, std::size_t width, std::size_t height, std::size_t stride)
{
auto g = [](auto x) { return x + 1; };
for (std::size_t y = 0; y < height; ++y)
{
auto lineptr = buffer + y * stride;
for (std::size_t x = 0; x < width; ++x)
lineptr[x] = g(lineptr[x]);
}
}
2. Pour chaque pixels, pour chaque voisin, effectuer une action. Par example,
un filtre moyenneur 3x3 s'écrira (on omet les problèmes d'arrondis et de bordure)::
void boxcar(const int* __restrict inbuffer, int* __restrict outbuffer, std::size_t width, std::size_t height, std::size_t stride)
{
for (std::size_t y = 0; y < height; ++y)
{
auto lineptr = inbuffer + y * stride;
auto outptr = outbuffer + y * stride;
for (std::size_t x = 0; x < width; ++x)
{
int sum = 0;
auto ptr = lineptr + x;
for (int dy = -1; dy <= 1; ++dy)
for (int dx = -1; dx <= 1; ++dx)
sum += ptr[dy * stride + dx];
outptr[x] = sum / 9;
}
}
}
Pour le premier code, tous les compilateurs (GCC, CLANG, MSVC, ICC) sont
capables de générer du code optimisé (vectorisé). Pour le second, il faut aider
le compilateur pour lui dire que les pointeurs ne sont pas aliasés. Mais même
alors, MSVC échoue à vectoriser. Si on déroule la boucle manuellement, alors il
y arrive.::
void boxcar_unrolled(const int* __restrict inbuffer, int* __restrict outbuffer, std::size_t width, std::size_t height, std::size_t stride)
{
for (std::size_t y = 0; y < height; ++y)
{
auto prevlineptr = inbuffer + (y-1) * stride;
auto nextlineptr = inbuffer + (y+1) * stride;
auto lineptr = inbuffer + y * stride;
auto outptr = outbuffer + y * stride;
for (std::size_t x = 0; x < width; ++x)
{
auto sum1 = prevlineptr[x-1] + prevlineptr[x] + prevlineptr[x+1];
auto sum2 = lineptr[x-1] + lineptr[x] + lineptr[x+1];
auto sum3 = nextlineptr[x-1] + nextlineptr[x] + nextlineptr[x+1];
outptr[x] = (sum1 + sum2 + sum3) / 9;
}
}
}
On se rend compte alors que cette écriture n'est pas ce que l'on veut:
* elle est non-générique (image 3D ?)
* elle est difficile à écrire et sujette à l'erreur
En revanche:
* elle permet d'obtenir un code rapide
Les itérateurs ont été introduits pour permettre une abstraction. Le code
suivant est le même pour le 2D et le 3D::
void change_constrat(image2d& f)
{
for (auto p : f.domain())
f(p) = f(p) + 1;
}
void boxcar(const image2d& f, image2d& g)
{
for (auto p : f.domain())
{
int sum = 0;
for (auto n : win3x3(p))
sum += f(n);
g(p) = sum / 9;
}
}
Le code est élégant mais l'asm produit horrible.
Design du code actuel et problèmes
==================================
L'itérateur sur pixel
---------------------
Les positions sont recalculées à chaque fois. `f(p) === f.buffer[p.y * f.stride + p.x]`.
* lors d'une itération externe, les pointeurs/indexes sont de simplement décalés
* lors d'une itération interne, les pointeurs/indexes sont des offsets par
rapport à la position courante.
(noter que pour le compilateur `f.buffer[?]` peut aliaser `f.stride` s'ils sont
du même type).
Les itérateurs sur pixels ont été introduits dans cet objectif. Ils fournissent
à la fois la valeur et la position dans l'image (`px.point()`, `px.val()`)::
void change_constrat(image2d& f)
{
for (auto px : f.pixels())
px.val() = px.val() + 1;
}
void boxcar(const image2d& f, image2d& g)
{
for (auto [pxin, pxout] : ziprange(f.pixels(), g.pixels()))
{
int sum = 0;
for (auto nx : win3x3(pxin))
sum += nx.val();
nxout = sum / 9;
}
}
Çà aurait pu être beau, mais attention à ce que fait `win3x3(pxin)`, il utilise
la position de `pxin` et créé un tableau d'offset. Si le compilateur n'est pas
capable de produire une "Loop Invariant Code Motion" de cette opération, il n'y
a aucun bénéfice. Aucun des compilateurs n'a été capable de le produire à
l'époque. On déplace donc la création des itérateurs à l'extérieur de la
boucle avec des itérateurs qui se *bind* sur d'autres itérateurs.::
void boxcar(const image2d& f, image2d& g)
{
auto&& inrng = f.pixels();
auto&& outrng = g.pixels();
auto pxin = std::begin(inrng);
auto pxend = std::end(inrng);
auto pxout = std::begin(outrng);
auto nbh = win3x3(pxin); // <- bounded to pxin
for(; pxin != pxend; ++pxin, ++pxend)
{
int sum = 0;
for (auto&& nx : nbh)
sum += nx.val();
pxout->val() = sum / 9;
}
}
Et le code devient lourd. C'est dû au fait que les itérateurs en C++ sont
flexibles mais lourd dans l'écriture. Les itérateurs ne connaissent ni leur
début, ni leur fin. On a donc pallié à cela::
class iterator
{
void init();
void next();
bool is_finished() const;
T dereference() const;
}
void boxcar(const image2d& f, image2d& g)
{
mln_iter(pxin, f.pixels());
mln_iter(pxout, g.pixels());
mln_iter(nx, win3x3(pxin));
mln_forall(pxin, pxout)
{
int sum = 0;
mln_forall(nx)
sum += nx->val();
pxout->val() = sum / 9;
}
}
Mais ce n'est pas rose:
* Incompatibilité avec la STL
* Des macros qui remplacent les "for"
* Et toujours des perfs manquantes... (whattt ??? )
Cas de la simple boucle
-----------------------
Pour permettre d'écrire ce code, des boucles imbriquées doivent être réécrites
en simple boucle. Un objet `nditerator` permet de faire un *flatten*.::
void change_constrat(image2d& f)
{
mln_foreach(auto px, f.pixels())
px.val() = px.val() + 1;
}
Avec `mln_foreach(X, RANGE)`. ::
{
auto it = RANGE.iter()
for (it.init(); !it.is_finished(); it.next())
{
X = it.dereference();
...
}
}
Cependant `it.next()` n'est pas trivial. C'est à cause de ça que les compilateurs ont du
mal à vectoriser. Typiquement::
void iterator2d::next()
{
x++;
if (x == xend)
{
y++;
x = 0;
}
}
C'est pourquoi une première optimisation a été que `mln_forall` et `mln_foreach`
soient deux boucles par défaut. `mln_foreach(X, RANGE)`::
{
auto it = RANGE.iter()
for (it.outer_init(); !it.outer_is_finished(); it.outer_next())
for (it.inner_init(); !it.inner_is_finished(); it.inner_next())
{
X = it.dereference();
...
}
}
Grâce à cela, on arrive à récupérer une performance proche du C mais le coût
est élevé:
* Non conforme au standard
* Perte de la compatibilité avec la STL
* Beaucoup de code non-visible du point de vue *end-user*
* Beaucoup de code tricky => trés mauvais pour la maintenance !!
Cas de la boucle imbriquée
--------------------------
On est loin de la performance du C. Le fait que le voisinage tienne une
back-référence vers le pixel itérateur empêche la SROA sur certains
compilateurs.
De plus, les problèmes d'aliasing d'images font que le code n'est pas
vectorisé. Le `__restrict` sur les pointeurs dans les `struct` ne sont pas
toujours compris.
Solution possible
=================
Les co-routines permettent de créer des générateurs n-dimentionnel à moindre
coût avec une SROA assuré. On espère:
* Un gain en perf
* Une réécriture conforme aux concepts du standard
* Du code grandement simplifié
References
==========
[1] Segmented Iterators and Hierarchical Algorithms, Austern 1998
(http://lafstern.org/matt/segmented.pdf)
[2]
About Neighborhood / Structuring Element / Windows
##################################################
Concept:
The three objects represents the same idea: given a location we can access near-by location.
.. topic:: Some mathematical properties:
* Flat or Non-flat (weighted):
* Flat: erosion, dilation...
* Non-flat (weighted): e.g. convolution, non-flat erosion
* Separable. Square -> 2 lines, Diamond
* Adaptative (they depend on the localisation)
We would like to iterate over:
* a *point* or a *point iterator*: it gives a *point range*::
auto p = point or point iterator;
for (auto q : nbh(p))
(void) q;
* a *pixel* or a *pixel iterator*: it gives a *pixel range*::
auto px = pixel or pixel iterator;
for (auto qx : nbh(px))
(void) qx;
Optimization issue:
For pixel iterator, we precompute offsets. It has to be declared outside the
outerloop::
auto px = mln_pixter(f);
auto neighbors = nbh(px); // Precompute offsets & binds to px
mln_forall(px) // px advances
for (auto nx : neighbors)
(void) nx;
or::
auto px = mln_pixter(f);
auto neighbors = nbh(px); // Precompute offsets & binds to px
mln_forall(px) // px advances
mln_forall(nx)
(void) nx;
or::
auto neighbors = nbh.with(f); // Return a new neighborhood with offsets computed
for (auto px : f.pixels())
for (auto nx : neighbors(px))
.. note:: With have to choose a syntax.
Algorithms that run with neighborhood
=====================================
* Erosion / Dilation: with an implementation specific for lines, and for
separable SE and incremental SE
* Hit or miss: Use two SEs with empty intersection
* Mean filter: Implementation with integral images or incremental updates
* Median / Rank filter: Implementation with incremental updates
* Convolution: can be separable
.. uml::
@startuml
abstract class WeightedNeighborhood {
}
abstract class Neighborhood {
}
WeightedNeighborhood <|-- Neighborhood
Neighborhood <|-- c4
Neighborhood <|-- c8
Neighborhood <|-- hline2d
Neighborhood <|-- vline2d
Neighborhood <|-- periodic_line2d
Neighborhood <|-- rectangle2d
Neighborhood <|-- ball2d
Neighborhood <|-- diamond2d_1
Neighborhood <|-- diamond2d
Neighborhood <|-- neighb2d
@enduml
Non-flat Structuring Element / Weighted Neighborhood Concept
============================================================
.. cpp:concept:: template <typename N> WeightedNeighborhood
Concept for generic weighted neighborbood
* `SE`: A model of :cpp:concept:`WeightedNeighborhood`
* `se`: An *constant* instance of `SE`
.. rubric:: `Type definition`
:class: concept-typedefs
+---------------------+--------------------------------+----------------------------------------------+
|Type | Definition |Requirements |
+=====================+================================+==============================================+
|`SE::category` | |Convertible to `adaptative_neighborhood_tag` |
+---------------------+--------------------------------+----------------------------------------------+
|`SE::is_incremental` | Either `std::true_type` or | |
| | `std:false_type` | |
+---------------------+--------------------------------+----------------------------------------------+
|`SE::is_decomposable`| Either `std::true_type` or | |
| | `std:false_type` | |
+---------------------+--------------------------------+----------------------------------------------+
|`SE::is_flat` | Either `std::true_type` or | |
| | `std:false_type` | |
+---------------------+--------------------------------+----------------------------------------------+
.. rubric:: `Valid expression`
:class: concept-expr
+----------------------+------------------------+----------------------------------------------------------+
|Expression | Return Type | Sementics |
+======================+========================+==========================================================+
| ``se(p)`` | `Range<WeightedPoint>` |Return a :concept:`Forward Range` of points centered in |
| | |the point `p` or given by `p` if `p` is a point iterator. |
| | | |
+----------------------+------------------------+----------------------------------------------------------+
| ``se(px)`` | `Range<WeightedPixel>` |Return a :concept:`Forward Range` of pixels centered in |
| | |the pixel `px` or given by `px` if `px` is a pixel |
| | |iterator. |
+----------------------+------------------------+----------------------------------------------------------+
| ``se.has(p)`` | `boolean` |True if `p` is in the se. |
+----------------------+------------------------+----------------------------------------------------------+
| ``se.is_centered()`` | `boolean` |Equivalent to ``se.has(literal::zero)`` |
+----------------------+------------------------+----------------------------------------------------------+
Structuring Element / Neighborhhod Properties
=============================================
Flat vs Weighted Property
*************************
Structuring Element can be weighted or not (there are so called Flat Structuring
Element). In that case, all there weights are implicitely 1 and `is_flat`
property is set. *Flat Structuring Elements* (resp. *Non-Weighted Neighborhood*) are
called simply named *Structing Elements* (resp. *Neighborhood*).
.. cpp:concept:: template <typename N> Neighborhood
It refines the :cpp:concept:`WeightedNeighborhood` having the `is_flat` property. It corresponds to flat
structuring element where all weights are 1.
.. rubric:: `Type definition`
:class: concept-typedefs
+----------------------------+--------------------+
|Type | Definition |
| | |
+============================+====================+
|`SE::is_flat` | `std::true_type` |
+----------------------------+--------------------+
Categories
**********
+-------------------------------+-------------------------------------------------------------------+
|Category |Operations |
+===============================+===================================================================+
|Adaptative SE |Iterate over SE whose elements depend on the current point |
| | |
| | |
| | |
+----+--------------------------+-------------------------------------------------------------------+
| | Constant SE |The structuring element is said *constant* if it knows size of the |
| | | SE. |
| | | |
| | | |
| | | |
+----+----+---------------------+-------------------------------------------------------------------+
| | | Static SE |The structuring element is said *static* if it knows size of the SE|
| | | |at compile time. |
| | | | |
| | | | |
| | | | |
+----+----+----+----------------+-------------------------------------------------------------------+
The categories are used generally for optimisation purpose. When the numbers of
elements is known, we can precompute the memory offsets and store them. Knowing
size at compile time allows to avoid dynamic allocation for temporay data.
Constant Structuring Element Property
-------------------------------------
.. rubric:: `Type definition`
:class: concept-typedefs
+----------------------------+----------------------------------------------+
|Type |Requirements |
| | |
+============================+==============================================+
| `SE::category` |Convertible to `constant_neighborhood_tag` |
+----------------------------+----------------------------------------------+
.. rubric:: `Valid expression`
:class: concept-expr
+--------------------------+------------------------------+------------------------------------------------------+
|Expression | Return Type | Sementics |
+==========================+==============================+======================================================+
| ``se.size()`` | `unsigned` | The number of elements in the SE. |
+--------------------------+------------------------------+------------------------------------------------------+
| ``se.offsets()`` | `Range<Point>` | The elements of the SE. |
+--------------------------+------------------------------+------------------------------------------------------+
| ``se.before(p)`` | `Range<WeightedPoint>` | The elements of the SE before the anchor centered in |
| | | `p`. |
+--------------------------+------------------------------+------------------------------------------------------+
| ``se.before(px)`` | `Range<WeightedPixel>` | The elements of the SE before the anchor centered in |
| | | `px`. |
+--------------------------+------------------------------+------------------------------------------------------+
| ``se.after(p)`` | `Range<WeightedPoint>` | The elements of the SE after the anchor centered in |
| | | `p`. |
+--------------------------+------------------------------+------------------------------------------------------+
| ``se.after(px)`` | `Range<WeightedPixel>` | The elements of the SE after the anchor centered in |
| | | `px`. |
+--------------------------+------------------------------+------------------------------------------------------+
Static Structuring Element Property
-----------------------------------
.. rubric:: `Type definition`
:class: concept-typedefs
+--------------------------------+----------------------------------------------+
|Type |Requirements |
| | |
+================================+==============================================+
| `SE::category` |Convertible to `static_neighborhood_tag` |
+--------------------------------+----------------------------------------------+
.. rubric:: `Valid expression`
:class: concept-expr
+--------------------+-----------------+------------------------------------------------------------------+
|Expression | Return Type | Sementics |
+====================+=================+==================================================================+
| ``se.size()`` | `unsigned` |The number of elements in the SE as a *constexpr expression* |
+--------------------+-----------------+------------------------------------------------------------------+
Incremental Structuring Element Property
****************************************
A SE is said to be *incremental*, if it enables to give the points
that are added or removed to the range given a *basic deplacement* of
the point, e.g. for `point2d`, the basic deplacement is `(0,1)`. This
is usually used to compute attributes over a sliding SE in linear
time.
.. rubric:: `Type definition`
:class: concept-typedefs
+--------------------------------+--------------------+----------------------------------------------+
|Type | Definition |Requirements |
| | | |
+================================+====================+==============================================+
|`SE::is_incremental` | `std::true_type` | A |
+--------------------------------+--------------------+----------------------------------------------+
|`SE::dec_type` | | A model of :cpp:concept:`SE` |
+--------------------------------+--------------------+----------------------------------------------+
|`SE::inc_type` | | A model of |
| | | :cpp:concept:`WeightedNeighborhood` |
+--------------------------------+--------------------+----------------------------------------------+
.. rubric:: `Valid expression`
:class: concept-expr
+----------------------+-------------------------+---------------------------------------------------------+
|Expression | Return Type | Sementics |
+======================+=========================+=========================================================+
| ``se.inc()`` | `SE::inc_type` |A SE equivalent to :math:`\Delta\mathcal{B}^+(p) = |
| | |\mathcal{B}(p) \setminus (\mathcal{B}(p) \cap |
| | |\mathcal{B}(\mathrm{prev}))` |
+----------------------+-------------------------+---------------------------------------------------------+
| ``se.dec()`` | `SE::dec_type` |A SE `s` equivalent to :math:`\Delta\mathcal{B}^-(p) = |
| | |\mathcal{B}(\mathrm{prev}) \setminus (\mathcal{B}(p) \cap|
| | |\mathcal{B}(\mathrm{prev}))` |
+----------------------+-------------------------+---------------------------------------------------------+
Decomposable Structring Element Property
****************************************
A SE is said to be *decomposable*, if it enables to provides a SE
decomposition. E.g. a square of size n×m is decomposable in a vertical line
of of length *m* and an horizontal line of line *m*.
.. rubric:: `Type definition`
:class: concept-typedefs
+--------------------------------+--------------------+
|Type | Definition |
| | |
+================================+====================+
|`SE::is_decomposable` | `std::true_type` |
+--------------------------------+--------------------+
.. rubric:: `Valid expression`
:class: concept-expr
+--------------------------+-------------------------------+------------------------------------------------------+
|Expression | Return Type | Sementics |
+==========================+===============================+======================================================+
| ``se.decomposable()`` | `boolean` | Check dynamically if the SE is decomposable. |
| | | |
+--------------------------+-------------------------------+------------------------------------------------------+
| ``se.decompose()`` | `Range<WeightedNeighborhood>` | Decompose a SE into a set of non-separable SEs. If |
| | | `se` is not separable, it returns a range of one |
| | | element: itself. |
+--------------------------+-------------------------------+------------------------------------------------------+
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment