/*----------------------------------------------------------------------------
Implementation of Canny/Devernay's sub-pixel edge detector. This code is part
of the following publication and was subject to peer review:
"A Sub-Pixel Edge Detector: an Implementation of the Canny/Devernay
Algorithm" by Rafael Grompone von Gioi and Gregory Randall,
Image Processing On Line, 2017. DOI:10.5201/ipol.2017.216
http://dx.doi.org/10.5201/ipol.2017.216
Copyright (c) 2016-2017 rafael grompone von gioi <grompone@gmail.com>,
Gregory Randall <randall@fing.edu.uy>
This program 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 <http://www.gnu.org/licenses/>.
----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
/*----------------------------------------------------------------------------*/
#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 */
/*----------------------------------------------------------------------------*/
/* 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 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<n; i++)
{
val = ( (double) i - mean ) / sigma;
kernel[i] = exp( -0.5 * val * val );
sum += kernel[i];
}
/* normalization */
if( sum > 0.0 ) for(i=0; i<n; i++) kernel[i] /= sum;
}
/*----------------------------------------------------------------------------*/
/* filter an image with a Gaussian kernel of parameter sigma. return a pointer
to a newly allocated filtered image, of the same size as the input image.
*/
static double * gaussian_filter(double * image, int X, int Y, double sigma)
{
int x,y,offset,i,j,nx2,ny2,n;
double * kernel;
double * tmp;
double * out;
double val,prec;
/* check input */
if( sigma <= 0.0 ) error("gaussian_filter: sigma must be positive");
if( image == NULL || X < 1 || Y < 1 ) error("gaussian_filter: invalid image");
/* get memory */
tmp = (double *) xmalloc( X * Y * sizeof(double) );
out = (double *) xmalloc( X * Y * sizeof(double) );
/* compute gaussian kernel */
/*
The size of the kernel is selected to guarantee that the first discarded
term is at least 10^prec times smaller than the central value. For that,
the half size of the kernel must be larger than x, with
e^(-x^2/2sigma^2) = 1/10^prec
Then,
x = sigma * sqrt( 2 * prec * ln(10) )
*/
prec = 3.0;
offset = (int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
n = 1 + 2 * offset; /* kernel size */
kernel = (double *) xmalloc( n * sizeof(double) );
gaussian_kernel(kernel, n, sigma, (double) offset);
/* auxiliary variables for the double of the image size */
nx2 = 2*X;
ny2 = 2*Y;
/* x axis convolution */
for(x=0; x<X; x++)
for(y=0; y<Y; y++)
{
val = 0.0;
for(i=0; i<n; i++)
{
j = x - offset + i;
/* symmetry boundary condition */
while(j<0) j += nx2;
while(j>=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<X; x++)
for(y=0; y<Y; y++)
{
val = 0.0;
for(i=0; i<n; i++)
{
j = y - offset + i;
/* symmetry boundary condition */
while(j<0) j += ny2;
while(j>=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;
}
/*----------------------------------------------------------------------------*/
/* 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, not a valid chaining */
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, not a valid chaining */
/* 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; /* incompatible gradient angles, not a valid chaining */
/* 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.
*/
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<X*Y; i++) Ex[i] = Ey[i] = -1.0;
/* explore pixels inside a 2 pixel margin (so modG[x,y +/- 1,1] is defined) */
for(x=2; x<(X-2); x++)
for(y=2; y<(Y-2); y++)
{
int Dx = 0; /* interpolation will be along Dx,Dy */
int Dy = 0; /* which will be selected below */
double mod = modG[x+y*X]; /* modG at pixel */
double L = modG[x-1+y*X]; /* modG at pixel on the left */
double R = modG[x+1+y*X]; /* modG at pixel on the right */
double U = modG[x+(y+1)*X]; /* modG at pixel up */
double D = modG[x+(y-1)*X]; /* modG at pixel below */
double gx = fabs( Gx[x+y*X] ); /* absolute value of Gx */
double gy = fabs( Gy[x+y*X] ); /* absolute value of Gy */
/* when local horizontal maxima of the gradient modulus and the
gradient direction is more horizontal (|Gx| >= |Gy|),
=> a "horizontal" (H) edge found
else, if local vertical maxima of the gradient modulus and the
gradient direction is more vertical (|Gx| <= |Gy|),
=> a "vertical" (V) edge found */
/* it can happen that two neighbor pixels have equal value and are both
maxima, for example when the edge is exactly between both pixels. in
such cases, as an arbitrary convention, the edge is marked on the
left one when an horizontal max or below when a vertical max. for
this the conditions are L < mod >= 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 */
if( greater(mod,L) && !greater(R,mod) && gx >= gy ) Dx = 1; /* H */
else if( greater(mod,D) && !greater(U,mod) && gx <= gy ) Dy = 1; /* V */
/* 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 in 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<X*Y; i++) next[i] = prev[i] = -1;
/* try each point to make local chains */
for(x=2; x<(X-2); x++) /* 2 pixel margin to include the tested neighbors */
for(y=2; y<(Y-2); y++)
if( Ex[x+y*X] >= 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.
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)
or equivalently ||a-b|| > ||b-c|| > ||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 the chaining 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 exists already 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 two chainings: 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 */
}
}
}
/*----------------------------------------------------------------------------*/
/* apply Canny thresholding with hysteresis
next and prev contain the number of next and previous edge points in the
chain or -1 when not chained. modG is modulus of the image gradient. X,Y is
the image size. th_h and th_l are the high and low thresholds, respectively.
this function modifies next and prev, removing chains not satisfying the
thresholds.
*/
static void thresholds_with_hysteresis( int * next, int * prev,
double * modG, int X, int Y,
double th_h, double th_l )
{
int * valid;
int i,j,k;
/* check input */
if( next == NULL || prev == NULL || modG == NULL )
error("thresholds_with_hysteresis: invalid input");
/* get memory */
valid = (int *) xmalloc( X * Y * sizeof(int) );
for(i=0; i<X*Y; i++) valid[i] = FALSE;
/* validate all edge points over th_h or connected to them and over th_l */
for(i=0; i<X*Y; i++) /* prev[i]>=0 or next[i]>=0 implies an edge point */
if( (prev[i] >= 0 || next[i] >= 0) && !valid[i] && modG[i] >= th_h )
{
valid[i] = TRUE; /* mark as valid the new point */
/* follow the chain of edge points forwards */
for(j=i; j>=0 && (k=next[j])>=0 && !valid[k]; j=next[j])
if( modG[k] < th_l )
{
next[j] = -1; /* cut the chain when the point is below th_l */
prev[k] = -1; /* j must be assigned to next[j] and not k,
so the loop is chained in this case */
}
else
valid[k] = TRUE; /* otherwise mark the new point as valid */
/* follow the chain of edge points backwards */
for(j=i; j>=0 && (k=prev[j])>=0 && !valid[k]; j=prev[j])
if( modG[k] < th_l )
{
prev[j] = -1; /* cut the chain when the point is below th_l */
next[k] = -1; /* j must be assigned to prev[j] and not k,
so the loop is chained in this case */
}
else
valid[k] = TRUE; /* otherwise mark the new point as valid */
}
/* remove any remaining non-valid chained point */
for(i=0; i<X*Y; i++) /* prev[i]>=0 or next[i]>=0 implies edge point */
if( (prev[i] >= 0 || next[i] >= 0) && !valid[i] )
prev[i] = next[i] = -1;
/* free memory */
free( (void *) valid );
}
/*----------------------------------------------------------------------------*/
/* 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<N) store the sub-pixel coordinates of the edge points.
curve_limits[j] (0<=j<=M) stores the limits of each chain in lists x and y.
x, y, and curve_limits will be allocated.
example:
curve number k (0<=k<M) consists of the edge points x[i],y[i]
for i determined by curve_limits[k] <= i < curve_limits[k+1].
curve k is closed if x[curve_limits[k]] == x[curve_limits[k+1] - 1] and
y[curve_limits[k]] == y[curve_limits[k+1] - 1].
*/
static void list_chained_edge_points( double ** x, double ** y, int * N,
int ** curve_limits, int * M,
int * next, int * prev,
double * Ex, double * Ey, int X, int Y )
{
int i,k,n;
/* initialize output: x, y, curve_limits, N, and M
there cannot be more than X*Y edge points to be put in the output list:
edge points must be local maxima of gradient modulus, so at most half of
the pixels could be so. when a closed curve is found, one edge point will
be put twice to the output. even if all possible edge points (half of the
pixels in the image) would form one pixel closed curves (which is not
possible) that would lead to output X*Y edge points.
for the same reason, there cannot be more than X*Y curves: the worst case
is when all possible edge points (half of the pixels in the image) would
form one pixel chains. in that case (which is not possible) one would need
a size for curve_limits of X*Y/2+1. so X*Y is enough.
(curve_limits requires one more item than the number of curves.
a simplest example is when only one chain of length 3 is present:
curve_limits[0] = 0, curve_limits[1] = 3.)
*/
*x = (double *) xmalloc( X * Y * sizeof(double) );
*y = (double *) xmalloc( X * Y * sizeof(double) );
*curve_limits = (int *) xmalloc( X * Y * sizeof(int) );
*N = 0;
*M = 0;
/* copy chained edge points to output */
for(i=0; i<X*Y; i++) /* prev[i]>=0 or next[i]>=0 implies an 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 beginning 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.
input:
image : the input image
X,Y : the size of the input image
sigma : standard deviation sigma for the Gaussian filtering
(if sigma=0 no filtering is performed)
th_h : high gradient threshold in Canny's hysteresis
th_l : low gradient threshold in Canny's hysteresis
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<N) store the sub-pixel coordinates of the edge points.
curve_limits[j] (0<=j<=M) stores the limits of each chain in lists x and y.
example:
curve number k (0<=k<M) consists of the edge points x[i],y[i]
for i determined by curve_limits[k] <= i < curve_limits[k+1].
curve k is closed if x[curve_limits[k]] == x[curve_limits[k+1] - 1] and
y[curve_limits[k]] == y[curve_limits[k+1] - 1].
*/
void devernay( double ** x, double ** y, int * N, int ** curve_limits, int * M,
double * image, int X, int Y,
double sigma, double th_h, double th_l )
{
double * Gx = (double *) xmalloc( X * Y * sizeof(double) ); /* grad_x */
double * Gy = (double *) xmalloc( X * Y * sizeof(double) ); /* grad_y */
double * modG = (double *) xmalloc( X * Y * sizeof(double) ); /* |grad| */
double * Ex = (double *) xmalloc( X * Y * sizeof(double) ); /* edge_x */
double * Ey = (double *) xmalloc( X * Y * sizeof(double) ); /* edge_y */
int * next = (int *) xmalloc( X * Y * sizeof(int) ); /* next point in chain */
int * prev = (int *) xmalloc( X * Y * sizeof(int) ); /* prev point in chain */
double * gauss = NULL;
if( sigma == 0.0 ) compute_gradient(Gx, Gy, modG, image, X, Y);
else
{
gauss = gaussian_filter(image, X, Y, sigma);
compute_gradient(Gx, Gy, modG, gauss, X, Y);
}
compute_edge_points(Ex, Ey, modG, Gx, Gy, X, Y);
chain_edge_points(next, prev, Ex, Ey, Gx, Gy, X, Y);
thresholds_with_hysteresis(next, prev, modG, X, Y, th_h, th_l);
list_chained_edge_points(x, y, N, curve_limits, M, next, prev, Ex, Ey, X, Y);
/* free memory */
if( gauss != NULL ) free( (void *) gauss );
free( (void *) Gx );
free( (void *) Gy );
free( (void *) modG );
free( (void *) Ex );
free( (void *) Ey );
free( (void *) next );
free( (void *) prev );
}
/*----------------------------------------------------------------------------*/
