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