component_tree.hpp 12.3 KB
Newer Older
1
2
3
4
#pragma once


#include <mln/accu/accumulator.hpp>
5
#include <mln/core/algorithm/clone.hpp>
6
7
#include <mln/core/algorithm/fill.hpp>
#include <mln/core/image/ndimage.hpp>
8
#include <mln/core/image/view/zip.hpp>
9
#include <mln/core/neighborhood/c4.hpp>
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include <mln/core/range/foreach.hpp>
#include <mln/core/trace.hpp>

#include <vector>

#include <range/v3/view/span.hpp>

namespace mln::morpho
{

  enum ct_filtering
  {
    CT_FILTER_DIRECT,
    CT_FILTER_MIN,
    CT_FILTER_MAX,
    // CT_FILTER_SUBTRACTIVE (not yet implemented)
  };

28
29
30
31
32
33
34
  template <typename P, typename W>
  struct edge_t
  {
    P p;
    P q;
    W w;
  };
35
36
37
38
39
40
41
42
43
44

  template <class T = void>
  class component_tree;


  template <>
  class component_tree<void>
  {
  public:
    using node_id_type = int;
45
    using node_map_t   = int;
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92


    /// \brief Filter the tree given a predicate that removes some nodes according to the selected strategy.
    ///
    /// The parent array is updated, as well as the node_map if provided. It does not invalidate the node id, i.e. node
    /// that remain keeps the same position.
    ///
    /// The strategy is one of DIRECT, MIN, MAX, SUBTRACTIVE
    ///
    /// * Min
    /// A node is remove if it does not pass the predicate and
    /// all its descendant are removed as well.
    /// Formally, a component Γ is kept if 𝓒(Γ) with
    /// 𝓒(Γ) = ⋁ { 𝓟(Γ'), Γ ⊆ Γ' }.
    ///
    /// * Max
    /// A node is remove if every node in its childhood does not pass
    /// the predicate. All its descendant are removed as well.  Formally, a
    /// component Γ is kept if 𝓒(Γ) with 𝓒(Γ) = ⋀ { 𝓟(Γ'), Γ' ⊆ Γ }.
    ///
    /// * Direct
    /// A node is remove if it does not pass the predicate, the others
    /// remain.
    ///
    ///
    /// \param strategy The filtering strategy
    /// \param pred A boolean predicate for each node id
    template <class F>
    void filter(ct_filtering strategy, F pred);

    template <class I, class F>
    void filter(ct_filtering strategy, I node_map, F pred);


    /// \brief Compute the depth attribute over a tree
    std::vector<int> compute_depth() const;


    /// \brief Compute attribute on values
    ///
    /// Compute attribute with values given from an image
    ///
    /// \param node_map Image point -> node_id mapping
    /// \param values Image point -> value mapping
    /// \param acc Accumulator to apply on values
    template <class I, class J, class Accu>
    std::vector<typename accu::result_of<Accu, image_value_t<J>>::type> //
93
    compute_attribute_on_values(I node_map, J values, Accu acc) const;
94
95
96
97
98
99
100
101
102
103

    /// \brief Compute attribute on values
    ///
    /// Compute attribute on points
    ///
    /// \param node_map Image point -> node_id mapping
    /// \param values Image point -> value mapping
    /// \param acc Accumulator to apply on points
    template <class I, class Accu>
    std::vector<typename accu::result_of<Accu, image_point_t<I>>::type> //
104
    compute_attribute_on_points(I node_map, Accu acc) const;
105
106
107
108
109
110
111
112
113
114
115


    /// \brief Compute attribute on pixels
    ///
    /// Compute attribute with pixels (point + value) given from an image
    ///
    /// \param node_map Image point -> node_id mapping
    /// \param values Image point -> value mapping
    /// \param acc Accumulator to apply on values
    template <class I, class J, class Accu>
    std::vector<typename accu::result_of<Accu, image_pixel_t<J>>::type> //
116
    compute_attribute_on_pixels(I node_map, J values, Accu acc) const;
117

118
119
120
121
122
123
124
125
    /// \brief Compute the horizontal cut of a hierarchie at level `threshold` and return a nodemap
    /// valued with the node indices of the lowest nodes satisfying levels[n] > threshold
    ///
    /// \param threshold The threshold of the cut
    /// \param nodemap Image point -> node_id mapping
    /// \param levels Altitude of each node in the tree
    template <class T, class I, class V>
    I horizontal_cut_from_levels(const T threshold, I nodemap, ::ranges::span<V> levels) const;
126
127
128
129
130
131
132
133
134


    /// \brief Reconstruct an image from an attribute map
    ///
    /// \param node_map Image point -> node_id mapping
    /// \param values node_id -> value mapping
    template <class I, class V>
    image_ch_value_t<I, V> reconstruct_from(I node_map, ::ranges::span<V> values) const;

135
    using node_t = int;
136

137
    /// \brief Produce a visualization of the given Component Tree using the Khalimsky grid of the saliency map
138
    ///        The component_tree must be built on a 2D image with a 4-connectivity.
139
    ///
140
    /// \param node_map Image point -> node_id mapping
141
    mln::image2d<double> saliency(mln::image2d<int> node_map, ::ranges::span<double> values) const;
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

    std::vector<node_t> parent;

  private:
    template <class I, class F>
    void update_node_map(I node_map, F pred) const;


    void filter_direct(const std::vector<bool>& pred);

    template <class F>
    void filter_direct_T(F pred);
  };


  template <class T>
  class component_tree : public component_tree<void>
  {
  public:
161
162
163
164
165
    template <class I>
    I horizontal_cut(const T threshold, I nodemap) const
    {
      return this->horizontal_cut_from_levels(threshold, nodemap, ::ranges::make_span(values.data(), values.size()));
    }
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195

    template <class I>
    image_ch_value_t<std::remove_reference_t<I>, T> reconstruct(I&& node_map)
    {
      return this->reconstruct_from(std::forward<I>(node_map), ::ranges::make_span(values.data(), values.size()));
    }


    std::vector<T> values;
  };


  /******************************************/
  /****          Implementation          ****/
  /******************************************/

  template <class F>
  void component_tree<void>::filter_direct_T(F pred)
  {
    std::size_t n = parent.size();

    for (std::size_t i = 1; i < n; ++i)
      if (parent[i] > 0 && !pred(parent[i])) // If the parent[i] > 0 <=> not the root
        parent[i] = parent[parent[i]];
  }


  template <class I, class F>
  void component_tree<void>::update_node_map(I node_map, F pred) const
  {
196
    mln_foreach (auto& id, node_map.values())
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
    {
      if (id > 0 && !pred(id))
        id = this->parent[id];
    }
  }


  template <class F>
  void component_tree<void>::filter(ct_filtering strategy, F pred)
  {
    mln_entering("mln::morpho::component_tree::filter");

    std::size_t n = parent.size();

    // Node pass status
    std::vector<bool> pass;

    if (strategy == CT_FILTER_DIRECT)
    {
      this->filter_direct_T(pred);
    }
    else if (strategy == CT_FILTER_MIN)
    {
      // A node is removed if one ancestror is remove
      pass.resize(n);
      pass[0] = true;


      // Propagate downward
      for (std::size_t i = 1; i < n; ++i)
      {
        pass[i] = pred(i) && pass[parent[i]];
        if (!pass[parent[i]])
          parent[i] = parent[parent[i]];
      }
    }
    else if (strategy == CT_FILTER_MAX)
    {
      // A node is removed if all its descandants are removed
      pass.assign(n, false);

      // Propagate upward
      for (std::size_t i = n - 1; i > 0; --i)
      {
241
        pass[i]         = pass[i] || pred(i);
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
        pass[parent[i]] = pass[parent[i]] || pass[i];
      }

      this->filter_direct(pass);
    }
  }

  template <class I, class F>
  void component_tree<void>::filter(ct_filtering strategy, I node_map, F pred)
  {
    mln_entering("mln::morpho::component_tree::filter");

    std::size_t n = parent.size();

    // Node pass status
    std::vector<bool> pass;

    if (strategy == CT_FILTER_DIRECT)
    {
      this->filter_direct_T(pred);
      this->update_node_map(node_map, pred);
      return;
    }
    else if (strategy == CT_FILTER_MIN)
    {
      // A node is removed if one ancestror is removed
      pass.resize(n);
      pass[0] = true;


      // Propagate downward
      for (std::size_t i = 1; i < n; ++i)
      {
        pass[i] = pred(static_cast<int>(i)) && pass[parent[i]];
        if (!pass[parent[i]])
          parent[i] = parent[parent[i]];
      }
    }
    else if (strategy == CT_FILTER_MAX)
    {
      // A node is removed if all its descandants are removed
      pass.assign(n, false);

      // Propagate upward
      for (std::size_t i = n - 1; i > 0; --i)
      {
288
        pass[i]         = pass[i] || pred(static_cast<int>(i));
289
290
291
292
293
294
295
296
297
298
299
        pass[parent[i]] = pass[parent[i]] || pass[i];
      }
      this->filter_direct(pass);
    }

    this->update_node_map(node_map, [&pass](int x) { return pass[x]; });
  }


  template <class I, class Accu>
  std::vector<typename accu::result_of<Accu, image_point_t<I>>::type>
300
  component_tree<void>::compute_attribute_on_points(I node_map, Accu acc) const
301
302
303
304
305
306
307
308
309
310
311
312
  {
    mln_entering("mln::morpho::component_tree::compute_attribute_on_points");

    static_assert(mln::is_a<I, mln::details::Image>());
    static_assert(mln::is_a<Accu, mln::AccumulatorLike>());
    using R = typename accu::result_of<Accu, image_point_t<I>>::type;

    auto a = accu::make_accumulator(acc, image_point_t<I>());

    std::vector<decltype(a)> attr(parent.size(), a);

    // Accumulate for each point
313
    mln_foreach (auto px, node_map.pixels())
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
      attr[px.val()].take(px.point());


    // Propgate to parent
    std::size_t n = parent.size();
    for (std::size_t i = n - 1; i > 0; --i)
      attr[parent[i]].take(attr[i]);


    // Extract values
    std::vector<R> out(n);
    for (std::size_t i = 0; i < n; ++i)
      out[i] = attr[i].to_result();

    return out;
  }

  template <class I, class J, class Accu>
  std::vector<typename accu::result_of<Accu, image_value_t<J>>::type> //
333
  component_tree<void>::compute_attribute_on_values(I node_map, J input, Accu acc) const
334
335
336
337
338
339
340
341
342
343
344
345
346
  {
    mln_entering("mln::morpho::component_tree::compute_attribute_on_values");

    static_assert(mln::is_a<I, mln::details::Image>());
    static_assert(mln::is_a<Accu, mln::AccumulatorLike>());
    using R = typename accu::result_of<Accu, image_value_t<J>>::type;

    auto a = accu::make_accumulator(acc, image_value_t<J>());

    std::vector<decltype(a)> attr(parent.size(), a);

    // Accumulate for each point
    auto zz = mln::view::zip(node_map, input);
347
    mln_foreach ((auto [node_id, val]), zz.values())
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
      attr[node_id].take(val);


    // Propgate to parent
    std::size_t n = parent.size();
    for (std::size_t i = n - 1; i > 0; --i)
      attr[parent[i]].take(attr[i]);


    // Extract values
    std::vector<R> out(n);
    for (std::size_t i = 0; i < n; ++i)
      out[i] = attr[i].to_result();

    return out;
  }

  template <class I, class J, class Accu>
  std::vector<typename accu::result_of<Accu, image_pixel_t<J>>::type> //
367
  component_tree<void>::compute_attribute_on_pixels(I node_map, J values, Accu acc) const
368
369
370
371
372
373
374
375
376
377
378
379
  {
    mln_entering("mln::morpho::component_tree::compute_attribute_on_pixels");

    static_assert(mln::is_a<I, mln::details::Image>());
    static_assert(mln::is_a<Accu, mln::AccumulatorLike>());
    using R = typename accu::result_of<Accu, image_pixel_t<J>>::type;

    auto a = accu::make_accumulator(acc, image_pixel_t<J>());

    std::vector<decltype(a)> attr(parent.size(), a);

    // Accumulate for each point
Edwin Carlinet's avatar
Edwin Carlinet committed
380
    mln_foreach (auto px, values.pixels())
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
      attr[node_map(px.point())].take(px);


    // Propgate to parent
    std::size_t n = parent.size();
    for (std::size_t i = n - 1; i > 0; --i)
      attr[parent[i]].take(attr[i]);


    // Extract values
    std::vector<R> out(n);
    for (std::size_t i = 0; i < n; ++i)
      out[i] = attr[i].to_result();

    return out;
  }

398
399
400
401
  template <class T, class I, class V>
  I component_tree<void>::horizontal_cut_from_levels(const T threshold, I nodemap, ::ranges::span<V> levels) const
  {
    mln_entering("mln::morpho::component_tree::horizontal_cut_from_levels");
402

403
404
405
406
    auto root_cut_cc = std::vector<int>(parent.size());
    for (std::size_t node = 0; node < parent.size(); ++node)
    {
      int parent_node   = parent[node];
Baptiste Esteban's avatar
Baptiste Esteban committed
407
      root_cut_cc[node] = levels[parent_node] > threshold ? static_cast<int>(node) : root_cut_cc[parent_node];
408
409
410
411
412
413
414
    }

    auto out = mln::clone(nodemap);
    mln_foreach (auto px, out.pixels())
      out(px.point()) = root_cut_cc[px.val()];
    return out;
  }
415
416
417
418
419
420
421
422
423

  template <class I, class V>
  image_ch_value_t<I, V> component_tree<void>::reconstruct_from(I node_map, ::ranges::span<V> values) const
  {
    mln_entering("mln::morpho::component_tree::reconstruction");

    image_ch_value_t<I, V> out = imchvalue<V>(node_map);

    auto zz = mln::view::zip(node_map, out);
Edwin Carlinet's avatar
Edwin Carlinet committed
424
    mln_foreach ((auto&& [node_id, val]), zz.values())
425
426
427
428
      val = values[node_id];

    return out;
  }
429
} // namespace mln::morpho