www.pudn.com > imageFeatures_Ver3_0_source.zip > canny.cpp
//////////////////////////////////////////////////////////////////
// Canny.cpp
// Original Code from:
// Algorithms for Image Processing and Computer Vision
// By J.R. Parker
// Modified for VisSDK by: Travis A Udd 11/30/00
//////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include
#include
#include "canny.h"
#include "profile.h"
/* Scale floating point magnitudes and angles to 8 bits */
#define ORI_SCALE 40.0
#define MAG_SCALE 20.0
#define PI_f 3.14159f
/* Biggest possible filter mask */
#define MAX_MASK_SIZE 20
//////////////////////////
// Travis A Udd 11/30/00
// Revised 6/02/03 Tyler Folsom
// Canny Constructor
//////////////////////////
Canny::Canny(CDefaultImage *imageSrc)
{
ratio = 0.1f;
WIDTH = 0;
m_pImage = imageSrc;
}
//////////////////////////////////////////////////////////////////
// CannyEdgeDetect()
// Travis A Udd 11/30/00
// Function: Fills 2D buffer with image
// and runs canny filter function on it.
// result is edges found in image.
//////////////////////////////////////////////////////////////////
void Canny::CannyEdgeDetect(CVisRGBAByteImage *featureImage)
{
PROFILE("CannyEdgeDetect");
// There is a memory leak somewhere! 06/02/03 Tyler Folsom
imagec *im;
int iWidth = m_pImage->Width(),
iHeight = m_pImage->Height(),
i, j, total;
CVisRGBA temp;
im=newimagec (iHeight,iWidth);
//Convert the RGB image to grayscale and load it into buffer
for(i=0 ; iPixel(j,i);
total = (temp.R()+temp.G()+temp.B()) / 3;
im->data[i][j] = total;
}
// Do canny filter on image
CannyFilter(im, HighTh,
LowTh,
Gaussian);
// Draw image
for(i=0 ; idata[i][j])
featureImage->Pixel(j,i) = m_pImage->Pixel(j,i);
else
featureImage->Pixel(j,i).SetRGB(255,0,0);
}
freeimagec (im);
}
float ** Canny::f2d (int nr, int nc)
{
float **x;
int i;
x = (float **)calloc ( nr, sizeof (float *) );
if (x == 0)
{
fprintf (stderr, "Out of storage: F2D.\n");
exit (1);
}
for (i=0; iinfo = (struct header *)malloc( sizeof(struct header) );
if (!(x->info)) {
printf ("Out of storage in NEWIMAGE.\n");
return 0;
}
x->info->nr = nr; x->info->nc = nc;
x->info->oi = x->info->oj = 0;
/* Allocate the pixel array */
x->data = (unsigned char **)malloc(sizeof(unsigned char *)*nr);
/* Pointers to rows */
if (!(x->data)) {
printf ("Out of storage in NEWIMAGE.\n");
return 0;
}
for (i=0; idata[i] = (unsigned char *)malloc(nc*sizeof(unsigned char));
if (x->data[i]==0)
{
printf ("Out of storage. Newimage - row %d\n", i);
exit(1);
}
}
return x;
}
void Canny::freeimagec (struct imagec *z)
{
/* Free the storage associated with the image Z */
int i;
if (z != 0)
{
for (i=0; iinfo->nr; i++)
free (z->data[i]);
free (z->info);
free (z->data);
free (z);
}
}
//////////////////////////////////////////////////////////////////
// CannyFilter()
// Original Code from:
// Algorithms for Image Processing and Computer Vision
// By J.R. Parker
// Modified for VisSDK by: Travis A Udd 11/30/00
// Function: Gets original image im and returns im canny filtered
// for display.
//////////////////////////////////////////////////////////////////
void Canny::CannyFilter(imagec *im,int low,int high,double gaussian)
{
PROFILE("CannyFilter");
int i,j;
float s;
imagec *magim, *oriim;
s=(float)gaussian;
magim = newimagec (im->info->nr, im->info->nc);
if (magim == NULL)
{
printf ("Out of storage: Magnitude\n");
exit (1);
}
oriim = newimagec (im->info->nr, im->info->nc);
if (oriim == NULL)
{
printf ("Out of storage: Orientation\n");
exit (1);
}
/* Apply the filter */
canny (s, im, magim, oriim);
/* Hysteresis thresholding of edge pixels */
hysteresis (high, low, im, magim, oriim);
for (i=0; iinfo->nc; j++)
im->data[i][j] = 255;
for (i=im->info->nr-1; i>im->info->nr-1-WIDTH; i--)
for (j=0; jinfo->nc; j++)
im->data[i][j] = 255;
for (i=0; iinfo->nr; i++)
for (j=0; jdata[i][j] = 255;
for (i=0; iinfo->nr; i++)
for (j=im->info->nc-WIDTH-1; jinfo->nc; j++)
im->data[i][j] = 255;
freeimagec(magim);
freeimagec(oriim);
}
float Canny::norm (float x, float y)
{
return (float) sqrt ( (double)(x*x + y*y) );
}
void Canny::canny (float s, imagec* im, imagec* mag, imagec* ori)
{
int width;
float **smx,**smy;
float **dx,**dy;
int i,j,n;
float gau[MAX_MASK_SIZE], dgau[MAX_MASK_SIZE], z;
/* Create a Gaussian and a derivative of Gaussian filter mask */
for(i=0; iinfo->nr, im->info->nc);
smy = f2d (im->info->nr, im->info->nc);
/* Convolution of source image with a Gaussian in X and Y directions */
seperable_convolution (im, gau, width, smx, smy);
/* Now convolve smoothed data with a derivative */
printf ("Convolution with the derivative of a Gaussian...\n");
dx = f2d (im->info->nr, im->info->nc);
dxy_seperable_convolution (smx, im->info->nr, im->info->nc,
dgau, width, dx, 1);
free(smx[0]); free(smx);
dy = f2d (im->info->nr, im->info->nc);
dxy_seperable_convolution (smy, im->info->nr, im->info->nc,
dgau, width, dy, 0);
free(smy[0]); free(smy);
/* Create an image of the norm of dx,dy */
for (i=0; iinfo->nr; i++)
for (j=0; jinfo->nc; j++)
{
z = norm (dx[i][j], dy[i][j]);
mag->data[i][j] = (unsigned char)(z*MAG_SCALE);
}
/* Non-maximum suppression - edge pixels should be a local max */
nonmax_suppress (dx, dy, (int)im->info->nr, (int)im->info->nc, mag, ori);
free(dx[0]); free(dx);
free(dy[0]); free(dy);
}
/* Gaussian */
float Canny::gauss(float x, float sigma)
{
float xx;
if (sigma == 0) return 0.0;
xx = (float)exp((double) ((-x*x)/(2*sigma*sigma)));
return xx;
}
float Canny::meanGauss (float x, float sigma)
{
float z;
z = (gauss(x,sigma)+gauss(x+0.5,sigma)+gauss(x-0.5,sigma))/3.0;
z = z/(PI_f*2.0*sigma*sigma);
return z;
}
/* First derivative of Gaussian */
float Canny::dGauss (float x, float sigma)
{
return -x/(sigma*sigma) * gauss(x, sigma);
}
/* HYSTERESIS thersholding of edge pixels. Starting at pixels with a
value greater than the HIGH threshold, trace a connected sequence
of pixels that have a value greater than the LOW threhsold. */
void Canny::hysteresis (int high, int low, imagec* im, imagec* mag, imagec* oriim)
{
int i,j;
printf ("Beginning hysteresis thresholding...\n");
for (i=0; iinfo->nr; i++)
for (j=0; jinfo->nc; j++)
im->data[i][j] = 0;
if (highinfo->nr; i++)
for (j=0; jinfo->nc; j++)
if (mag->data[i][j] >= high)
trace (i, j, low, im, mag, oriim);
/* Make the edge black (to be the same as the other methods) */
for (i=0; iinfo->nr; i++)
for (j=0; jinfo->nc; j++)
if (im->data[i][j] == 0) im->data[i][j] = 255;
else im->data[i][j] = 0;
}
/* Check that a pixel index is in range. Return TRUE(1) if so. */
int Canny::range (imagec* im, int i, int j)
{
if ((i<0) || (i>=im->info->nr)) return 0;
if ((j<0) || (j>=im->info->nc)) return 0;
return 1;
}
/* TRACE - recursively trace edge pixels that have a threshold > the low
edge threshold, continuing from the pixel at (i,j). */
int Canny::trace (int i, int j, int low, imagec* im,imagec* mag, imagec* ori)
{
int n,m;
char flag = 0;
if (im->data[i][j] == 0)
{
im->data[i][j] = 255;
flag=0;
for (n= -1; n<=1; n++)
{
for(m= -1; m<=1; m++)
{
if (i==0 && m==0) continue;
if (range(mag, i+n, j+m) && mag->data[i+n][j+m] >= low)
if (trace(i+n, j+m, low, im, mag, ori))
{
flag=1;
break;
}
}
if (flag) break;
}
return(1);
}
return(0);
}
void Canny::seperable_convolution (imagec* im, float *gau, int width,
float **smx, float **smy)
{
int i,j,k, I1, I2, nr, nc;
float x, y;
nr = im->info->nr;
nc = im->info->nc;
for (i=0; idata[i][j]; y = gau[0] * im->data[i][j];
for (k=1; kdata[I1][j] + gau[k]*im->data[I2][j];
I1 = (j+k)%nc; I2 = (j-k+nc)%nc;
x += gau[k]*im->data[i][I1] + gau[k]*im->data[i][I2];
}
smx[i][j] = x; smy[i][j] = y;
}
}
void Canny::dxy_seperable_convolution (float** im, int nr, int nc, float *gau,
int width, float **sm, int which)
{
int i,j,k, I1, I2;
float x;
for (i=0; iinfo->nr-1; i++)
{
for (j=1; jinfo->nc-1; j++)
{
mag->data[i][j] = 0;
/* Treat the x and y derivatives as components of a vector */
xc = dx[i][j];
yc = dy[i][j];
if (fabs(xc)<0.01 && fabs(yc)<0.01) continue;
g = norm (xc, yc);
/* Follow the gradient direction, as indicated by the direction of
the vector (xc, yc); retain pixels that are a local maximum. */
if (fabs(yc) > fabs(xc))
{
/* The Y component is biggest, so gradient direction is basically UP/DOWN */
xx = fabs(xc)/fabs(yc);
yy = 1.0;
g2 = norm (dx[i-1][j], dy[i-1][j]);
g4 = norm (dx[i+1][j], dy[i+1][j]);
if (xc*yc > 0.0)
{
g3 = norm (dx[i+1][j+1], dy[i+1][j+1]);
g1 = norm (dx[i-1][j-1], dy[i-1][j-1]);
} else
{
g3 = norm (dx[i+1][j-1], dy[i+1][j-1]);
g1 = norm (dx[i-1][j+1], dy[i-1][j+1]);
}
} else
{
/* The X component is biggest, so gradient direction is basically LEFT/RIGHT */
xx = fabs(yc)/fabs(xc);
yy = 1.0;
g2 = norm (dx[i][j+1], dy[i][j+1]);
g4 = norm (dx[i][j-1], dy[i][j-1]);
if (xc*yc > 0.0)
{
g3 = norm (dx[i-1][j-1], dy[i-1][j-1]);
g1 = norm (dx[i+1][j+1], dy[i+1][j+1]);
}
else
{
g1 = norm (dx[i-1][j+1], dy[i-1][j+1]);
g3 = norm (dx[i+1][j-1], dy[i+1][j-1]);
}
}
/* Compute the interpolated value of the gradient magnitude */
if ( (g > (xx*g1 + (yy-xx)*g2)) &&
(g > (xx*g3 + (yy-xx)*g4)) )
{
if (g*MAG_SCALE <= 255)
mag->data[i][j] = (unsigned char)(g*MAG_SCALE);
else
mag->data[i][j] = 255;
ori->data[i][j] = atan2 (yc, xc) * ORI_SCALE;
} else
{
mag->data[i][j] = 0;
ori->data[i][j] = 0;
}
}
}
}
void Canny::estimate_thresh (imagec* mag, int *hi, int *low)
{
int i,j,k, hist[256], count;
/* Build a histogram of the magnitude image. */
for (k=0; k<256; k++) hist[k] = 0;
for (i=WIDTH; iinfo->nr-WIDTH; i++)
for (j=WIDTH; jinfo->nc-WIDTH; j++)
hist[mag->data[i][j]]++;
/* The high threshold should be > 80 or 90% of the pixels
j = (int)(ratio*mag->info->nr*mag->info->nc);
*/
j = mag->info->nr;
if (jinfo->nc) j = mag->info->nc;
j = (int)(0.9*j);
k = 255;
count = hist[255];
while (count < j)
{
k--;
if (k<0) break;
count += hist[k];
}
*hi = k;
i=0;
while (hist[i]==0) i++;
*low = (*hi+i)/2.0;
}