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