component_tree.hpp 12.7 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


    /// \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
91
    /// \param propagate Boolean to propagate the values to the parent
92
93
    template <class I, class J, class Accu>
    std::vector<typename accu::result_of<Accu, image_value_t<J>>::type> //
94
    compute_attribute_on_values(I node_map, J values, Accu acc, bool propagate = true) const;
95
96
97
98
99
100
101
102

    /// \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
103
    /// \param propagate Boolean to propagate the values to the parent
104
105
    template <class I, class Accu>
    std::vector<typename accu::result_of<Accu, image_point_t<I>>::type> //
106
    compute_attribute_on_points(I node_map, Accu acc, bool propagate = true) const;
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
116
    /// \param propagate Boolean to propagate the values to the parent
117
118
    template <class I, class J, class Accu>
    std::vector<typename accu::result_of<Accu, image_pixel_t<J>>::type> //
119
    compute_attribute_on_pixels(I node_map, J values, Accu acc, bool propagate = true) const;
120

121
122
123
124
125
126
127
128
    /// \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;
129
130
131
132
133
134
135
136
137


    /// \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;

138
    using node_t = int;
139

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

    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:
164
165
166
167
168
    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()));
    }
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
196
197
198

    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
  {
199
    mln_foreach (auto& id, node_map.values())
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
241
242
243
    {
      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)
      {
244
        pass[i]         = pass[i] || pred(i);
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
288
289
290
        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)
      {
291
        pass[i]         = pass[i] || pred(static_cast<int>(i));
292
293
294
295
296
297
298
299
300
301
302
        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>
303
  component_tree<void>::compute_attribute_on_points(I node_map, Accu acc, bool propagate) const
304
305
306
307
308
309
310
311
312
313
314
315
  {
    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
316
    mln_foreach (auto px, node_map.pixels())
317
318
      attr[px.val()].take(px.point());

319
320
321
322
323
324
325
    const std::size_t n = parent.size();
    if (propagate)
    {
      // Propgate to parent
      for (std::size_t i = n - 1; i > 0; --i)
        attr[parent[i]].take(attr[i]);
    }
326
327
328
329
330
331
332
333
334
335
336

    // 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> //
337
  component_tree<void>::compute_attribute_on_values(I node_map, J input, Accu acc, bool propagate) const
338
339
340
341
342
343
344
345
346
347
348
349
350
  {
    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);
351
    mln_foreach ((auto [node_id, val]), zz.values())
352
353
      attr[node_id].take(val);

354
355
356
357
358
359
360
    const std::size_t n = parent.size();
    if (propagate)
    {
      // Propgate to parent
      for (std::size_t i = n - 1; i > 0; --i)
        attr[parent[i]].take(attr[i]);
    }
361
362
363
364
365
366
367
368
369
370
371

    // 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> //
372
  component_tree<void>::compute_attribute_on_pixels(I node_map, J values, Accu acc, bool propagate) const
373
374
375
376
377
378
379
380
381
382
383
384
  {
    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
385
    mln_foreach (auto px, values.pixels())
386
387
      attr[node_map(px.point())].take(px);

388
389
390
391
392
393
394
    const std::size_t n = parent.size();
    if (propagate)
    {
      // Propgate to parent
      for (std::size_t i = n - 1; i > 0; --i)
        attr[parent[i]].take(attr[i]);
    }
395
396
397
398
399
400
401
402
403

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

    return out;
  }

404
405
406
407
  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");
408

409
410
411
412
    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
413
      root_cut_cc[node] = levels[parent_node] > threshold ? static_cast<int>(node) : root_cut_cc[parent_node];
414
415
416
417
418
419
420
    }

    auto out = mln::clone(nodemap);
    mln_foreach (auto px, out.pixels())
      out(px.point()) = root_cut_cc[px.val()];
    return out;
  }
421
422
423
424
425
426
427
428
429

  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
430
    mln_foreach ((auto&& [node_id, val]), zz.values())
431
432
433
434
      val = values[node_id];

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