1 // Copyright (C) 2015 Carnë Draug <carandraug@octave.org>
2 //
3 // This program is free software; you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation; either version 3 of the License, or
6 // (at your option) any later version.
7 //
8 // This program is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU 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 <http://www.gnu.org/licenses/>.
15
16 #include <queue>
17
18 #include <octave/oct.h>
19 #include <octave/error.h>
20 #include <octave/parse.h>
21 #include <octave/Cell.h>
22 #include <octave/oct-map.h>
23
24 #define WANTS_FEVAL 1
25 #define WANTS_OCTAVE_IMAGE_VALUE 1
26 #include "octave-wrappers.h"
27
28 #include "connectivity.h"
29 using namespace octave::image;
30
31 template<class T>
32 static boolNDArray
imregionalmin(const T & im,const connectivity & conn)33 imregionalmin (const T& im, const connectivity& conn)
34 {
35 octave_value_list args (2);
36 args(0) = im;
37 args(1) = conn.mask;
38 const octave_value regional_min = octave_image::feval ("imregionalmin",
39 args)(0);
40 return regional_min.bool_array_value ();
41 }
42
43 static NDArray
bwlabeln(const boolNDArray & bw,const connectivity & conn)44 bwlabeln (const boolNDArray& bw, const connectivity& conn)
45 {
46 octave_value_list args (2);
47 args(0) = bw;
48 args(1) = conn.mask;
49 const octave_value label = octave_image::feval ("bwlabeln", args)(0);
50 return label.array_value ();
51 }
52
53 // Implements watershed in a quite naïve way. From the wikipedia, named
54 // "Meyer's flooding algorithm" (but I could not find the actual paper
55 // that reports it). There are faster (and also nicer results) algorithms,
56 // but this is the only one I found that matches Matlab results.
57 //
58 // 1. A set of markers, pixels where the flooding shall start, are chosen.
59 // Each is given a different label.
60 // 2. The neighboring pixels of each marked area are inserted into a
61 // priority queue with a priority level corresponding to the gray level
62 // of the pixel.
63 // 3. The pixel with the lowest priority level is extracted from the
64 // priority queue. If the neighbors of the extracted pixel that have
65 // already been labeled all have the same label, then the pixel is
66 // labeled with their label. All non-marked neighbors that are not yet
67 // in the priority queue are put into the priority queue.
68 // 4. Redo step 3 until the priority queue is empty.
69 //
70 // There is a detail missing on the description above. On step 3, if the
71 // labeled neighbours do *not* have the same label, should the non-labeled
72 // neighbours be added to the queue? Apparently not.
73
74 template<class P>
75 class
76 Voxel
77 {
78 public:
79 P val;
80 octave_idx_type idx;
81 // We need this to sort elements with the same priority. We need them
82 // to come out in the same order they went in.
83 octave_idx_type pos;
84
Voxel(const P val,const octave_idx_type idx,const octave_idx_type pos)85 Voxel (const P val, const octave_idx_type idx, const octave_idx_type pos)
86 : val (val), idx (idx), pos (pos)
87 { }
88
89 inline bool
operator >(const Voxel & rhs) const90 operator>(const Voxel& rhs) const
91 {
92 if (val == rhs.val)
93 return pos > rhs.pos;
94 else
95 return val > rhs.val;
96 }
97 };
98
99 // As part of this algorithm, we will check the neighbourhood for existing
100 // labels. We don't know in advance the number of labeled neighbours, or
101 // where the first label will be. But we do know the length of the
102 // neighbourhood.
103 template<class T>
104 class
105 Collection
106 {
107 public:
Collection(const octave_idx_type n)108 explicit Collection (const octave_idx_type n) : data (new T[n])
109 { }
110
~Collection(void)111 ~Collection (void)
112 { delete [] data; }
113
114 inline octave_idx_type
numel(void) const115 numel (void) const
116 { return count; }
117
118 inline void
push_back(const T val)119 push_back (const T val)
120 { data[count++] = val; }
121
122 inline void
reset(void)123 reset (void)
124 { count = 0; }
125
126 protected:
127 T* data = NULL;
128 octave_idx_type count = 0;
129
130 private:
131 // Disable default and copy constructor and assignment
132 Collection (void);
133 Collection (Collection const& other);
134 Collection& operator = (Collection const& other);
135 };
136
137
138 class
139 LabelsCollection : public Collection<double>
140 {
141 public:
142 using Collection<double>::Collection;
143
144 inline double
label(void) const145 label (void) const
146 { return *data; }
147
148 inline bool
all_equal(void) const149 all_equal (void) const
150 {
151 for (octave_idx_type i = 0; i < count; i++)
152 if (data[0] != data[i])
153 return false;
154 return true;
155 }
156 };
157
158 class
159 IdxCollection : public Collection<octave_idx_type>
160 {
161 public:
162 using Collection<octave_idx_type>::Collection;
163
164 inline octave_idx_type
operator [](octave_idx_type i) const165 operator [] (octave_idx_type i) const
166 { return data[i]; }
167 };
168
169
170 template<class T>
171 NDArray
watershed(const T & im,const connectivity & conn)172 watershed (const T& im, const connectivity& conn)
173 {
174 typedef typename T::element_type P;
175
176 // 1. A set of markers, pixels where the flooding shall start, are chosen.
177 // Each is given a different label.
178 const boolNDArray markers = imregionalmin (im, conn);
179 boolNDArray padded_markers = conn.create_padded (markers, false);
180 NDArray label_array = bwlabeln (padded_markers, conn);
181 double* label = label_array.fortran_vec ();
182
183 const T padded_im_array = conn.create_padded (im, 0);
184 const P* padded_im = padded_im_array.fortran_vec ();
185
186 const Array<octave_idx_type> neighbours_array
187 = conn.deleted_neighbourhood (padded_im_array.dims ());
188 const octave_idx_type* neighbours = neighbours_array.fortran_vec ();
189 const octave_idx_type n_neighbours = neighbours_array.numel ();
190
191 // We need two flags per voxel for this implementation:
192 // 1. Whether a voxel has been labelled or not. (TODO profile this later,
193 // maybe it's enough to do label > 0)
194 // 2. Whether a voxel can go into the queue. Reasons to not go into
195 // the queue are: it's a padding voxel, it's already in the queue,
196 // it's already been labelled.
197
198 bool* label_flag = padded_markers.fortran_vec ();
199
200 boolNDArray queue_flag_array (padded_markers);
201 connectivity::set_padding (markers.dims (), padded_markers.dims (),
202 queue_flag_array, true);
203 bool* queue_flag = queue_flag_array.fortran_vec ();
204
205 const octave_idx_type n = padded_im_array.numel ();
206 octave_idx_type pos = 0;
207
208 // 2. The neighboring pixels of each marked area are inserted into a
209 // priority queue with a priority level corresponding to the gray level
210 // of the pixel.
211 std::priority_queue<Voxel<P>, std::vector<Voxel<P>>, std::greater<Voxel<P>>> q;
212 for (octave_idx_type i = 0; i < n; i++)
213 if (label_flag[i])
214 for (octave_idx_type j = 0; j < n_neighbours; j++)
215 {
216 const octave_idx_type ij = i + neighbours[j];
217 if (! queue_flag[ij])
218 {
219 queue_flag[ij] = true;
220 q.push (Voxel<P> (padded_im[ij], ij, pos++));
221 }
222 }
223
224 // 3. The pixel with the lowest priority level is extracted from the
225 // priority queue. If the neighbors of the extracted pixel that have
226 // already been labeled all have the same label, then the pixel is
227 // labeled with their label. All non-marked neighbors that are not yet
228 // in the priority queue are put into the priority queue.
229 // 4. Redo step 3 until the priority queue is empty.
230 //
231 // There is a detail missing on the description above. On step 3, if the
232 // labeled neighbours do *not* have the same label, should the non-labeled
233 // neighbours be added to the queue? Apparently not.
234 LabelsCollection lc (n_neighbours);
235 IdxCollection ic (n_neighbours);
236 while (! q.empty ())
237 {
238 Voxel<P> v = q.top ();
239 q.pop ();
240
241 lc.reset ();
242 ic.reset ();
243 for (octave_idx_type j = 0; j < n_neighbours; j++)
244 {
245 const octave_idx_type ij = v.idx + neighbours[j];
246 if (label_flag[ij])
247 lc.push_back(label[ij]);
248 else if (! queue_flag[ij])
249 ic.push_back(ij);
250 }
251 if (lc.numel () > 0 && lc.all_equal ())
252 {
253 label[v.idx] = lc.label ();
254 label_flag[v.idx] = true;
255 for (octave_idx_type i = 0; i < ic.numel (); i++)
256 {
257 const octave_idx_type ij = ic[i];
258 queue_flag[ij] = true;
259 q.push (Voxel<P> (padded_im[ij], ij, pos++));
260 }
261 }
262 }
263
264 conn.unpad (label_array);
265 return label_array;
266 }
267
268
269 DEFUN_DLD(watershed, args, , "\
270 -*- texinfo -*-\n\
271 @deftypefn {Function File} {} watershed (@var{im})\n\
272 @deftypefnx {Function File} {} watershed (@var{im}, @var{conn})\n\
273 Compute watershed transform.\n\
274 \n\
275 Computes by immersion\n\
276 \n\
277 Element connectivity @var{conn}, to define the size of objects, can be\n\
278 specified with a numeric scalar (number of elements in the neighborhood):\n\
279 \n\
280 @table @samp\n\
281 @item 4 or 8\n\
282 for 2 dimensional matrices;\n\
283 @item 6, 18 or 26\n\
284 for 3 dimensional matrices;\n\
285 @end table\n\
286 \n\
287 or with a binary matrix representing a connectivity array. Defaults to\n\
288 @code{conndef (ndims (@var{bw}), \"maximal\")} which is equivalent to\n\
289 @var{conn} of 8 and 26 for 2 and 3 dimensional matrices respectively.\n\
290 \n\
291 @seealso{bwdist, bwlabeln, regionprops}\n\
292 @end deftypefn")
293 {
294 const octave_idx_type nargin = args.length ();
295 if (nargin < 1 || nargin > 2)
296 print_usage ();
297
298 connectivity conn;
299 if (nargin > 1)
300 conn = octave::image::conndef (args(1));
301 else
302 {
303 try
304 {
305 conn = connectivity (args(0).ndims (), "maximal");
306 }
307 catch (invalid_connectivity& e)
308 {
309 error ("bwconncomp: failed to create MASK (%s)", e.what ());
310 }
311 }
312
313 const octave_image::value im (args(0));
314 #define IF_TYPE(IS_TYPE, VALUE_TYPE) \
315 if (im.is ## IS_TYPE ()) \
316 return octave_value (watershed (im. VALUE_TYPE ## array_value (), \
317 conn)); \
318
319 // My guess is that uint8, uint16, and double will be the most common types.
320 IF_TYPE(_uint8_type, uint8_)
321 else IF_TYPE(_uint16_type, uint16_)
322 else if (im.isfloat ())
323 {
324 if (im.iscomplex ())
325 {
326 IF_TYPE(_double_type, complex_)
327 else IF_TYPE(_single_type, float_complex_)
328 }
329 else
330 {
331 IF_TYPE(_double_type, )
332 else IF_TYPE(_single_type, float_)
333 }
334 }
335 else IF_TYPE(_uint32_type, uint32_)
336 else IF_TYPE(_uint64_type, uint64_)
337 else IF_TYPE(_int8_type, int8_)
338 else IF_TYPE(_int16_type, int16_)
339 else IF_TYPE(_int32_type, int32_)
340 else IF_TYPE(_int64_type, int64_)
341 else IF_TYPE(_uint8_type, uint8_)
342 else IF_TYPE(logical, bool_)
343
344 // default case if all other above fail.
345 error ("watershed: IM of unsupported class `%s'",
346 im.class_name ().c_str ());
347
348 #undef IF_TYPE
349 }
350
351 /*
352 ## Some simple tests that will check the multiple ways to measure
353 ## distances (comes to light on plateus)
354 %!test
355 %! ex = tril (ones (50), -1) + triu (repmat (2, [50 50]), 2);
356 %! ex(1, 1) = 1;
357 %! ex(end, end) = 1;
358 %!
359 %! in = ones (50);
360 %! in(end,1) = 0;
361 %! in(1,end) = 0;
362 %! assert (watershed (in), ex)
363
364 %!test
365 %! ex = tril (ones (49), -1) + triu (repmat (2, [49 49]), 2);
366 %! ex(1, 1) = 1;
367 %! ex(end, end) = 1;
368 %!
369 %! in = ones (49);
370 %! in(end,1) = 0;
371 %! in(1,end) = 0;
372 %! assert (watershed (in), ex)
373 %!
374 %! c = (fspecial ('disk', 5) > 0) + 1;
375 %! in(20:30,20:30) = c;
376 %! c = (fspecial ('disk', 4) > 0) + 2;
377 %! in(21:29,21:29) = c;
378 %! assert (watershed (in), ex)
379
380 %!test
381 %! ex = tril (ones (49), -1) + triu (repmat (2, [49 49]), 2);
382 %! ex(1:28,1:28) = (tril (ones (28) ,7) + triu (repmat (2, [28 28]), 10));
383 %! ex(1,9) = 1;
384 %! ex(end,end) = 1;
385 %! ex(20:29, 29) = 0;
386 %!
387 %! in = ones (49);
388 %! in(end,1) = 0;
389 %! in(1,end) = 0;
390 %! c = (fspecial ("disk", 5) > 0) + 1;
391 %! in(1:11,38:48) = c;
392 %!
393 %! assert (watershed (in), ex)
394
395 ## See http://perso.esiee.fr/~info/tw/index.html for a page on topological
396 ## watershed. The following test cases were taken from a powerpoint
397 ## presentation there http://perso.esiee.fr/~info/tw/isis03b.ppt
398 ## "A topological approach to watersheds". Presentation made by Gilles Bertrand
399 ## at the ISIS Workshop on Mathematical Morphology in Paris, France, 2003.
400 ##
401 ## From that presentation, the algorithm we must implement for Matlab
402 ## compatibility is named "Meyer".
403
404 %!test
405 %! im = [
406 %! 3 4 5 6 0
407 %! 2 3 4 5 6
408 %! 1 2 3 4 5
409 %! 0 1 2 3 4
410 %! 1 0 1 2 3];
411 %!
412 %! labeled8 = [
413 %! 1 1 1 0 2
414 %! 1 1 1 0 0
415 %! 1 1 1 1 1
416 %! 1 1 1 1 1
417 %! 1 1 1 1 1];
418 %! labeled4 = [
419 %! 1 1 1 0 3
420 %! 1 1 1 0 0
421 %! 1 1 0 2 2
422 %! 1 0 2 2 2
423 %! 0 2 2 2 2];
424 %! labeled_weird = [
425 %! 1 1 1 0 2
426 %! 1 1 1 1 0
427 %! 1 1 1 1 1
428 %! 1 1 1 1 1
429 %! 1 1 1 1 1];
430 %!
431 %! assert (watershed (im), labeled8);
432 %! assert (watershed (im, 8), labeled8);
433 %! assert (watershed (im, 4), labeled4);
434 %! assert (watershed (im, [1 1 0; 1 1 1; 0 1 1]), labeled_weird);
435
436 %!test
437 %! im = [
438 %! 2 3 30 2
439 %! 3 30 3 30
440 %! 255 31 30 4
441 %! 2 255 31 30
442 %! 1 2 255 5];
443 %!
444 %! labeled4 = [
445 %! 1 1 0 4
446 %! 1 0 3 0
447 %! 0 2 0 5
448 %! 2 2 2 0
449 %! 2 2 0 6];
450 %! labeled_weird = [
451 %! 1 1 0 3
452 %! 1 1 1 0
453 %! 0 1 1 1
454 %! 2 0 0 0
455 %! 2 2 0 4];
456 %!
457 %! assert (watershed (im, 4), labeled4);
458 %! assert (watershed (im, [1 1 0; 1 1 1; 0 1 1]), labeled_weird);
459
460 %!xtest
461 %! ## The following test is required for Matlab compatibility. There must be
462 %! ## something specific about their implementation that causes it to return
463 %! ## this value. Even when solving it on paper, we get different results.
464 %! im = [
465 %! 2 3 30 2
466 %! 3 30 3 30
467 %! 255 31 30 4
468 %! 2 255 31 30
469 %! 1 2 255 5];
470 %!
471 %! labeled8 = [
472 %! 1 1 0 3
473 %! 1 1 0 3
474 %! 0 0 0 0
475 %! 2 2 0 4
476 %! 2 2 0 4];
477 %! assert (watershed (im), labeled8);
478 %! assert (watershed (im, 8), labeled8);
479
480 %!test
481 %! im = [
482 %! 2 2 2 2 2 2 2
483 %! 2 2 30 30 30 2 2
484 %! 2 30 20 20 20 30 2
485 %! 40 40 20 20 20 40 40
486 %! 1 40 20 20 20 40 0
487 %! 1 1 40 20 40 0 0
488 %! 1 1 1 20 0 0 0];
489 %!
490 %! labeled8 = [
491 %! 1 1 1 1 1 1 1
492 %! 1 1 1 1 1 1 1
493 %! 1 1 1 1 1 1 1
494 %! 0 0 0 0 0 0 0
495 %! 2 2 2 0 3 3 3
496 %! 2 2 2 0 3 3 3
497 %! 2 2 2 0 3 3 3];
498 %! labeled4 = [
499 %! 1 1 1 1 1 1 1
500 %! 1 1 1 1 1 1 1
501 %! 1 1 1 1 1 1 1
502 %! 0 1 1 1 1 1 0
503 %! 2 0 1 1 1 0 3
504 %! 2 2 0 1 0 3 3
505 %! 2 2 2 0 3 3 3];
506 %! labeled_weird = [
507 %! 1 1 1 1 1 1 1
508 %! 1 1 1 1 1 1 1
509 %! 1 1 1 1 1 1 1
510 %! 0 1 1 0 0 0 0
511 %! 2 0 0 0 3 3 3
512 %! 2 2 0 3 3 3 3
513 %! 2 2 2 0 3 3 3];
514 %!
515 %! assert (watershed (im), labeled8);
516 %! assert (watershed (im, 8), labeled8);
517 %! assert (watershed (im, 4), labeled4);
518 %! assert (watershed (im, [1 1 0; 1 1 1; 0 1 1]), labeled_weird);
519
520 %!test
521 %! im = [
522 %! 40 40 40 40 40 40 40 40 40 40 40 40 40
523 %! 40 3 3 5 5 5 10 10 10 10 15 20 40
524 %! 40 3 3 5 5 30 30 30 10 15 15 20 40
525 %! 40 3 3 5 30 20 20 20 30 15 15 20 40
526 %! 40 40 40 40 40 20 20 20 40 40 40 40 40
527 %! 40 10 10 10 40 20 20 20 40 10 10 10 40
528 %! 40 5 5 5 10 40 20 40 10 10 5 5 40
529 %! 40 1 3 5 10 15 20 15 10 5 1 0 40
530 %! 40 1 3 5 10 15 20 15 10 5 1 0 40
531 %! 40 40 40 40 40 40 40 40 40 40 40 40 40];
532 %!
533 %! labeled8 = [
534 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
535 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
536 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
537 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
538 %! 0 0 0 0 0 0 0 0 0 0 0 0 0
539 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
540 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
541 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
542 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
543 %! 2 2 2 2 2 2 0 3 3 3 3 3 3];
544 %! labeled4 = [
545 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
546 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
547 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
548 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
549 %! 0 0 0 0 1 1 1 1 1 0 0 0 0
550 %! 2 2 2 2 0 1 1 1 0 3 3 3 3
551 %! 2 2 2 2 2 0 1 0 3 3 3 3 3
552 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
553 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
554 %! 2 2 2 2 2 2 0 3 3 3 3 3 3];
555 %! labeled_weird = [
556 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
557 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
558 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
559 %! 1 1 1 1 1 1 1 1 1 1 1 1 1
560 %! 0 0 0 0 1 1 0 0 0 0 0 0 0
561 %! 2 2 2 2 0 0 0 3 3 3 3 3 3
562 %! 2 2 2 2 2 0 3 3 3 3 3 3 3
563 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
564 %! 2 2 2 2 2 2 0 3 3 3 3 3 3
565 %! 2 2 2 2 2 2 0 3 3 3 3 3 3];
566 %!
567 %! assert (watershed (im), labeled8);
568 %! assert (watershed (im, 8), labeled8);
569 %! assert (watershed (im, 4), labeled4);
570 %! assert (watershed (im, [1 1 0; 1 1 1; 0 1 1]), labeled_weird);
571
572 %!xtest
573 %! ## This test is failing for Matlab compatibility
574 %! im_full = [
575 %! 1 2 10 3 8 7 5
576 %! 3 2 5 10 8 1 4
577 %! 1 8 2 3 8 3 6];
578 %!
579 %! matlab_result_full = [
580 %! 1 1 0 3 0 4 4
581 %! 0 0 0 0 0 4 4
582 %! 2 2 2 0 4 4 4];
583 %!
584 %! assert (watershed (im_full), matlab_result_full);
585 %!
586 %! im_crop = [
587 %! 2 10 3 8 7 5
588 %! 2 5 10 8 1 4
589 %! 8 2 3 8 3 6];
590 %!
591 %! matlab_result_crop = [
592 %! 1 0 2 0 3 3
593 %! 1 0 0 0 3 3
594 %! 1 1 1 0 3 3];
595 %!
596 %! assert (watershed (im_crop), matlab_result_crop);
597 */
598