1 // Copyright (C) 2002 Jeffrey E. Boyd <boyd@cpsc.ucalgary.ca>
2 // Copyright (C) 2011-2012 Jordi Gutiérrez Hermoso <jordigh@octave.org>
3 // Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
4 //
5 // This program is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU General Public License as
7 // published by the Free Software Foundation; either version 3 of the
8 // License, or (at your option) any later version.
9 //
10 // This program is distributed in the hope that it will be useful, but
11 // WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with this program; if not, see
17 // <http://www.gnu.org/licenses/>.
18
19 // Copyright
20 // Jeffrey E. Boyd and Carnë Draug for bwlabel_2d
21 // Jordi Gutiérrez Hermoso for bwlabel_nd
22
23 #include <vector>
24 #include <algorithm>
25 #include <unordered_map>
26
27 #include <octave/oct.h>
28
29 #include "union-find.h"
30 #include "connectivity.h"
31
32 #define WANTS_OCTAVE_IMAGE_VALUE 1
33 #include "octave-wrappers.h"
34
35 using namespace octave::image;
36
37 static union_find
pre_label(NDArray & L,const connectivity & conn)38 pre_label (NDArray& L, const connectivity& conn)
39 {
40 double* L_vec = L.fortran_vec ();
41 const octave_idx_type numel = L.numel ();
42
43 const Array<octave_idx_type> neighbours
44 = conn.negative_neighbourhood (L.dims ());
45 const octave_idx_type* nbr = neighbours.fortran_vec ();
46 const octave_idx_type nbr_numel = neighbours.numel ();
47
48 union_find u_f (numel);
49 for (octave_idx_type Lidx = 0; Lidx < numel; Lidx++)
50 {
51 // The boundary is always zero, so we'll always skip it, so
52 // we're never considering the neighbours of the boundary. Thus,
53 // there is no possibility of out-of-bounds error below.
54 if (L_vec[Lidx])
55 {
56 //Insert this one into its group
57 u_f.add (Lidx);
58
59 for (octave_idx_type i = 0; i < nbr_numel; i++)
60 {
61 octave_idx_type n = *nbr++ + Lidx;
62 if (L_vec[n])
63 u_f.unite (n, Lidx);
64 }
65 nbr -= nbr_numel;
66 }
67 }
68 return u_f;
69 }
70
71 static octave_idx_type
paint_labels(NDArray & L,union_find & u_f)72 paint_labels (NDArray& L, union_find& u_f)
73 {
74 double* L_vec = L.fortran_vec ();
75
76 std::unordered_map<octave_idx_type, octave_idx_type> ids_to_label;
77 octave_idx_type next_label = 1;
78
79 std::vector<octave_idx_type> idxs = u_f.get_ids (L);
80 for (auto idx = idxs.begin (); idx != idxs.end (); idx++)
81 {
82 octave_idx_type label;
83 octave_idx_type id = u_f.find (*idx);
84 auto try_label = ids_to_label.find (id);
85 if (try_label == ids_to_label.end ())
86 {
87 label = next_label++;
88 ids_to_label[id] = label;
89 }
90 else
91 label = try_label->second;
92
93 L_vec[*idx] = label;
94 }
95 return ids_to_label.size ();
96 }
97
98 static octave_value_list
bwlabel_nd(const boolNDArray & BW,const connectivity & conn)99 bwlabel_nd (const boolNDArray& BW, const connectivity& conn)
100 {
101 boolNDArray conn_mask = conn.mask;
102
103 const dim_vector size_vec = BW.dims ();
104
105 // Use temporary array with borders padded with zeros. Labels will
106 // also go in here eventually.
107 NDArray L = conn.create_padded (BW, 0);
108
109 union_find u_f = pre_label (L, conn);
110 octave_idx_type n_labels = paint_labels (L, u_f);
111
112 // Remove the zero padding...
113 conn.unpad (L);
114
115 octave_value_list rval;
116 rval(0) = L;
117 rval(1) = n_labels;
118 return rval;
119 }
120
121 static octave_idx_type
find(std::vector<octave_idx_type> & lset,octave_idx_type x)122 find (std::vector<octave_idx_type>& lset, octave_idx_type x)
123 {
124 // Follow lset until we find a value that points to itself
125 while (lset[x] != x)
126 x = lset[x];
127 return x;
128 }
129
130 static octave_value_list
bwlabel_2d(const boolMatrix & BW,const octave_idx_type & n)131 bwlabel_2d (const boolMatrix& BW, const octave_idx_type& n)
132 {
133 // This algorithm was derived from BKP Horn, Robot Vision, MIT Press,
134 // 1986, p 65 - 89 by Jeffrey E. Boyd in 2002. Some smaller changes
135 // were then introduced by Carnë Draug in 2013 to speed up by iterating
136 // down a column, and what values to use when connecting two labels
137 // to increase chance of getting them in the right order in the end.
138
139 const octave_idx_type nr = BW.rows ();
140 const octave_idx_type nc = BW.columns ();
141
142 // The labelled image
143 Matrix L (nr, nc);
144
145 std::vector<octave_idx_type> lset (nc*nr); // label table/tree
146
147 octave_idx_type ntable = 0; // number of elements in the component table/tree
148 octave_idx_type ind = 0; // linear index
149
150 bool n4, n6, n8;
151 n4 = n6 = n8 = false;
152 if (n == 4)
153 n4 = true;
154 else if (n == 6)
155 n6 = true;
156 else if (n == 8)
157 n8 = true;
158
159 const bool* BW_vec = BW.data ();
160 double* L_vec = L.fortran_vec ();
161
162 for (octave_idx_type c = 0; c < nc; c++)
163 {
164 for (octave_idx_type r = 0; r < nr; r++, ind++)
165 {
166 if (BW_vec[ind]) // if A is an object
167 {
168 octave_idx_type stride = ind - nr;
169 // Get the neighboring pixels B, C, D, and E
170 //
171 // D B
172 // C A <-- ind is linear index to A
173 // E
174 //
175 // C and B will always be needed so we get them here, but
176 // D is only needed when n is 6 or 8, and E when n is 8.
177
178 octave_idx_type B, C;
179 if (c == 0)
180 C = 0;
181 else
182 C = find (lset, L_vec[stride]);
183
184 if (r == 0)
185 B = 0;
186 else
187 B = find (lset, L_vec[ind -1]);
188
189 if (n4)
190 {
191 // apply 4 connectedness
192 if (B && C) // B and C are labeled
193 {
194 if (B != C)
195 lset[B] = C;
196
197 L_vec[ind] = C;
198 }
199 else if (B) // B is object but C is not
200 L_vec[ind] = B;
201 else if (C) // C is object but B is not
202 L_vec[ind] = C;
203 else // B, C not object - new object
204 {
205 // label and put into table
206 ntable++;
207 L_vec[ind] = lset[ntable] = ntable;
208 }
209 }
210 else if (n6)
211 {
212 // Apply 6 connectedness. Seem there's more than one
213 // possible way to do this for 2D images but for some
214 // reason, the most common seems to be the top left pixel
215 // and the bottom right
216 // See http://en.wikipedia.org/wiki/Pixel_connectivity
217
218 octave_idx_type D;
219 // D is only required for n6 and n8
220 if (r == 0 || c == 0)
221 D = 0;
222 else
223 D = find (lset, L_vec[stride -1]);
224
225 if (D) // D object, copy label and move on
226 L_vec[ind] = D;
227 else if (B && C) // B and C are labeled
228 {
229 if (B == C)
230 L_vec[ind] = B;
231 else
232 {
233 octave_idx_type tlabel = std::min (B, C);
234 lset[B] = tlabel;
235 lset[C] = tlabel;
236 L_vec[ind] = tlabel;
237 }
238 }
239 else if (B) // B is object but C is not
240 L_vec[ind] = B;
241 else if (C) // C is object but B is not
242 L_vec[ind] = C;
243 else // B, C, D not object - new object
244 {
245 // label and put into table
246 ntable++;
247 L_vec[ind] = lset[ntable] = ntable;
248 }
249 }
250 else if (n8)
251 {
252 octave_idx_type D, E;
253 // D is only required for n6 and n8
254 if (r == 0 || c == 0)
255 D = 0;
256 else
257 D = find (lset, L_vec[stride -1]);
258
259 // E is only required for n8
260 if (c == 0 || r == nr -1)
261 E = 0;
262 else
263 E = find (lset, L_vec[stride +1]);
264
265 // apply 8 connectedness
266 if (B || C || D || E)
267 {
268 octave_idx_type tlabel = D;
269 if (D)
270 ; // do nothing (tlabel is already D)
271 else if (C)
272 tlabel = C;
273 else if (E)
274 tlabel = E;
275 else if (B)
276 tlabel = B;
277
278 L_vec[ind] = tlabel;
279
280 if (B && B != tlabel)
281 lset[B] = tlabel;
282 if (C && C != tlabel)
283 lset[C] = tlabel;
284 if (D)
285 // we don't check if B != tlabel since if B
286 // is true, tlabel == B
287 lset[D] = tlabel;
288 if (E && E != tlabel)
289 lset[E] = tlabel;
290 }
291 else
292 {
293 // label and put into table
294 ntable++; // run image through the look-up table
295 L_vec[ind] = lset[ntable] = ntable;
296 }
297 }
298 }
299 else
300 L_vec[ind] = 0; // A is not an object so leave it
301 }
302 }
303
304 const octave_idx_type numel = BW.numel ();
305
306 // consolidate component table
307 for (octave_idx_type i = 0; i <= ntable; i++)
308 lset[i] = find (lset, i);
309
310 // run image through the look-up table
311 for (octave_idx_type ind = 0; ind < numel; ind++)
312 L_vec[ind] = lset[L_vec[ind]];
313
314 // count up the objects in the image
315 for (octave_idx_type i = 0; i <= ntable; i++)
316 lset[i] = 0;
317
318 for (octave_idx_type ind = 0; ind < numel; ind++)
319 lset[L_vec[ind]]++;
320
321 // number the objects from 1 through n objects
322 octave_idx_type nobj = 0;
323 lset[0] = 0;
324 for (octave_idx_type i = 1; i <= ntable; i++)
325 if (lset[i] > 0)
326 lset[i] = ++nobj;
327
328 // Run through the look-up table again, so that their numbers
329 // match the number of labels
330 for (octave_idx_type ind = 0; ind < numel; ind++)
331 L_vec[ind] = lset[L_vec[ind]];
332
333 octave_value_list rval;
334 rval(0) = L;
335 rval(1) = double (nobj);
336 return rval;
337 }
338
339 DEFUN_DLD(bwlabeln, args, , "\
340 -*- texinfo -*-\n\
341 @deftypefn {Loadable Function} {[@var{l}, @var{num}] =} bwlabeln (@var{bw})\n\
342 @deftypefnx {Loadable Function} {[@var{l}, @var{num}] =} bwlabeln (@var{bw}, @var{n})\n\
343 Label foreground objects in the n-dimensional binary image @var{bw}.\n\
344 \n\
345 The optional argument @var{n} sets the connectivity and defaults 26,\n\
346 for 26-connectivity in 3-D images. Other possible values are 18 and 6\n\
347 for 3-D images, 4 and 8 for 2-D images, or an arbitrary N-dimensional\n\
348 binary connectivity mask where each dimension is of size 3.\n\
349 \n\
350 The output @var{l} is an Nd-array where 0 indicates a background\n\
351 pixel, 1 indicates that the pixel belongs to object number 1, 2 that\n\
352 the pixel belongs to object number 2, etc. The total number of objects\n\
353 is @var{num}.\n\
354 \n\
355 The algorithm used is a disjoint-set data structure, a.k.a. union-find.\n\
356 See, for example, http://en.wikipedia.org/wiki/Union-find\n\
357 \n\
358 @seealso{bwconncomp, bwlabel, regionprops}\n\
359 @end deftypefn\n\
360 ")
361 {
362 octave_value_list rval;
363
364 const octave_idx_type nargin = args.length ();
365 if (nargin < 1 || nargin > 2)
366 print_usage ();
367
368 octave_image::value bw_value (args(0));
369 if (! bw_value.isnumeric () && ! bw_value.islogical ())
370 error ("bwlabeln: BW must be a numeric or logical matrix");
371
372 boolNDArray BW = bw_value.bool_array_value ();
373 dim_vector size_vec = BW.dims ();
374
375 connectivity conn;
376 if (nargin == 2)
377 conn = conndef (args(1));
378 else
379 {
380 try
381 {
382 conn = connectivity (BW.ndims (), "maximal");
383 }
384 catch (invalid_connectivity& e)
385 {
386 error ("bwlabeln: faild to create MASK (%s)", e.what ());
387 }
388 }
389
390 // The implementation in bwlabel_2d is faster so use it if we can
391 const octave_idx_type ndims = BW.ndims ();
392 if (ndims == 2 && boolMatrix (conn.mask) == connectivity (4).mask)
393 rval = bwlabel_2d (BW, 4);
394 else if (ndims == 2 && boolMatrix (conn.mask) == connectivity (8).mask)
395 rval = bwlabel_2d (BW, 8);
396 else
397 rval = bwlabel_nd (BW, conn);
398
399 return rval;
400 }
401
402 /*
403 %!shared a2d, a3d
404 %! a2d = [1 0 0 0 0 0 1 0 0 1
405 %! 1 0 0 1 0 1 0 1 0 1
406 %! 1 0 1 0 0 0 0 0 0 0
407 %! 0 0 0 0 0 0 0 0 0 0
408 %! 0 1 0 0 0 0 0 0 0 0
409 %! 1 1 0 1 1 1 0 0 0 0
410 %! 1 1 0 1 0 0 0 1 0 0
411 %! 1 1 0 0 0 0 1 0 1 0
412 %! 1 1 0 0 0 0 0 0 0 0
413 %! 1 1 0 0 0 1 1 0 0 1];
414 %!
415 %! a3d = a2d;
416 %! a3d(:,:,2) = [
417 %! 0 0 0 0 0 0 0 0 0 0
418 %! 1 0 0 1 1 0 0 1 0 0
419 %! 0 0 0 1 0 0 0 0 0 0
420 %! 0 0 0 0 0 0 0 0 0 0
421 %! 0 1 0 0 0 0 0 0 0 0
422 %! 1 1 0 0 1 1 0 0 0 0
423 %! 1 1 0 1 0 0 0 0 0 0
424 %! 1 0 0 0 0 0 1 0 0 0
425 %! 0 1 0 0 0 0 0 0 0 1
426 %! 1 1 0 0 0 0 1 0 0 0];
427 %!
428 %! a3d(:,:,3) = [
429 %! 1 0 0 0 0 0 0 0 0 0
430 %! 0 1 0 1 1 0 0 1 0 0
431 %! 0 0 0 1 0 0 0 0 0 0
432 %! 0 0 0 0 0 0 0 0 0 0
433 %! 0 0 0 0 0 0 0 0 0 0
434 %! 0 0 0 1 1 1 0 0 0 0
435 %! 0 0 0 0 0 0 0 0 0 0
436 %! 1 0 0 0 0 0 0 0 0 0
437 %! 1 1 0 0 0 0 0 0 0 1
438 %! 1 1 0 0 0 0 0 0 0 0];
439
440 %!test
441 %! label2dc4 = [
442 %! 1 0 0 0 0 0 8 0 0 13
443 %! 1 0 0 4 0 6 0 10 0 13
444 %! 1 0 3 0 0 0 0 0 0 0
445 %! 0 0 0 0 0 0 0 0 0 0
446 %! 0 2 0 0 0 0 0 0 0 0
447 %! 2 2 0 5 5 5 0 0 0 0
448 %! 2 2 0 5 0 0 0 11 0 0
449 %! 2 2 0 0 0 0 9 0 12 0
450 %! 2 2 0 0 0 0 0 0 0 0
451 %! 2 2 0 0 0 7 7 0 0 14];
452 %! assert (bwlabeln (a2d, 4), label2dc4)
453 %! assert (bwlabeln (a2d, [0 1 0; 1 1 1; 0 1 0]), label2dc4)
454 %! assert (bwlabeln (a2d, conndef (2, "minimal")), label2dc4)
455 %! assert (bwlabeln (a2d, conndef (3, "minimal")), label2dc4)
456
457 %!test
458 %! label2dc8 = [
459 %! 1 0 0 0 0 0 5 0 0 8
460 %! 1 0 0 3 0 5 0 5 0 8
461 %! 1 0 3 0 0 0 0 0 0 0
462 %! 0 0 0 0 0 0 0 0 0 0
463 %! 0 2 0 0 0 0 0 0 0 0
464 %! 2 2 0 4 4 4 0 0 0 0
465 %! 2 2 0 4 0 0 0 7 0 0
466 %! 2 2 0 0 0 0 7 0 7 0
467 %! 2 2 0 0 0 0 0 0 0 0
468 %! 2 2 0 0 0 6 6 0 0 9];
469 %! assert (bwlabeln (a2d, 8), label2dc8)
470 %! assert (bwlabeln (a2d, ones (3)), label2dc8)
471 %! assert (bwlabeln (a2d, conndef (2, "maximal")), label2dc8)
472 %! assert (bwlabeln (a2d, conndef (3, "maximal")), label2dc8)
473
474 %!test
475 %! label3dc8 = [
476 %! 1 0 0 0 0 0 5 0 0 8
477 %! 1 0 0 3 0 5 0 5 0 8
478 %! 1 0 3 0 0 0 0 0 0 0
479 %! 0 0 0 0 0 0 0 0 0 0
480 %! 0 2 0 0 0 0 0 0 0 0
481 %! 2 2 0 4 4 4 0 0 0 0
482 %! 2 2 0 4 0 0 0 7 0 0
483 %! 2 2 0 0 0 0 7 0 7 0
484 %! 2 2 0 0 0 0 0 0 0 0
485 %! 2 2 0 0 0 6 6 0 0 9];
486 %! label3dc8(:,:,2) = [
487 %! 0 0 0 0 0 0 0 0 0 0
488 %! 10 0 0 12 12 0 0 16 0 0
489 %! 0 0 0 12 0 0 0 0 0 0
490 %! 0 0 0 0 0 0 0 0 0 0
491 %! 0 11 0 0 0 0 0 0 0 0
492 %! 11 11 0 0 13 13 0 0 0 0
493 %! 11 11 0 13 0 0 0 0 0 0
494 %! 11 0 0 0 0 0 14 0 0 0
495 %! 0 11 0 0 0 0 0 0 0 17
496 %! 11 11 0 0 0 0 15 0 0 0];
497 %! label3dc8(:,:,3) = [
498 %! 18 0 0 0 0 0 0 0 0 0
499 %! 0 18 0 20 20 0 0 22 0 0
500 %! 0 0 0 20 0 0 0 0 0 0
501 %! 0 0 0 0 0 0 0 0 0 0
502 %! 0 0 0 0 0 0 0 0 0 0
503 %! 0 0 0 21 21 21 0 0 0 0
504 %! 0 0 0 0 0 0 0 0 0 0
505 %! 19 0 0 0 0 0 0 0 0 0
506 %! 19 19 0 0 0 0 0 0 0 23
507 %! 19 19 0 0 0 0 0 0 0 0];
508 %! assert (bwlabeln (a3d, 8), label3dc8)
509 %! assert (bwlabeln (a3d, ones (3, 3)), label3dc8)
510 %! assert (bwlabeln (a3d, conndef (2, "maximal")), label3dc8)
511
512 %!test
513 %! label3dc26 = [
514 %! 1 0 0 0 0 0 3 0 0 7
515 %! 1 0 0 3 0 3 0 3 0 7
516 %! 1 0 3 0 0 0 0 0 0 0
517 %! 0 0 0 0 0 0 0 0 0 0
518 %! 0 2 0 0 0 0 0 0 0 0
519 %! 2 2 0 4 4 4 0 0 0 0
520 %! 2 2 0 4 0 0 0 6 0 0
521 %! 2 2 0 0 0 0 6 0 6 0
522 %! 2 2 0 0 0 0 0 0 0 0
523 %! 2 2 0 0 0 5 5 0 0 6];
524 %! label3dc26(:,:,2) = [
525 %! 0 0 0 0 0 0 0 0 0 0
526 %! 1 0 0 3 3 0 0 3 0 0
527 %! 0 0 0 3 0 0 0 0 0 0
528 %! 0 0 0 0 0 0 0 0 0 0
529 %! 0 2 0 0 0 0 0 0 0 0
530 %! 2 2 0 0 4 4 0 0 0 0
531 %! 2 2 0 4 0 0 0 0 0 0
532 %! 2 0 0 0 0 0 6 0 0 0
533 %! 0 2 0 0 0 0 0 0 0 6
534 %! 2 2 0 0 0 0 5 0 0 0];
535 %! label3dc26(:,:,3) = [
536 %! 1 0 0 0 0 0 0 0 0 0
537 %! 0 1 0 3 3 0 0 3 0 0
538 %! 0 0 0 3 0 0 0 0 0 0
539 %! 0 0 0 0 0 0 0 0 0 0
540 %! 0 0 0 0 0 0 0 0 0 0
541 %! 0 0 0 4 4 4 0 0 0 0
542 %! 0 0 0 0 0 0 0 0 0 0
543 %! 2 0 0 0 0 0 0 0 0 0
544 %! 2 2 0 0 0 0 0 0 0 6
545 %! 2 2 0 0 0 0 0 0 0 0];
546 %! assert (bwlabeln (a3d, 26), label3dc26)
547 %! assert (bwlabeln (a3d, ones (3, 3, 3)), label3dc26)
548 %! assert (bwlabeln (a3d, conndef (3, "maximal")), label3dc26)
549
550 %!test
551 %! label3dc18 = [
552 %! 1 0 0 0 0 0 3 0 0 7
553 %! 1 0 0 3 0 3 0 3 0 7
554 %! 1 0 3 0 0 0 0 0 0 0
555 %! 0 0 0 0 0 0 0 0 0 0
556 %! 0 2 0 0 0 0 0 0 0 0
557 %! 2 2 0 4 4 4 0 0 0 0
558 %! 2 2 0 4 0 0 0 6 0 0
559 %! 2 2 0 0 0 0 6 0 6 0
560 %! 2 2 0 0 0 0 0 0 0 0
561 %! 2 2 0 0 0 5 5 0 0 8];
562 %! label3dc18(:,:,2) = [
563 %! 0 0 0 0 0 0 0 0 0 0
564 %! 1 0 0 3 3 0 0 3 0 0
565 %! 0 0 0 3 0 0 0 0 0 0
566 %! 0 0 0 0 0 0 0 0 0 0
567 %! 0 2 0 0 0 0 0 0 0 0
568 %! 2 2 0 0 4 4 0 0 0 0
569 %! 2 2 0 4 0 0 0 0 0 0
570 %! 2 0 0 0 0 0 6 0 0 0
571 %! 0 2 0 0 0 0 0 0 0 8
572 %! 2 2 0 0 0 0 5 0 0 0];
573 %! label3dc18(:,:,3) = [
574 %! 1 0 0 0 0 0 0 0 0 0
575 %! 0 1 0 3 3 0 0 3 0 0
576 %! 0 0 0 3 0 0 0 0 0 0
577 %! 0 0 0 0 0 0 0 0 0 0
578 %! 0 0 0 0 0 0 0 0 0 0
579 %! 0 0 0 4 4 4 0 0 0 0
580 %! 0 0 0 0 0 0 0 0 0 0
581 %! 2 0 0 0 0 0 0 0 0 0
582 %! 2 2 0 0 0 0 0 0 0 8
583 %! 2 2 0 0 0 0 0 0 0 0];
584 %! assert (bwlabeln (a3d, 18), label3dc18)
585
586 %!test
587 %! label2dc3 = [
588 %! 1 0 0 0 0 0 11 0 0 17
589 %! 1 0 0 5 0 8 0 14 0 17
590 %! 1 0 4 0 0 0 0 0 0 0
591 %! 0 0 0 0 0 0 0 0 0 0
592 %! 0 3 0 0 0 0 0 0 0 0
593 %! 2 3 0 6 7 9 0 0 0 0
594 %! 2 3 0 6 0 0 0 15 0 0
595 %! 2 3 0 0 0 0 12 0 16 0
596 %! 2 3 0 0 0 0 0 0 0 0
597 %! 2 3 0 0 0 10 13 0 0 18];
598 %! assert (bwlabeln (a2d, [1 1 1]'), label2dc3)
599 %!
600 %! label3dc3 = label2dc3;
601 %! label3dc3(:,:,2) = [
602 %! 0 0 0 0 0 0 0 0 0 0
603 %! 19 0 0 24 26 0 0 31 0 0
604 %! 0 0 0 24 0 0 0 0 0 0
605 %! 0 0 0 0 0 0 0 0 0 0
606 %! 0 22 0 0 0 0 0 0 0 0
607 %! 20 22 0 0 27 28 0 0 0 0
608 %! 20 22 0 25 0 0 0 0 0 0
609 %! 20 0 0 0 0 0 29 0 0 0
610 %! 0 23 0 0 0 0 0 0 0 32
611 %! 21 23 0 0 0 0 30 0 0 0];
612 %! label3dc3(:,:,3) = [
613 %! 33 0 0 0 0 0 0 0 0 0
614 %! 0 35 0 37 39 0 0 42 0 0
615 %! 0 0 0 37 0 0 0 0 0 0
616 %! 0 0 0 0 0 0 0 0 0 0
617 %! 0 0 0 0 0 0 0 0 0 0
618 %! 0 0 0 38 40 41 0 0 0 0
619 %! 0 0 0 0 0 0 0 0 0 0
620 %! 34 0 0 0 0 0 0 0 0 0
621 %! 34 36 0 0 0 0 0 0 0 43
622 %! 34 36 0 0 0 0 0 0 0 0];
623 %! assert (bwlabeln (a3d, [1 1 1]'), label3dc3)
624
625 %!test
626 %! label2dc1 = zeros (size (a2d));
627 %! label2dc1(a2d != 0) = 1:nnz (a2d);
628 %! assert (bwlabeln (a2d, [1]), label2dc1);
629 %! assert (bwlabeln (a2d, [0 1 0]'), label2dc1);
630 %!
631 %! label3dc1 = zeros (size (a3d));
632 %! label3dc1(a3d != 0) = 1:nnz (a3d);
633 %! assert (bwlabeln (a3d, [1]), label3dc1);
634 %! assert (bwlabeln (a3d, [0 1 0]'), label3dc1);
635 */
636
637 // PKG_ADD: autoload ("bwlabel", which ("bwlabeln"));
638 // PKG_DEL: autoload ("bwlabel", which ("bwlabeln"), "remove");
639 DEFUN_DLD(bwlabel, args, , "\
640 -*- texinfo -*-\n\
641 @deftypefn {Loadable Function} {[@var{l}, @var{num}] =} bwlabel(@var{BW})\n\
642 @deftypefnx {Loadable Function} {[@var{l}, @var{num}] =} bwlabel(@var{BW}, @var{n})\n\
643 Label binary 2 dimensional image.\n\
644 \n\
645 Labels foreground objects in the binary image @var{bw}.\n\
646 The output @var{l} is a matrix where 0 indicates a background pixel,\n\
647 1 indicates that the pixel belongs to object number 1, 2 that the pixel\n\
648 belongs to object number 2, etc.\n\
649 The total number of objects is @var{num}.\n\
650 \n\
651 Two pixels belong to the same object if they are neighbors. By default\n\
652 the algorithm uses 8-connectivity to define a neighborhood, but this\n\
653 can be changed through the argument @var{n}, which can be either 4, 6, or 8.\n\
654 \n\
655 @seealso{bwconncomp, bwlabeln, regionprops}\n\
656 @end deftypefn\n\
657 ")
658 {
659 octave_value_list rval;
660
661 const octave_idx_type nargin = args.length ();
662 if (nargin < 1 || nargin > 2)
663 print_usage ();
664
665 // We do not check error state after conversion to boolMatrix
666 // because what we want is to actually get a boolean matrix
667 // with all non-zero elements as true (Matlab compatibility).
668 octave_image::value bw_value (args(0));
669 if ((! bw_value.isnumeric () && ! bw_value.islogical ()) ||
670 bw_value.ndims () != 2)
671 error ("bwlabel: BW must be a 2D matrix");
672
673 // For some reason, we can't use bool_matrix_value() to get a
674 // a boolMatrix since it will error if there's values other
675 // than 0 and 1 (whatever bool_array_value() does, bool_matrix_value()
676 // does not).
677 const boolMatrix BW = bw_value.bool_array_value ();
678
679 // N-hood connectivity
680 const octave_idx_type n = nargin < 2 ? 8 : args(1).idx_type_value ();
681 if (n != 4 && n!= 6 && n != 8)
682 error ("bwlabel: BW must be a 2 dimensional matrix");
683
684 return bwlabel_2d (BW, n);
685 }
686
687 /*
688 %!shared in
689 %! in = rand (10) > 0.8;
690 %!assert (bwlabel (in, 4), bwlabeln (in, 4));
691 %!assert (bwlabel (in, 4), bwlabeln (in, [0 1 0; 1 1 1; 0 1 0]));
692 %!assert (bwlabel (in, 8), bwlabeln (in, 8));
693 %!assert (bwlabel (in, 8), bwlabeln (in, [1 1 1; 1 1 1; 1 1 1]));
694
695 %!assert (bwlabel (logical ([0 1 0; 0 0 0; 1 0 1])), [0 2 0; 0 0 0; 1 0 3]);
696 %!assert (bwlabel ([0 1 0; 0 0 0; 1 0 1]), [0 2 0; 0 0 0; 1 0 3]);
697
698 ## Support any type of real non-zero value
699 %!assert (bwlabel ([0 -1 0; 0 0 0; 5 0 0.2]), [0 2 0; 0 0 0; 1 0 3]);
700
701 %!shared in, out
702 %!
703 %! in = [ 0 1 1 0 0 1 0 0 0 0
704 %! 0 0 0 1 0 0 0 0 0 1
705 %! 0 1 1 0 0 0 0 0 1 1
706 %! 1 0 0 0 0 0 0 1 0 0
707 %! 0 0 0 0 0 1 1 0 0 0
708 %! 0 0 0 0 0 0 0 0 0 0
709 %! 0 0 0 1 0 0 0 0 0 0
710 %! 0 0 0 0 1 1 0 1 0 0
711 %! 0 0 0 1 0 1 0 1 0 1
712 %! 1 1 0 0 0 0 0 1 1 0];
713 %!
714 %! out = [ 0 3 3 0 0 9 0 0 0 0
715 %! 0 0 0 5 0 0 0 0 0 13
716 %! 0 4 4 0 0 0 0 0 13 13
717 %! 1 0 0 0 0 0 0 11 0 0
718 %! 0 0 0 0 0 10 10 0 0 0
719 %! 0 0 0 0 0 0 0 0 0 0
720 %! 0 0 0 6 0 0 0 0 0 0
721 %! 0 0 0 0 8 8 0 12 0 0
722 %! 0 0 0 7 0 8 0 12 0 14
723 %! 2 2 0 0 0 0 0 12 12 0];
724 %!assert (nthargout ([1 2], @bwlabel, in, 4), {out, 14});
725 %!assert (nthargout ([1 2], @bwlabel, logical (in), 4), {out, 14});
726 %!
727 %! out = [ 0 3 3 0 0 7 0 0 0 0
728 %! 0 0 0 3 0 0 0 0 0 11
729 %! 0 4 4 0 0 0 0 0 11 11
730 %! 1 0 0 0 0 0 0 9 0 0
731 %! 0 0 0 0 0 8 8 0 0 0
732 %! 0 0 0 0 0 0 0 0 0 0
733 %! 0 0 0 5 0 0 0 0 0 0
734 %! 0 0 0 0 5 5 0 10 0 0
735 %! 0 0 0 6 0 5 0 10 0 12
736 %! 2 2 0 0 0 0 0 10 10 0];
737 %!assert (nthargout ([1 2], @bwlabel, in, 6), {out, 12});
738 %!assert (nthargout ([1 2], @bwlabel, logical (in), 6), {out, 12});
739 %!
740 %! ## The labeled image is not the same as Matlab, but they are
741 %! ## labeled correctly. Do we really need to get them properly
742 %! ## ordered? (the algorithm in bwlabeln does it)
743 %! mout = [0 1 1 0 0 4 0 0 0 0
744 %! 0 0 0 1 0 0 0 0 0 5
745 %! 0 1 1 0 0 0 0 0 5 5
746 %! 1 0 0 0 0 0 0 5 0 0
747 %! 0 0 0 0 0 5 5 0 0 0
748 %! 0 0 0 0 0 0 0 0 0 0
749 %! 0 0 0 3 0 0 0 0 0 0
750 %! 0 0 0 0 3 3 0 6 0 0
751 %! 0 0 0 3 0 3 0 6 0 6
752 %! 2 2 0 0 0 0 0 6 6 0];
753 %!
754 %! out = [ 0 2 2 0 0 4 0 0 0 0
755 %! 0 0 0 2 0 0 0 0 0 5
756 %! 0 2 2 0 0 0 0 0 5 5
757 %! 2 0 0 0 0 0 0 5 0 0
758 %! 0 0 0 0 0 5 5 0 0 0
759 %! 0 0 0 0 0 0 0 0 0 0
760 %! 0 0 0 3 0 0 0 0 0 0
761 %! 0 0 0 0 3 3 0 6 0 0
762 %! 0 0 0 3 0 3 0 6 0 6
763 %! 1 1 0 0 0 0 0 6 6 0];
764 %!assert (nthargout ([1 2], @bwlabel, in, 8), {out, 6});
765 %!assert (nthargout ([1 2], @bwlabel, logical (in), 8), {out, 6});
766 %!
767 %!error bwlabel (rand (10, 10, 10) > 0.8, 4)
768 %!error bwlabel (rand (10) > 0.8, "text")
769 %!error bwlabel ("text", 6)
770 */
771