| /* | |
| * edtaa3() | |
| * | |
| * Sweep-and-update Euclidean distance transform of an | |
| * image. Positive pixels are treated as object pixels, | |
| * zero or negative pixels are treated as background. | |
| * An attempt is made to treat antialiased edges correctly. | |
| * The input image must have pixels in the range [0,1], | |
| * and the antialiased image should be a box-filter | |
| * sampling of the ideal, crisp edge. | |
| * If the antialias region is more than 1 pixel wide, | |
| * the result from this transform will be inaccurate. | |
| * | |
| * By Stefan Gustavson (stefan.gustavson@gmail.com). | |
| * | |
| * Originally written in 1994, based on a verbal | |
| * description of the SSED8 algorithm published in the | |
| * PhD dissertation of Ingemar Ragnemalm. This is his | |
| * algorithm, I only implemented it in C. | |
| * | |
| * Updated in 2004 to treat border pixels correctly, | |
| * and cleaned up the code to improve readability. | |
| * | |
| * Updated in 2009 to handle anti-aliased edges. | |
| * | |
| * Updated in 2011 to avoid a corner case infinite loop. | |
| * | |
| * Updated 2012 to change license from LGPL to MIT. | |
| */ | |
| /* | |
| Copyright (C) 2009-2012 Stefan Gustavson (stefan.gustavson@gmail.com) | |
| The code in this file is distributed under the MIT license: | |
| Permission is hereby granted, free of charge, to any person obtaining a copy | |
| of this software and associated documentation files (the "Software"), to deal | |
| in the Software without restriction, including without limitation the rights | |
| to use, copy, modify, merge, publish, distribute, sublicense, and/or sell | |
| copies of the Software, and to permit persons to whom the Software is | |
| furnished to do so, subject to the following conditions: | |
| The above copyright notice and this permission notice shall be included in | |
| all copies or substantial portions of the Software. | |
| THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR | |
| IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
| FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | |
| AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
| LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | |
| OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN | |
| THE SOFTWARE. | |
| */ | |
| #include "edtaa3.h" | |
| #include <math.h> | |
| #if EDTAA_UNSIGNED_CHAR_INPUT | |
| #define IMG(i) ((double)(img[i] & 0xff)/256.0) | |
| #else | |
| #define IMG(i) (img[i]) | |
| #endif | |
| namespace EDTAA { | |
| /* | |
| * Compute the local gradient at edge pixels using convolution filters. | |
| * The gradient is computed only at edge pixels. At other places in the | |
| * image, it is never used, and it's mostly zero anyway. | |
| */ | |
| void computegradient(EdtaaImageType *img, int w, int h, double *gx, double *gy) | |
| { | |
| int i,j,k; | |
| double glength; | |
| #define SQRT2 1.4142136 | |
| for(i = 1; i < h-1; i++) { // Avoid edges where the kernels would spill over | |
| for(j = 1; j < w-1; j++) { | |
| k = i*w + j; | |
| if((IMG(k)>0.0) && (IMG(k)<1.0)) { // Compute gradient for edge pixels only | |
| 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); | |
| 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); | |
| glength = gx[k]*gx[k] + gy[k]*gy[k]; | |
| if(glength > 0.0) { // Avoid division by zero | |
| glength = sqrt(glength); | |
| gx[k]=gx[k]/glength; | |
| gy[k]=gy[k]/glength; | |
| } | |
| } | |
| } | |
| } | |
| // TODO: Compute reasonable values for gx, gy also around the image edges. | |
| // (These are zero now, which reduces the accuracy for a 1-pixel wide region | |
| // around the image edge.) 2x2 kernels would be suitable for this. | |
| } | |
| /* | |
| * A somewhat tricky function to approximate the distance to an edge in a | |
| * certain pixel, with consideration to either the local gradient (gx,gy) | |
| * or the direction to the pixel (dx,dy) and the pixel greyscale value a. | |
| * The latter alternative, using (dx,dy), is the metric used by edtaa2(). | |
| * Using a local estimate of the edge gradient (gx,gy) yields much better | |
| * accuracy at and near edges, and reduces the error even at distant pixels | |
| * provided that the gradient direction is accurately estimated. | |
| */ | |
| static double edgedf(double gx, double gy, double a) | |
| { | |
| double df, glength, temp, a1; | |
| if ((gx == 0) || (gy == 0)) { // Either A) gu or gv are zero, or B) both | |
| df = 0.5-a; // Linear approximation is A) correct or B) a fair guess | |
| } else { | |
| glength = sqrt(gx*gx + gy*gy); | |
| if(glength>0) { | |
| gx = gx/glength; | |
| gy = gy/glength; | |
| } | |
| /* Everything is symmetric wrt sign and transposition, | |
| * so move to first octant (gx>=0, gy>=0, gx>=gy) to | |
| * avoid handling all possible edge directions. | |
| */ | |
| gx = fabs(gx); | |
| gy = fabs(gy); | |
| if(gx<gy) { | |
| temp = gx; | |
| gx = gy; | |
| gy = temp; | |
| } | |
| a1 = 0.5*gy/gx; | |
| if (a < a1) { // 0 <= a < a1 | |
| df = 0.5*(gx + gy) - sqrt(2.0*gx*gy*a); | |
| } else if (a < (1.0-a1)) { // a1 <= a <= 1-a1 | |
| df = (0.5-a)*gx; | |
| } else { // 1-a1 < a <= 1 | |
| df = -0.5*(gx + gy) + sqrt(2.0*gx*gy*(1.0-a)); | |
| } | |
| } | |
| return df; | |
| } | |
| static double distaa3(EdtaaImageType *img, double *gximg, double *gyimg, int w, int c, int xc, int yc, int xi, int yi) | |
| { | |
| double di, df, dx, dy, gx, gy, a; | |
| int closest; | |
| closest = c-xc-yc*w; // Index to the edge pixel pointed to from c | |
| a = IMG(closest); // Grayscale value at the edge pixel | |
| gx = gximg[closest]; // X gradient component at the edge pixel | |
| gy = gyimg[closest]; // Y gradient component at the edge pixel | |
| if(a > 1.0) a = 1.0; | |
| if(a < 0.0) a = 0.0; // Clip grayscale values outside the range [0,1] | |
| if(a == 0.0) return 1000000.0; // Not an object pixel, return "very far" ("don't know yet") | |
| dx = (double)xi; | |
| dy = (double)yi; | |
| di = sqrt(dx*dx + dy*dy); // Length of integer vector, like a traditional EDT | |
| if(di==0) { // Use local gradient only at edges | |
| // Estimate based on local gradient only | |
| df = edgedf(gx, gy, a); | |
| } else { | |
| // Estimate gradient based on direction to edge (accurate for large di) | |
| df = edgedf(dx, dy, a); | |
| } | |
| return di + df; // Same metric as edtaa2, except at edges (where di=0) | |
| } | |
| // Shorthand macro: add ubiquitous parameters dist, gx, gy, img and w and call distaa3() | |
| #define DISTAA(c,xc,yc,xi,yi) (distaa3(img, gx, gy, w, c, xc, yc, xi, yi)) | |
| void edtaa3(EdtaaImageType *img, double *gx, double *gy, int w, int h, short *distx, short *disty, double *dist) | |
| { | |
| int x, y, i, c; | |
| int offset_u, offset_ur, offset_r, offset_rd, | |
| offset_d, offset_dl, offset_l, offset_lu; | |
| double olddist, newdist; | |
| int cdistx, cdisty, newdistx, newdisty; | |
| int changed; | |
| double epsilon = 1e-3; | |
| double a; | |
| /* Initialize index offsets for the current image width */ | |
| offset_u = -w; | |
| offset_ur = -w+1; | |
| offset_r = 1; | |
| offset_rd = w+1; | |
| offset_d = w; | |
| offset_dl = w-1; | |
| offset_l = -1; | |
| offset_lu = -w-1; | |
| /* Initialize the distance images */ | |
| for(i=0; i<w*h; i++) { | |
| distx[i] = 0; // At first, all pixels point to | |
| disty[i] = 0; // themselves as the closest known. | |
| a = IMG(i); | |
| if(a <= 0.0) | |
| { | |
| dist[i]= 1000000.0; // Big value, means "not set yet" | |
| } | |
| else if (a<1.0) { | |
| dist[i] = edgedf(gx[i], gy[i], a); // Gradient-assisted estimate | |
| } | |
| else { | |
| dist[i]= 0.0; // Inside the object | |
| } | |
| } | |
| /* Perform the transformation */ | |
| do | |
| { | |
| changed = 0; | |
| /* Scan rows, except first row */ | |
| for(y=1; y<h; y++) | |
| { | |
| /* move index to leftmost pixel of current row */ | |
| i = y*w; | |
| /* scan right, propagate distances from above & left */ | |
| /* Leftmost pixel is special, has no left neighbors */ | |
| olddist = dist[i]; | |
| if(olddist > 0) // If non-zero distance or not set yet | |
| { | |
| c = i + offset_u; // Index of candidate for testing | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_ur; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| i++; | |
| /* Middle pixels have all neighbors */ | |
| for(x=1; x<w-1; x++, i++) | |
| { | |
| olddist = dist[i]; | |
| if(olddist <= 0) continue; // No need to update further | |
| c = i+offset_l; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_lu; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_u; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_ur; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| /* Rightmost pixel of row is special, has no right neighbors */ | |
| olddist = dist[i]; | |
| if(olddist > 0) // If not already zero distance | |
| { | |
| c = i+offset_l; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_lu; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_u; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty+1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| /* Move index to second rightmost pixel of current row. */ | |
| /* Rightmost pixel is skipped, it has no right neighbor. */ | |
| i = y*w + w-2; | |
| /* scan left, propagate distance from right */ | |
| for(x=w-2; x>=0; x--, i--) | |
| { | |
| olddist = dist[i]; | |
| if(olddist <= 0) continue; // Already zero distance | |
| c = i+offset_r; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| } | |
| /* Scan rows in reverse order, except last row */ | |
| for(y=h-2; y>=0; y--) | |
| { | |
| /* move index to rightmost pixel of current row */ | |
| i = y*w + w-1; | |
| /* Scan left, propagate distances from below & right */ | |
| /* Rightmost pixel is special, has no right neighbors */ | |
| olddist = dist[i]; | |
| if(olddist > 0) // If not already zero distance | |
| { | |
| c = i+offset_d; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_dl; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| i--; | |
| /* Middle pixels have all neighbors */ | |
| for(x=w-2; x>0; x--, i--) | |
| { | |
| olddist = dist[i]; | |
| if(olddist <= 0) continue; // Already zero distance | |
| c = i+offset_r; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_rd; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_d; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_dl; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| /* Leftmost pixel is special, has no left neighbors */ | |
| olddist = dist[i]; | |
| if(olddist > 0) // If not already zero distance | |
| { | |
| c = i+offset_r; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_rd; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx-1; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| olddist=newdist; | |
| changed = 1; | |
| } | |
| c = i+offset_d; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx; | |
| newdisty = cdisty-1; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| /* Move index to second leftmost pixel of current row. */ | |
| /* Leftmost pixel is skipped, it has no left neighbor. */ | |
| i = y*w + 1; | |
| for(x=1; x<w; x++, i++) | |
| { | |
| /* scan right, propagate distance from left */ | |
| olddist = dist[i]; | |
| if(olddist <= 0) continue; // Already zero distance | |
| c = i+offset_l; | |
| cdistx = distx[c]; | |
| cdisty = disty[c]; | |
| newdistx = cdistx+1; | |
| newdisty = cdisty; | |
| newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty); | |
| if(newdist < olddist-epsilon) | |
| { | |
| distx[i]=newdistx; | |
| disty[i]=newdisty; | |
| dist[i]=newdist; | |
| changed = 1; | |
| } | |
| } | |
| } | |
| } | |
| while(changed); // Sweep until no more updates are made | |
| /* The transformation is completed. */ | |
| } | |
| } // namespace EDTAA |