Merge branch 'development/se-decomposition' into 'dev'

SE decomposition of disc and rect2d for dilation based algorithms.

See merge request !11

.gitlab-ci.yml

@@ -3,10 +3,11 @@ stages:
     - test
     - deploy
 
+image: carlinet/buildimage
+
 
 .job-build-linux-base: &distcheck-linux-base
     stage: test
-    image: carlinet/buildimage
     script:
         - mkdir build && cd build
         - cmake .. -DCMAKE_BUILD_TYPE=$PYLENE_CONFIGURATION
@@ -42,7 +43,6 @@ distcheck-linux-clang-debug:
         CC: "clang"
 
 build-doc-images:
-    image: carlinet/buildimage
     stage: build
     script:
         - mkdir build && cd build
@@ -58,7 +58,6 @@ build-doc-images:
 #################
 
 .documentation-base:  &documentation-base
-    image: carlinet/doc-image
     stage: test
     script:
         - cd doc

bench/BMDilation.cpp

+#include <mln/core/extension/fill.hpp>
+#include <mln/core/image/image2d.hpp>
+#include <mln/io/imread.hpp>
+#include <mln/io/imsave.hpp>
+
+#include <benchmark/benchmark.h>
+
+#include <mln/core/se/disc.hpp>
+#include <mln/core/se/rect2d.hpp>
+#include <mln/morpho/structural/dilate.hpp>
+
+using namespace mln;
+
+
+class BMDilation : public benchmark::Fixture
+{
+public:
+  BMDilation()
+  {
+    const char* filename = MLN_IMG_PATH "/lena.pgm";
+    io::imread(filename, m_input);
+    int nr = m_input.nrows();
+    int nc = m_input.ncols();
+    resize(m_output, m_input);
+    m_bytes = nr * nc * sizeof(uint8);
+  }
+
+  template <class SE, class Compare = mln::productorder_less<mln::uint8>>
+  void run(benchmark::State& st, const SE& se, Compare cmp = Compare())
+  {
+    for (auto _ : st)
+      mln::morpho::structural::dilate(m_input, se, m_output, cmp);
+    st.SetBytesProcessed(st.iterations() * m_bytes);
+  }
+
+protected:
+  image2d<uint8> m_input;
+  image2d<uint8> m_output;
+  std::size_t m_bytes;
+};
+
+BENCHMARK_DEFINE_F(BMDilation, EuclideanDisc_naive)(benchmark::State& st)
+{
+  int radius = st.range(0);
+  mln::se::disc se(radius, 0);
+  this->run(st, se, std::less<mln::uint8>());
+}
+
+BENCHMARK_DEFINE_F(BMDilation, EuclideanDisc_incremental)(benchmark::State& st)
+{
+  int radius = st.range(0);
+  mln::se::disc se(radius, 0);
+  this->run(st, se);
+}
+
+BENCHMARK_DEFINE_F(BMDilation, ApproximatedDisc)(benchmark::State& st)
+{
+  int radius = st.range(0);
+  mln::se::disc se(radius);
+  this->run(st, se);
+}
+
+
+BENCHMARK_DEFINE_F(BMDilation, Square)(benchmark::State& st)
+{
+  int width = 2 * st.range(0) + 1;
+  mln::se::rect2d se(width, width);
+  this->run(st, se);
+}
+
+constexpr int max_range = 2 << 6;
+
+BENCHMARK_REGISTER_F(BMDilation, ApproximatedDisc)->RangeMultiplier(2)->Range(2, max_range);
+BENCHMARK_REGISTER_F(BMDilation, EuclideanDisc_naive)->RangeMultiplier(2)->Range(2, max_range);
+BENCHMARK_REGISTER_F(BMDilation, EuclideanDisc_incremental)->RangeMultiplier(2)->Range(2, max_range);
+BENCHMARK_REGISTER_F(BMDilation, Square)->RangeMultiplier(2)->Range(2, max_range);
+
+
+BENCHMARK_MAIN();

bench/CMakeLists.txt

@@ -35,7 +35,9 @@ link_libraries(benchmark::benchmark)
 add_executable(BMErosion        BMErosion.cpp BMErosion_main.cpp)
 add_executable(BMNeighborhood   BMNeighborhood.cpp BMNeighborhood_main.cpp)
 add_executable(BMRotation       BMRotation.cpp)
+add_executable(BMDilation       BMDilation.cpp)
 add_executable(BMMorphers       BMMorphers.cpp   BMMorphers_main.cpp)
 target_link_libraries(BMRotation ${FreeImage_LIBRARIES})
 target_link_libraries(BMMorphers ${FreeImage_LIBRARIES})
 target_link_libraries(BMErosion ${FreeImage_LIBRARIES})
+target_link_libraries(BMDilation ${FreeImage_LIBRARIES})

bench/plotbench2.py

+#! /usr/bin/python3
+
+import json
+import sys
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+data = json.load(open(sys.argv[1]))
+entries = data["benchmarks"]
+
+for entry in entries:
+    _, algo, size = entry["name"].split("/")
+    entry["algo"] = algo
+    entry["size"] = int(size)
+    entry["bytes_per_second"] = entry["bytes_per_second"] / (2 ** 20)
+    del entry["name"]
+    del entry["cpu_time"]
+    del entry["iterations"]
+    del entry["real_time"]
+    del entry["time_unit"]
+
+df = pd.DataFrame(entries)
+
+
+df.set_index(["size", "algo"], inplace=True)
+df = df.unstack()["bytes_per_second"]
+df.plot(logy=True)
+
+
+plt.xlabel("SE Size")
+plt.ylabel("Performance (MB/s)")
+plt.show()

doc/source/conf.py

@@ -28,7 +28,8 @@ import pylene_lexer
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.ifconfig', 'breathe', "ext.conceptdoc"]
+extensions = ['sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.ifconfig', 'breathe', "ext.conceptdoc",
+              'matplotlib.sphinxext.plot_directive']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -254,6 +255,10 @@ breathe_default_project = "Pylene"
 breathe_path = os.path.realpath("../doxyxml/")
 
 
+# Matplotlib
+plot_include_source = False
+plot_html_show_formats = False
+plot_formats = [ ("png", 150) ]
 
 # Substitution
 rst_prolog = """

doc/source/core.rst

+Core Module
+==============================
+
+Structural Elements
+*******************
+
+.. toctree::
+   :maxdepth: 1
+
+   core/se/disc
+   core/se/rectangle
+   core/se/periodic_lines

doc/source/core/se/disc.rst

+Disc
+====
+
+.. doxygenstruct:: mln::se::disc
+   :members:
+
+
+
+Approximated discs
+------------------
+
+.. rubric:: Quality of the approximation
+
+The approximated disc 𝔇ₐ of the true euclidean disc 𝔇 uses a radial
+decomposition in 8 periodic lines :math:`k_i.L_\theta`
+of orientations θ ∈ \{0°, 30°, 45°, 60°, 90°, 120°, 135°, 150°\} and length kᵢ.
+
+
+Such approximation introduces an error:
+
+.. math::
+
+   err = \frac{|\mathcal{D} \Delta \mathcal{D}_a|}{r}
+
+Below is the error of our approximation, the 8-approximation of matlab and the
+best approximation possible using 8 periodic lines (computed by exhaustive
+search).
+
+.. |euc| image:: ./disc_euclidean.gif
+                 :width: 200px
+
+.. |approx| image:: ./disc_approximated.gif
+                    :width: 200px
+
+.. |matlab| image:: ./disc_matlab.gif
+                    :width: 200px
+
+.. table::
+   :widths: 33 33 33
+
+   +---------------------+---------------------------+-------------------------+
+   | |euc|               |   |matlab|                |   |approx|              |
+   +---------------------+---------------------------+-------------------------+
+   |   Euclidean disc    | Matlab 8-appoximated disc |  Our appoximated disc   |
+   +---------------------+---------------------------+-------------------------+
+
+
+
+
+
+.. figure:: ./disc_error.png
+   :width: 100%
+
+   Error of our disc approximation (labeled custom above), matlab approximation
+   and the best appromximation relative to the euclidean disc.
+
+
+.. rubric:: Performance
+
+Using approximated disc can speed up the running of some morphological
+operations. Below is the running time of the *dilation* by the euclidean disc
+*vs* the approximated disc.
+
+.. plot::
+   :width: 100%
+
+   from pyplots import plotbysize
+   plotbysize("disc_dilation.json")
+
+
+Dilation by a square is given as reference. The running time of the dilation by
+the approximated disc does not depend on the radius (like for a square) because
+it uses a decomposition in periodic lines (the SE is decomposable) [Adam93]_
+[JoSo96a]_ [JoSo96b]_. The dilation by the euclidean disc is :math:`O(r.n)`
+because of the optimization of incremental SEs by [VaTa96]_. It contrasts with
+the naive implementation which is :math:`O(r^2.n)` for disc.
+
+
+
+References
+----------
+
+.. [Adam93] Adams, R. (1993), ‘Radial decomposition of discs and spheres.’, Computer
+            Vision, Graphics, and Image Processing: Graphical Models and Image Processing
+            55(5), 325–332.
+
+.. [JoSo96a] Jones, R., & Soille, P. (1996). Periodic lines: Definition, cascades, and
+             application to granulometries. Pattern Recognition Letters, 17(10), 1057-1063.
+
+.. [JoSo96b] Jones, R., & Soille, P. (1996). Periodic lines and their application to
+             granulometries. In Mathematical Morphology and its Applications to Image and
+             Signal Processing (pp. 263-272). Springer, Boston, MA.
+
+.. [VaTa96] Van Droogenbroeck, M., & Talbot, H. (1996). Fast computation of
+            morphological operations with arbitrary structuring elements. Pattern
+            recognition letters, 17(14), 1451-1460.

doc/source/core/se/disc_dilation.json

+{
+  "context": {
+    "date": "2018-01-26 17:05:14",
+    "num_cpus": 4,
+    "mhz_per_cpu": 3500,
+    "cpu_scaling_enabled": true,
+    "caches": [
+      {
+        "type": "Data",
+        "level": 1,
+        "size": 32000000,
+        "num_sharing": 2
+      },
+      {
+        "type": "Instruction",
+        "level": 1,
+        "size": 32000000,
+        "num_sharing": 2
+      },
+      {
+        "type": "Unified",
+        "level": 2,
+        "size": 256000000,
+        "num_sharing": 2
+      },
+      {
+        "type": "Unified",
+        "level": 3,
+        "size": 4096000000,
+        "num_sharing": 4
+      }
+    ],
+    "library_build_type": "release"
+  },
+  "benchmarks": [
+    {
+      "name": "BMDilation/Approximated Disc/2",
+      "iterations": 311,
+      "real_time": 2.2145544437396228e+06,
+      "cpu_time": 2.2122267459807075e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1849779886997448e+08
+    },
+    {
+      "name": "BMDilation/Approximated Disc/4",
+      "iterations": 152,
+      "real_time": 4.6295076842254372e+06,
+      "cpu_time": 4.6245967171052638e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 5.6684726482288241e+07
+    },
+    {
+      "name": "BMDilation/Approximated Disc/8",
+      "iterations": 96,
+      "real_time": 7.0571836354247360e+06,
+      "cpu_time": 7.0493722187500028e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 3.7186857476861037e+07
+    },
+    {
+      "name": "BMDilation/Approximated Disc/16",
+      "iterations": 101,
+      "real_time": 6.7710320990144201e+06,
+      "cpu_time": 6.7638077722772276e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 3.8756867259659246e+07
+    },
+    {
+      "name": "BMDilation/Approximated Disc/32",
+      "iterations": 77,
+      "real_time": 8.8636344285491351e+06,
+      "cpu_time": 8.8537914415584449e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 2.9608106507855281e+07
+    },
+    {
+      "name": "BMDilation/Approximated Disc/64",
+      "iterations": 79,
+      "real_time": 8.7361794683579020e+06,
+      "cpu_time": 8.7265214050632920e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 3.0039919440053072e+07
+    },
+    {
+      "name": "BMDilation/Approximated Disc/128",
+      "iterations": 77,
+      "real_time": 8.8617318051783275e+06,
+      "cpu_time": 8.8521802857142892e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 2.9613495380686026e+07
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/2",
+      "iterations": 84,
+      "real_time": 8.2215158571281275e+06,
+      "cpu_time": 8.2140779166666707e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 3.1913989940136805e+07
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/4",
+      "iterations": 18,
+      "real_time": 3.8401868888791390e+07,
+      "cpu_time": 3.8368693055555545e+07,
+      "time_unit": "ns",
+      "bytes_per_second": 6.8322368869961603e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/8",
+      "iterations": 6,
+      "real_time": 1.2213701449945801e+08,
+      "cpu_time": 1.2203581916666669e+08,
+      "time_unit": "ns",
+      "bytes_per_second": 2.1480906326525728e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/16",
+      "iterations": 2,
+      "real_time": 4.3583392350046778e+08,
+      "cpu_time": 4.3536064950000066e+08,
+      "time_unit": "ns",
+      "bytes_per_second": 6.0213067097604007e+05
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/32",
+      "iterations": 1,
+      "real_time": 1.6635573169987764e+09,
+      "cpu_time": 1.6620425710000007e+09,
+      "time_unit": "ns",
+      "bytes_per_second": 1.5772399851483703e+05
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/64",
+      "iterations": 1,
+      "real_time": 6.4490619220014200e+09,
+      "cpu_time": 6.4423174940000010e+09,
+      "time_unit": "ns",
+      "bytes_per_second": 4.0690947045709196e+04
+    },
+    {
+      "name": "BMDilation/Euclidean Disc (naive)/128",
+      "iterations": 1,
+      "real_time": 2.5237285107999924e+10,
+      "cpu_time": 2.5207719910000000e+10,
+      "time_unit": "ns",
+      "bytes_per_second": 1.0399353885870751e+04
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/2",
+      "iterations": 37,
+      "real_time": 1.8754427108089317e+07,
+      "cpu_time": 1.8732530891891714e+07,
+      "time_unit": "ns",
+      "bytes_per_second": 1.3994051391820619e+07
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/4",
+      "iterations": 24,
+      "real_time": 2.9199526999946102e+07,
+      "cpu_time": 2.9164850375000089e+07,
+      "time_unit": "ns",
+      "bytes_per_second": 8.9883540161998589e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/8",
+      "iterations": 18,
+      "real_time": 3.8627775000147969e+07,
+      "cpu_time": 3.8580057722222470e+07,
+      "time_unit": "ns",
+      "bytes_per_second": 6.7948058006404340e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/16",
+      "iterations": 12,
+      "real_time": 5.7734777083472483e+07,
+      "cpu_time": 5.7669801250000320e+07,
+      "time_unit": "ns",
+      "bytes_per_second": 4.5456026259497218e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/32",
+      "iterations": 7,
+      "real_time": 1.0060752242861782e+08,
+      "cpu_time": 1.0049030228571354e+08,
+      "time_unit": "ns",
+      "bytes_per_second": 2.6086497307439027e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/64",
+      "iterations": 4,
+      "real_time": 1.9608837674968526e+08,
+      "cpu_time": 1.9586752849999911e+08,
+      "time_unit": "ns",
+      "bytes_per_second": 1.3383739612561720e+06
+    },
+    {
+      "name": "BMDilation/Euclidean Disc/128",
+      "iterations": 2,
+      "real_time": 4.0451671699884170e+08,
+      "cpu_time": 4.0401250300000060e+08,
+      "time_unit": "ns",
+      "bytes_per_second": 6.4885120646872558e+05
+    },
+    {
+      "name": "BMDilation/Square/2",
+      "iterations": 282,
+      "real_time": 2.4817838581600506e+06,
+      "cpu_time": 2.4786692056737524e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.0575997773319007e+08
+    },
+    {
+      "name": "BMDilation/Square/4",
+      "iterations": 297,
+      "real_time": 2.3502951144752908e+06,
+      "cpu_time": 2.3473787710437682e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1167520267018388e+08
+    },
+    {
+      "name": "BMDilation/Square/8",
+      "iterations": 306,
+      "real_time": 2.2860967222192390e+06,
+      "cpu_time": 2.2829784019607869e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1482544021215960e+08
+    },
+    {
+      "name": "BMDilation/Square/16",
+      "iterations": 309,
+      "real_time": 2.2427470776694939e+06,
+      "cpu_time": 2.2398339223301108e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1703724878284290e+08
+    },
+    {
+      "name": "BMDilation/Square/32",
+      "iterations": 307,
+      "real_time": 2.2185559315960179e+06,
+      "cpu_time": 2.2156340130293099e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1831556947511630e+08
+    },
+    {
+      "name": "BMDilation/Square/64",
+      "iterations": 311,
+      "real_time": 2.2275530482323049e+06,
+      "cpu_time": 2.2245858778135087e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1783946064498754e+08
+    },
+    {
+      "name": "BMDilation/Square/128",
+      "iterations": 310,
+      "real_time": 2.2337310741990693e+06,
+      "cpu_time": 2.2308704354838571e+06,
+      "time_unit": "ns",
+      "bytes_per_second": 1.1750749654949960e+08
+    }
+  ]
+}

doc/source/core/se/periodic_lines.rst

+Periodic Line
+=============
+
+.. doxygenstruct:: mln::se::periodic_line2d
+   :members:
+
+
+.. rubric:: Performance
+
+The speed of a dilation (or erosion) by a periodic line does not depend on its
+length.
+
+
+

doc/source/core/se/rectangle.rst

+Rectangle
+=========
+
+.. doxygenstruct:: mln::se::rect2d
+   :members:
+
+
+.. rubric:: Performance
+
+The rectangle is decomposable in two lines, so that the speed of a dilation (or
+erosion) by a rectangle does not depend on the size of the rectangle.
+
+
+Example 1. Dilation by rectangle
+--------------------------------
+
+.. code-block:: pylenecpp
+
+   #include <mln/core/se/rect2d.hpp>
+
+   mln::se::rect2d (21, 21);
+   auto output = mln::morpho::structural::dilate(input, rect);
+
+
+.. image:: /images/lena_gray.jpg
+           :width: 49%
+
+.. image:: /images/morpho_dilation_1.png
+           :width: 49%

doc/source/index.rst

@@ -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`
 
 
+

doc/source/morpho/dilation.rst

@@ -59,7 +59,7 @@ Example 1 : Dilation by a square on a gray-level image
 
    // Define a square SE of size 21x21
    auto input = ...;
-   auto rect = mln::make_rectangle2d(21,21);
+   mln::se::rect2d rect(21,21);
    auto output = mln::morpho::structural::dilate(input, rect);
 
 

doc/source/pyplots/__init__.py

+#! /usr/bin/python3
+import json
+import sys
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+def plotbysize(filename):
+    data = json.load(open(filename))
+    entries = data["benchmarks"]
+
+    for entry in entries:
+        _, algo, size = entry["name"].split("/")
+        entry["algo"] = algo
+        entry["size"] = int(size)
+        entry["bytes_per_second"] = entry["bytes_per_second"] / (2 ** 20)
+        del entry["name"]
+        del entry["cpu_time"]
+        del entry["iterations"]
+        del entry["real_time"]
+        del entry["time_unit"]
+
+    df = pd.DataFrame(entries)
+
+
+    df.set_index(["size", "algo"], inplace=True)
+    df = df.unstack()["bytes_per_second"]
+    df.plot(logy=True)
+
+
+    plt.xlabel("SE Size")
+    plt.ylabel("Performance (MB/s)")
+    plt.show()

doc/source/reference.rst

@@ -5,6 +5,7 @@ Reference
 .. toctree::
    concept
    images
+   core
    algorithms
    morpho
 

doc/source/ruminations.rst

+Ruminations
+###########