cannyedge.c

Canny Edge sample good luck

         result = 0;          // reset
      }

   // smoothing vertical direction (Y)
   t = 0;
   for(j= 0; j < x; j++, t=j)       // columns
      for(i = 0; i < y; i++, t+=x)  // rows
      {
         for(k = -center; k <= center; k++)
            if(((i+k) >= 0) && ((i+k) < y))
               result += out[t+(k*x)] * kernel[k+center] + 0.5;
            else
               result += out[t] * kernel[k+center] + 0.5;
        
         out[t] = result;
         result = 0;          // reset
      }
}




///////////////////////////////////////////////////////////////////////////////
// smooth image with Gaussian filter
///////////////////////////////////////////////////////////////////////////////
void gaussian(unsigned char *inBuf, int x, int y, float sigma,
              unsigned char *outBuf)
{
   int kernelSize;
   float *kernel;

   // determine size of kernel (odd #)
   // 0.0 <= sigma < 0.5 : 3
   // 0.5 <= sigma < 1.0 : 5
   // 1.0 <= sigma < 1.5 : 7
   // 1.5 <= sigma < 2.0 : 9
   // 2.0 <= sigma < 2.5 : 11
   // 2.5 <= sigma < 3.0 : 13 ...
   kernelSize = 2 * int(2*sigma) + 3;
  
   // allocate memory space for kernel
   kernel = (float *)malloc(kernelSize*sizeof(float));
   if(kernel == NULL){
      printf("ERROR: Memory Allocation Failed..\n");
      exit(-1);
   }

   // generate 1-D gaussian kernel
   makeGaussianKernel(sigma, kernel, kernelSize);

   // smoothing image with Gaussian filter
   // it is convolution of image array and kernel
   // Use 1D convolution because Gaussian Filter is seperable
   convolve(inBuf, x, y, kernel, kernelSize, outBuf);

   // free up memory space for kernel
   free(kernel);

#if 1
   printf("sigma: %3.2f,  Kernel Size: %d\n", sigma, kernelSize);
#endif
}

   
///////////////////////////////////////////////////////////////////////////////
// generate 1D Gaussian kernel
// kernel size should be odd number (3, 5, 7, 9, ...)
///////////////////////////////////////////////////////////////////////////////
void makeGaussianKernel(float sigma, float *kernel, int kernelSize)
{
   //const double PI = 3.14159265;       // PI
   
   int i, center;
   float sum = 0;                      // used for normalization
   double result = 0;                  // result of gaussian func

   // compute kernel elements normal distribution equation(Gaussian)
   // do only half(positive area) and mirror to negative side
   // because Gaussian is even function, symmetric to Y-axis.
   center = kernelSize / 2;   // center value of n-array(0 ~ n-1)

   if(sigma == 0){
      for(i = 0; i <= center; i++)
         kernel[center+i] = kernel[center-i] = 0;
      kernel[center] = 1.0;
   }

   else{
      for(i = 0; i <= center; i++){
         //result = exp(-(i*i)/(double)(2*sigma*sigma)) / (sqrt(2*PI)*sigma);
         // dividing (sqrt(2*PI)*sigma) is not needed because normalizing result later
         result = exp(-(i*i)/(double)(2*sigma*sigma));
         kernel[center+i] = kernel[center-i] = (float)result;
         sum += (float)result;
         if(i != 0) sum += (float)result;
      }

      // normalize kernel
      // make sum of all elements in kernel to 1
      for(i = 0; i <= center; i++)
         kernel[center+i] = kernel[center-i] /= sum;
   }

#if 0 // print result for debug
   printf("1-D Gaussian Kernel\n");
   printf("===================\n");
   sum = 0;
   for(i = 0; i < kernelSize; i++){
      printf("%3d: %10.9f\n", i, kernel[i]);
      sum += kernel[i];
   }
   printf("KERNEL SUM: %f\n", sum);
#endif
}
      

///////////////////////////////////////////////////////////////////////////////
// compute gradient (first order derivative x and y)
///////////////////////////////////////////////////////////////////////////////
void gradient(unsigned char *pix, int x, int y, short *deltaX, short *deltaY)
{
   int i, j;
   int t;

   // compute delta X ***************************
   // deltaX = f(x+1) - f(x-1)
   for(i = 0; i < y; i++)
   {
      t = i * x;  // current position X per line 

      // gradient at the first pixel of each line
      // note: the edge,pix[t-1] is NOT exsit
      deltaX[t] = pix[t+1] - pix[t];

      // gradients where NOT edge
      for(j = 1; j < x-1; j++){
         t++;
         deltaX[t] = pix[t+1] - pix[t-1];
      //printf("%d", deltaX[t]);
      }

      // gradient at the last pixel of each line
      t++;
      deltaX[t] = pix[t] - pix[t-1];

      //printf("%d", deltaX[t]);
   }

   // compute delta Y ***************************
   // deltaY = f(y+1) - f(y-1)
   for(j = 0; j < x; j++)
   {
      t = j;   // current Y position per column

      // gradient at the first pixel
      deltaY[t] = pix[t+x] - pix[t];

      // gradients for NOT edge pixels
      for(i = 1; i < y-1; i++){
         t += x;
         deltaY[t] = pix[t+x] - pix[t-x];
      }

      // gradient at the last pixel of each column
      t += x;
      deltaY[t] = pix[t] - pix[t-x];
   }

}


///////////////////////////////////////////////////////////////////////////////
// compute magnitude of gradient(deltaX & deltaY) per pixel
///////////////////////////////////////////////////////////////////////////////
void magnitude(short *deltaX, short *deltaY, int x, int y, unsigned short *mag)
{
   int i, j, t;

#if 1 // accurate computation (SLOW)
   t = 0;
   for(i = 0; i < y; i++)
      for(j = 0; j < x; j++, t++)
      {
         mag[t] = (unsigned short)(sqrt((double)deltaX[t]*deltaX[t] + 
                               (double)deltaY[t]*deltaY[t]) + 0.5) ;
      }
#endif

#if 0 // faster routine (use absolute value)
   t = 0;
   for(i = 0; i < y; i++)
      for(j = 0; j < x; j++, t++)
         mag[t] = (unsigned short)(abs(deltaX[t]) + abs(deltaY[t]));
#endif
}


///////////////////////////////////////////////////////////////////////////////
// compute direction of edges for each pixel (angle of 1st derivative of image)
// quantize the normal directions into 16 plus 0(dx=dy=0)
///////////////////////////////////////////////////////////////////////////////
void direction(short *deltaX, short *deltaY, int x, int y, unsigned char *orient)
{
   int i, j, t;

#if 1
   // compute direction into 16 ways plus 0(dx=dy=0)
   // it is symmetry on origin where theta > 180 
   // 0 : no direction, when dx = dy = 0
   // 1 :       theta = 0,   when  dy = 0
   // 2 :   0 < theta < 45,  when  dx > dy > 0
   // 3 :       theta = 45,  when  dx = dy > 0
   // 4 :  45 < theta < 90,  when  dy > dx > 0 
   // 5 :       theta = 90,  when  dx = 0, dy > 0
   // 6 :  90 < theta < 135, when  dy > -dx > 0
   // 7 :       theta = 135, when  dy = -dx > 0
   // 8 : 135 < theta < 180, when -dx > dy > 0
   // 9 :       theta = 180, when  dy = 0, dx > 0
   //10 : 180 < theta < 225, when  -dx > -dy > 0
   //11 :       theta = 225, when  dx = dy < 0
   //12 : 225 < theta < 270, when  dy < dx < 0
   //13 :       theta = 270, when  dx = 0, dy < 0
   //14 : 270 < theta < 315, when -dy > dx > 0
   //15 :       theta = 315, when  dx = -dy > 0
   //16 : 315 < theta < 360, when  dx > -dy > 0
   t = 0;
   for(j = 0; j < y; j++)
      for(i = 0; i < x; i++)
      {
         if(deltaX[t] == 0) // all axis directions
            if(deltaY[t] == 0) orient[t] = 0;
            else if(deltaY[t] > 0) orient[t] = 5;
            else orient[t] = 13;

         else if(deltaX[t] > 0)
            if(deltaY[t] == 0) orient[t] = 1;
            else if(deltaY[t] > 0)
               if(deltaX[t] - deltaY[t] == 0) orient[t] = 3;
               else if(deltaX[t] - deltaY[t] > 0) orient[t] = 2;
               else orient[t] = 4;
            else
               if(deltaX[t] + deltaY[t] == 0) orient[t] = 15;
               else if(deltaX[t] + deltaY[t] > 0) orient[t] = 16;
               else orient[t] = 14;

         else
            if(deltaY[t] == 0) orient[t] = 9;
            else if(deltaY[t] > 0)
               if(deltaY[t] + deltaX[t] == 0) orient[t] = 7;
               else if(deltaY[t] + deltaX[t] > 0) orient[t] = 6;
               else orient[t] = 8;
            else
               if(deltaY[t] - deltaX[t] == 0) orient[t] = 11;
               else if(deltaY[t] - deltaX[t] > 0) orient[t] = 10;
               else orient[t] = 12;

         //if(orient[t] == 16) printf("%d,%d, %d    ", deltaX[t], deltaY[t],orient[t]);
         t++;
      }
#endif


#if (0)
   t = 0;
   for(j = 0; j < y; j++)
      for(i = 0; i < x; i++)
      {
         if(deltaX[t] == 0 && deltaY[t] == 0)
            orient[t] = 0;
         else if(deltaX[t] == 0
         else if(deltaX[t] >= 0 && deltaY[t] >= 0) // quadrant I
            deltaX[t] - deltaY[t] >= 0 ? orient[t] = 1 : orient[t] = 2;
         else if(deltaX[t] < 0 && deltaY[t] >= 0) // quadrant II
            deltaY[t] + deltaX[t] >= 0 ? orient[t] = 3 : orient[t] = 4;
         else if(deltaX[t] < 0 && deltaY[t] < 0) // quadrant III
            deltaY[t] - deltaX[t] >= 0 ? orient[t] = 1 : orient[t] = 2;
         else if(deltaX[t] >= 0 && deltaY[t] < 0) // quadrant IV
            deltaX[t] + deltaY[t] >= 0 ? orient[t] = 4 : orient[t] = 3;

         printf("%d, %d %d    ", deltaX[t], deltaY[t],orient[t]);
         t++; // next pixel
      }
#endif
}


///////////////////////////////////////////////////////////////////////////////
// Non Maximal Suppression
// If the centre pixel is not greater than neighboured pixels in the direction,
// then the center pixel is set to zero.
// This process results in one pixel wide ridges.
///////////////////////////////////////////////////////////////////////////////
void nonMaxSupp(unsigned short *mag, short *deltaX, short *deltaY, int x, int y,
                unsigned char *nms)
{
   int i, j, t;
   float alpha;
   float mag1, mag2;
   const unsigned char SUPPRESSED = 0;

   // put zero all boundaries of image
   // TOP edge line of the image
   for(j = 0; j < x; j++)
      nms[j] = 0;

   // BOTTOM edge line of image
   t = (y-1)*x;
   for(j = 0; j < x; j++, t++)
      nms[t] = 0;

   // LEFT & RIGHT edge line
   t = x;
   for(i = 1; i < y; i++, t+=x){
      nms[t] = 0;
      nms[t+x-1] = 0;
   }


   t = x + 1;  // skip boundaries of image
   // start and stop 1 pixel inner pixels from boundaries
   for(i = 1; i < y-1; i++, t+=2){
      for(j = 1; j < x-1; j++, t++){
         // if magnitude = 0, no edge
         if(mag[t] == 0) nms[t] = SUPPRESSED;
         else{
            // looking for 8 directions (clockwise)
            // NOTE: positive dy means down(south) direction on screen(image)
            // because dy is defined as f(x,y+1)-f(x,y-1).
            // 1: dy/dx <= 1,  dx > 0, dy > 0
            // 2: dy/dx > 1,   dx > 0, dy > 0
            // 3: dy/dx >= -1, dx < 0, dy > 0
            // 4: dy/dx < -1,  dx < 0, dy > 0
            // 5: dy/dx <= 1,  dx < 0, dy < 0
            // 6: dy/dx > 1,   dx < 0, dy < 0
            // 7: dy/dx <= -1, dx > 0, dy < 0
            // 8: dy/dx > -1,  dx > 0, dy < 0
            // After determine direction,
            // compute magnitudes of nieghbour pixels using linear interpolation
            if(deltaX[t] >= 0){
               if(deltaY[t] >= 0){  // dx >= 0, dy >= 0
                  if((deltaX - deltaY[t]) >= 0){      // direction 1 (SEE, South-East-East)
                     alpha = (float)deltaY[t] / deltaX[t];
                     mag1 = (1-alpha)*mag[t+1] + alpha*mag[t+x+1];
                     mag2 = (1-alpha)*mag[t-1] + alpha*mag[t-x-1];

                  }else{                              // direction 2 (SSE)
                     alpha = (float)deltaX[t] / deltaY[t];
                     mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x+1];
                     mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x-1];
                  }
      
               }else{ // dx >= 0, dy < 0
                  if((deltaX[t] + deltaY[t]) >= 0){   // direction 8 (NEE)
                     alpha = (float)-deltaY[t] / deltaX[t];
                     mag1 = (1-alpha)*mag[t+1] + alpha*mag[t-x+1];
                     mag2 = (1-alpha)*mag[t-1] + alpha*mag[t+x-1];

                  }else{                              // direction 7 (NNE)
                     alpha = (float)deltaX[t] / -deltaY[t];
                     mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x-1];
                     mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x+1];
                  }
               }

            }else{
               if(deltaY[t] >= 0){ // dx < 0, dy >= 0
                  if((deltaX[t] + deltaY[t]) >= 0){   // direction 3 (SSW)
                     alpha = (float)-deltaX[t] / deltaY[t];
                     mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x-1];
                     mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x+1];

                  }else{                              // direction 4 (SWW)
                     alpha = (float)deltaY[t] / -deltaX[t];
                     mag1 = (1-alpha)*mag[t-1] + alpha*mag[t+x-1];
                     mag2 = (1-alpha)*mag[t+1] + alpha*mag[t-x+1];
                  }

               }else{ // dx < 0, dy < 0