1 /*
2 * edtaa3()
3 *
4 * Sweep-and-update Euclidean distance transform of an
5 * image. Positive pixels are treated as object pixels,
6 * zero or negative pixels are treated as background.
7 * An attempt is made to treat antialiased edges correctly.
8 * The input image must have pixels in the range [0,1],
9 * and the antialiased image should be a box-filter
10 * sampling of the ideal, crisp edge.
11 * If the antialias region is more than 1 pixel wide,
12 * the result from this transform will be inaccurate.
13 *
14 * By Stefan Gustavson (stefan.gustavson@gmail.com).
15 *
16 * Originally written in 1994, based on a verbal
17 * description of the SSED8 algorithm published in the
18 * PhD dissertation of Ingemar Ragnemalm. This is his
19 * algorithm, I only implemented it in C.
20 *
21 * Updated in 2004 to treat border pixels correctly,
22 * and cleaned up the code to improve readability.
23 *
24 * Updated in 2009 to handle anti-aliased edges.
25 *
26 * Updated in 2011 to avoid a corner case infinite loop.
27 *
28 * Updated 2012 to change license from LGPL to MIT.
29 *
30 * Updated 2014 to fix a bug with the 'gy' gradient computation.
31 *
32 */
33
34 /*
35 Copyright (C) 2009-2012 Stefan Gustavson (stefan.gustavson@gmail.com)
36 The code in this file is distributed under the MIT license:
37
38 Permission is hereby granted, free of charge, to any person obtaining a copy
39 of this software and associated documentation files (the "Software"), to deal
40 in the Software without restriction, including without limitation the rights
41 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
42 copies of the Software, and to permit persons to whom the Software is
43 furnished to do so, subject to the following conditions:
44
45 The above copyright notice and this permission notice shall be included in
46 all copies or substantial portions of the Software.
47
48 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
49 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
50 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
51 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
52 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
53 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
54 THE SOFTWARE.
55 */
56
57 #include <math.h>
58 #include "edtaa3func.h"
59
60 /*
61 * Compute the local gradient at edge pixels using convolution filters.
62 * The gradient is computed only at edge pixels. At other places in the
63 * image, it is never used, and it's mostly zero anyway.
64 */
computegradient(double * img,int w,int h,double * gx,double * gy)65 void computegradient(double *img, int w, int h, double *gx, double *gy)
66 {
67 int i,j,k,p,q;
68 double glength, phi, phiscaled, ascaled, errsign, pfrac, qfrac, err0, err1, err;
69 #define SQRT2 1.4142136
70 for(i = 1; i < h-1; i++) { // Avoid edges where the kernels would spill over
71 for(j = 1; j < w-1; j++) {
72 k = i*w + j;
73 if((img[k]>0.0) && (img[k]<1.0)) { // Compute gradient for edge pixels only
74 gx[k] = -img[k-w-1] - SQRT2*img[k-1] - img[k+w-1] + img[k-w+1] + SQRT2*img[k+1] + img[k+w+1];
75 gy[k] = -img[k-w-1] - SQRT2*img[k-w] - img[k-w+1] + img[k+w-1] + SQRT2*img[k+w] + img[k+w+1];
76 glength = gx[k]*gx[k] + gy[k]*gy[k];
77 if(glength > 0.0) { // Avoid division by zero
78 glength = sqrt(glength);
79 gx[k]=gx[k]/glength;
80 gy[k]=gy[k]/glength;
81 }
82 }
83 }
84 }
85 // TODO: Compute reasonable values for gx, gy also around the image edges.
86 // (These are zero now, which reduces the accuracy for a 1-pixel wide region
87 // around the image edge.) 2x2 kernels would be suitable for this.
88 }
89
90 /*
91 * A somewhat tricky function to approximate the distance to an edge in a
92 * certain pixel, with consideration to either the local gradient (gx,gy)
93 * or the direction to the pixel (dx,dy) and the pixel greyscale value a.
94 * The latter alternative, using (dx,dy), is the metric used by edtaa2().
95 * Using a local estimate of the edge gradient (gx,gy) yields much better
96 * accuracy at and near edges, and reduces the error even at distant pixels
97 * provided that the gradient direction is accurately estimated.
98 */
edgedf(double gx,double gy,double a)99 double edgedf(double gx, double gy, double a)
100 {
101 double df, glength, temp, a1;
102
103 if ((gx == 0) || (gy == 0)) { // Either A) gu or gv are zero, or B) both
104 df = 0.5-a; // Linear approximation is A) correct or B) a fair guess
105 } else {
106 glength = sqrt(gx*gx + gy*gy);
107 if(glength>0) {
108 gx = gx/glength;
109 gy = gy/glength;
110 }
111 /* Everything is symmetric wrt sign and transposition,
112 * so move to first octant (gx>=0, gy>=0, gx>=gy) to
113 * avoid handling all possible edge directions.
114 */
115 gx = fabs(gx);
116 gy = fabs(gy);
117 if(gx<gy) {
118 temp = gx;
119 gx = gy;
120 gy = temp;
121 }
122 a1 = 0.5*gy/gx;
123 if (a < a1) { // 0 <= a < a1
124 df = 0.5*(gx + gy) - sqrt(2.0*gx*gy*a);
125 } else if (a < (1.0-a1)) { // a1 <= a <= 1-a1
126 df = (0.5-a)*gx;
127 } else { // 1-a1 < a <= 1
128 df = -0.5*(gx + gy) + sqrt(2.0*gx*gy*(1.0-a));
129 }
130 }
131 return df;
132 }
133
distaa3(double * img,double * gximg,double * gyimg,int w,int c,int xc,int yc,int xi,int yi)134 double distaa3(double *img, double *gximg, double *gyimg, int w, int c, int xc, int yc, int xi, int yi)
135 {
136 double di, df, dx, dy, gx, gy, a;
137 int closest;
138
139 closest = c-xc-yc*w; // Index to the edge pixel pointed to from c
140 a = img[closest]; // Grayscale value at the edge pixel
141 gx = gximg[closest]; // X gradient component at the edge pixel
142 gy = gyimg[closest]; // Y gradient component at the edge pixel
143
144 if(a > 1.0) a = 1.0;
145 if(a < 0.0) a = 0.0; // Clip grayscale values outside the range [0,1]
146 if(a == 0.0) return 1000000.0; // Not an object pixel, return "very far" ("don't know yet")
147
148 dx = (double)xi;
149 dy = (double)yi;
150 di = sqrt(dx*dx + dy*dy); // Length of integer vector, like a traditional EDT
151 if(di==0) { // Use local gradient only at edges
152 // Estimate based on local gradient only
153 df = edgedf(gx, gy, a);
154 } else {
155 // Estimate gradient based on direction to edge (accurate for large di)
156 df = edgedf(dx, dy, a);
157 }
158 return di + df; // Same metric as edtaa2, except at edges (where di=0)
159 }
160
161 // Shorthand macro: add ubiquitous parameters dist, gx, gy, img and w and call distaa3()
162 #define DISTAA(c,xc,yc,xi,yi) (distaa3(img, gx, gy, w, c, xc, yc, xi, yi))
163
edtaa3(double * img,double * gx,double * gy,int w,int h,short * distx,short * disty,double * dist)164 void edtaa3(double *img, double *gx, double *gy, int w, int h, short *distx, short *disty, double *dist)
165 {
166 int x, y, i, c;
167 int offset_u, offset_ur, offset_r, offset_rd,
168 offset_d, offset_dl, offset_l, offset_lu;
169 double olddist, newdist;
170 int cdistx, cdisty, newdistx, newdisty;
171 int changed;
172 double epsilon = 1e-3;
173
174 /* Initialize index offsets for the current image width */
175 offset_u = -w;
176 offset_ur = -w+1;
177 offset_r = 1;
178 offset_rd = w+1;
179 offset_d = w;
180 offset_dl = w-1;
181 offset_l = -1;
182 offset_lu = -w-1;
183
184 /* Initialize the distance images */
185 for(i=0; i<w*h; i++) {
186 distx[i] = 0; // At first, all pixels point to
187 disty[i] = 0; // themselves as the closest known.
188 if(img[i] <= 0.0)
189 {
190 dist[i]= 1000000.0; // Big value, means "not set yet"
191 }
192 else if (img[i]<1.0) {
193 dist[i] = edgedf(gx[i], gy[i], img[i]); // Gradient-assisted estimate
194 }
195 else {
196 dist[i]= 0.0; // Inside the object
197 }
198 }
199
200 /* Perform the transformation */
201 do
202 {
203 changed = 0;
204
205 /* Scan rows, except first row */
206 for(y=1; y<h; y++)
207 {
208
209 /* move index to leftmost pixel of current row */
210 i = y*w;
211
212 /* scan right, propagate distances from above & left */
213
214 /* Leftmost pixel is special, has no left neighbors */
215 olddist = dist[i];
216 if(olddist > 0) // If non-zero distance or not set yet
217 {
218 c = i + offset_u; // Index of candidate for testing
219 cdistx = distx[c];
220 cdisty = disty[c];
221 newdistx = cdistx;
222 newdisty = cdisty+1;
223 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
224 if(newdist < olddist-epsilon)
225 {
226 distx[i]=newdistx;
227 disty[i]=newdisty;
228 dist[i]=newdist;
229 olddist=newdist;
230 changed = 1;
231 }
232
233 c = i+offset_ur;
234 cdistx = distx[c];
235 cdisty = disty[c];
236 newdistx = cdistx-1;
237 newdisty = cdisty+1;
238 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
239 if(newdist < olddist-epsilon)
240 {
241 distx[i]=newdistx;
242 disty[i]=newdisty;
243 dist[i]=newdist;
244 changed = 1;
245 }
246 }
247 i++;
248
249 /* Middle pixels have all neighbors */
250 for(x=1; x<w-1; x++, i++)
251 {
252 olddist = dist[i];
253 if(olddist <= 0) continue; // No need to update further
254
255 c = i+offset_l;
256 cdistx = distx[c];
257 cdisty = disty[c];
258 newdistx = cdistx+1;
259 newdisty = cdisty;
260 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
261 if(newdist < olddist-epsilon)
262 {
263 distx[i]=newdistx;
264 disty[i]=newdisty;
265 dist[i]=newdist;
266 olddist=newdist;
267 changed = 1;
268 }
269
270 c = i+offset_lu;
271 cdistx = distx[c];
272 cdisty = disty[c];
273 newdistx = cdistx+1;
274 newdisty = cdisty+1;
275 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
276 if(newdist < olddist-epsilon)
277 {
278 distx[i]=newdistx;
279 disty[i]=newdisty;
280 dist[i]=newdist;
281 olddist=newdist;
282 changed = 1;
283 }
284
285 c = i+offset_u;
286 cdistx = distx[c];
287 cdisty = disty[c];
288 newdistx = cdistx;
289 newdisty = cdisty+1;
290 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
291 if(newdist < olddist-epsilon)
292 {
293 distx[i]=newdistx;
294 disty[i]=newdisty;
295 dist[i]=newdist;
296 olddist=newdist;
297 changed = 1;
298 }
299
300 c = i+offset_ur;
301 cdistx = distx[c];
302 cdisty = disty[c];
303 newdistx = cdistx-1;
304 newdisty = cdisty+1;
305 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
306 if(newdist < olddist-epsilon)
307 {
308 distx[i]=newdistx;
309 disty[i]=newdisty;
310 dist[i]=newdist;
311 changed = 1;
312 }
313 }
314
315 /* Rightmost pixel of row is special, has no right neighbors */
316 olddist = dist[i];
317 if(olddist > 0) // If not already zero distance
318 {
319 c = i+offset_l;
320 cdistx = distx[c];
321 cdisty = disty[c];
322 newdistx = cdistx+1;
323 newdisty = cdisty;
324 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
325 if(newdist < olddist-epsilon)
326 {
327 distx[i]=newdistx;
328 disty[i]=newdisty;
329 dist[i]=newdist;
330 olddist=newdist;
331 changed = 1;
332 }
333
334 c = i+offset_lu;
335 cdistx = distx[c];
336 cdisty = disty[c];
337 newdistx = cdistx+1;
338 newdisty = cdisty+1;
339 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
340 if(newdist < olddist-epsilon)
341 {
342 distx[i]=newdistx;
343 disty[i]=newdisty;
344 dist[i]=newdist;
345 olddist=newdist;
346 changed = 1;
347 }
348
349 c = i+offset_u;
350 cdistx = distx[c];
351 cdisty = disty[c];
352 newdistx = cdistx;
353 newdisty = cdisty+1;
354 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
355 if(newdist < olddist-epsilon)
356 {
357 distx[i]=newdistx;
358 disty[i]=newdisty;
359 dist[i]=newdist;
360 changed = 1;
361 }
362 }
363
364 /* Move index to second rightmost pixel of current row. */
365 /* Rightmost pixel is skipped, it has no right neighbor. */
366 i = y*w + w-2;
367
368 /* scan left, propagate distance from right */
369 for(x=w-2; x>=0; x--, i--)
370 {
371 olddist = dist[i];
372 if(olddist <= 0) continue; // Already zero distance
373
374 c = i+offset_r;
375 cdistx = distx[c];
376 cdisty = disty[c];
377 newdistx = cdistx-1;
378 newdisty = cdisty;
379 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
380 if(newdist < olddist-epsilon)
381 {
382 distx[i]=newdistx;
383 disty[i]=newdisty;
384 dist[i]=newdist;
385 changed = 1;
386 }
387 }
388 }
389
390 /* Scan rows in reverse order, except last row */
391 for(y=h-2; y>=0; y--)
392 {
393 /* move index to rightmost pixel of current row */
394 i = y*w + w-1;
395
396 /* Scan left, propagate distances from below & right */
397
398 /* Rightmost pixel is special, has no right neighbors */
399 olddist = dist[i];
400 if(olddist > 0) // If not already zero distance
401 {
402 c = i+offset_d;
403 cdistx = distx[c];
404 cdisty = disty[c];
405 newdistx = cdistx;
406 newdisty = cdisty-1;
407 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
408 if(newdist < olddist-epsilon)
409 {
410 distx[i]=newdistx;
411 disty[i]=newdisty;
412 dist[i]=newdist;
413 olddist=newdist;
414 changed = 1;
415 }
416
417 c = i+offset_dl;
418 cdistx = distx[c];
419 cdisty = disty[c];
420 newdistx = cdistx+1;
421 newdisty = cdisty-1;
422 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
423 if(newdist < olddist-epsilon)
424 {
425 distx[i]=newdistx;
426 disty[i]=newdisty;
427 dist[i]=newdist;
428 changed = 1;
429 }
430 }
431 i--;
432
433 /* Middle pixels have all neighbors */
434 for(x=w-2; x>0; x--, i--)
435 {
436 olddist = dist[i];
437 if(olddist <= 0) continue; // Already zero distance
438
439 c = i+offset_r;
440 cdistx = distx[c];
441 cdisty = disty[c];
442 newdistx = cdistx-1;
443 newdisty = cdisty;
444 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
445 if(newdist < olddist-epsilon)
446 {
447 distx[i]=newdistx;
448 disty[i]=newdisty;
449 dist[i]=newdist;
450 olddist=newdist;
451 changed = 1;
452 }
453
454 c = i+offset_rd;
455 cdistx = distx[c];
456 cdisty = disty[c];
457 newdistx = cdistx-1;
458 newdisty = cdisty-1;
459 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
460 if(newdist < olddist-epsilon)
461 {
462 distx[i]=newdistx;
463 disty[i]=newdisty;
464 dist[i]=newdist;
465 olddist=newdist;
466 changed = 1;
467 }
468
469 c = i+offset_d;
470 cdistx = distx[c];
471 cdisty = disty[c];
472 newdistx = cdistx;
473 newdisty = cdisty-1;
474 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
475 if(newdist < olddist-epsilon)
476 {
477 distx[i]=newdistx;
478 disty[i]=newdisty;
479 dist[i]=newdist;
480 olddist=newdist;
481 changed = 1;
482 }
483
484 c = i+offset_dl;
485 cdistx = distx[c];
486 cdisty = disty[c];
487 newdistx = cdistx+1;
488 newdisty = cdisty-1;
489 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
490 if(newdist < olddist-epsilon)
491 {
492 distx[i]=newdistx;
493 disty[i]=newdisty;
494 dist[i]=newdist;
495 changed = 1;
496 }
497 }
498 /* Leftmost pixel is special, has no left neighbors */
499 olddist = dist[i];
500 if(olddist > 0) // If not already zero distance
501 {
502 c = i+offset_r;
503 cdistx = distx[c];
504 cdisty = disty[c];
505 newdistx = cdistx-1;
506 newdisty = cdisty;
507 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
508 if(newdist < olddist-epsilon)
509 {
510 distx[i]=newdistx;
511 disty[i]=newdisty;
512 dist[i]=newdist;
513 olddist=newdist;
514 changed = 1;
515 }
516
517 c = i+offset_rd;
518 cdistx = distx[c];
519 cdisty = disty[c];
520 newdistx = cdistx-1;
521 newdisty = cdisty-1;
522 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
523 if(newdist < olddist-epsilon)
524 {
525 distx[i]=newdistx;
526 disty[i]=newdisty;
527 dist[i]=newdist;
528 olddist=newdist;
529 changed = 1;
530 }
531
532 c = i+offset_d;
533 cdistx = distx[c];
534 cdisty = disty[c];
535 newdistx = cdistx;
536 newdisty = cdisty-1;
537 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
538 if(newdist < olddist-epsilon)
539 {
540 distx[i]=newdistx;
541 disty[i]=newdisty;
542 dist[i]=newdist;
543 changed = 1;
544 }
545 }
546
547 /* Move index to second leftmost pixel of current row. */
548 /* Leftmost pixel is skipped, it has no left neighbor. */
549 i = y*w + 1;
550 for(x=1; x<w; x++, i++)
551 {
552 /* scan right, propagate distance from left */
553 olddist = dist[i];
554 if(olddist <= 0) continue; // Already zero distance
555
556 c = i+offset_l;
557 cdistx = distx[c];
558 cdisty = disty[c];
559 newdistx = cdistx+1;
560 newdisty = cdisty;
561 newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
562 if(newdist < olddist-epsilon)
563 {
564 distx[i]=newdistx;
565 disty[i]=newdisty;
566 dist[i]=newdist;
567 changed = 1;
568 }
569 }
570 }
571 }
572 while(changed); // Sweep until no more updates are made
573
574 /* The transformation is completed. */
575
576 }
577