canny.cpp

A tutorial and open source code for finding edges and corners based on the filters used in primary v

//////////////////////////////////////////////////////////////////
// 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 <math.h>
#include <stdio.h>
#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<unsigned char> temp;
	im=newimagec (iHeight,iWidth);
	//Convert the RGB image to grayscale and load it into buffer
	for(i=0 ; i<iHeight ; i++)
		for(j=0 ; j<iWidth ; j++)
		{
            temp = m_pImage->Pixel(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 ; i<iHeight ; i++)
		for(j=0 ; j<iWidth ; j++)
		{
			if(im->data[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; i<nr; i++)
	{  
	  x[i] = (float *) calloc ( nc, sizeof (float)  );
	  if (x[i] == 0)
	  {
	    fprintf (stderr, "Out of storage: F2D %d.\n", i);
	    exit (1);
	  }
	}
	return x;
}


struct imagec  * Canny::newimagec (int nr, int nc)
{
	struct imagec  *x;                /* New image */
//	unsigned char *ptr;             /* new pixel array */
	int i;

	if (nr < 0 || nc < 0) {
		printf ("Error: Bad image size (%d,%d)\n", nr, nc);
		return 0;
	}

/*      Allocate the image structure    */
	x = (struct imagec  *) malloc( sizeof (struct imagec) );
	if (!x) {
		printf ("Out of storage in NEWIMAGE.\n");
		return 0;
	}

/*      Allocate and initialize the header      */

	x->info = (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; i<nr; i++) 
	{
	  x->data[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; i<z->info->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; i<WIDTH; i++)
	  for (j=0; j<im->info->nc; j++)
	    im->data[i][j] = 255;

	for (i=im->info->nr-1; i>im->info->nr-1-WIDTH; i--)
	  for (j=0; j<im->info->nc; j++)
	    im->data[i][j] = 255;

	for (i=0; i<im->info->nr; i++)
	  for (j=0; j<WIDTH; j++)
	    im->data[i][j] = 255;

	for (i=0; i<im->info->nr; i++)
	  for (j=im->info->nc-WIDTH-1; j<im->info->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; i<MAX_MASK_SIZE; i++)
	{
	  gau[i] = meanGauss ((float)i, s);
	  if (gau[i] < 0.005)
	  {
		width = i;
		break;
	  }
	  dgau[i] = dGauss ((float)i, s);
	}

	n = width+width + 1;
	WIDTH = width/2;
	printf ("Smoothing with a Gaussian (width = %d) ...\n", n);

	smx = f2d (im->info->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; i<im->info->nr; i++)
	  for (j=0; j<im->info->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; i<im->info->nr; i++)
	  for (j=0; j<im->info->nc; j++)
	    im->data[i][j] = 0;

	if (high<low)
	{
	  estimate_thresh (mag, &high, &low);
	  printf ("Hysteresis thresholds (from image): HI %d LOW %D\n",
			high, low);
	}
/* For each edge with a magnitude above the high threshold, begin
   tracing edge pixels that are above the low threshold.                */

	for (i=0; i<im->info->nr; i++)
	  for (j=0; j<im->info->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; i<im->info->nr; i++)
	  for (j=0; j<im->info->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; i<nr; i++)
	  for (j=0; j<nc; j++)
	  {
	    x = gau[0] * im->data[i][j]; y = gau[0] * im->data[i][j];
	    for (k=1; k<width; k++)
	    {
	      I1 = (i+k)%nr; I2 = (i-k+nr)%nr;
	      y += gau[k]*im->data[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; i<nr; i++)
	  for (j=0; j<nc; j++)
	  {
	    x = 0.0;
	    for (k=1; k<width; k++)
	    {
	      if (which == 0)
	      {
		I1 = (i+k)%nr; I2 = (i-k+nr)%nr;
		x += -gau[k]*im[I1][j] + gau[k]*im[I2][j];
	      }
	      else
	      {
		I1 = (j+k)%nc; I2 = (j-k+nc)%nc;
		x += -gau[k]*im[i][I1] + gau[k]*im[i][I2];
	      }
	    }
	    sm[i][j] = x;
	  }
}

void Canny::nonmax_suppress (float **dx, float **dy, int nr, int nc, 
		imagec* mag, imagec* ori)
{
	int i,j;
	float xx, yy, g2, g1, g3, g4, g, xc, yc;

	for (i=1; i<mag->info->nr-1; i++)
	{
	  for (j=1; j<mag->info->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; i<mag->info->nr-WIDTH; i++)
	  for (j=WIDTH; j<mag->info->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 (j<mag->info->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;
}