www.pudn.com > roadextr.rar > canny.cpp
/* ---------------------------------------------------------------------------------
"Canny" edge detector code:
---------------------------
This text file contains the source code for a "Canny" edge detector. It
was written by Mike Heath (heath@csee.usf.edu) using some pieces of a
Canny edge detector originally written by someone at Michigan State
University.
There are three 'C' source code files in this text file. They are named
"canny_edge.c", "hysteresis.c" and "pgm_io.c". They were written and compiled
under SunOS 4.1.3. Since then they have also been compiled under Solaris.
To make an executable program: (1) Separate this file into three files with
the previously specified names, and then (2) compile the code using
gcc -o canny_edge canny_edge.c hysteresis.c pgm_io.c -lm
(Note: You can also use optimization such as -O3)
The resulting program, canny_edge, will process images in the PGM format.
Parameter selection is left up to the user. A broad range of parameters to
use as a starting point are: sigma 0.60-2.40, tlow 0.20-0.50 and,
thigh 0.60-0.90.
If you are using a Unix system, PGM file format conversion tools can be found
at ftp://wuarchive.wustl.edu/graphics/graphics/packages/pbmplus/.
Otherwise, it would be easy for anyone to rewrite the image I/O procedures
because they are listed in the separate file pgm_io.c.
If you want to check your compiled code, you can download grey-scale and edge
images from http://marathon.csee.usf.edu/edge/edge_detection.html. You can use
the parameters given in the edge filenames and check whether the edges that
are output from your program match the edge images posted at that address.
Mike Heath
(10/29/96)
<------------------------- begin canny_edge.c ------------------------->
*/
/*******************************************************************************
* PROGRAM: canny_edge
* PURPOSE: This program implements a "Canny" edge detector. The processing
* steps are as follows:
*
* 1) Convolve the image with a separable gaussian filter.
* 2) Take the dx and dy the first derivatives using [-1,0,1] and [1,0,-1]'.
* 3) Compute the magnitude: sqrt(dx*dx+dy*dy).
* 4) Perform non-maximal suppression.
* 5) Perform hysteresis.
*
* The user must input three parameters. These are as follows:
*
* sigma = The standard deviation of the gaussian smoothing filter.
* tlow = Specifies the low value to use in hysteresis. This is a
* fraction (0-1) of the computed high threshold edge strength value.
* thigh = Specifies the high value to use in hysteresis. This fraction (0-1)
* specifies the percentage point in a histogram of the gradient of
* the magnitude. Magnitude values of zero are not counted in the
* histogram.
*
* NAME: Mike Heath
* Computer Vision Laboratory
* University of South Floeida
* heath@csee.usf.edu
*
* DATE: 2/15/96
*
* Modified: 5/17/96 - To write out a floating point RAW headerless file of
* the edge gradient "up the edge" where the angle is
* defined in radians counterclockwise from the x direction.
* (Mike Heath)
******************************************************************************
*/
#include "stdafx.h"
#include
#include
#include
#include
#include "canny.h"
#define VERBOSE 0
#define BOOSTBLURFACTOR 90.0
#define VERBOSE 0
#define NOEDGE 255
#define POSSIBLE_EDGE 128
#define EDGE 0
#define M_PI 3.1415926
/*******************************************************************************
* PROCEDURE: canny
* PURPOSE: To perform canny edge detection.
* NAME: Mike Heath
* DATE: 2/15/96
*******************************************************************************/
void canny(unsigned char *image, int rows, int cols, float sigma,
float tlow, float thigh, unsigned char *edge)
//, char *fname)
{
FILE *fpdir=NULL; ///* File to write the gradient image to.
unsigned char *nms; ///* Points that are local maximal magnitude.
short int *smoothedim, ///* The image after gaussian smoothing.
*delta_x, ///* The first devivative image, x-direction.
*delta_y, ////* The first derivative image, y-direction.
*magnitude; ////* The magnitude of the gadient image.
int r, c, pos;
float *dir_radians=NULL; ////* Gradient direction image.
//* Perform gaussian smoothing on the image using the input standard
//* deviation.
gaussian_smooth(image, rows, cols, sigma, &smoothedim);
/// * Compute the first derivative in the x and y directions.
derrivative_x_y(smoothedim, rows, cols, &delta_x, &delta_y);
///* This option to write out the direction of the edge gradient was added
///* to make the information available for computing an edge quality figure
///* of merit.
if(1) //fname != NULL
{
///* Compute the direction up the gradient, in radians that are
///* specified counteclockwise from the positive x-axis.
radian_direction(delta_x, delta_y, rows, cols, &dir_radians, -1, -1);
free(dir_radians);
}
///// Compute the magnitude of the gradient.
magnitude_x_y(delta_x, delta_y, rows, cols, &magnitude);
/////////// for finding ridge
// for(r=0;r= 0){
if(y >= 0) return(ang);
else return(2*M_PI - ang);
}
else{
if(y >= 0) return(M_PI - ang);
else return(M_PI + ang);
}
}
/*******************************************************************************
* PROCEDURE: magnitude_x_y
* PURPOSE: Compute the magnitude of the gradient. This is the square root of
* the sum of the squared derivative values.
* NAME: Mike Heath
* DATE: 2/15/96
*******************************************************************************/
void magnitude_x_y(short int *delta_x, short int *delta_y, int rows, int cols,
short int **magnitude)
{
int r, c, pos, sq1, sq2;
/****************************************************************************
* Allocate an image to store the magnitude of the gradient.
****************************************************************************/
if((*magnitude = (short *) calloc(rows*cols, sizeof(short))) == NULL){
fprintf(stderr, "Error allocating the magnitude image.\n");
exit(1);
}
for(r=0,pos=0;r= 0) && ((c+cc) < cols)){
dot += (float)image[r*cols+(c+cc)] * kernel[center+cc];
sum += kernel[center+cc];
}
}
tempim[r*cols+c] = dot/sum;
}
}
/****************************************************************************
* Blur in the y - direction.
****************************************************************************/
// if(VERBOSE) printf(" Bluring the image in the Y-direction.\n");
for(c=0;c= 0) && ((r+rr) < rows)){
dot += tempim[(r+rr)*cols+c] * kernel[center+rr];
sum += kernel[center+rr];
}
}
(*smoothedim)[r*cols+c] = (short int)(dot*BOOSTBLURFACTOR/sum + 0.5);
}
}
free(tempim);
free(kernel);
}
/*******************************************************************************
* PROCEDURE: gaussian_smooth
* PURPOSE: Blur an image with a gaussian filter.
* NAME: Mike Heath
* DATE: 2/15/96
*******************************************************************************/
void gaussian_smoothImg(unsigned char *image, int rows, int cols, float sigma)
{
int r, c, rr, cc, /* Counter variables. */
windowsize, /* Dimension of the gaussian kernel. */
center; /* Half of the windowsize. */
float *tempim, /* Buffer for separable filter gaussian smoothing. */
*kernel, /* A one dimensional gaussian kernel. */
dot, /* Dot product summing variable. */
sum; /* Sum of the kernel weights variable. */
/****************************************************************************
* Create a 1-dimensional gaussian smoothing kernel.
****************************************************************************/
make_gaussian_kernel(sigma, &kernel, &windowsize);
center = windowsize / 2;
/****************************************************************************
* Allocate a temporary buffer image and the smoothed image.
****************************************************************************/
if((tempim = (float *) calloc(rows*cols, sizeof(float))) == NULL)
{
fprintf(stderr, "Error allocating the buffer image.\n");
exit(1);
}
/****************************************************************************
* Blur in the x - direction.
****************************************************************************/
// if(VERBOSE) printf(" Bluring the image in the X-direction.\n");
for(r=0;r= 0) && ((c+cc) < cols))
{
dot += (float)image[r*cols+(c+cc)] * kernel[center+cc];
sum += kernel[center+cc];
}
}
tempim[r*cols+c] = dot/sum;
}
}
/****************************************************************************
* Blur in the y - direction.
****************************************************************************/
// if(VERBOSE) printf(" Bluring the image in the Y-direction.\n");
for(c=0;c= 0) && ((r+rr) < rows))
{
dot += tempim[(r+rr)*cols+c] * kernel[center+rr];
sum += kernel[center+rr];
}
}
(image)[r*cols+c] = (BYTE)(dot/sum + 0.5);
}
}
free(tempim);
free(kernel);
}
/****************************************************************************
* PROCEDURE: make_gaussian_kernel
* PURPOSE: Create a one dimensional gaussian kernel.
* NAME: Mike Heath
* DATE: 2/15/96
*****************************************************************************/
void make_gaussian_kernel(float sigma, float **kernel, int *windowsize)
{
int i, center;
float x, fx, sum=0.0;
*windowsize = 1 + 2 * ceil(2.5 * sigma);
center = (*windowsize) / 2;
if(VERBOSE) printf(" The kernel has %d elements.\n", *windowsize);
if((*kernel = (float *) calloc((*windowsize), sizeof(float))) == NULL){
fprintf(stderr, "Error callocing the gaussian kernel array.\n");
exit(1);
}
for(i=0;i<(*windowsize);i++){
x = (float)(i - center);
fx = pow(2.71828, -0.5*x*x/(sigma*sigma)) / (sigma * sqrt(6.2831853));
(*kernel)[i] = fx;
sum += fx;
}
for(i=0;i<(*windowsize);i++) (*kernel)[i] /= sum;
if(VERBOSE){
printf("The filter coefficients are:\n");
for(i=0;i<(*windowsize);i++)
printf("kernel[%d] = %f\n", i, (*kernel)[i]);
}
}
//// ------------------------- end canny_edge.c ------------------------->
///// ------------------------- begin hysteresis.c ------------------------->
/*******************************************************************************
* FILE: hysteresis.c
* This code was re-written by Mike Heath from original code obtained indirectly
* from Michigan State University. heath@csee.usf.edu (Re-written in 1996).
*******************************************************************************/
/*******************************************************************************
* PROCEDURE: follow_edges
* PURPOSE: This procedure edges is a recursive routine that traces edgs along
* all paths whose magnitude values remain above some specifyable lower
* threshhold.
* NAME: Mike Heath
* DATE: 2/15/96
*******************************************************************************/
follow_edges(unsigned char *edgemapptr, short *edgemagptr, short lowval,
int cols)
{
short *tempmagptr;
unsigned char *tempmapptr;
int i;
float thethresh;
int x[8] = {1,1,0,-1,-1,-1,0,1},
y[8] = {0,1,1,1,0,-1,-1,-1};
for(i=0;i<8;i++){
tempmapptr = edgemapptr - y[i]*cols + x[i];
tempmagptr = edgemagptr - y[i]*cols + x[i];
if((*tempmapptr == POSSIBLE_EDGE) && (*tempmagptr > lowval)){
*tempmapptr = (unsigned char) EDGE;
follow_edges(tempmapptr,tempmagptr, lowval, cols);
}
}
}
/*******************************************************************************
* PROCEDURE: apply_hysteresis
* PURPOSE: This routine finds edges that are above some high threshhold or
* are connected to a high pixel by a path of pixels greater than a low
* threshold.
* NAME: Mike Heath
* DATE: 2/15/96
*******************************************************************************/
void apply_hysteresis(short int *mag, unsigned char *nms, int rows, int cols,
float tlow, float thigh, unsigned char *edge)
{
int r, c, pos, numedges, lowcount, highcount, lowthreshold, highthreshold,
i, hist[32768], rr, cc;
short int maximum_mag, sumpix;
/****************************************************************************
* Initialize the edge map to possible edges everywhere the non-maximal
* suppression suggested there could be an edge except for the border. At
* the border we say there can not be an edge because it makes the
* follow_edges algorithm more efficient to not worry about tracking an
* edge off the side of the image.
****************************************************************************/
for(r=0,pos=0;r= highthreshold)){
edge[pos] = EDGE;
follow_edges((edge+pos), (mag+pos), lowthreshold, cols);
}
}
}
/****************************************************************************
* Set all the remaining possible edges to non-edges.
****************************************************************************/
for(r=0,pos=0;r= 0){
if(gy >= 0){
if (gx >= gy)
{
/* 111 */
/* Left point */
z1 = *(magptr - 1);
z2 = *(magptr - ncols - 1);
mag1 = (m00 - z1)*xperp + (z2 - z1)*yperp;
/* Right point */
z1 = *(magptr + 1);
z2 = *(magptr + ncols + 1);
mag2 = (m00 - z1)*xperp + (z2 - z1)*yperp;
}
else
{
/* 110 */
/* Left point */
z1 = *(magptr - ncols);
z2 = *(magptr - ncols - 1);
mag1 = (z1 - z2)*xperp + (z1 - m00)*yperp;
/* Right point */
z1 = *(magptr + ncols);
z2 = *(magptr + ncols + 1);
mag2 = (z1 - z2)*xperp + (z1 - m00)*yperp;
}
}
else
{
if (gx >= -gy)
{
/* 101 */
/* Left point */
z1 = *(magptr - 1);
z2 = *(magptr + ncols - 1);
mag1 = (m00 - z1)*xperp + (z1 - z2)*yperp;
/* Right point */
z1 = *(magptr + 1);
z2 = *(magptr - ncols + 1);
mag2 = (m00 - z1)*xperp + (z1 - z2)*yperp;
}
else
{
/* 100 */
/* Left point */
z1 = *(magptr + ncols);
z2 = *(magptr + ncols - 1);
mag1 = (z1 - z2)*xperp + (m00 - z1)*yperp;
/* Right point */
z1 = *(magptr - ncols);
z2 = *(magptr - ncols + 1);
mag2 = (z1 - z2)*xperp + (m00 - z1)*yperp;
}
}
}
else
{
if ((gy = *gyptr) >= 0)
{
if (-gx >= gy)
{
/* 011 */
/* Left point */
z1 = *(magptr + 1);
z2 = *(magptr - ncols + 1);
mag1 = (z1 - m00)*xperp + (z2 - z1)*yperp;
/* Right point */
z1 = *(magptr - 1);
z2 = *(magptr + ncols - 1);
mag2 = (z1 - m00)*xperp + (z2 - z1)*yperp;
}
else
{
/* 010 */
/* Left point */
z1 = *(magptr - ncols);
z2 = *(magptr - ncols + 1);
mag1 = (z2 - z1)*xperp + (z1 - m00)*yperp;
/* Right point */
z1 = *(magptr + ncols);
z2 = *(magptr + ncols - 1);
mag2 = (z2 - z1)*xperp + (z1 - m00)*yperp;
}
}
else
{
if (-gx > -gy)
{
/* 001 */
/* Left point */
z1 = *(magptr + 1);
z2 = *(magptr + ncols + 1);
mag1 = (z1 - m00)*xperp + (z1 - z2)*yperp;
/* Right point */
z1 = *(magptr - 1);
z2 = *(magptr - ncols - 1);
mag2 = (z1 - m00)*xperp + (z1 - z2)*yperp;
}
else
{
/* 000 */
/* Left point */
z1 = *(magptr + ncols);
z2 = *(magptr + ncols + 1);
mag1 = (z2 - z1)*xperp + (m00 - z1)*yperp;
/* Right point */
z1 = *(magptr - ncols);
z2 = *(magptr - ncols - 1);
mag2 = (z2 - z1)*xperp + (m00 - z1)*yperp;
}
}
}
/* Now determine if the current point is a maximum point */
if ((mag1 > 0.0) || (mag2 > 0.0))
{
*resultptr = (unsigned char) NOEDGE;
}
else
{
if (mag2 == 0.0)
*resultptr = (unsigned char) NOEDGE;
else
*resultptr = (unsigned char) POSSIBLE_EDGE;
}
}
}
}
//// ------------------------- end hysteresis.c ------------------------->