1 // Copyright (C) 2009 Stefan Gustavson <stefan.gustavson@gmail.com>
2 // Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
3 //
4 // This program is free software; you can redistribute it and/or modify it under
5 // the terms of the GNU General Public License as published by the Free Software
6 // Foundation; either version 3 of the License, or (at your option) any later
7 // version.
8 //
9 // This program is distributed in the hope that it will be useful, but WITHOUT
10 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 // details.
13 //
14 // You should have received a copy of the GNU General Public License along with
15 // this program; if not, see <http://www.gnu.org/licenses/>.
16
17 #include <vector>
18
19 #include <octave/oct.h>
20
21 /*
22 edtfunc - Euclidean distance transform of a binary image
23
24 This is a sweep-and-update Euclidean distance transform of
25 a binary image. All positive pixels are considered object
26 pixels, zero or negative pixels are treated as background.
27
28 By Stefan Gustavson (stefan.gustavson@gmail.com).
29
30 Originally written in 1994, based on paper-only descriptions
31 of the SSED8 algorithm, invented by Per-Erik Danielsson
32 and improved by Ingemar Ragnemalm. This is a classic algorithm
33 with roots in the 1980s, still very good for the 2D case.
34
35 Updated in 2004 to treat pixels at image edges correctly,
36 and to improve code readability.
37
38 Edited in 2009 to form the foundation for Octave BWDIST:
39 added #define-configurable distance measure and function name
40
41 Edited in 2013 for C++, removed the #define stuff, and other
42 fixes for matlab compatibility.
43 */
44
edtfunc(float (* func)(short int,short int),const boolNDArray & img,std::vector<short> & distx,std::vector<short> & disty)45 void edtfunc (float (*func)(short int, short int),
46 const boolNDArray &img,
47 std::vector<short>& distx,
48 std::vector<short>& disty)
49 {
50 const int w = img.cols ();
51 const int h = img.rows ();
52 const int numel = img.numel ();
53
54 // Initialize the distance images to be all large values
55 const bool* elem = img.fortran_vec ();
56 for (octave_idx_type i = 0; i < numel; i++)
57 if(! elem[i])
58 {
59 // Large but still representable in a short, and 32000^2 +
60 // 32000^2 does not overflow an int
61 distx[i] = 32000;
62 disty[i] = 32000;
63 }
64
65 if (h == 1 || w == 1)
66 {
67 // special treatment for 1D input
68 for (octave_idx_type i = 0; i < numel; i++)
69 disty[i] = 0;
70 for (octave_idx_type i = 1; i < numel; i++)
71 if (distx[i] != 0)
72 distx[i] = distx[i-1] +1;
73
74 for (octave_idx_type i = numel-2; i >= 0; i--)
75 {
76 if (distx[i] != 0)
77 {
78 int tmp = distx[i+1] -1;
79 // signed distance, to allow index calculation
80 if (abs (tmp) < distx[i])
81 distx[i] = tmp;
82 }
83 }
84
85 return;
86 }
87
88 double olddist2, newdist2, newdistx, newdisty;
89 bool changed;
90
91 // Initialize index offsets for the current image width
92 const int offset_u = -h;
93 const int offset_ur = -h+1;
94 const int offset_r = 1;
95 const int offset_rd = h+1;
96 const int offset_d = h;
97 const int offset_dl = h-1;
98 const int offset_l = -1;
99 const int offset_lu = -h-1;
100
101 // Perform the transformation
102 int x, y, i;
103 do
104 {
105 changed = false;
106
107 // Scan rows, except first row
108 for (y = 1; y < w; y++)
109 {
110 OCTAVE_QUIT;
111 // move index to leftmost pixel of current row
112 octave_idx_type i = y*h;
113
114 /* scan right, propagate distances from above & left */
115
116 /* Leftmost pixel is special, has no left neighbors */
117 olddist2 = (*func)(distx[i], disty[i]);
118 if(olddist2 > 0) // If not already zero distance
119 {
120 newdistx = distx[i+offset_u];
121 newdisty = disty[i+offset_u]+1;
122 newdist2 = (*func)(newdistx, newdisty);
123 if(newdist2 < olddist2)
124 {
125 distx[i]=newdistx;
126 disty[i]=newdisty;
127 olddist2=newdist2;
128 changed = true;
129 }
130
131 newdistx = distx[i+offset_ur]-1;
132 newdisty = disty[i+offset_ur]+1;
133 newdist2 = (*func)(newdistx, newdisty);
134 if(newdist2 < olddist2)
135 {
136 distx[i]=newdistx;
137 disty[i]=newdisty;
138 changed = true;
139 }
140 }
141 i++;
142
143 /* Middle pixels have all neighbors */
144 for(x=1; x<h-1; x++, i++)
145 {
146 olddist2 = (*func)(distx[i], disty[i]);
147 if(olddist2 == 0) continue; // Already zero distance
148
149 newdistx = distx[i+offset_l]+1;
150 newdisty = disty[i+offset_l];
151 newdist2 = (*func)(newdistx, newdisty);
152 if(newdist2 < olddist2)
153 {
154 distx[i]=newdistx;
155 disty[i]=newdisty;
156 olddist2=newdist2;
157 changed = true;
158 }
159
160 newdistx = distx[i+offset_lu]+1;
161 newdisty = disty[i+offset_lu]+1;
162 newdist2 = (*func)(newdistx, newdisty);
163 if(newdist2 < olddist2)
164 {
165 distx[i]=newdistx;
166 disty[i]=newdisty;
167 olddist2=newdist2;
168 changed = true;
169 }
170
171 newdistx = distx[i+offset_u];
172 newdisty = disty[i+offset_u]+1;
173 newdist2 = (*func)(newdistx, newdisty);
174 if(newdist2 < olddist2)
175 {
176 distx[i]=newdistx;
177 disty[i]=newdisty;
178 olddist2=newdist2;
179 changed = true;
180 }
181
182 newdistx = distx[i+offset_ur]-1;
183 newdisty = disty[i+offset_ur]+1;
184 newdist2 = (*func)(newdistx, newdisty);
185 if(newdist2 < olddist2)
186 {
187 distx[i]=newdistx;
188 disty[i]=newdisty;
189 changed = true;
190 }
191 }
192
193 /* Rightmost pixel of row is special, has no right neighbors */
194 olddist2 = (*func)(distx[i], disty[i]);
195 if(olddist2 > 0) // If not already zero distance
196 {
197 newdistx = distx[i+offset_l]+1;
198 newdisty = disty[i+offset_l];
199 newdist2 = (*func)(newdistx, newdisty);
200 if(newdist2 < olddist2)
201 {
202 distx[i]=newdistx;
203 disty[i]=newdisty;
204 olddist2=newdist2;
205 changed = true;
206 }
207
208 newdistx = distx[i+offset_lu]+1;
209 newdisty = disty[i+offset_lu]+1;
210 newdist2 = (*func)(newdistx, newdisty);
211 if(newdist2 < olddist2)
212 {
213 distx[i]=newdistx;
214 disty[i]=newdisty;
215 olddist2=newdist2;
216 changed = true;
217 }
218
219 newdistx = distx[i+offset_u];
220 newdisty = disty[i+offset_u]+1;
221 newdist2 = (*func)(newdistx, newdisty);
222 if(newdist2 < olddist2)
223 {
224 distx[i]=newdistx;
225 disty[i]=newdisty;
226 olddist2=newdist2;
227 changed = true;
228 }
229 }
230
231 /* Move index to second rightmost pixel of current row. */
232 /* Rightmost pixel is skipped, it has no right neighbor. */
233 i = y*h + h-2;
234
235 /* scan left, propagate distance from right */
236 for(x=h-2; x>=0; x--, i--)
237 {
238 olddist2 = (*func)(distx[i], disty[i]);
239 if(olddist2 == 0) continue; // Already zero distance
240
241 newdistx = distx[i+offset_r]-1;
242 newdisty = disty[i+offset_r];
243 newdist2 = (*func)(newdistx, newdisty);
244 if(newdist2 < olddist2)
245 {
246 distx[i]=newdistx;
247 disty[i]=newdisty;
248 changed = true;
249 }
250 }
251 }
252
253 /* Scan rows in reverse order, except last row */
254 for(y=w-2; y>=0; y--)
255 {
256 OCTAVE_QUIT;
257 /* move index to rightmost pixel of current row */
258 i = y*h + h-1;
259
260 /* Scan left, propagate distances from below & right */
261
262 /* Rightmost pixel is special, has no right neighbors */
263 olddist2 = (*func)(distx[i], disty[i]);
264 if(olddist2 > 0) // If not already zero distance
265 {
266 newdistx = distx[i+offset_d];
267 newdisty = disty[i+offset_d]-1;
268 newdist2 = (*func)(newdistx, newdisty);
269 if(newdist2 < olddist2)
270 {
271 distx[i]=newdistx;
272 disty[i]=newdisty;
273 olddist2=newdist2;
274 changed = true;
275 }
276
277 newdistx = distx[i+offset_dl]+1;
278 newdisty = disty[i+offset_dl]-1;
279 newdist2 = (*func)(newdistx, newdisty);
280 if(newdist2 < olddist2)
281 {
282 distx[i]=newdistx;
283 disty[i]=newdisty;
284 changed = true;
285 }
286 }
287 i--;
288
289 /* Middle pixels have all neighbors */
290 for(x=h-2; x>0; x--, i--)
291 {
292 olddist2 = (*func)(distx[i], disty[i]);
293 if(olddist2 == 0) continue; // Already zero distance
294
295 newdistx = distx[i+offset_r]-1;
296 newdisty = disty[i+offset_r];
297 newdist2 = (*func)(newdistx, newdisty);
298 if(newdist2 < olddist2)
299 {
300 distx[i]=newdistx;
301 disty[i]=newdisty;
302 olddist2=newdist2;
303 changed = true;
304 }
305
306 newdistx = distx[i+offset_rd]-1;
307 newdisty = disty[i+offset_rd]-1;
308 newdist2 = (*func)(newdistx, newdisty);
309 if(newdist2 < olddist2)
310 {
311 distx[i]=newdistx;
312 disty[i]=newdisty;
313 olddist2=newdist2;
314 changed = true;
315 }
316
317 newdistx = distx[i+offset_d];
318 newdisty = disty[i+offset_d]-1;
319 newdist2 = (*func)(newdistx, newdisty);
320 if(newdist2 < olddist2)
321 {
322 distx[i]=newdistx;
323 disty[i]=newdisty;
324 olddist2=newdist2;
325 changed = true;
326 }
327
328 newdistx = distx[i+offset_dl]+1;
329 newdisty = disty[i+offset_dl]-1;
330 newdist2 = (*func)(newdistx, newdisty);
331 if(newdist2 < olddist2)
332 {
333 distx[i]=newdistx;
334 disty[i]=newdisty;
335 changed = true;
336 }
337 }
338 /* Leftmost pixel is special, has no left neighbors */
339 olddist2 = (*func)(distx[i], disty[i]);
340 if(olddist2 > 0) // If not already zero distance
341 {
342 newdistx = distx[i+offset_r]-1;
343 newdisty = disty[i+offset_r];
344 newdist2 = (*func)(newdistx, newdisty);
345 if(newdist2 < olddist2)
346 {
347 distx[i]=newdistx;
348 disty[i]=newdisty;
349 olddist2=newdist2;
350 changed = true;
351 }
352
353 newdistx = distx[i+offset_rd]-1;
354 newdisty = disty[i+offset_rd]-1;
355 newdist2 = (*func)(newdistx, newdisty);
356 if(newdist2 < olddist2)
357 {
358 distx[i]=newdistx;
359 disty[i]=newdisty;
360 olddist2=newdist2;
361 changed = true;
362 }
363
364 newdistx = distx[i+offset_d];
365 newdisty = disty[i+offset_d]-1;
366 newdist2 = (*func)(newdistx, newdisty);
367 if(newdist2 < olddist2)
368 {
369 distx[i]=newdistx;
370 disty[i]=newdisty;
371 olddist2=newdist2;
372 changed = true;
373 }
374 }
375
376 /* Move index to second leftmost pixel of current row. */
377 /* Leftmost pixel is skipped, it has no left neighbor. */
378 i = y*h + 1;
379 for(x=1; x<h; x++, i++)
380 {
381 /* scan right, propagate distance from left */
382 olddist2 = (*func)(distx[i], disty[i]);
383 if(olddist2 == 0) continue; // Already zero distance
384
385 newdistx = distx[i+offset_l]+1;
386 newdisty = disty[i+offset_l];
387 newdist2 = (*func)(newdistx, newdisty);
388 if(newdist2 < olddist2)
389 {
390 distx[i]=newdistx;
391 disty[i]=newdisty;
392 changed = true;
393 }
394 }
395 }
396 }
397 while (changed); // Sweep until no more updates are made
398 // The transformation is completed
399 }
400
401 // The different functions used to calculate distance, as a
402 // class so its typename can be used for edtfunc template
403 static float
euclidean(short x,short y)404 euclidean (short x, short y)
405 {
406 // the actual euclidean distance, is the square root of this. But
407 // squaring does not change the order of the distances, so we can
408 // do it in the end and only in the values that matter
409 return ((int)(x)*(x) + (y)*(y));
410 }
411
412 static float
chessboard(short x,short y)413 chessboard (short x, short y)
414 { return std::max (abs (y), abs (x)); }
415
416 static float
cityblock(short x,short y)417 cityblock (short x, short y)
418 { return abs (x) + abs (y); }
419
420 static float
quasi_euclidean(short x,short y)421 quasi_euclidean (short x, short y)
422 {
423 static const float sqrt2_1 = sqrt (2) - 1;
424 return abs(x)>abs(y) ? (abs(x) + sqrt2_1 * abs(y)) :
425 (sqrt2_1 * abs(x) + abs(y)) ;
426 }
427
428 static FloatMatrix
calc_distances(float (* func)(short,short),const boolNDArray & bw,std::vector<short> & xdist,std::vector<short> & ydist)429 calc_distances (float (*func)(short, short), const boolNDArray& bw,
430 std::vector<short>& xdist, std::vector<short>& ydist)
431 {
432 FloatMatrix dist (bw.dims ());
433 edtfunc (func, bw, xdist, ydist);
434 const int numel = dist.numel ();
435 float* dist_vec = dist.fortran_vec ();
436 for (int i = 0; i < numel; i++)
437 dist_vec[i] = (*func)(xdist[i], ydist[i]);
438 return dist;
439 }
440
441 template <class T>
calc_index(const boolNDArray & bw,const std::vector<short> & xdist,const std::vector<short> & ydist)442 T calc_index (const boolNDArray& bw, const std::vector<short>& xdist,
443 const std::vector<short>& ydist)
444 {
445 typedef typename T::element_type P;
446
447 T idx (bw.dims ());
448 const int numel = bw.numel ();
449 const int rows = bw.rows ();
450
451 P* idx_vec = idx.fortran_vec ();
452 for(int i = 0; i < numel; i++)
453 idx_vec[i] = i+1 - xdist[i] - ydist[i]*rows;
454
455 return idx;
456 }
457
458 DEFUN_DLD (bwdist, args, nargout,
459 "-*- texinfo -*-\n\
460 @deftypefn {Loadable Function} {@var{dist} =} bwdist (@var{bw})\n\
461 @deftypefnx {Loadable Function} {@var{dist} =} bwdist (@var{bw}, @var{method})\n\
462 @deftypefnx {Loadable Function} {[@var{dist}, @var{idx}] =} bwdist (@dots{})\n\
463 Compute distance transform in binary image.\n\
464 \n\
465 The image @var{bw} must be a binary matrix For @sc{matlab} compatibility, no\n\
466 check is performed, all non-zero values are considered object pixels.\n\
467 The return value @var{dist}, is the distance of each background pixel to the\n\
468 closest object pixel in a matrix of class @code{single}.\n\
469 \n\
470 @var{idx} is the linear index for the closest object, used to calculate the\n\
471 distance for each of the pixels. Its class is dependent on the number of\n\
472 elements in @var{bw}, @code{uint64} if less than 2^32 elements, @code{uint32}\n\
473 otherwise.\n\
474 \n\
475 The distance can be measured through different @var{method}s:\n\
476 \n\
477 @table @asis\n\
478 @item euclidean (default)\n\
479 \n\
480 @item chessboard\n\
481 \n\
482 @item cityblock\n\
483 \n\
484 @item quasi-euclidean\n\
485 \n\
486 @end table\n\
487 \n\
488 Currently, only 2D images are supported.\n\
489 \n\
490 @end deftypefn")
491 {
492 octave_value_list retval;
493
494 const int nargin = args.length ();
495 if (nargin < 1 || nargin > 2)
496 print_usage ();
497
498 // for matlab compatibility, we do not actually check if the values are all
499 // 0 and 1, any non-zero value is considered true
500 const boolNDArray bw = args (0).bool_array_value ();
501
502 std::string method = (nargin > 1) ? args (1).string_value () : "euclidean";
503 for (octave_idx_type q = 0; q < octave_idx_type (method.length ()); q++)
504 method[q] = tolower (method[q]);
505
506
507 // Special case of there being no foreground element (bug #50874).
508 // Because of the way the function was originally structured,
509 // handling this case as part of the rest gets a bit hairy.
510 if (! (static_cast<boolNDArray>(bw.as_column())).any()(0))
511 {
512 FloatMatrix dist (bw.dims (), std::numeric_limits<float>::infinity ());
513 retval(0) = dist;
514
515 // Compute optional 'index to closest object pixel', only if
516 // requested.
517 if (nargout > 1)
518 {
519 if (bw.numel () >= pow (2, 32))
520 retval(1) = uint64NDArray (bw.dims (), 0);
521 else
522 retval(1) = uint32NDArray (bw.dims (), 0);
523 }
524 }
525 else
526 {
527 // Allocate two arrays for temporary output values
528 const int numel = bw.numel ();
529 std::vector<short> xdist (numel, 0);
530 std::vector<short> ydist (numel, 0);
531
532 FloatMatrix dist;
533 if (method == "euclidean")
534 {
535 dist = calc_distances (euclidean, bw, xdist, ydist);
536 const Array<octave_idx_type> positions = (!bw).find ();
537 const int zpos = positions.numel();
538 const octave_idx_type* pos_vec = positions.fortran_vec ();
539 float* dist_vec = dist.fortran_vec ();
540 for (int i = 0; i < zpos; i++)
541 dist_vec[pos_vec[i]] = sqrt (dist_vec[pos_vec[i]]);
542 }
543 else if (method == "chessboard")
544 dist = calc_distances (chessboard, bw, xdist, ydist);
545 else if (method == "cityblock")
546 dist = calc_distances (cityblock, bw, xdist, ydist);
547 else if (method == "quasi-euclidean")
548 dist = calc_distances (quasi_euclidean, bw, xdist, ydist);
549 else
550 error ("bwdist: unknown METHOD '%s'", method.c_str ());
551
552 retval(0) = dist;
553
554 // Compute optional 'index to closest object pixel', only if
555 // requested
556 if (nargout > 1)
557 {
558 if (numel >= pow (2, 32))
559 retval(1) = calc_index<uint64NDArray> (bw, xdist, ydist);
560 else
561 retval(1) = calc_index<uint32NDArray> (bw, xdist, ydist);
562 }
563 }
564
565 return retval;
566 }
567
568 /*
569 %!shared bw
570 %!
571 %! bw = [0 1 0 1 0 1 1 0
572 %! 0 0 0 1 1 0 0 0
573 %! 0 0 0 1 1 0 0 0
574 %! 0 0 0 1 1 0 0 0
575 %! 0 0 1 1 1 1 1 1
576 %! 1 1 1 1 0 0 0 1
577 %! 1 1 1 0 0 0 1 0
578 %! 0 0 1 0 0 0 1 1];
579
580 %!test
581 %! out = [ 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 1.00000
582 %! 1.41421 1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 1.41421
583 %! 2.23607 2.00000 1.00000 0.00000 0.00000 1.00000 2.00000 2.00000
584 %! 2.00000 1.41421 1.00000 0.00000 0.00000 1.00000 1.00000 1.00000
585 %! 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
586 %! 0.00000 0.00000 0.00000 0.00000 1.00000 1.00000 1.00000 0.00000
587 %! 0.00000 0.00000 0.00000 1.00000 1.41421 1.00000 0.00000 1.00000
588 %! 1.00000 1.00000 0.00000 1.00000 2.00000 1.00000 0.00000 0.00000];
589 %! out = single (out);
590 %!
591 %! assert (bwdist (bw), out, 0.0001); # default is euclidean
592 %! assert (bwdist (bw, "euclidean"), out, 0.0001);
593 %! assert (bwdist (logical (bw), "euclidean"), out, 0.0001);
594
595 %!test
596 %! out = [ 1 0 1 0 1 0 0 1
597 %! 1 1 1 0 0 1 1 1
598 %! 2 2 1 0 0 1 2 2
599 %! 2 1 1 0 0 1 1 1
600 %! 1 1 0 0 0 0 0 0
601 %! 0 0 0 0 1 1 1 0
602 %! 0 0 0 1 1 1 0 1
603 %! 1 1 0 1 2 1 0 0];
604 %! out = single (out);
605 %!
606 %! assert (bwdist (bw, "chessboard"), out);
607
608 %!test
609 %! out = [ 1 0 1 0 1 0 0 1
610 %! 2 1 1 0 0 1 1 2
611 %! 3 2 1 0 0 1 2 2
612 %! 2 2 1 0 0 1 1 1
613 %! 1 1 0 0 0 0 0 0
614 %! 0 0 0 0 1 1 1 0
615 %! 0 0 0 1 2 1 0 1
616 %! 1 1 0 1 2 1 0 0];
617 %! out = single (out);
618 %!
619 %! assert (bwdist (bw, "cityblock"), out);
620
621 %!test
622 %! out = [ 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 1.00000
623 %! 1.41421 1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 1.41421
624 %! 2.41421 2.00000 1.00000 0.00000 0.00000 1.00000 2.00000 2.00000
625 %! 2.00000 1.41421 1.00000 0.00000 0.00000 1.00000 1.00000 1.00000
626 %! 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
627 %! 0.00000 0.00000 0.00000 0.00000 1.00000 1.00000 1.00000 0.00000
628 %! 0.00000 0.00000 0.00000 1.00000 1.41421 1.00000 0.00000 1.00000
629 %! 1.00000 1.00000 0.00000 1.00000 2.00000 1.00000 0.00000 0.00000];
630 %! out = single (out);
631 %!
632 %! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
633 %!
634 %! bw(logical (bw)) = 3; # there is no actual check if matrix is binary or 0 and 1
635 %! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
636 %!
637 %! bw(logical (bw)) = -2; # anything non-zero is considered object
638 %! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
639
640 %!test
641 %! bw = [ 1 1 1 1 0 1 1 1 1
642 %! 1 1 1 1 0 1 1 1 1
643 %! 1 1 0 1 1 1 1 1 1
644 %! 0 1 1 1 1 1 1 1 1];
645 %!
646 %! dist = [ 0 0 0 0 1 0 0 0 0
647 %! 0 0 0 0 1 0 0 0 0
648 %! 0 0 1 0 0 0 0 0 0
649 %! 1 0 0 0 0 0 0 0 0];
650 %! dist = single (dist);
651 %!
652 %! c = [ 1 5 9 13 13 21 25 29 33
653 %! 2 6 10 14 14 22 26 30 34
654 %! 3 7 10 15 19 23 27 31 35
655 %! 8 8 12 16 20 24 28 32 36];
656 %! c = uint32 (c);
657 %!
658 %! [dout, cout] = bwdist (bw, "euclidean");
659 %! assert (dout, dist)
660 %! assert (cout, c)
661
662 ## The quasi-euclidean method is apparently sensitive to a machine precision
663 ## error that happens in x86 systems only. This test will cause an endless
664 ## loop in case of a regression.
665 %!test
666 %! bw = [ 0 1 1 0 0 0 1 0
667 %! 0 0 0 0 0 0 0 0
668 %! 1 1 0 0 0 0 0 0
669 %! 0 0 0 0 0 0 1 0
670 %! 0 0 0 0 1 0 0 1
671 %! 0 0 0 0 0 0 0 0
672 %! 1 0 0 0 0 0 0 0
673 %! 0 0 1 0 0 1 1 0];
674 %! out = single ([
675 %! 1.00000 0.00000 0.00000 1.00000 2.00000 1.00000 0.00000 1.00000
676 %! 1.00000 1.00000 1.00000 sqrt(2) sqrt(2)+1 sqrt(2) 1.00000 sqrt(2)
677 %! 0.00000 0.00000 1.00000 2.00000 2.00000 sqrt(2) 1.00000 sqrt(2)
678 %! 1.00000 1.00000 sqrt(2) sqrt(2) 1.00000 1.00000 0.00000 1.00000
679 %! 2.00000 2.00000 2.00000 1.00000 0.00000 1.00000 1.00000 0.00000
680 %! 1.00000 sqrt(2) 2.00000 sqrt(2) 1.00000 sqrt(2) sqrt(2) 1.00000
681 %! 0.00000 1.00000 1.00000 sqrt(2) sqrt(2) 1.00000 1.00000 sqrt(2)
682 %! 1.00000 1.00000 0.00000 1.00000 1.00000 0.00000 0.00000 1.00000
683 %! ]);
684 %! assert (bwdist (bw, "quasi-euclidean"), out);
685
686 %!error <unknown METHOD> bwdist (bw, "not a valid method");
687
688 %!test
689 %! ## Special case of there being no foreground element (bug #50874)
690 %! expected_dist = single (Inf (2, 2));
691 %! expected_idx = uint32 ([0 0; 0 0]);
692 %!
693 %! [dist, idx] = bwdist (false (2, 2));
694 %! assert (dist, expected_dist)
695 %! assert (idx, expected_idx)
696 %!
697 %! [dist, idx] = bwdist (zeros (2, 2));
698 %! assert (dist, expected_dist)
699 %! assert (idx, expected_idx)
700
701 %!test
702 %! ## Special case of 1D input (bug #50874)
703 %! assert (bwdist ([1 0]), single ([0 1]))
704 %! assert (bwdist ([1 0]'), single ([0 1]'))
705 %! assert (bwdist ([0 1 0 0 0 0 1 1]), single ([1 0 1 2 2 1 0 0]))
706 %! assert (bwdist ([1 1 0 0 0 0 1 1]'), single ([0 0 1 2 2 1 0 0])')
707 %! assert (bwdist ([1 0], "euclidean"), single ([0 1]))
708 %! assert (bwdist ([1 0], "chessboard"), single ([0 1]))
709 %! assert (bwdist ([1 0], "cityblock"), single ([0 1]))
710 %! assert (bwdist ([1 0], "quasi-euclidean"), single ([0 1]))
711
712 %!test
713 %! ## test 1D input with 2nd output argument (indices) (bug #50874)
714 %! expected_dist = single ([1 0 1]);
715 %! expected_idx = uint32 ([2 2 2]);
716 %!
717 %! [dist, idx] = bwdist ([0 1 0]);
718 %! assert (dist, expected_dist)
719 %! assert (idx, expected_idx)
720 %!
721 %! [dist, idx] = bwdist ([0 1 0]');
722 %! assert (dist, expected_dist')
723 %! assert (idx, expected_idx')
724 %!
725 %! expected_dist = single ([0 0 1 0 0]);
726 %! expected_idx = uint32 ([1 2 2 4 5]);
727 %! [dist, idx] = bwdist ([1 1 0 1 1]);
728 %! assert (dist, expected_dist)
729 %! assert (idx, expected_idx)
730 %!
731 %! expected_dist = single ([1 0 1 2 1 0 0 0 1 1 0 0 0 0 1 2 3 4]);
732 %! expected_idx = uint32 ([2 2 2 2 6 6 7 8 8 11 11 12 13 14 14 14 14 14]);
733 %! [dist, idx] = bwdist ([0 1 0 0 0 1 1 1 0 0 1 1 1 1 0 0 0 0]);
734 %! assert (dist, expected_dist)
735 %! assert (idx, expected_idx)
736 %!
737 %! expected_dist = single ([0 0 1 2 1 0 0 0 1 1 0 0 0 0 1 2 1 0]);
738 %! expected_idx = uint32 ([1 2 2 2 6 6 7 8 8 11 11 12 13 14 14 14 18 18]);
739 %! [dist, idx] = bwdist ([1 1 0 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1]);
740 %! assert (dist, expected_dist)
741 %! assert (idx, expected_idx)
742
743 %!test
744 %! assert (bwdist ([0 0]), single ([Inf, Inf]))
745 %! assert (bwdist ([0 0]'), single ([Inf, Inf]'))
746
747 %!xtest
748 %! ## This is Matlab incompatible because the bottom right corners is
749 %! ## equally distant to the top right and bottom left corners. However,
750 %! ## both are correct answers, and the returned value is just
751 %! ## implementation dependent.
752 %! bw = logical ([
753 %! 0 0 1
754 %! 0 0 0
755 %! 1 0 0
756 %! ]);
757 %! expected_dist = single ([
758 %! 2.0 1.0 0.0
759 %! 1.0 sqrt(2) 1.0
760 %! 0.0 1.0 2.0
761 %! ]);
762 %! expected_idx = uint32 ([
763 %! 3 7 7
764 %! 3 3 7
765 %! 3 3 3
766 %! ]);
767 %! [dist, idx] = bwdist (bw);
768 %! assert (dist, expected_dist)
769 %! assert (idx, expected_idx)
770
771 */
772