mtos_example.cpp 4.25 KB
Newer Older
1
#include <mln/accu/accumulators/count.hpp>
2
3
#include <mln/accu/accumulators/max.hpp>
#include <mln/core/algorithm/accumulate.hpp>
4
#include <mln/core/colors.hpp>
5
#include <mln/core/extension/padding.hpp>
6
7
8
9
10
11
12
#include <mln/core/image/ndimage.hpp>
#include <mln/core/image/view/channel.hpp>
#include <mln/core/image/view/transform.hpp>
#include <mln/io/imread.hpp>
#include <mln/io/imsave.hpp>
#include <mln/morpho/tos.hpp>

13
#include <mln/morpho/private/satmaxtree.hpp>
14
15
16
17
18
19
20
21
22
#include <mln/morpho/private/trees_fusion.hpp>

#include <iostream>
#include <vector>

#include "lut.hpp"

namespace
{
23
24
25
  struct mean_node_accu : mln::Accumulator<mean_node_accu>
  {
    using result_type = decltype(std::declval<mln::rgb8>() + std::declval<mln::rgb8>());
26

27
  public:
Baptiste Esteban's avatar
Baptiste Esteban committed
28
    void take(const mln::rgb8& v)
29
    {
Baptiste Esteban's avatar
Baptiste Esteban committed
30
      m_sum += v;
31
32
33
34
35
36
37
38
      m_count++;
    }

    void take(const mean_node_accu&) {}

    result_type to_result() const { return m_count > 1 ? static_cast<result_type>(m_sum / m_count) : m_sum; }

  private:
39
40
    result_type m_sum   = {0, 0, 0};
    int         m_count = 0;
41
42
  };

43
  mln::rgb8 get_median_on_border(mln::image2d<mln::rgb8> ima)
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
  {
    std::vector<mln::rgb8> border;
    border.reserve(2 * ima.width() + 2 * ima.height() - 4);

    for (int i = 0; i < ima.width(); i++)
    {
      border.push_back(ima(mln::point2d{i, 0}));
      border.push_back(ima(mln::point2d{i, ima.height() - 1}));
    }

    for (int i = 1; i < ima.height() - 1; i++)
    {
      border.push_back(ima(mln::point2d{0, i}));
      border.push_back(ima(mln::point2d{ima.width() - 1, i}));
    }

    std::sort(border.begin(), border.end(), [](const mln::rgb8& a, const mln::rgb8& b) {
61
62
63
64
      int i = 0;
      while (i < 3 && a[i] == b[i])
        i++;
      return i < 3 && a[i] < b[i];
65
    });
66
    return border[border.size() / 2];
67
68
  }

69
70
71
  /// \brief Reduce the size of a nodemap by a factor 2
  mln::image2d<int> reduce_nodemap(mln::image2d<int> n)
  {
72
73
74
75
    // mln::image2d<int> res((n.width() + 3) / 4, (n.height() + 3) / 4);
    const mln::point2d pmin(n.domain().tl());
    const mln::point2d pmax(n.domain().br());
    mln::image2d<int>  res(mln::box2d{{(pmin[0] - 3) / 4, (pmin[1] - 3) / 4}, {(pmax[0] + 3) / 4, (pmax[1] + 3) / 4}});
76

77
    mln_foreach (auto p, n.domain())
78
    {
79
      if (p[0] % 4 == 0 && p[1] % 4 == 0)
80
      {
81
        res(mln::point2d{p[0] / 4, p[1] / 4}) = n(p);
82
83
84
85
86
      }
    }

    return res;
  }
87
88
89
90
} // namespace

int main(int argc, char* argv[])
{
91
  if (argc < 4)
92
  {
93
94
    std::cerr << "Invalid number of arguments\nUsage: " << argv[0]
              << " input_filename depth_map_filename rec_filename\n";
95
96
97
98
99
    return 1;
  }

  mln::image2d<mln::rgb8> ima;
  mln::io::imread(argv[1], ima);
100
101
102
103
104

  const auto med = get_median_on_border(ima);
  ima.inflate_domain(1);
  constexpr int borders[2][2] = {{1, 1}, {1, 1}};
  mln::pad(ima, mln::PAD_CONSTANT, borders, med);
105
106
107
108
109
110
111

  mln::morpho::component_tree<> trees[3];
  mln::image2d<int>             nodemaps[3];
  std::vector<int>              depths[3];

  for (int c = 0; c < 3; c++)
  {
112
    std::tie(trees[c], nodemaps[c]) = mln::morpho::tos(mln::view::channel(ima, c), {-1, -1});
113
114
115
116
117
    depths[c]                       = trees[c].compute_depth();
  }

  const auto [gos, tree_to_graph] = mln::morpho::details::compute_inclusion_graph(trees, nodemaps, depths, 3);
  auto depth_map                  = mln::morpho::details::compute_depth_map(gos, tree_to_graph, nodemaps);
118
  {
119
    std::uint16_t max_depth = mln::accumulate(depth_map, mln::accu::accumulators::max<std::uint16_t>());
120
121
122
123
124
    auto          normalized_depth =
        mln::view::transform(depth_map, [&max_depth](std::uint16_t a) -> float { return (float)a / (float)max_depth; });
    auto heat_depth = mln::view::transform(normalized_depth, heat_lut);
    mln::io::imsave(heat_depth, argv[2]);
  }
125

126
127
128
129
130
  auto [t, nm] = mln::morpho::details::satmaxtree(depth_map, {-1, -1});
  nm           = reduce_nodemap(nm);
  ima.inflate_domain(-1);
  nm.inflate_domain(-1);

131
  auto area = t.compute_attribute_on_points(nm, mln::accu::accumulators::count<int>());
Baptiste Esteban's avatar
Baptiste Esteban committed
132
  t.filter(mln::morpho::CT_FILTER_DIRECT, nm, [&area](int n) { return area[n] >= 100; });
Baptiste Esteban's avatar
Baptiste Esteban committed
133
  auto mean = t.compute_attribute_on_values(nm, ima, mean_node_accu());
134
  auto rec  = t.reconstruct_from(nm, ranges::make_span(mean.data(), mean.size()));
135
  mln::io::imsave(mln::view::cast<mln::rgb8>(rec), argv[3]);
136
137
138

  return 0;
}