1 // Copyright (C) 2014 Carnë Draug <carandraug@octave.org>
2 //
3 // This program is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU General Public License as
5 // published by the Free Software Foundation; either version 3 of the
6 // License, or (at your option) any later version.
7 //
8 // This program is distributed in the hope that it will be useful, but
9 // WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with this program; if not, see
15 // <http://www.gnu.org/licenses/>.
16
17 // This file implements imreconstruct as described on "Vincent, L. (1993).
18 // Morphological grayscale reconstruction in image analysis: applications
19 // and efficient algorithms. Image Processing, IEEE Transactions on, 2(2),
20 // 176-201."
21 //
22 // Our strategy to handle elements in the border is to simply pad it with
23 // the lowest value for the type, which will be ignored on the comparisons.
24 // This should still be more efficient than using subscript indices to find
25 // when we are on the border.
26
27 #include <functional>
28 #include <queue>
29
30 #include <octave/lo-mappers.h>
31
32 #include <octave/defun-dld.h>
33 #include <octave/defun-int.h>
34 #include <octave/error.h>
35 #include <octave/ovl.h>
36
37 #include "connectivity.h"
38
39 using namespace octave::image;
40
41 #define WANTS_OCTAVE_IMAGE_VALUE 1
42 #include "octave-wrappers.h"
43
44 /*
45 ## A dirty implementation of the fast hybrid reconstruction as m file
46 ## for testing purposes.
47
48 function marker = fast_hybrid_reconstruction (marker, mask)
49 ## works for 10x10 matrices, padded to 12x12 with zeros, when
50 ## connectivity is ones (3, 3)
51
52 offsets = [-13 -12 -11 -1 1 11 12 13];
53 pos_offsets = [0 1 11 12 13]; # don't forget the zero
54
55 neg_offsets = - [pos_offsets];
56
57 ## the raster scan
58 for c = 2:(columns(marker) -1)
59 for r = 2:(rows(marker) -1)
60 i = sub2ind (size (marker), r, c);
61 marker(r,c) = min (mask(r,c), max (marker(i + [neg_offsets])));
62 endfor
63 endfor
64
65 ## the antiraster scan
66 fifo = [];
67 for c = (columns(marker) -1):-1:2
68 for r = (rows(marker) -1):-1:2
69 i = sub2ind (size (marker), r, c);
70 offs = marker(i + [pos_offsets]);
71 marker(r,c) = min (mask(r,c), max (offs));
72
73 offs(1) = []; #remove itself
74 picks = offs < marker(i) & offs < mask(i + [pos_offsets(2:end)]);
75 if (any (picks))
76 fifo(end+1) = i;
77 endif
78 endfor
79 endfor
80
81 ## the propagation step
82 while (numel (fifo) != 0)
83 p = fifo(1);
84 fifo(1) = [];
85
86 for i = offsets;
87 if (marker(p +i) < marker(p) && mask(p+i) != marker(p+i))
88 marker(p +i) = min (marker(p), mask(p+i));
89 fifo(end+1) = p+i;
90 endif
91 endfor
92 endwhile
93
94 endfunction
95 */
96
97 template<class T>
98 static void
scan_raster_order(T & padded_marker,const T & padded_mask,const dim_vector & original_size,const Array<octave_idx_type> & padding_lengths,const Array<octave_idx_type> & raster_neighbours)99 scan_raster_order (T& padded_marker, const T& padded_mask,
100 const dim_vector& original_size,
101 const Array<octave_idx_type>& padding_lengths,
102 const Array<octave_idx_type>& raster_neighbours)
103 {
104 typedef typename T::element_type P;
105
106 P* J = padded_marker.fortran_vec ();
107 const P* I = padded_mask.fortran_vec ();
108 const octave_idx_type* pads = padding_lengths.fortran_vec ();
109 const octave_idx_type* neighbours = raster_neighbours.fortran_vec ();
110
111 const octave_idx_type n_neighbours = raster_neighbours.numel ();
112
113 // We probably should not be using this but converting to Array
114 // just to have fortran_vec seems a bit too much.
115 const octave_idx_type* s = original_size.to_jit ();
116
117 std::function<void(const octave_idx_type)> scan;
118 scan = [&] (const octave_idx_type dim) -> void
119 {
120 J += pads[dim];
121 I += pads[dim];
122
123 if (dim == 0)
124 {
125 for (octave_idx_type k = 0; k < s[0]; k++, J++, I++)
126 {
127 for (octave_idx_type i = 0; i < n_neighbours; i++)
128 if (*J < J[neighbours[i]])
129 *J = J[neighbours[i]];
130
131 if (*J > *I)
132 *J = *I;
133 }
134 }
135 else
136 for (octave_idx_type i = 0; i < s[dim]; i++)
137 scan (dim-1);
138
139 J += pads[dim];
140 I += pads[dim];
141 return;
142 };
143 scan (original_size.length () -1);
144 return;
145 }
146
147 template<class T>
148 static std::queue<octave_idx_type>
scan_antiraster_order(T & padded_marker,const T & padded_mask,const dim_vector & original_size,const Array<octave_idx_type> & padding_lengths,const Array<octave_idx_type> & antiraster_neighbours)149 scan_antiraster_order (T& padded_marker, const T& padded_mask,
150 const dim_vector& original_size,
151 const Array<octave_idx_type>& padding_lengths,
152 const Array<octave_idx_type>& antiraster_neighbours)
153 {
154 typedef typename T::element_type P;
155 std::queue<octave_idx_type> unfinished;
156
157 P* J = padded_marker.fortran_vec ();
158 const P* I = padded_mask.fortran_vec ();
159 const octave_idx_type* pads = padding_lengths.fortran_vec ();
160 const octave_idx_type* neighbours = antiraster_neighbours.fortran_vec ();
161
162 const octave_idx_type n_neighbours = antiraster_neighbours.numel ();
163
164 // We probably should not be using this but converting to Array
165 // just to have fortran_vec seems a bit too much.
166 const octave_idx_type* s = original_size.to_jit ();
167
168 J += padded_marker.numel () -1;
169 I += padded_marker.numel () -1;
170
171 octave_idx_type ind = padded_marker.numel () -1;
172 std::function<void(const octave_idx_type)> scan;
173 scan = [&] (const octave_idx_type dim) -> void
174 {
175 J -= pads[dim];
176 I -= pads[dim];
177 ind -= pads[dim];
178
179 if (dim == 0)
180 {
181 for (octave_idx_type k = 0; k < s[0]; k++, J--, I--, ind--)
182 {
183 for (octave_idx_type i = 0; i < n_neighbours; i++)
184 if (*J < J[neighbours[i]])
185 *J = J[neighbours[i]];
186
187 if (*J > *I)
188 *J = *I;
189
190 for (octave_idx_type i = 0; i < n_neighbours; i++)
191 if (J[neighbours[i]] < *J && J[neighbours[i]] < I[neighbours[i]])
192 unfinished.push (ind);
193 }
194 }
195 else
196 for (octave_idx_type i = 0; i < s[dim]; i++)
197 scan (dim-1);
198
199 J -= pads[dim];
200 I -= pads[dim];
201 ind -= pads[dim];
202
203 return;
204 };
205 scan (original_size.length () -1);
206 return unfinished;
207 }
208
209 template<class T>
210 static void
propagation_step(T & padded_marker,const T & padded_mask,std::queue<octave_idx_type> & unfinished,const Array<octave_idx_type> & deleted_neighbours)211 propagation_step (T& padded_marker, const T& padded_mask,
212 std::queue<octave_idx_type>& unfinished,
213 const Array<octave_idx_type>& deleted_neighbours)
214 {
215 typedef typename T::element_type P;
216
217 P* J = padded_marker.fortran_vec ();
218 const P* I = padded_mask.fortran_vec ();
219 const octave_idx_type* neighbours = deleted_neighbours.fortran_vec ();
220
221 const octave_idx_type n_neighbours = deleted_neighbours.numel ();
222
223 while (! unfinished.empty ())
224 {
225 octave_idx_type p = unfinished.front ();
226 unfinished.pop ();
227 for (octave_idx_type k = 0; k < n_neighbours; k++)
228 {
229 octave_idx_type q = p + neighbours[k];
230 if (J[q] < J[p] && I[q] != J[q])
231 {
232 J[q] = octave::math::min (J[p], I[q]);
233 unfinished.push (q);
234 }
235 }
236 OCTAVE_QUIT;
237 }
238 return;
239 }
240
241 template<class T>
242 static T
fast_hybrid_reconstruction(const T & marker,const T & mask,const connectivity & conn)243 fast_hybrid_reconstruction (const T& marker, const T& mask,
244 const connectivity& conn)
245 {
246 typedef typename T::element_type P;
247
248 const dim_vector original_size = marker.dims ();
249
250 T padded_marker = conn.create_padded (marker, connectivity::min_value<P> ());
251 const T padded_mask = conn.create_padded (mask, connectivity::min_value<P> ());
252
253 const dim_vector padded_size = padded_marker.dims ();
254 const Array<octave_idx_type> padding_lengths
255 = connectivity::padding_lengths (original_size, padded_size);
256
257 scan_raster_order (padded_marker, padded_mask, original_size,
258 padding_lengths,
259 conn.negative_neighbourhood (padded_size));
260
261 OCTAVE_QUIT;
262
263 std::queue<octave_idx_type> unfinished =
264 scan_antiraster_order (padded_marker, padded_mask, original_size,
265 padding_lengths,
266 conn.positive_neighbourhood (padded_size));
267
268 OCTAVE_QUIT;
269
270 propagation_step (padded_marker, padded_mask, unfinished,
271 conn.deleted_neighbourhood (padded_size));
272
273 conn.unpad (padded_marker);
274 return padded_marker;
275 }
276
277 template<class T>
278 static T
reconstruct(const T & marker,const T & mask,const connectivity & conn)279 reconstruct (const T& marker, const T& mask, const connectivity& conn)
280 {
281 return fast_hybrid_reconstruction (marker, mask, conn);
282 }
283
284 // TODO implement the following by reusing the code in bwlabeln
285 //static boolNDArray
286 //reconstruct (const boolNDArray& marker, const boolNDArray& mask,
287 // const connectivity& conn)
288 //{
289 // /*
290 // 1. Label the connected components of the mask image, i.e., each of these
291 // components is assigned a unique number. Note that this step can itself
292 // be implemented very efficiently by using algorithms based on chain an
293 // loops [16] or queues of pixels [23, 26].
294 // 2. Determine the labels of the connected components which contain at
295 // least a pixel of the marker image.
296 // 3. Remove all the connected components whose label is not one of the
297 // previous ones.
298 // */
299 // return boolNDArray ();
300 //}
301
302 DEFUN_DLD(imreconstruct, args, , "\
303 -*- texinfo -*-\n\
304 @deftypefn {Loadable Function} {} imreconstruct (@var{marker}, @var{mask})\n\
305 @deftypefnx {Loadable Function} {} imreconstruct (@var{marker}, @var{mask}, @var{conn})\n\
306 \n\
307 @seealso{imclearborder, imdilate, imerode}\n\
308 @end deftypefn")
309 {
310 const octave_idx_type nargin = args.length ();
311
312 if (nargin < 2 || nargin > 3)
313 print_usage ();
314 if (args(0).class_name () != args(1).class_name ())
315 error ("imreconstruct: MARKER and MASK must be of same class");
316
317 connectivity conn;
318 if (nargin > 2)
319 conn = conndef (args(2));
320 else
321 {
322 try
323 {
324 conn = connectivity (args(0).ndims (), "maximal");
325 }
326 catch (invalid_connectivity& e)
327 {
328 error ("imreconstruct: unable to create connectivity (%s)", e.what ());
329 }
330 }
331 octave_image::value marker (args(0));
332
333 #define RECONSTRUCT(TYPE) \
334 ret = reconstruct (marker.TYPE ## _array_value (), \
335 args(1).TYPE ## _array_value (), conn);
336
337 #define IF_TYPE(TYPE) \
338 if (marker.is_ ## TYPE ## _type ()) \
339 RECONSTRUCT (TYPE)
340
341 #define INT_BRANCH(TYPE) \
342 IF_TYPE(u ## TYPE) \
343 else IF_TYPE(TYPE)
344
345 #define FLOAT_BRANCH(CR) \
346 if (marker.is_single_type ()) \
347 ret = reconstruct (marker.float_ ## CR ## array_value (), \
348 args(1).float_ ## CR ## array_value (), conn); \
349 else \
350 ret = reconstruct (marker.CR ## array_value (), \
351 args(1).CR ## array_value (), conn);
352
353 octave_value ret;
354 if (marker.islogical ())
355 RECONSTRUCT(bool)
356 else INT_BRANCH (int8)
357 else INT_BRANCH (int16)
358 else INT_BRANCH (int32)
359 else INT_BRANCH (int64)
360 else if (marker.isreal ())
361 {
362 FLOAT_BRANCH()
363 }
364 else if (marker.iscomplex ())
365 {
366 FLOAT_BRANCH(complex_)
367 }
368 else
369 error ("imreconstruct: unsupported class %s for MARKER",
370 marker.class_name ().c_str ());
371
372 #undef IF_TYPE
373 #undef INT_BRANCH
374 #undef FLOAT_BRANCH
375
376 return ret;
377 }
378
379 /*
380 ## When using the fast hybrid reconstruction (and specially with random
381 ## images), and if the images are small, it is often finished after the
382 ## antiraster scan and before the propagation step. Using larger images
383 ## makes sure we get in the propagation step and that we catch bugs in there.
384
385 ## This function does exactly what imreconstruct is meant to but is, in
386 ## the words of Luc Vicent 1993, and I can attest to it, "[...] not suited
387 ## to conventional computers, where its execution time is often of several
388 ## minutes."
389 %!function recon = parallel_reconstruction (marker, mask,
390 %! conn = conndef (ndims (marker), "maximal"))
391 %! do
392 %! previous = marker;
393 %! marker = imdilate (marker, conn);
394 %! ## FIXME https://savannah.gnu.org/bugs/index.php?43712
395 %! if (strcmp (class (marker), "logical"))
396 %! marker = marker & mask;
397 %! else
398 %! marker = min (marker, mask);
399 %! endif
400 %! until (all ((marker == previous)(:)))
401 %! recon = marker;
402 %!endfunction
403
404 %!test
405 %! for cl = {"int8", "uint8", "int16", "uint16", "int32", "uint32"}
406 %! cl = cl{1};
407 %! a = randi ([intmin(cl) intmax(cl)-30], 100, 100, cl);
408 %! b = a + randi (20, 100, 100, cl);
409 %! assert (imreconstruct (a, b), parallel_reconstruction (a, b))
410 %! endfor
411 %! for cl = {"double", "single"}
412 %! cl = cl{1};
413 %! a = (rand (100, 100, cl) - 0.5) .* 1000;
414 %! b = a + rand (100, 100, cl) * 100;
415 %! assert (imreconstruct (a, b), parallel_reconstruction (a, b))
416 %! endfor
417
418 %!test
419 %! for cl = {"int8", "uint8", "int16", "uint16", "int32", "uint32"}
420 %! cl = cl{1};
421 %! a = randi ([intmin(cl) intmax(cl)-30], 100, 100, cl);
422 %! b = a + randi (20, 100, 100, cl);
423 %! c = [0 1 0; 1 1 1; 0 1 0];
424 %! assert (imreconstruct (a, b, c), parallel_reconstruction (a, b, c))
425 %! endfor
426
427 %!test
428 %! a = randi (210, 100, 100);
429 %! b = a + randi (20, 100, 100);
430 %! c = ones (3, 1);
431 %! assert (imreconstruct (a, b, c), parallel_reconstruction (a, b, c))
432
433 %!test
434 %! a = randi (210, 500, 500, 10, 4);
435 %! b = a + randi (20, 500, 500, 10, 4);
436 %! c = ones (3, 3, 3);
437 %! assert (imreconstruct (a, b, c), parallel_reconstruction (a, b, c))
438
439 %!test
440 %! a = randi (210, 500, 500, 10, 4);
441 %! b = a + randi (20, 500, 500, 10, 4);
442 %! c = conndef (4, "minimal");
443 %! assert (imreconstruct (a, b, c), parallel_reconstruction (a, b, c))
444
445 %!test
446 %! a = [ 0 0 0 0 0 0 0 1 0 0
447 %! 0 0 0 0 0 0 0 1 0 0
448 %! 1 0 0 0 0 0 0 0 0 0
449 %! 0 0 0 0 0 0 0 0 0 0
450 %! 0 0 0 0 0 0 0 1 0 0
451 %! 0 0 0 0 0 0 1 0 0 0
452 %! 0 0 0 0 0 0 0 0 0 0
453 %! 0 0 0 0 0 0 0 0 0 0
454 %! 0 0 0 0 1 0 0 0 0 0
455 %! 0 0 0 0 0 0 0 1 0 0];
456 %!
457 %! b = [ 0 1 0 0 0 0 0 1 1 0
458 %! 1 1 0 0 0 1 0 1 1 0
459 %! 1 1 0 0 1 0 0 0 0 0
460 %! 1 1 0 0 0 1 1 0 0 0
461 %! 1 0 0 0 0 0 1 1 0 0
462 %! 0 1 0 0 0 0 1 1 0 0
463 %! 0 0 0 1 0 0 0 0 0 0
464 %! 0 0 0 0 1 1 0 0 0 0
465 %! 0 0 0 1 1 0 0 0 0 0
466 %! 1 0 0 0 1 0 0 1 0 1];
467 %!
468 %! c = [ 0 1 0 0 0 0 0 1 1 0
469 %! 1 1 0 0 0 1 0 1 1 0
470 %! 1 1 0 0 1 0 0 0 0 0
471 %! 1 1 0 0 0 1 1 0 0 0
472 %! 1 0 0 0 0 0 1 1 0 0
473 %! 0 1 0 0 0 0 1 1 0 0
474 %! 0 0 0 1 0 0 0 0 0 0
475 %! 0 0 0 0 1 1 0 0 0 0
476 %! 0 0 0 1 1 0 0 0 0 0
477 %! 0 0 0 0 1 0 0 1 0 0];
478 %! assert (imreconstruct (logical (a), logical (b)), logical (c));
479 %!
480 %! c = [ 0 1 0 0 0 0 0 1 1 0
481 %! 1 1 0 0 0 0 0 1 1 0
482 %! 1 1 0 0 0 0 0 0 0 0
483 %! 1 1 0 0 0 1 1 0 0 0
484 %! 1 0 0 0 0 0 1 1 0 0
485 %! 0 0 0 0 0 0 1 1 0 0
486 %! 0 0 0 0 0 0 0 0 0 0
487 %! 0 0 0 0 1 1 0 0 0 0
488 %! 0 0 0 1 1 0 0 0 0 0
489 %! 0 0 0 0 1 0 0 1 0 0];
490 %! assert (imreconstruct (logical (a), logical (b), [0 1 0; 1 1 1; 0 1 0]),
491 %! logical (c));
492
493 %!test
494 %! do
495 %! b = rand (100, 100, 100) > 0.98;
496 %! until (nnz (b) > 4)
497 %! b = imdilate (b, ones (5, 5, 5));
498 %! a = false (size (b));
499 %! f = find (b);
500 %! a(f(randi (numel (f), 6, 1))) = true;
501 %! assert (imreconstruct (a, b), parallel_reconstruction (a, b))
502
503 ## we try to be smart about the padding so make sure this works. There
504 ## was a nasty bug during development which this test brings up.
505 %!test
506 %! a = randi (200, 100,100, 10, 10);
507 %! b = a + randi (20, 100,100, 10, 10);
508 %! c1 = ones (3, 3, 3);
509 %! c2 = zeros (3, 3, 3, 3);
510 %! c2(:,:,:,2) = c1;
511 %! assert (imreconstruct (a, b, c1), imreconstruct (a, b, c2))
512
513 %!test
514 %! ## Values in MARKER above MASK should be clipped (bug #48794)
515 %! ## (well, treated internally as if they were clipped)
516 %! mask = logical ([1 1 1; 1 0 1; 1 1 1]);
517 %! assert (imreconstruct (true (3, 3), mask), mask)
518 %!
519 %! mask = ones (5, 5);
520 %! mask(2:4,2:4) = 0;
521 %! assert (imreconstruct (ones (5, 5), mask), mask)
522 %!
523 %! mask = ones (5, 5);
524 %! mask(2:4,2:4) = 0;
525 %! assert (imreconstruct (repmat (2, [5, 5]), mask), mask)
526 %!
527 %! mask = ones (5, 5);
528 %! mask(2:4,2:4) = 0;
529 %! assert (imreconstruct (repmat (2, [5, 5]), mask), mask)
530 %!
531 %! marker = ones (3, 3, 3, 3);
532 %! mask = marker;
533 %! mask(2, 2, 2, 2) = 0;
534 %! assert (imreconstruct (marker, mask), mask)
535 %!
536 %! marker = randi (210, 100, 100);
537 %! assert (imreconstruct (marker +1, marker), marker)
538 %! assert (imreconstruct (marker +1, marker), imreconstruct (marker, marker))
539 */
540