Index: third_party/edtaa/edtaa3func.cpp |
diff --git a/third_party/edtaa/edtaa3func.cpp b/third_party/edtaa/edtaa3func.cpp |
new file mode 100755 |
index 0000000000000000000000000000000000000000..beed19831efd91df08d3f0d5d497504bebcd30b5 |
--- /dev/null |
+++ b/third_party/edtaa/edtaa3func.cpp |
@@ -0,0 +1,586 @@ |
+/* |
+ * 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 |