alphatree.hpp 10.8 KB
Newer Older
1
2
3
4
5
6
7
#pragma once


#include <mln/core/concepts/image.hpp>
#include <mln/core/concepts/neighborhood.hpp>

#include <mln/core/algorithm/for_each.hpp>
8
#include <mln/core/functional_ops.hpp>
9
#include <mln/morpho/canvas/unionfind.hpp>
10
11
12
#include <mln/morpho/component_tree.hpp>

#include <mln/morpho/private/directional_hqueue.hpp>
13
14

#include <algorithm>
15
#include <type_traits>
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#include <range/v3/functional/concepts.hpp>


namespace mln::morpho
{

  /// Compute the alphatree of an image
  ///
  /// \param input The input image
  /// \param neighborhood The neighborhood relation
  /// \param distance Distance function
  template <class I, class N, class F = mln::functional::l2dist_t<>>
  std::pair<component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>, image_ch_value_t<I, int>> //
  alphatree(I input, N nbh, F distance = F{});


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


  namespace internal
  {
40
41
    /// \brief Canvas for the edges in the alphatree. Using different data
    /// structures related to the type of the edges.
42
    template <typename P, typename N, typename W>
43
    class alphatree_edges;
44

45
46
    template <typename P, typename N, typename W>
    requires(std::is_integral_v<W>&& std::is_unsigned_v<W> && sizeof(W) <= 2) class alphatree_edges<P, N, W>
47
    {
48
    public:
49
      void                push(int dir, W w, P p) { m_cont.insert(dir, w, p); }
50
51
52
53
54
55
      std::tuple<P, P, W> pop() { return m_cont.pop(); }
      W                   top() const { return m_cont.current_level(); }
      bool                empty() const { return m_cont.empty(); }
      void                on_finish_insert() {}

    private:
56
      details::directional_hqueue<P, N, W> m_cont;
57
58
    };

Baptiste Esteban's avatar
Baptiste Esteban committed
59
60
61
62
63
64
65
66
    template <typename P, typename W>
    struct edge_t
    {
      P p;
      P q;
      W w;
    };

67
    template <typename P, typename N, typename W>
68
69
70
    class alphatree_edges
    {
    public:
71
      void                push(int dir, W w, P p) { m_cont.push_back({p, p + cn.after_offsets()[dir], w}); }
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
      std::tuple<P, P, W> pop()
      {
        assert(m_current < m_cont.size());
        const auto e = m_cont[m_current++];
        return {std::move(e.p), std::move(e.q), std::move(e.w)};
      }

      W top() const
      {
        assert(m_current < m_cont.size());
        return m_cont[m_current].w;
      }
      bool empty() const { return m_cont.size() == m_current; }
      void on_finish_insert()
      {
Baptiste Esteban's avatar
Baptiste Esteban committed
87
        std::sort(m_cont.begin(), m_cont.end(), [](const edge_t<P, W>& a, const edge_t<P, W>& b) { return a.w < b.w; });
88
89
90
      }

    private:
Baptiste Esteban's avatar
Baptiste Esteban committed
91
92
      static constexpr auto     cn = N();
      std::vector<edge_t<P, W>> m_cont;
93
      std::size_t               m_current = 0;
94
95
96
97
98
    };


    // Compute a node_id for each flat zone
    template <class I, class J>
99
    [[gnu::noinline]] std::size_t alphatree_create_nodemap(I node_map, J zpar)
100
101
    {
      std::size_t node_count = 0;
102
      mln_foreach (auto px, node_map.pixels())
103
104
105
106
107
108
109
110
111
112
113
114
115
      {
        auto p  = px.point();
        auto rp = canvas::impl::zfindroot(zpar, p);
        if (node_map(rp) < 0)
          node_map(rp) = static_cast<int>(node_count++);
        node_map(p) = node_map(rp);
      }
      return node_count;
    }

    // Compute flat zone of the image using union-find structures
    // It returns the number of edges processed
    template <class E, class J>
116
    void alphatree_compute_flatzones(E& edges, J zpar)
117
118
119
    {
      canvas::impl::union_find_init_par(zpar);

120
      while (!edges.empty() && edges.top() == 0)
121
      {
122
123
124
        const auto [p, q, w] = edges.pop();
        auto rp              = canvas::impl::zfindroot(zpar, p);
        auto rq              = canvas::impl::zfindroot(zpar, q);
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
        if (rp != rq) // Merge set
          zpar(rq) = rp;
      }
    }

    // Inverse the node array, updating the parent relation
    template <class W>
    [[gnu::noinline]] void alphatree_reorder_nodes(int* parent, W* levels, std::size_t n)
    {
      std::size_t i = 0;
      std::size_t j = n - 1;

      while (i < j)
      {
        std::swap(levels[i], levels[j]);
        std::swap(parent[i], parent[j]);
        parent[i] = static_cast<int>(n) - parent[i] - 1;
        parent[j] = static_cast<int>(n) - parent[j] - 1;
        i++;
        j--;
      }
      if (i == j)
        parent[i] = static_cast<int>(n) - parent[i] - 1;
    }


151
152
    template <class I, class N, class F, class C>
    void alphatree_compute_edges(I input, N nbh, F distance, C& edges)
153
    {
154
      static_assert(is_a_v<I, mln::details::Image>);
155

156
      auto dom = input.domain();
Edwin Carlinet's avatar
Edwin Carlinet committed
157
      mln_foreach (auto p, dom)
158
      {
159
        int i = 0;
160
        for (auto n : nbh.after(p))
161
        {
162
          if (dom.has(n))
163
164
165
            edges.push(i, distance(input(p), input(n)), p);
          ++i;
        }
166
      }
167
      edges.on_finish_insert();
168
169
170
171
    }

    //
    //
172
    template <class E, class I, class W, class M>
173
174
175
    std::size_t alphatree_compute_hierarchy(E& edges, I node_map,         //
                                            std::size_t       node_count, //
                                            std::vector<int>& par,        //
176
177
                                            std::vector<W>&   levels,     //
                                            std::vector<M>*   mst)
178
179
180
    {
      static_assert(mln::is_a<I, mln::details::Image>());

181
182
      std::vector<int> zpar(node_count);  // Parent of ufind structure
      std::vector<int> links(node_count); // Node index in the alpha tree
183
184
185
186

      std::iota(std::begin(zpar), std::end(zpar), 0);
      std::iota(std::begin(links), std::end(links), 0);

187
      while (!edges.empty())
188
      {
189
        auto [p, q, w] = edges.pop();
190
191
192
193
194
195
196
197
198
199
200
201
202
        auto rp        = canvas::impl::zfindroot(zpar.data(), node_map(p));
        auto rq        = canvas::impl::zfindroot(zpar.data(), node_map(q));

        if (rp != rq)
        {
          // Merge set

          assert(levels[rp] <= w && levels[rq] <= w);
          zpar[rq] = rp;

          // Do we need to create a new node
          int rp_root = links[rp];
          int rq_root = links[rq];
203
204
205
206

          int max_root = std::max(rp_root, rq_root);
          int min_root = std::min(rp_root, rq_root);

207
          int new_root_id;
208
          if (levels[max_root] == w)
209
          {
210
            new_root_id = max_root;
211
          }
212
          else if (levels[min_root] == w)
213
          {
214
            new_root_id = min_root;
215
216
217
218
219
220
221
222
          }
          else
          {
            // Create a new_node
            new_root_id = static_cast<int>(node_count++);
            par.push_back(new_root_id);
            levels.push_back(w);
          }
223

224
225
226
          par[rp_root] = new_root_id;
          par[rq_root] = new_root_id;
          links[rp]    = new_root_id;
227
228
229

          if (mst != nullptr)
            mst->push_back({p, q, w});
230
231
232
233
234
        }
      }
      return node_count;
    }

Quentin Kaci's avatar
Quentin Kaci committed
235
236
237
238
    template <class W>
    std::pair<std::vector<int>, std::vector<W>> canonize_component_tree(const std::vector<int>& par,    //
                                                                        const std::vector<W>&   levels, //
                                                                        std::size_t             node_count)
Quentin Kaci's avatar
Quentin Kaci committed
239
    {
Quentin Kaci's avatar
Quentin Kaci committed
240
241
      std::vector<int> canonized_par;
      std::vector<W>   canonized_levels;
242

243
      // Root initialization
Quentin Kaci's avatar
Quentin Kaci committed
244
245
      canonized_par.push_back(0);
      canonized_levels.push_back(levels[0]);
246

Quentin Kaci's avatar
Quentin Kaci committed
247
      std::vector<int> translation_map(node_count);
248

Quentin Kaci's avatar
Quentin Kaci committed
249
      translation_map[0] = 0;
250
      int count          = 1;
251

Quentin Kaci's avatar
Quentin Kaci committed
252
253
      // Build canonized component tree
      for (std::size_t i = 1; i < node_count; ++i)
Quentin Kaci's avatar
Quentin Kaci committed
254
      {
Quentin Kaci's avatar
Quentin Kaci committed
255
256
257
258
259
260
261
262
        if (levels[i] != levels[par[i]]) // Keep the node: Update tree
        {
          translation_map[i] = count++;
          canonized_par.push_back(translation_map[par[i]]);
          canonized_levels.push_back(levels[i]);
        }
        else // Deleted node: retrieve the parent translation
          translation_map[i] = translation_map[par[i]];
Quentin Kaci's avatar
Quentin Kaci committed
263
      }
264

Quentin Kaci's avatar
Quentin Kaci committed
265
      return {canonized_par, canonized_levels};
266
267
    }

268
269
270
271
    template <class I, class N, class F,
              class M = edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>>
    std::pair<component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>, image_ch_value_t<I, int>> //
    __alphatree(I input, N nbh, F distance, std::vector<M>* mst = nullptr, bool canonize_tree = true)
272
    {
273
274
275
      static_assert(mln::is_a<I, mln::details::Image>());
      static_assert(mln::is_a<N, mln::details::Neighborhood>());
      static_assert(::ranges::cpp20::invocable<F, image_value_t<I>, image_value_t<I>>);
276

277
278
279
      using V = image_value_t<I>;
      using P = image_point_t<I>;
      using W = std::invoke_result_t<F, V, V>;
280

281
282
      static_assert(::concepts::totally_ordered<W>);
      static_assert(std::is_same<M, edge_t<P, W>>());
283

284
      // 1. Get the list of edges
285
286
      auto edges = alphatree_edges<P, N, W>();
      internal::alphatree_compute_edges(std::move(input), std::move(nbh), std::move(distance), edges);
287

288
289
290
291
292
293
      std::size_t              flatzones_count;
      image_ch_value_t<I, int> node_map = imchvalue<int>(input).set_init_value(-1);
      {
        image_ch_value_t<I, P> zpar = imchvalue<P>(input);
        // 2. Compute flat zone of the image
        internal::alphatree_compute_flatzones(edges, zpar);
294

295
296
297
        // 3. Compute a node_id for each flat zone
        flatzones_count = internal::alphatree_create_nodemap(node_map, zpar);
      }
298

299
      std::size_t node_count = flatzones_count;
300

301
302
303
304
305
306
307
      std::vector<int> par(node_count);
      std::vector<W>   levels(node_count, 0);
      // 4. Compute the hierarchy
      {
        std::iota(std::begin(par), std::end(par), 0);
        node_count = internal::alphatree_compute_hierarchy(edges, node_map, node_count, par, levels, mst);
      }
308

309
310
      // 5. Parent / levels are ordered from leaves to root, we need to reverse
      internal::alphatree_reorder_nodes(par.data(), levels.data(), node_count);
311

312
313
314
315
316
317
318
319
      // Optional tree canonization: remove useless nodes
      if (canonize_tree)
      {
        auto [canonized_par, canonized_levels] = internal::canonize_component_tree(par, levels, node_count);
        par                                    = canonized_par;
        levels                                 = canonized_levels;
        node_count                             = par.size();
      }
Quentin Kaci's avatar
Quentin Kaci committed
320

321
322
      // 6. Update the node_map
      mln::for_each(node_map, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });
323

324
325
326
      component_tree<W> t;
      t.parent = std::move(par);
      t.values = std::move(levels);
327

328
329
330
331
332
333
334
335
336
      return {std::move(t), std::move(node_map)};
    }
  } // namespace internal

  template <class I, class N, class F>
  std::pair<component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>, image_ch_value_t<I, int>> //
  alphatree(I input, N nbh, F distance)
  {
    return internal::__alphatree(input, nbh, distance);
337
338
  }

339
} // namespace mln::morpho