cannyedge.c

Canny Edge sample good luck

// cannyedge.cxx
// =============
// Canny Edge detection method
// Processing steps are;
// 1) smoothing image using gaussian filter
// 2) compute gradient deltaX, deltaY of smoothed image
// 3) compute magnitude of deltaX & deltaY
// 4) non-maximal suppression
// 5) hysteresis threshold
//
// Runtime Keyboard option:
// ESC: exit program
// 
// COMPILE: CC -O3 -o cannyedge cannyedge.cxx -lglut -lGLU -lGL -lX11 -lXmu -lm
// 
//  AUTHOR: Song Ho Ahn ([e-mail])
// CREATED: 2002.03.14
// UPDATED: 2002.03.27
///////////////////////////////////////////////////////////////////////////////

// header files
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <GL/glut.h>


// image data structure for openGL
struct imageData{
   GLint  x;         // resolution X
   GLint  y;         // resolution Y
   GLenum format;    // openGL data format (RGB or INDEX..)
   GLenum type;      // openGL data type (8bit, 16bit or 32bit..)
   GLvoid *buf;      // image pixel bits
};



// CALLBACK functions /////////////////////////////////////////////////////////
void displayCB(void);
void displaySubWin1CB(void);
void reshapeCB(int w, int h);
void keyboardHandlerCB(unsigned char key, int x, int y);

// fuctions ///////////////////////////////////////////////////////////////////
void  init(void);
int   initSharedMem(void);
int   loadImage(char *filename);
void  cannyEdge(imageData *img, float sigma, float tLow, float tHigh,
                unsigned char *result); 
void  convolve(unsigned char *in, int x, int y, float *kernel, int kernelSize,
               unsigned char *out);
void  gaussian(unsigned char *in, int x, int y, float sigma,
               unsigned char *result);
void  makeGaussianKernel(float sigma, float *kernel, int kernelSize);
void  gradient(unsigned char *in, int x, int y, short *deltaX, short *deltaY);
void  magnitude(short *deltaX, short *deltaY, int x, int y,
                unsigned short *magnitude);
//void  direction(short *deltaX, short *deltaY, int x, int y, unsigned char *orient);
void  nonMaxSupp(unsigned short *mag, short *deltaX, short *deltaY, int x, int y,
                 unsigned char *nms);
void  hysteresis(unsigned char *nms, int x, int y, float tLow, float tHigh,
                 unsigned char *out);
void  traceEdge(unsigned char *buf, int t, int cols, int lowThreshold,
                unsigned char *out);
void  histogram(unsigned char *buf, int x, int y, unsigned short *hosto,
                int *min, int *max);
void  histogram(unsigned short *buf, int x, int y, unsigned int *hosto,
                int *min, int *max);
float linearInterp(float p1, float p2, float t);
float bilinearInterp(float p1, float p2, float p3, float p4, float tx, float ty);
void  threshold(unsigned char *buf, int x, int y, float p);
void drawString(char *str, int x, int y, void *font);



// CONSTANTS //////////////////////////////////////////////////////////////////
#define SIZE_X 256                     // image width
#define SIZE_Y 256                     // image height


// global variables ///////////////////////////////////////////////////////////
imageData      *imgData;
unsigned char  *inBuf;                          // input image pixels
unsigned char  *outBuf;                         // output image pixels
char           *fileName;                       // image file name
int            mainWin, subWin1, subWin2;       // GLUT window ID
int            font = (int)GLUT_BITMAP_8_BY_13; // GLUT font type
int            fontWidth = 8;                   // GLUT font width
int            fontHeight = 13;                 // GLUT font height
float          sigma = 1.0;                     // standard deviation(positive)


///////////////////////////////////////////////////////////////////////////////
void main(int argc, char **argv)
{

   // use default image file if not specified
   if(argc == 2)
      fileName = argv[1];
   else{
      printf("Usage: %s <frm-file>\n", argv[0]); 
      fileName = "data/lena.frm";
      printf("\nUse default image \"%s\"\n", fileName);
   }

   if(initSharedMem() == -1) return;

   // read 256 grey-level raw image file (256x256) 
   if(loadImage(fileName) == -1)   // exit program if failed to load image
      return;

   printf("\n");
   printf("ESC: Exit program.\n");


   // GLUT stuff for windowing ////////////////////////////////////////////////
   // initialization openGL window.
   // it is called before any other GLUT routine
   glutInit(&argc, argv);

   //glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE); // display mode
   glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE); // display mode
   glutInitWindowSize(2*imgData->x, imgData->y); // window size
   glutInitWindowPosition(100, 100); // window location

   // finally, create a window with openGL context
   // Window will not displayed until glutMainLoop() is called
   // it returns a unique ID
   mainWin = glutCreateWindow(argv[0]); // param is the title of window

   // initialize OpenGL stuff
   init();

   // register GLUT callback functions for main window
   glutDisplayFunc(displayCB);
   glutReshapeFunc(reshapeCB);
   glutKeyboardFunc(keyboardHandlerCB);

   // sub-window //////////////////////////////////////////////////////////////
   // each sub-windows has its own openGL context, callbacks
   // subWin1: dispaly for original image
   subWin1 = glutCreateSubWindow(mainWin, SIZE_X, 0, SIZE_X, SIZE_Y);
   glutDisplayFunc(displaySubWin1CB);
   glutKeyboardFunc(keyboardHandlerCB);
   init();




   // Canny Edge detection
   clock_t start = clock();
   double elapsed = 0;
   cannyEdge(imgData, sigma, 0.5, 0.92, outBuf);
   //cannyEdge(imgData, 1.0, 0.5, 0.96, outBuf);
   elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;
   printf("Canny Edge Elapsed Time: %4.2f sec\n", elapsed);



   // the last GLUT call (LOOP) ///////////////////////////////////////////////
   // window will be shown and display callback is triggered by events
   // NOTE: this call never return main().
   glutMainLoop(); /* Start GLUT event-processing loop */

} // END OF MAIN //////////////////////////////////////////////////////////////



///////////////////////////////////////////////////////////////////////////////
// openGL initialization 
// disable unnecessary OpenGL operations
///////////////////////////////////////////////////////////////////////////////
void init(void)
{
   glClearColor(0, 0, 0, 0);
   glShadeModel(GL_FLAT);  // shading mathod: GL_SMOOTH or GL_FLAT
   //glPixelStorei(GL_UNPACK_ALIGNMENT, 1);

   // TUNNING IMAGING PERFORMANCE
   // turn off  all pixel path and per-fragment operation
   // which slow down OpenGL imaging operation (glDrawPixels glReadPixels...).
   glDisable(GL_ALPHA_TEST);
   glDisable(GL_BLEND);
   glDisable(GL_DEPTH_TEST);
   glDisable(GL_DITHER);
   glDisable(GL_FOG);
   glDisable(GL_LIGHTING);
   glDisable(GL_LOGIC_OP);
   glDisable(GL_STENCIL_TEST);
   glDisable(GL_TEXTURE_1D);
   glDisable(GL_TEXTURE_2D);

   glPixelTransferi(GL_MAP_COLOR, GL_FALSE);
   glPixelTransferi(GL_RED_SCALE, 1);
   glPixelTransferi(GL_RED_BIAS, 0);
   glPixelTransferi(GL_GREEN_SCALE, 1);
   glPixelTransferi(GL_GREEN_BIAS, 0);
   glPixelTransferi(GL_BLUE_SCALE, 1);
   glPixelTransferi(GL_BLUE_BIAS, 0);
   glPixelTransferi(GL_ALPHA_SCALE, 1);
   glPixelTransferi(GL_ALPHA_BIAS, 0);
   //glPixelZoom(1.0, 1.0);

   // disable extensions that could slow down glDrawPixels.
   const GLubyte* extString = glGetString(GL_EXTENSIONS);
   if(extString != NULL){
      if(strstr((char*) extString, "GL_EXT_convolution") != NULL){
         glDisable(GL_CONVOLUTION_1D_EXT);
         glDisable(GL_CONVOLUTION_2D_EXT);
         glDisable(GL_SEPARABLE_2D_EXT);
      }
      if(strstr((char*) extString, "GL_EXT_histogram") != NULL){
         glDisable(GL_HISTOGRAM_EXT);
         glDisable(GL_MINMAX_EXT);
      }
      if(strstr((char*) extString, "GL_EXT_texture3D") != NULL){
         glDisable(GL_TEXTURE_3D_EXT);
      }
   }
}


///////////////////////////////////////////////////////////////////////////////
// initialize global variables
///////////////////////////////////////////////////////////////////////////////
int initSharedMem(void)
{
   int returnVal = 0;

   imgData = (imageData *)malloc(sizeof(imageData));
   if(imgData == NULL)
   {
      printf("ERROR: Memory Allocation Failed..\n");
      returnVal = -1;
   }
   return returnVal;
}


///////////////////////////////////////////////////////////////////////////////
// read input image data and store it to system memory
///////////////////////////////////////////////////////////////////////////////
int loadImage(char *fileName)
{
   FILE *fp;   // input file pointer

   // open the image file
   if((fp = fopen(fileName, "r")) == NULL){
      printf("Cannot open %s\n", fileName);
      return -1;
   }

   // get memory space based on image data size
   size_t bufsize = SIZE_X * SIZE_Y;
   inBuf = (unsigned char *)malloc(bufsize);
   outBuf = (unsigned char *)malloc(bufsize);

   if(inBuf == NULL){
      printf("ERROR: Memory Allocation Failed..\n");
      return -1;
   }
   if(outBuf == NULL){
      printf("ERROR: Memory Allocation Failed..\n");
      return -1;
   }

   // read image pixels
   fread(inBuf, 1, bufsize, fp);
   fclose(fp);

   // copy inBuf to outBuf
   memcpy(outBuf, inBuf, bufsize);

   // get openGL image data structure
   imgData->x = SIZE_X;
   imgData->y = SIZE_Y;
   imgData->format = GL_LUMINANCE;
   imgData->type = GL_UNSIGNED_BYTE;
   imgData->buf = (GLvoid*)inBuf;

   return 0;
}


///////////////////////////////////////////////////////////////////////////////
// Canny Edge Detection (Thinning edges)
// 1) smooth image using Gaussian filter
// 2) compute gradient X & Y
// 3) compute magnitude of gradient X & Y
// 4) non-maximal suppression (thinning ridge)
// 5) hysteresis thrsholding
///////////////////////////////////////////////////////////////////////////////
void cannyEdge(imageData *img, float sigma, float tLow, float tHigh, unsigned char *out)
{
   int x, y;
   short *deltaX, *deltaY;
   unsigned short *mag;
   unsigned char *imgBuf, *nms;

   imgBuf = (unsigned char*)img->buf;  // input image buffer
   x = img->x;                         // resolution X
   y = img->y;                         // rsolution Y

   // allocate memory space
   if((deltaX = (short*)malloc(x*y*sizeof(short))) == NULL){
      printf("ERROR: Memory Allocation Failed.\n");
      exit(0);
   }
   if((deltaY = (short*)malloc(x*y*sizeof(short))) == NULL){
      printf("ERROR: Memory Allocation Failed.\n");
      exit(0);
   }
   if((mag = (unsigned short*)malloc(x*y*sizeof(unsigned short))) == NULL){
      printf("ERROR: Memory Allocation Failed.\n");
      exit(0);
   }
   if((nms = (unsigned char*)malloc(x*y)) == NULL){
      printf("ERROR: Memory Allocation Failed.\n");
      exit(0);
   }


   // 1) smoothing
   if(sigma != 0){   // skip smoothing when sigma = 0.
      gaussian(imgBuf, x, y, sigma, out);
      imgBuf = out;  // assign gaussian result as source of later processes
   }

   // 2) compute derivative X & Y
   gradient(imgBuf, x, y, deltaX, deltaY);
         
   // 3) compute magnitude of dx & dy
   magnitude(deltaX, deltaY, x, y, mag);

   // 4) non-maximal suppression
   nonMaxSupp(mag, deltaX, deltaY, x, y, nms);

   // 5) hysteresis thresholding
   // RECOMMENDATION: 0.8 <= tHigh <= 0.9,  0.3 <= tLow <= 0.5
   hysteresis(nms, x, y, tLow, tHigh, out);


   // free up memory
   free(deltaX);
   free(deltaY);
   free(mag);
   free(nms);
}



///////////////////////////////////////////////////////////////////////////////
// 1D convolution of image data and kernel(filter) matrix
///////////////////////////////////////////////////////////////////////////////
void convolve(unsigned char *buf, int  x, int y, float *kernel, int kernelSize,
              unsigned char *out)
{
   int i, j, k, t;
   int center;
   unsigned char result = 0;

   center = kernelSize / 2;

   // smoothing horizontal direction (X)
   t = 0;
   for(i= 0; i < y; i++)            // rows
      for(j = 0; j < x; j++, t++)   // columns
      {
         for(k = -center; k <= center; k++)
            if(((j+k) >= 0) && ((j+k) < x))
               result += buf[t+k] * kernel[k+center] + 0.5;
            else  // if out of bound, assume same as center
               result += buf[t] * kernel[k+center] + 0.5;

         out[t] = result;