/*---------------------------------------------------------------------------- Smooth Contours: an unsupervised method for detecting smooth contours in digital images. This code is part of the following publication and was subject to peer review: "Unsupervised Smooth Contour Detection" by Rafael Grompone von Gioi and Gregory Randall, Image Processing On Line, 2016. http://dx.doi.org/10.5201/ipol.2016.175 Copyright (c) 2016 rafael grompone von gioi , Gregory Randall Smooth Contours is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details. You should have received a copy of the GNU Affero General Public License along with this program. If not, see . ----------------------------------------------------------------------------*/ #include #include #include #include /*----------------------------------------------------------------------------*/ #ifndef FALSE #define FALSE 0 #endif /* !FALSE */ #ifndef TRUE #define TRUE 1 #endif /* !TRUE */ /*----------------------------------------------------------------------------*/ /* PI */ #ifndef M_PI #define M_PI 3.14159265358979323846 #endif /* !M_PI */ /*----------------------------------------------------------------------------*/ #define MIN(a,b) ( (a)<(b) ? (a) : (b) ) /*----------------------------------------------------------------------------*/ /* structure to store an element of the lateral regions of pixels to an arc operator. for each pixel it stores its value, the lateral distance to the arc defining the operator, and to which of the two lateral regions belongs */ struct region { double val; /* pixel value */ double w; /* absolute value of the lateral distance to the arc */ int reg; /* number of the lateral region to which it belongs: 1 or 2 */ }; /*----------------------------------------------------------------------------*/ /* structure to store a parameterization of an arc of circle */ struct arc_of_circle { int is_line_segment; double a,b,c; /* a*x+b*y+c = 0 -> line through the arc endpoints */ double d; /* b*x-a*y+d = 0 -> orthogonal to previous through midpoint */ double len; /* arc length */ double xc,yc,radius; /* center and radius of circle containing the arc */ double ang_ref,ang_span; int dir; int bbx0,bby0,bbx1,bby1; /* bounding box */ }; /*----------------------------------------------------------------------------*/ /* fatal error, print a message to standard error and exit */ static void error(char * msg) { fprintf(stderr,"error: %s\n",msg); exit(EXIT_FAILURE); } /*----------------------------------------------------------------------------*/ /* memory allocation, print an error and exit if fail */ static void * xmalloc(size_t size) { void * p; if( size == 0 ) error("xmalloc: zero size"); p = malloc(size); if( p == NULL ) error("xmalloc: out of memory"); return p; } /*----------------------------------------------------------------------------*/ /* compute a > b considering the rounding errors due to the representation of double numbers */ static int greater(double a, double b) { if( a <= b ) return FALSE; /* trivial case, return as soon as possible */ if( (a-b) < 1000 * DBL_EPSILON ) return FALSE; return TRUE; /* greater */ } /*----------------------------------------------------------------------------*/ /* Euclidean distance between x1,y1 and x2,y2 */ static double dist(double x1, double y1, double x2, double y2) { return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ); } /*----------------------------------------------------------------------------*/ /* compute the error function erf using Winitzki's approximation. the formula is given in "A handy approximation for the error function and its inverse" by Sergei Winitzki, February 6, 2008, http://sites.google.com/site/winitzki/sergei-winitzkis-files/erf-approx.pdf the error function is defined as erf(x) = 2/sqrt(pi) \int_0^x e^{-t^2} dt and the approximation is erf(x) ~ (1 - exp(-x^2 * (4/pi + axˆ2) / (1 + ax^2) ) )^1/2 with a = 8/3pi * (pi - 3) / (4 - pi) */ static double erf_winitzki(double x) { static double a = 8.0 / 3.0 / M_PI * (M_PI - 3.0) / (4.0 - M_PI); if( x < 0.0 ) return -erf_winitzki(-x); return sqrt( 1.0 - exp(-x*x * (4.0/M_PI + a*x*x) / (1.0 + a*x*x)) ); } /*----------------------------------------------------------------------------*/ /* compute a Gaussian kernel of length n, standard deviation sigma, and centered at value mean. for example, if mean=0.5, the Gaussian will be centered in the middle point between values kernel[0] and kernel[1]. kernel must be allocated to a size n. */ static void gaussian_kernel(double * kernel, int n, double sigma, double mean) { double sum = 0.0; double val; int i; /* check input */ if( kernel == NULL ) error("gaussian_kernel: kernel not allocated"); if( sigma <= 0.0 ) error("gaussian_kernel: sigma must be positive"); /* compute Gaussian kernel */ for(i=0; i 0.0 ) for(i=0; i=nx2) j -= nx2; if( j >= X ) j = nx2-1-j; val += image[j+y*X] * kernel[i]; } tmp[x+y*X] = val; } /* y axis convolution */ for(x=0; x=ny2) j -= ny2; if( j >= Y ) j = ny2-1-j; val += tmp[x+j*X] * kernel[i]; } out[x+y*X] = val; } /* free memory */ free( (void *) kernel ); free( (void *) tmp ); return out; } /*----------------------------------------------------------------------------*/ /* non oriented angle difference, returns a value in [0,pi] */ static double diff_0_pi(double a, double b) { a -= b; while( a <= -M_PI ) a += 2.0 * M_PI; while( a > M_PI ) a -= 2.0 * M_PI; if( a < 0.0 ) a = -a; return a; } /*----------------------------------------------------------------------------*/ /* return a score for chaining pixels 'from' to 'to', favoring closet point: = 0.0 invalid chaining > 0.0 valid forward chaining; the larger the value, the better the chaining < 0.0 valid backward chaining; the smaller the value, the better the chaining input: from, to the two pixel IDs to evaluate their potential chaining Ex[i], Ey[i] the sub-pixel position of point i, if i is an edge point; they take values -1,-1 if i is not an edge point Gx[i], Gy[i] the image gradient at pixel i X, Y the size of the image */ static double chain( int from, int to, double * Ex, double * Ey, double * Gx, double * Gy, int X, int Y ) { double dx,dy; /* check input */ if( Ex == NULL || Ey == NULL || Gx == NULL || Gy == NULL ) error("chain: invalid input"); if( from < 0 || to < 0 || from >= X*Y || to >= X*Y ) error("chain: one of the points is out the image"); /* check that the points are different and valid edge points, otherwise return invalid chaining */ if( from == to ) return 0.0; /* same pixel */ if( Ex[from] < 0.0 || Ey[from] < 0.0 || Ex[to] < 0.0 || Ey[to] < 0.0 ) return 0.0; /* one of them is not an edge point */ /* in a good chaining, the gradient should be roughly orthogonal to the line joining the two points to be chained: Gx,Gy | ------> dx,dy | thus from x-------x to ---> Gy,-Gx (orthogonal to the gradient) when Gy * dx - Gx * dy > 0, it corresponds to a forward chaining, when Gy * dx - Gx * dy < 0, it corresponds to a backward chaining. (this choice is arbitrary) first check that the gradient at both points to be chained agree in one direction, otherwise return invalid chaining. */ dx = Ex[to] - Ex[from]; dy = Ey[to] - Ey[from]; if( (Gy[from] * dx - Gx[from] * dy) * (Gy[to] * dx - Gx[to] * dy) <= 0.0 ) return 0.0; /* return the chaining score: positive for forward chaining, negative for backwards. the score is the inverse of the distance to the chaining point, to give preference to closer points */ if( (Gy[from] * dx - Gx[from] * dy) >= 0.0 ) return 1.0 / dist(Ex[from],Ey[from],Ex[to],Ey[to]); /* forward chaining */ else return -1.0 / dist(Ex[from],Ey[from],Ex[to],Ey[to]); /* backward chaining */ } /*----------------------------------------------------------------------------*/ /* compute the image gradient, giving its x and y components as well as the modulus. Gx, Gy, and modG must be already allocated. */ static void compute_gradient( double * Gx, double * Gy, double * modG, double * image, int X, int Y ) { int x,y; /* check input */ if( Gx == NULL || Gy == NULL || modG == NULL || image == NULL ) error("compute_gradient: invalid input"); /* approximate image gradient using centered differences */ for(x=1; x<(X-1); x++) for(y=1; y<(Y-1); y++) { Gx[x+y*X] = image[(x+1)+y*X] - image[(x-1)+y*X]; Gy[x+y*X] = image[x+(y+1)*X] - image[x+(y-1)*X]; modG[x+y*X] = sqrt( Gx[x+y*X] * Gx[x+y*X] + Gy[x+y*X] * Gy[x+y*X] ); } } /*----------------------------------------------------------------------------*/ /* compute sub-pixel edge points using adapted Canny and Devernay methods. input: Gx, Gy, and modG are the x and y components and modulus of the image gradient, respectively. X,Y is the image size. output: Ex and Ey will have the x and y sub-pixel coordinates of the edge points found, or -1 and -1 when not an edge point. Ex and Ey must be already allocated. a modified Canny non maximal suppression [1] is used to select edge points, and a modified Devernay sub-pixel correction [2] is used to improve the position accuracy. in both cases, the modification boils down to using only vertical or horizontal non maximal suppression and sub-pixel correction. no threshold is used on the gradient. [1] J.F. Canny, "A computational approach to edge detection", IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.8, no.6, pp.679-698, 1986. [2] F. Devernay, "A Non-Maxima Suppression Method for Edge Detection with Sub-Pixel Accuracy", Rapport de recherche 2724, INRIA, Nov. 1995. the reason for this modification is that Devernay correction is inconsistent for some configurations at 45, 225, -45 or -225 degree. in edges that should go exactly in the middle of a pixels like (5 pixels drawn): ___ | ___| | ___| | ___| | ___| the correction terms of both sides of the perfect edge are not compatible, leading to edge points with "oscillations" like: . . . . . . but the Devernay correction works very well and is very consistent when used to interpolate only along horizontal or vertical direction. this modified version requires that a pixel, to be an edge point, must be a local maximum horizontally or vertically, depending on the gradient orientation: if the x component of the gradient is larger than the y component, Gx > Gy, this means that the gradient is roughly horizontal and a horizontal maximum is required to be an edge point. when the pixel is both, a horizontal and vertical maximum, the direction with largest contrast is used for the interpolation. this usually leads to a more accurate estimation of the position. using only horizontal or vertical interpolation may lead occasionally to some missing edge points for edges at 45, 225, -45 or -225 degree. these does not happen very often and seems to be an acceptable price to pay for removing the "oscillations" in slanted edges. moreover, as it is very rare to have two consecutive missing points, using a chaining strategy that looks slightly farther away for neighbors, the chains of edge points are recovered well enough. this is done in the function chain_edge_points() below. */ static void compute_edge_points( double * Ex, double * Ey, double * modG, double * Gx, double * Gy, int X, int Y ) { int x,y,i; /* check input */ if( Ex == NULL || Ey == NULL || modG == NULL || Gx == NULL || Gy == NULL ) error("compute_edge_points: invalid input"); /* initialize Ex and Ey as non edge points for all pixels */ for(i=0; i= R and D < mod >= U, respectively. the comparisons are done using the function greater() instead of the operators > or >= so numbers differing only due to rounding errors are considered equal */ int lHm = greater(mod,L) && !greater(R,mod); /* local horizontal max */ int lVm = greater(mod,D) && !greater(U,mod); /* local vertical max */ /* if both horizontal and vertical max, interpolate along direction of max contrast; otherwise, requires max along approx gradient dir */ if( lHm && lVm && MIN(L,R) < MIN(U,D) ) Dx = 1; /* both */ else if( lHm && lVm ) Dy = 1; /* both but more contrast in vert. */ else if( lHm && gx >= gy ) Dx = 1; /* only horizontal maximum */ else if( lVm && gx <= gy ) Dy = 1; /* only vertical maximum */ /* Devernay sub-pixel correction [2] the edge point position is selected as the one of the maximum of a quadratic interpolation of the magnitude of the gradient along a unidimensional direction. the pixel must be a local maximum. so we have the values: . b a . | x = -1, |Gx| = a | | . c x = 0, |Gx| = b | | | x = 1, |Gx| = c ------------------> x -1 0 1 the x position of the maximum of the parabola passing through (-1,a), (0,b), and (1,c) is offset = (a - c) / 2(a - 2b + c) and because b >= a and b >= c, -0.5 <= offset <= 0.5 */ if( Dx > 0 || Dy > 0 ) { /* offset value is in [-0.5, 0.5] */ double a = modG[ x-Dx + (y-Dy) * X ]; double b = modG[ x + y * X ]; double c = modG[ x+Dx + (y+Dy) * X ]; double offset = 0.5 * (a - c) / (a - b - b + c); /* store edge point */ Ex[x+y*X] = x + offset * Dx; Ey[x+y*X] = y + offset * Dy; } } } /*----------------------------------------------------------------------------*/ /* chain edge points input: Ex and Ey are the sub-pixel coordinates when an edge point is present or -1,-1 otherwise. Gx, Gy and modG are the x and y components and the modulus of the image gradient, respectively. X,Y is the image size. output: next and prev will contain the number of next and previous edge points in the chain. when not chained on one of the directions, the corresponding value is set to -1. next and prev must be allocated before calling. */ static void chain_edge_points( int * next, int * prev, double * Ex, double * Ey, double * Gx, double * Gy, int X, int Y ) { int x,y,i,j,alt; /* check input */ if( next==NULL || prev==NULL || Ex==NULL || Ey==NULL || Gx==NULL || Gy==NULL ) error("chain_edge_points: invalid input"); /* initialize next and prev as non linked */ for(i=0; i= 0.0 && Ey[x+y*X] >= 0.0 ) /* must be an edge point */ { int from = x+y*X; /* edge point to be chained */ double fwd_s = 0.0; /* score of best forward chaining */ double bck_s = 0.0; /* score of best backward chaining */ int fwd = -1; /* edge point of best forward chaining */ int bck = -1; /* edge point of best backward chaining */ /* try all neighbors two pixels apart or less. the detection of edge points occasionally fails to detect some edge points at 45, 225, -45 or -225 degree. as explained before, this is small price to pay for getting accurate detections. looking for candidates for chaining two pixels apart, in most such cases, is enough to obtain good chains of edge points that accurately describes the edge. */ for(i=-2; i<=2; i++) for(j=-2; j<=2; j++) { int to = x+i + (y+j)*X; /* candidate edge point to be chained */ double s = chain(from,to,Ex,Ey,Gx,Gy,X,Y); /* score from-to */ if( s > fwd_s ) /* a better forward chaining found */ { fwd_s = s; /* set the new best forward chaining */ fwd = to; } if( s < bck_s ) /* a better backward chaining found */ { bck_s = s; /* set the new best backward chaining */ bck = to; } } /* before making the new chain, check whether the target was already chained and in that case, whether the alternative chaining is better than the proposed one. x alt x alt \ / \ / from x---------x fwd bck x---------x from we know that the best forward chain starting at from is from-fwd. but it is possible that there is an alternative chaining arriving at fwd that is better, such that alt-fwd is to be preferred to from-fwd. an analogous situation is possible in backward chaining, where an alternative link bck-alt may be better than bck-from. before making the new link, check if fwd/bck are already chained, and in such case compare the scores of the proposed chaining to the existing one, and keep only the best of the two. there is an undesirable aspect of this procedure: the result may depend on the order of exploration. consider the following configuration: a x-------x b / / c x---x d with score(a-b) < score(c-b) < score(c-d) let us consider two possible orders of exploration. order: a,b,c we will first chain a-b when exploring a. when analyzing the backward links of b, we will prefer c-b, and a-b will be unlinked. finally, when exploring c, c-d will be preferred and c-b will be unlinked. the result is just c-d. order: c,b,a we will first chain c-d when exploring c. then, when exploring the backward connections of b, c-b will be the preferred link; but because c-d was already done and has a better score, c-b cannot be linked. finally, when exploring a, the link a-b will be created because there is no better backward linking of b. the result is c-d and a-b. we did not found yet a simple algorithm to solve this problem. by simple, we mean an algorithm without two passes or the need to re-evaluate the chaining of points where one link is cut. for most edge points, there is only one possible chaining and this problem does not arise. but it does happen and a better solution is desirable. */ if( fwd >= 0 && next[from] != fwd && ((alt=prev[fwd]) < 0 || chain(alt,fwd,Ex,Ey,Gx,Gy,X,Y) < fwd_s) ) { if( next[from] >= 0 ) /* remove previous from-x link if one */ prev[next[from]] = -1; /* only prev requires explicit reset */ next[from] = fwd; /* set next of from-fwd link */ if( alt >= 0 ) /* remove alt-fwd link if one */ next[alt] = -1; /* only next requires explicit reset */ prev[fwd] = from; /* set prev of from-fwd link */ } if( bck >= 0 && prev[from] != bck && ((alt=next[bck]) < 0 || chain(alt,bck,Ex,Ey,Gx,Gy,X,Y) > bck_s ) ) { if( alt >= 0 ) /* remove bck-alt link if one */ prev[alt] = -1; /* only prev requires explicit reset */ next[bck] = from; /* set next of bck-from link */ if( prev[from] >= 0 ) /* remove previous x-from link if one */ next[prev[from]] = -1; /* only next requires explicit reset */ prev[from] = bck; /* set prev of bck-from link */ } } } /*----------------------------------------------------------------------------*/ /* simplify a chain i-j-k-l-m into i-j-l-m when the latter is more regular next and prev contain the number of next and previous edge points in the chain or -1 when not chained. Ex and Ey are the sub-pixel coordinates when an edge point is present or -1,-1 otherwise. Ex and Ey are the sub-pixel coordinates when an edge point is present or -1,-1 otherwise. this function modifies next and prev when simplifications are found. */ static void simplify_chains( int * next, int * prev, double * Ex, double * Ey, int X, int Y ) { int i,j,k,l,m; /* check input */ if( next == NULL || prev == NULL || Ex == NULL || Ey == NULL ) error("simplify_chains: invalid input"); /* look for chains of 5 consecutive chained edge points, i-j-k-l-m: x k / \ x---x x---x i j l m we will consider to link j directly to l, leaving k out of the chain. an additional requirement is that the distance of the potential link j-l is less than 2. the latter is to avoid a general regularization of the chain. in normal linking along edges, the mean distance between edge points is in the interval from 1 (when horizontal or vertical) to sqrt(2) (45 degree). so, the mean distance between edge points at two chaining steps is in [2, 2*sqrt(2)]. restricting the simplifications to cases where the two step distance is less than 2, avoid regularizing normal chaining and concentrates on its objective of simplifying complex corners. to state the conditions to simplify the chain, let us define the following angles: A is angle(i-j,l-m) x---x x---x i j l m x / B is angle(i-j,j-k) x---x i j x k \ C is angle(k-l,l-m) x---x l m the condition to simplify the chain is A < B and A < C or equivalently A < MIN(B,C). x k / \ to sum up, a chain x---x x---x becomes x---x---x---x when the latter i j l m i j l m is more regular. */ for(i=0; i=0 -> edge point, no explicit test needed */ if( (j=next[i])>=0 && (k=next[j])>=0 && (l=next[k])>=0 && (m=next[l])>=0 && dist(Ex[j],Ey[j],Ex[l],Ey[l]) < 2.0 ) { double a = atan2( Ey[j] - Ey[i], Ex[j] - Ex[i] ); /* angle(i-j) */ double b = atan2( Ey[k] - Ey[j], Ex[k] - Ex[j] ); /* angle(j-k) */ double c = atan2( Ey[l] - Ey[k], Ex[l] - Ex[k] ); /* angle(k-l) */ double d = atan2( Ey[m] - Ey[l], Ex[m] - Ex[l] ); /* angle(l-m) */ /* if A < MIN(B,C) => simplify: link j-l, unlink j-k and k-l */ if( diff_0_pi(a,d) < MIN( diff_0_pi(a,b), diff_0_pi(c,d) ) ) { next[j] = l; prev[l] = j; next[k] = -1; prev[k] = -1; } } } /*----------------------------------------------------------------------------*/ /* create a list of chained edge points composed of 3 lists x, y and curve_limits; it also computes N (the number of edge points) and M (the number of curves). x[i] and y[i] (0<=i=0 or next[i]>=0 implies edge point */ if( prev[i] >= 0 || next[i] >= 0 ) { /* a new chain found, set chain starting index to the current point and then increase the curve counter */ (*curve_limits)[*M] = *N; ++(*M); /* set k to the begining of the chain, or to i if closed curve */ for(k=i; (n=prev[k])>=0 && n!=i; k=n); /* follow the chain of edge points starting on k */ do { /* store the current point coordinates in the output lists */ (*x)[*N] = Ex[k]; (*y)[*N] = Ey[k]; ++(*N); n = next[k]; /* save the id of the next point in the chain */ next[k] = -1; /* unlink chains from k so it is not used again */ prev[k] = -1; /* for closed curves, the initial point is included again as the last point of the chain. actually, testing if the first and last points are equal is the only way to know that it is a closed curve. to understand that this code actually repeats the first point, consider a closed chain as follows: a--b | | d--c let us say that the algorithm starts by point a. it will store the coordinates of point a and then unlink a-b. then, will store point b and unlink b-c, and so on. but the link d-a is still there. (point a is no longer pointing backwards to d, because both links are removed at each step. but d is indeed still pointing to a.) so it will arrive at point a again and store its coordinates again as last point. there, it cannot continue because the link a-b was removed, there would be no next point, k would be -1 and the curve is finished. */ k = n; /* set the current point to the next in the chain */ } while( k >= 0 ); /* continue while there is a next point in the chain */ } (*curve_limits)[*M] = *N; /* store end of the last chain */ } /*----------------------------------------------------------------------------*/ /* chained, sub-pixel edge detector. based on a modified Canny non-maximal suppression and a modified Devernay sub-pixel correction. the input image is assumed to be blurred as desired. x, y, and curve_limits will be allocated. the output are the chained edge points given as 3 lists x, y and curve_limits, as well as the number N of edge points, and the number M of curves. x[i] and y[i] (0<=ixc = h; arc->yc = k; arc->radius = dist(x1,y1,h,k); } /*----------------------------------------------------------------------------*/ /* oriented difference of angles modulo 2pi */ static double diff_0_2pi(double a, double b) { a -= b; while( a < 0.0 ) a += 2.0*M_PI; while( a > 2.0*M_PI ) a -= 2.0*M_PI; return a; } /*----------------------------------------------------------------------------*/ /* initialize an arc structure from point i to point k in the list x[],y[] and with width max_w. return TRUE if the arc is smooth up to precision sigma, FALSE otherwise */ static int smooth_segment( struct arc_of_circle * arc, double * x, double * y, int i, int k, double sigma, double max_w, int X, int Y ) { int l; int j = i + (k-i) / 2; /* middle point */ double ang_i,ang_j,ang_k; /* full circles are not handled as one arc by the validation step. this is to avoid considering the special case of closed circles. in practice, the impact is minor, it only means that full circles will not be validated all together. but usually most of the partial arcs of the circle will be validated, so the result would be a curve following the full circle. the real risk is of losing small circles, which may have the limit size where removing even one pixel leads to a reject. */ if( x[i]==x[k] && y[i]==y[k] ) return FALSE; /* exactly the same x,y values are repeated in closed curves, so in this case there is no risk in using the operator == to compare doubles */ /* compute the parameters of the line through x[i],y[i] and x[k],y[k], the orthogonal line to the previous, and len as a line segment */ arc->len = dist(x[i],y[i],x[k],y[k]); /* if arc will be re-computed later */ arc->a = -(y[k] - y[i]) / arc->len; arc->b = (x[k] - x[i]) / arc->len; arc->c = -arc->a*x[i] - arc->b*y[i]; arc->d = -arc->b*0.5*(x[i]+x[k]) + arc->a*0.5*(y[i]+y[k]); /* initialize bounding box */ arc->bbx0 = arc->bbx1 = (int) x[i]; arc->bby0 = arc->bby1 = (int) y[i]; /* line segment or arc? decide evaluating the middle point */ if( fabs( arc->a*x[j] + arc->b*y[j] + arc->c ) < sigma ) { /* line segment */ arc->is_line_segment = TRUE; /* check that all points are not farther than sigma from the line */ for(l=i; l<=k; l++) { /* update bounding box */ if( x[l] < (double) arc->bbx0 ) arc->bbx0 = (int) x[l]; if( y[l] < (double) arc->bby0 ) arc->bby0 = (int) y[l]; if( x[l] > (double) arc->bbx1 ) arc->bbx1 = (int) x[l]; if( y[l] > (double) arc->bby1 ) arc->bby1 = (int) y[l]; if( fabs(arc->a*x[l] + arc->b*y[l] + arc->c) > sigma ) return FALSE; /* not smooth */ } } else { /* arc */ arc->is_line_segment = FALSE; /* compute circle center and radius */ circle_through(arc,x[i],y[i],x[j],y[j],x[k],y[k]); /* arc direction */ ang_i = atan2(y[i] - arc->yc, x[i] - arc->xc); ang_j = atan2(y[j] - arc->yc, x[j] - arc->xc); ang_k = atan2(y[k] - arc->yc, x[k] - arc->xc); if( diff_0_2pi(ang_j,ang_i) < diff_0_2pi(ang_k,ang_i) ) { arc->dir = 1; /* direction of turn of the arc */ arc->ang_span = diff_0_2pi(ang_k,ang_i); /* span */ arc->ang_ref = ang_i; } else { arc->dir = -1; /* direction of turn of the arc */ arc->ang_span = diff_0_2pi(ang_i,ang_k); /* span */ arc->ang_ref = ang_k; } /* compute arc length */ arc->len = arc->radius * arc->ang_span; /* check that all points are not farther than sigma from the arc */ for(l=i; l<=k; l++) { /* update bounding box */ if( x[l] < (double) arc->bbx0 ) arc->bbx0 = (int) x[l]; if( y[l] < (double) arc->bby0 ) arc->bby0 = (int) y[l]; if( x[l] > (double) arc->bbx1 ) arc->bbx1 = (int) x[l]; if( y[l] > (double) arc->bby1 ) arc->bby1 = (int) y[l]; if( fabs(dist(x[l],y[l],arc->xc,arc->yc) - arc->radius) > sigma ) return FALSE; /* not smooth */ } } /* correct bounding box */ arc->bbx0 -= max_w; /* add the width of the max lateral regions */ arc->bby0 -= max_w; arc->bbx1 += max_w + 1; arc->bby1 += max_w + 1; if( arc->bbx0 < 0 ) arc->bbx0 = 0; /* keep the bounding box inside image */ if( arc->bby0 < 0 ) arc->bby0 = 0; if( arc->bbx1 > X ) arc->bbx1 = X; if( arc->bby1 > Y ) arc->bby1 = Y; return TRUE; /* is smooth */ } /*----------------------------------------------------------------------------*/ /* get the lateral regions to an arc operator: the sub-case of a line segment input: pointer to an image of size X,Y, a pointer to the arc operator, and the width of lateral region w. output: the array pointed by reg will be filled with the elements of the regions. its size is returned in n. */ static void get_region_line( int * n, struct region * reg, double * image, int X, int Y, struct arc_of_circle * arc, double w ) { int x,y; /* check input */ if( image==NULL || X<=0 || Y<=0 ) error("get_region_line: invalid image"); if( arc==NULL ) error("get_region_line: invalid arc"); /* count points */ *n = 0; for(x=arc->bbx0; xbbx1; x++) for(y=arc->bby0; ybby1; y++) { double d_lat = arc->a*x + arc->b*y + arc->c; double d_lon = arc->b*x - arc->a*y + arc->d; if( fabs(d_lon) <= 0.5*arc->len && fabs(d_lat) <= w ) { if( d_lat < 0.0 ) reg[*n].reg = 1; else reg[*n].reg = 2; reg[*n].val = image[x+y*X]; reg[*n].w = fabs(d_lat); (*n)++; } } } /*----------------------------------------------------------------------------*/ /* get the lateral regions to an arc operator: the sub-case of an arc of circle input: pointer to an image of size X,Y, a pointer to the arc operator, and the width of lateral region w. output: the array pointed by reg will be filled with the elements of the regions. its size is returned in n. */ static void get_region_arc( int * n, struct region * reg, double * image, int X, int Y, struct arc_of_circle * arc, double w ) { int x,y; /* check input */ if( image==NULL || X<=0 || Y<=0 ) error("get_region_arc: invalid image"); if( arc==NULL || arc->is_line_segment ) error("get_region_arc: invalid arc"); /* count points */ *n = 0; for(x=arc->bbx0; xbbx1; x++) for(y=arc->bby0; ybby1; y++) { double r = dist(arc->xc, arc->yc, (double) x, (double) y); double offset = r - arc->radius; double ang = atan2(y - arc->yc, x - arc->xc); double ang_diff = diff_0_2pi(ang,arc->ang_ref); /* is the point in the angle sector? */ if( ang_diff >= 0.0 && ang_diff <= arc->ang_span && fabs(offset) <= w ) { if( (offset < 0.0 && arc->dir < 0) || (offset > 0.0 && arc->dir > 0) ) reg[*n].reg = 1; else reg[*n].reg = 2; reg[*n].val = image[x+y*X]; reg[*n].w = fabs(offset); (*n)++; } } } /*----------------------------------------------------------------------------*/ /* comparison function for type 'struct region' to be used with qsort */ static int comp_region(const void * a, const void * b) { if( ( (struct region *) a)->val < ( (struct region *) b)->val ) return -1; if( ( (struct region *) a)->val > ( (struct region *) b)->val ) return 1; return 0; } /*----------------------------------------------------------------------------*/ /* get the lateral regions to an arc operator, apply the correcting term for pixel quantization, and sort its elements according to pixel value input: pointer to an image of size X,Y, a pointer to the arc operator, the width of lateral region w, and the pixel quantization step Q. output: the array pointed by reg will be filled with the elements of the regions. its size returned in n. */ static void get_region( int * n, struct region * reg, double * image, int X, int Y, struct arc_of_circle * arc, double w, double Q ) { double q_offset = 0.616793 * Q; int i; /* collect pixel values */ if( arc->is_line_segment ) get_region_line(n,reg,image,X,Y,arc,w); else get_region_arc(n,reg,image,X,Y,arc,w); /* compensate pixel values quantization */ for(i=0; i<*n; i++) if( reg[i].reg == 1 ) reg[i].val += q_offset; /* sort region by value */ qsort( (void *) reg, (size_t) *n, sizeof(struct region), &comp_region ); } /*----------------------------------------------------------------------------*/ /* compute arc NFA using Mann-Whitney U test http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test input: X,Y is the image size, arc is a pointer to the arc, w is the width of the lateral regions of the current arc operator, reg is a pointer to the lateral regions of size n, W is the number of width used in the method, and gap is the size in pixel units of the space to separate both regions. return: log(NFA) */ static double arc_log_nfa( int X, int Y, struct arc_of_circle * arc, double w, int n, struct region * reg, int W, double gap ) { int n1,n2,i,rank,sum_tied_ranks,num_tied,num_tied_r2; double tie_val,adjusted_rank,sum_rank_r2,u,m,s,z,pvalue; /* Number of Tests: NT = sqrt(XY) * XY * (4 + pi^2)/3 * l^2 * num_width XY : number of centers sqrt(XY) : number of arc lengths pi*l : number of orientations of arcs of length l (4/pi + pi)/3 * l : number of arcs of a given length, center, orientation; this includes from full circles of perimeter l, to straight line segments of length l. num_width : number of arc width considered (counting with or without central gap as different) */ double logNT = 1.5 * log10( (double) X ) + 1.5 * log10( (double) Y ) + log10(4.6232) + 2.0 * log10( (double) arc->len ) + log10( (double) W ); /* compute Mann-Whitney U statistic. the elements in reg are assumed to be ordered by increasing values. the assignation of the initial rank, the assignation of the adjusted rank, and the sum of ranks in region 2, all is done in one single loop. for this there are variables to keep the value of the current tied pixel group, its size, and how many belong to region 2. each new value evaluated can have the same value as in the current tied group or a larger one. if it is the same, it belongs to the tied group and the variables are updated. if it is larger, that means that the previous tied group is completed, its adjusted rank can be assigned, and the corresponding quantity added to the sum of ranks in region 2. the variables are updated to reflect the new tied group including initially only the new pixel. */ n1 = n2 = 0; sum_rank_r2 = 0.0; rank = 0; tie_val = reg[0].val; /* set tied pixel group for the first pixel in reg */ sum_tied_ranks = num_tied = num_tied_r2 = 0; for(i=0; i 0.5*gap && reg[i].w <= w ) /* evaluate only pixels inside the current width */ { if( greater(reg[i].val, tie_val) ) /* a new tied pixel group */ { /* compute the adjusted rank and assign the rank values of region 2 in the last tied pixel group */ if( num_tied > 0 && num_tied_r2 > 0 ) { adjusted_rank = (double) sum_tied_ranks / (double) num_tied; sum_rank_r2 += (double) num_tied_r2 * adjusted_rank; } /* initialize new tied pixel group */ tie_val = reg[i].val; sum_tied_ranks = num_tied = num_tied_r2 = 0; } ++rank; /* rank in the ordering among the pixels in the region */ sum_tied_ranks += rank; ++num_tied; /* count pixels in region 1 and 2 */ if( reg[i].reg == 1 ) ++n1; else { ++n2; ++num_tied_r2; } } if( num_tied > 0 && num_tied_r2 > 0 ) /* assign ranks of last tied group */ sum_rank_r2 += (double) num_tied_r2 * sum_tied_ranks / num_tied; u = sum_rank_r2 - 0.5 * n2 * (n2 + 1.0); /* compute u value */ /* compute z, a version of u with standard normal distribution, N(0,1) */ m = 0.5 * n1 * n2; s = sqrt( n1 * n2 * (n1+n2+1.0) / 12.0 ); if( n1 > 0 && n2 > 0 && s > 0.0 ) z = (u - m) / s; else return logNT; /* one of the regions has no pixel => not meaningful */ /* compute the p-value using the standard error function */ pvalue = 0.5 * ( 1.0 - erf_winitzki(z/sqrt(2.0)) ); if( pvalue <= 0.0 ) /* the p-value should always be larger than zero. this condition reveals a numeric overflow, due to a very very small p-value. then the arc is meaningful => return a negative log10(NFA) */ return (double) DBL_MIN_10_EXP; /* minimal negative exponent in doubles */ /* return the log10(NFA) */ return logNT + log10(pvalue); } /*----------------------------------------------------------------------------*/ /* copy the meaningful curves from structures xx,yy,curve,NN,MM to structures x,y,curve_limits,N,M input: xx,yy : list of point coordinates curve : curve limits in lists xx,yy NN : number of points MM : number of curves meaningful : status of each point (TRUE=meaningful; FALSE otherwise) output: x,y : list of point coordinates curve_limits : curve limits in lists x,y N : number of points M : number of curves the output arrays x,y and curve_limits will be allocated. */ static void keep_meaningful_curves( double ** x, double ** y, int * N, int ** curve_limits, int * M, double * xx, double * yy, int NN, int * curve, int MM, int * meaningful ) { int c,i; *N = *M = 0; /* initialize N and M to zero point, zero curve */ /* if the input is empty, return NULL pointers */ if( NN <= 0 ) { *x = *y = NULL; *curve_limits = NULL; return; } /* allocate memory worst case: each point is kept (n points so 2n coordinates) and each one is one chain (n+1 limits) */ *x = (double *) xmalloc( NN * sizeof(double) ); *y = (double *) xmalloc( NN * sizeof(double) ); *curve_limits = (int *) xmalloc( (NN+1) * sizeof(int) ); /* create list of meaningful curves */ (*curve_limits)[0] = 0; /* this is always right but needs to be initialized */ for(c=0; c=log_eps; min_l++) { /* most meaningful arc when all points are in the right order. that is, u gets its maximal possible value of n*n. then, z = (n*n - 0.5*n*n) / sqrt(.) = 0.5*n*n / sqrt(.) */ double n = min_l * max_w; double z = 0.5*n*n / sqrt(n*n*(n+n+1.0)/12.0); /* Number of Tests: NT = sqrt(XY) * XY * (4 + pi^2)/3 * l^2 * num_width XY : number of centers sqrt(XY) : number of arc lengths pi*l : number of orientations of arcs of length l (4/pi + pi)/3 * l : num. of arcs of a given length, center, orientation this includes from full circles of perimeter l, to straight line segments of length l num_width : number of arc width considered (counting with or without central gap as different) */ double logNT = 1.5 * log10( (double) X ) + 1.5 * log10( (double) Y ) + log10(4.6232) + 2.0 * log10( (double) min_l ) + log10( (double) W ); l_nfa = logNT + log10( 0.5 * (1.0 - erf_winitzki(z/sqrt(2.0))) ); } return min_l; } /*----------------------------------------------------------------------------*/ /* Smooth Contours is an algorithm for detecting smooth contours on digital images. The output contours are given as chained sub-pixel edge points. Input: image : the input image X,Y : the size of the input image Q : the pixel quantization step Output: x,y : lists of sub-pixel coordinates of edge points curve_limits : the limits of each curve in lists x and y N : number of edge points M : number of curves The input is a XxY graylevel image given as a pointer to an array of doubles such that image[x+y*X] is the value at coordinates x,y (for 0 <= x < X and 0 <= y < Y). The output are the chained edge points given as 3 allocated lists: x, y and curve_limits. Also the numbers N (size of lists x and y) and M (number of curves). x[i] and y[i] (0<=i=min_l; k--) /* last point of curve segment */ if( !used[i] || !used[k] ) /* test segment only if one end not used */ if( smooth_segment(&arc,xx,yy,i,k,sigma,max_w,X,Y) ) { /* heuristic to speed up: mark as used center part of segment */ for(n=i+3; n<=(k-3); n++) used[n] = TRUE; /* get pixel values in lateral regions to the arc. to speed-up it is done only once, thus outside w-loop and arc_log_nfa() */ get_region(®_n, reg, diff, X, Y, &arc, max_w, Q); if( reg_n <= 0 ) continue; /* empty region */ /* lateral width loop. two gaps are tried for each: 0.0 and 1.0 */ for(w=min_w; w<=max_w; w*=fac_w) if( arc_log_nfa(X,Y,&arc,w,reg_n,reg,W,0.0) < log_eps || arc_log_nfa(X,Y,&arc,w,reg_n,reg,W,1.0) < log_eps ) { for(n=i; n<=k; n++) meaningful[n] = used[n] = TRUE; break; /* arc already meaningful, no need to test other w */ } } keep_meaningful_curves(x,y,N,curve_limits,M,xx,yy,NN,curve,MM,meaningful); /* free memory */ free( (void *) gauss ); free( (void *) diff ); free( (void *) meaningful ); free( (void *) used ); free( (void *) xx ); free( (void *) yy ); free( (void *) curve ); free( (void *) reg ); } /*----------------------------------------------------------------------------*/