https://doi.org/10.5201/ipol.2017.216
Raw File
Tip revision: 604958267cd62092c9872fbd9828949be88ae114 authored by Software Heritage on 10 October 2017, 00:00 UTC
ipol: Deposit 1314 in collection ipol
Tip revision: 6049582
devernay.c
/*----------------------------------------------------------------------------

  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 );
}
/*----------------------------------------------------------------------------*/
back to top