mesh-complex-max-curv-1-collapse.cc 12.7 KB
Newer Older
1
// Copyright (C) 2008-2013 EPITA Research and Development Laboratory (LRDE)
2
//
3
// This file is part of Olena.
4
//
5
6
7
8
9
// Olena is free software: you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation, version 2 of the License.
//
// Olena is distributed in the hope that it will be useful,
10
11
12
13
14
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
15
// along with Olena.  If not, see <http://www.gnu.org/licenses/>.
16
17
//
// As a special exception, you may use this file as part of a free
18
// software project without restriction.  Specifically, if other files
19
// instantiate templates or use macros or inline functions from this
20
21
22
23
24
// file, or you compile this file and link it with other files to produce
// an executable, this file does not by itself cause the resulting
// executable to be covered by the GNU General Public License.  This
// exception does not however invalidate any other reasons why the
// executable file might be covered by the GNU General Public License.
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

/// \file
/// \brief A program computing the maximal curvature values from the
/// surface of the (triangle) mesh of a statue, then performing a
/// 1-collapse of this surface, using a complex-based image.

// FIXME: Factor with mesh-complex-max-curv-skel.cc and
// mesh-complex-pinv-curv-skel.cc, mesh-complex-max-curv-2-collapse.cc.

#include <iostream>

#include <mln/core/image/complex_image.hh>
#include <mln/core/image/complex_neighborhoods.hh>

#include <mln/core/image/dmorph/image_if.hh>
#include <mln/core/image/dmorph/sub_image.hh>
41
#include <mln/core/routine/extend.hh>
42
43
44
45
46
47
48
#include <mln/core/routine/mutable_extend.hh>
#include <mln/data/paste.hh>

#include <mln/value/label_16.hh>

#include <mln/labeling/regional_minima.hh>
#include <mln/morpho/closing/area.hh>
49
#include <mln/morpho/dilation.hh>
50
51
52
53

#include <mln/topo/is_n_face.hh>
#include <mln/topo/is_simple_pair.hh>
#include <mln/topo/detach_pair.hh>
54
55
56
#include <mln/topo/skeleton/priority_driven_thinning.hh>

#include <mln/arith/revert.hh>
57
58
59
60

#include <mln/io/vtk/load.hh>
#include <mln/io/vtk/save.hh>

61
62
#include <mln/util/timer.hh>

63
#include "misc.hh"
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


int
main(int argc, char* argv[])
{
  if (argc != 4)
    {
      std::cerr << "usage: " << argv[0] << " input.vtk lambda output.vtk"
		<< std::endl;
      std::exit(1);
    }

  std::string input_filename = argv[1];
  unsigned lambda = atoi(argv[2]);
  std::string output_filename = argv[3];

  /*----------------.
  | Complex image.  |
  `----------------*/

  // Curvature image type.
  typedef mln::float_2complex_image3df float_ima_t;
  // Dimension of the image (and therefore of the complex).
  static const unsigned D = float_ima_t::dim;
  // Geometry of the image.
  typedef mln_geom_(float_ima_t) G;

  mln::bin_2complex_image3df bin_input;
  mln::io::vtk::load(bin_input, input_filename);
93
94
95
96
97
  // FIXME: mln::complex_image should provide shorcuts for this.
  const mln::topo::complex<D>& cplx = bin_input.domain().cplx();
  for (unsigned n = 0; n <= D; ++n)
    std::cout << cplx.nfaces_of_dim(n) << ' ' << n << "-faces" << std::endl;

98
99
100
101
102
103
104
105
106
107
108
109
110
  std::pair<float_ima_t, float_ima_t> curv =
    mln::geom::mesh_curvature(bin_input.domain());

  // Compute the max curvature at each vertex.
  float_ima_t float_ima(bin_input.domain());
  mln::p_n_faces_fwd_piter<D, G> v(float_ima.domain(), 0);
  for_all(v)
    {
      // Max curvature.
      float_ima(v) = mln::math::max(mln::math::sqr(curv.first(v)),
				    mln::math::sqr(curv.second(v)));
    }

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
  // Neighborhood type returning the set of (n-1)-faces adjacent to a
  // an n-face.
  typedef mln::complex_lower_neighborhood<D, G> lower_adj_nbh_t;
  lower_adj_nbh_t lower_adj_nbh;

  // Values on edges.
  /* FIXME: We could probably simplify this by using a
     convolution-like operator and morphers (see
     apps/graph-morpho).  */
  mln::p_n_faces_fwd_piter<D, G> e(float_ima.domain(), 1);
  // For each edge (1-face) E, iterate on the the set of vertices
  // (0-faces) adjacent to E.
  mln_niter_(lower_adj_nbh_t) adj_v(lower_adj_nbh, e);
  // Iterate on edges (1-faces).
  for_all(e)
  {
    float s = 0.0f;
    unsigned n = 0;
    // Iterate on vertices (0-faces).
    for_all(adj_v)
    {
      s += float_ima(adj_v);
      ++n;
    }
    float_ima(e) = s / n;
    // An edge should be adjacent to exactly two vertices.
    mln_invariant(n == 2);
  }

140
  // Values on triangles.
141
142
143
  /* FIXME: We could probably simplify this by using a
     convolution-like operator and morphers (see
     apps/graph-morpho).  */
144
  mln::p_n_faces_fwd_piter<D, G> t(float_ima.domain(), 2);
145
146
147
  // For each triangle (2-face) T, iterate on the the set of edges
  // (1-faces) adjacent to T.
  mln_niter_(lower_adj_nbh_t) adj_e(lower_adj_nbh, t);
148
149
150
151
152
  // Iterate on triangles (2-faces).
  for_all(t)
  {
    float s = 0.0f;
    unsigned n = 0;
153
154
    // Iterate on edges (1-faces).
    for_all(adj_e)
155
    {
156
      s += float_ima(adj_e);
157
158
159
      ++n;
    }
    float_ima(t) = s / n;
160
161
    // A triangle should be adjacent to exactly three edges.
    mln_invariant(n == 3);
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  }

  // Convert the float image into an unsigned image because some
  // algorithms do not handle float images well.
  /* FIXME: This is bad: float images should be handled as-is.  At
     least, we should use a decent conversion, using an optimal affine
     transformation (stretching/shrinking) or any other kind of
     interpolation.  */
  typedef mln::unsigned_2complex_image3df ima_t;
  ima_t ima(float_ima.domain());
  // Process only triangles, as edges and vertices are set afterwards
  // (see FIXME below).
  for_all(t)
    ima(t) = 1000 * float_ima(t);

  /*-----------------.
  | Simplification.  |
  `-----------------*/

  /// Adjacent triangles are connected by shared edges.
  typedef mln::complex_lower_dim_connected_n_face_neighborhood<D, G> nbh_t;
  nbh_t nbh;

185
186
187
188
189
190
191
192
193
  // Predicate type: is a face a triangle (2-face)?
  typedef mln::topo::is_n_face<mln_psite_(ima_t), 2> is_a_triangle_t;
  is_a_triangle_t is_a_triangle;

  // Consider only triangles.
  ima_t closed_ima = mln::duplicate(ima);
  mln::data::paste(mln::morpho::closing::area(ima | is_a_triangle,
					      nbh, lambda),
		   closed_ima);
194
195
196
197
198
199
200
201

  /*---------------.
  | Local minima.  |
  `---------------*/

  typedef mln::value::label_16 label_t;
  label_t nminima;

202
  // Consider only triangles.
203
  typedef mln_ch_value_(ima_t, label_t) label_ima_t;
204
205
206
207
208
  label_ima_t minima;
  mln::initialize(minima, closed_ima);
  mln::data::paste(mln::labeling::regional_minima(closed_ima | is_a_triangle,
						  nbh, nminima),
	     minima);
209
210
211
212
213
214
215
216
217
218
219
220
221
222

  /*-----------------------.
  | Initial binary image.  |
  `-----------------------*/

  /* Careful: creating ``holes'' in the surface obviously changes its
     topology, but it may also split a single connected component in
     two or more components, resulting in a disconnected skeleton.  We
     may want to improve this step either by forbidding any splitting,
     or by incrementally ``digging'' a regional minima as long as no
     splitting occurs.  */

  typedef mln_ch_value_(ima_t, bool) bin_ima_t;
  bin_ima_t surface(minima.domain());
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239

  // Predicate type: is a face an edge (1-face)?
  typedef mln::topo::is_n_face<mln_psite_(ima_t), 1> is_an_edge_t;
  is_an_edge_t is_an_edge;
  // Predicate type: is a face a vertex (0-face)?
  typedef mln::topo::is_n_face<mln_psite_(ima_t), 0> is_a_vertex_t;
  is_a_vertex_t is_a_vertex;

  // Neighborhood type returning the set of (n+1)-faces adjacent to a
  // an n-face.
  typedef mln::complex_higher_neighborhood<D, G> higher_adj_nbh_t;
  higher_adj_nbh_t higher_adj_nbh;

  mln::data::fill(surface, false);
  // Set non minima triangles to true;
  mln::data::fill
    ((surface |
240
      (mln::pw::value(minima) == mln::pw::cst(mln::literal::zero))).rw(),
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
     true);
  // Extend non minima values from triangles to edges.
  mln::data::paste (mln::morpho::dilation(mln::extend(surface | is_an_edge,
						      surface),
					  /* Dilations require windows,
					     not neighborhoods.  */
					  higher_adj_nbh.win()),
		    surface);
  // Extend non minima values from edges to vertices.
  mln::data::paste(mln::morpho::dilation(mln::extend(surface | is_a_vertex,
						     surface),
					  /* Dilations require windows,
					     not neighborhoods.  */
					 higher_adj_nbh.win()),
		   surface);
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
291
292
293
294
295
296

  /*-------------.
  | 2-collapse.  |
  `-------------*/

  // ------------------------------- //
  // Image restricted to triangles.  //
  // ------------------------------- //

  // Surface image type, of which domain is restricted to triangles.
  typedef mln::image_if<bin_ima_t, is_a_triangle_t> bin_triangle_only_ima_t;
  // Surface image type, of which iteration (not domain) is restricted
  // to triangles.
  typedef mln::mutable_extension_ima<bin_triangle_only_ima_t, bin_ima_t>
    bin_triangle_ima_t;

  // ------------------------ //
  // Simple point predicate.  //
  // ------------------------ //

  // Predicate type: is a triangle (2-face) simple?
  typedef mln::topo::is_simple_pair< bin_triangle_ima_t,
                                     lower_adj_nbh_t,
                                     higher_adj_nbh_t >
    is_simple_triangle_t;
  is_simple_triangle_t is_simple_triangle(lower_adj_nbh, higher_adj_nbh);

  // ------------------------------- //
  // Simple point detach procedure.  //
  // ------------------------------- //

  // Functor detaching a cell.
  typedef mln::topo::detach_pair< bin_triangle_ima_t,
                                  lower_adj_nbh_t,
                                  higher_adj_nbh_t > detach_triangle_t;
  detach_triangle_t detach_triangle(lower_adj_nbh, higher_adj_nbh);

  // ------------------------ //
  // Thinning by 2-collapse.  //
  // ------------------------ //

297
298
299
300
  // Create a priority function (actually, an image) using the inverse
  // of the curvature image.
  ima_t priority = mln::arith::revert(closed_ima);

301
302
303
  mln_concrete_(bin_ima_t) surface_2_collapse;
  mln::initialize(surface_2_collapse, surface);
  mln::data::paste
304
    (mln::topo::skeleton::priority_driven_thinning
305
306
307
     (mln::mutable_extend((surface | is_a_triangle).rw(), surface),
      nbh,
      is_simple_triangle,
308
309
      detach_triangle,
      priority)
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
     /* Before pasting the result of the computation into
        SURFACE_2_COLLAPSE, re-expand its domain to the initial site
        set, to ensure data from all faces (i.e., both the 2-faces,
        directly processed; and the 1-faces from the extension,
        undirectly processed).  */
     | surface.domain(),
     surface_2_collapse);

  /*-------------.
  | 1-collapse.  |
  `-------------*/

  // --------------------------- //
  // Image restricted to edges.  //
  // --------------------------- //

  // Surface image type, of which domain is restricted to edges.
  typedef mln::image_if<bin_ima_t, is_an_edge_t> bin_edge_only_ima_t;
  // Surface image type, of which iteration (not domain) is restricted
  // to edges.
  typedef mln::mutable_extension_ima<bin_edge_only_ima_t, bin_ima_t>
    bin_edge_ima_t;

  // ------------------------ //
  // Simple point predicate.  //
  // ------------------------ //

  // Predicate type: is a edge (1-face) simple?
  typedef mln::topo::is_simple_pair< bin_edge_ima_t,
                                     lower_adj_nbh_t,
                                     higher_adj_nbh_t >
    is_simple_edge_t;
  is_simple_edge_t is_simple_edge(lower_adj_nbh, higher_adj_nbh);

  // ------------------------------- //
  // Simple point detach procedure.  //
  // ------------------------------- //

  // Functor detaching a cell.
  typedef mln::topo::detach_pair< bin_edge_ima_t,
                                  lower_adj_nbh_t,
                                  higher_adj_nbh_t > detach_edge_t;
  detach_edge_t detach_edge(lower_adj_nbh, higher_adj_nbh);

  // ------------------------ //
  // Thinning by 1-collapse.  //
  // ------------------------ //

  mln_concrete_(bin_ima_t) surface_1_collapse;
  mln::initialize(surface_1_collapse, surface_2_collapse);
360
361
362

  mln::util::timer time;
  time.start();
363
  mln::data::paste
364
    (mln::topo::skeleton::priority_driven_thinning
365
366
367
368
     (mln::mutable_extend((surface_2_collapse | is_an_edge).rw(),
			  surface_2_collapse),
      nbh,
      is_simple_edge,
369
370
      detach_edge,
      priority)
371
372
373
374
375
376
377
     /* Likewise, before pasting the result of the computation into
        SURFACE_1_COLLAPSE, re-expand its domain to the initial site
        set, to ensure data from all faces (i.e., both the 1-faces,
        directly processed; and the 0-faces from the extension,
        undirectly processed).  */
     | surface_2_collapse.domain(),
     surface_1_collapse);
378
379
  time.stop();
  std::cout << time.read() << " s" << std::endl;
380
381
382
383
384
385
386

  /*---------.
  | Output.  |
  `---------*/

  mln::io::vtk::save(surface_1_collapse, output_filename);
}