Canny edge detector

Canny edge detector is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Task: Write a program that performs so-called canny edge detection on an image. The algorithm consists of the following steps:

1. Noise reduction. May be performed by Gaussian filter.
2. Compute intensity gradient (matrices $G_{x}$ and $G_{y}$) and its magnitude $G$.
      $G={\sqrt {G_{x}^{2}+G_{y}^{2}}}$
   May be performed by convolution of an image with Sobel operators.
3. Non-maximum suppression. For each pixel compute the orientation of intensity gradient vector: $\theta ={\rm {atan2}}\left(G_{y},\,G_{x}\right)$. Transform angle $\theta$ to one of four directions: 0, 45, 90, 135 degrees. Compute new array $N$: if
      $G\left(p_{a}\right)<G\left(p\right)<G\left(p_{b}\right)$
   where $p$ is the current pixel, $p_{a}$ and $p_{b}$ are the two neighbour pixels in the direction of gradient, then $N(p)=G(p)$, otherwise $N(p)=0$. Nonzero pixels in resulting array correspond to local maxima of $G$ in direction $\theta (p)$.
4. Tracing edges with hysteresis. At this stage two thresholds for the values of $G$ are introduced: $T_{min}$ and $T_{max}$. Starting from pixels with $N(p)\geqslant T_{max}$ find all paths of pixels with $N(p)\geqslant T_{min}$ and put them to the resulting image.

C

The following program reads an 8 bits per pixel grayscale BMP file and saves the result to `out.bmp'. Compile with `-lm'. <lang c>#include <stdint.h>

1. include <stdio.h>
2. include <stdlib.h>
3. include <float.h>
4. include <math.h>
5. include <string.h>

1. define MAX_BRIGHTNESS 255

/*

* Loading part taken from
* http://www.vbforums.com/showthread.php?t=261522
* BMP info:
* http://en.wikipedia.org/wiki/BMP_file_format
*/

/*

* Note: the magic number has been removed from the bmpfile_header structure
* since it causes alignment problems
*     struct bmpfile_magic should be written/read first
* followed by the
*     struct bmpfile_header
* [this avoids compiler-specific alignment pragmas etc.]
*/

struct bmpfile_magic {

   unsigned char magic[2];

};

struct bmpfile_header {

   uint32_t filesz;
   uint16_t creator1;
   uint16_t creator2;
   uint32_t bmp_offset;

};

typedef struct {

   uint32_t header_sz;
   int32_t width;
   int32_t height;
   uint16_t nplanes;
   uint16_t bitspp;
   uint32_t compress_type;
   uint32_t bmp_bytesz;
   int32_t hres;
   int32_t vres;
   uint32_t ncolors;
   uint32_t nimpcolors;

} BITMAPINFOHEADER;

typedef struct {

   uint8_t r;
   uint8_t g;
   uint8_t b;
   uint8_t null;

} rgb_t;

BITMAPINFOHEADER ih;

// Use int instead `unsigned char' so that we can store negative values. typedef int pixel_t;

pixel_t *load_bmp(char *filename, BITMAPINFOHEADER *bitmapInfoHeader) {

   FILE *filePtr; //our file pointer
   struct bmpfile_magic  mag;
   struct bmpfile_header bitmapFileHeader; //our bitmap file header
   pixel_t *bitmapImage;  //store image data
   int i;
   unsigned char c;

   filePtr = fopen(filename,"r");
   if (filePtr == NULL) {
       perror("fopen()");
       exit(1);
   }

   fread(&mag, sizeof(struct bmpfile_magic),1,filePtr);
   //verify that this is a bmp file by check bitmap id
   if ( *((uint16_t*) mag.magic) != 0x4D42 ) {
       fprintf(stderr, "Not a BMP file: magic=%c%c\n", mag.magic[0], mag.magic[1]);
       fclose(filePtr);
       return NULL;
   }

   //read the bitmap file header
   fread(&bitmapFileHeader, sizeof(struct bmpfile_header), 1, filePtr);

   //read the bitmap info header
   fread(bitmapInfoHeader, sizeof(BITMAPINFOHEADER), 1, filePtr);

   if( bitmapInfoHeader->compress_type != 0)
       fprintf(stderr, "Warning, compression is not supported.\n");

   //move file point to the beginning of bitmap data
   fseek(filePtr, bitmapFileHeader.bmp_offset, SEEK_SET);

   //allocate enough memory for the bitmap image data
   bitmapImage = (pixel_t *)malloc(bitmapInfoHeader->bmp_bytesz*(sizeof(pixel_t)));

   //verify memory allocation
   if (!bitmapImage) {
       free(bitmapImage);
       fclose(filePtr);
       return NULL;
   }

   //read in the bitmap image data
   for(i=0; i<bitmapInfoHeader->bmp_bytesz; i++){
       fread(&c, sizeof(unsigned char), 1, filePtr);
       bitmapImage[i] = (int) c;
   }

   // If we were using unsigned char as pixel_t, then:
   //fread(bitmapImage, 1, bitmapInfoHeader->bmp_bytesz, filePtr);

   //close file and return bitmap image data
   fclose(filePtr);
   return bitmapImage;

}

// Return: nonzero on error. int save_bmp(char *filename, BITMAPINFOHEADER *bmp_ih, pixel_t *data) {

   unsigned int offset =
       sizeof(struct bmpfile_magic)
       + sizeof(struct bmpfile_header)
       + sizeof(BITMAPINFOHEADER)
       + (1<<bmp_ih->bitspp)*4;
   struct bmpfile_header bmp_fh = {
       .filesz = offset + bmp_ih->bmp_bytesz,
       .creator1 = 0,
       .creator2 = 0,
       .bmp_offset = offset
   };
   struct bmpfile_magic mag = Template:0x42, 0x4d;
   rgb_t color = {0, 0, 0, 0};
   int i;
   FILE* fp = fopen(filename,"w");

   if (fp == NULL)
       return 1;

   fwrite(&mag, 1, sizeof(struct bmpfile_magic), fp);
   fwrite(&bmp_fh, 1, sizeof(struct bmpfile_header), fp);
   fwrite(bmp_ih, 1, sizeof(BITMAPINFOHEADER), fp);

   // Palette
   for(i=0; i < (1<<bmp_ih->bitspp); i++) {
       color.r = i;
       color.g = i;
       color.b = i;
       fwrite(&color, 1, sizeof(rgb_t), fp);
   }

   // We use int instead of uchar, so we can't write img in 1 call any more.
   //fwrite(data, 1, bmp_ih->bmp_bytesz, fp);
   for(i=0; i<bmp_ih->bmp_bytesz; i++) {
       unsigned char c = (unsigned char) data[i];

       fwrite(&c, sizeof(unsigned char), 1, fp);
   }

   fclose(fp);
   return 0;

}

// if norm==1, map pixels to range 0..MAX_BRIGHTNESS void convolution(pixel_t *in, pixel_t *out, float *kernel, int nx, int ny, int kn, int norm) {

   int i, j, m, n, c;
   int khalf = floor(kn/2.);
   float pixel, min=DBL_MAX, max=DBL_MIN;

   if(norm)
       for(m=khalf; m<nx-khalf; m++)
           for(n=khalf; n<ny-khalf; n++) {
               pixel = 0;
               c = 0;
               for(j=-khalf; j<=khalf; j++)
                   for(i=-khalf; i<=khalf; i++)
                       pixel += in[(n-j)*nx + m-i]*kernel[c++];
               if(pixel < min)
                   min = pixel;
               if(pixel > max)
                   max = pixel;
               }

   for(m=khalf; m<nx-khalf; m++)
       for(n=khalf; n<ny-khalf; n++) {
           pixel = 0;
           c = 0;
           for(j=-khalf; j<=khalf; j++)
               for(i=-khalf; i<=khalf; i++)
                   pixel += in[(n-j)*nx + m-i]*kernel[c++];

           if(norm)
               pixel = MAX_BRIGHTNESS*(pixel-min)/(max-min);
           out[n*nx + m] = (pixel_t) pixel;
       }

}

// http://www.songho.ca/dsp/cannyedge/cannyedge.html // 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;

void gaussian_filter(pixel_t *in, pixel_t *out, int nx, int ny, float sigma) {

   int i, j, c=0;
   const int n = 2*(int)(2*sigma) + 3;
   float mean = floor(n/2.);
   float kernel[n*n];

   fprintf(stderr, "gaussian_filter: kernel size %d, sigma=%g\n", n, sigma);
   for(i=0; i<n; i++)
       for(j=0; j<n; j++)
           kernel[c++]=exp(-0.5*( pow((i-mean)/sigma,2.0) +pow((j-mean)/sigma,2.0)))
               /(2*M_PI*sigma*sigma);
   convolution(in, out, kernel, nx, ny, n, 1);

}

/*

* Links:
* http://en.wikipedia.org/wiki/Canny_edge_detector
* http://www.tomgibara.com/computer-vision/CannyEdgeDetector.java
* http://fourier.eng.hmc.edu/e161/lectures/canny/node1.html
* http://www.songho.ca/dsp/cannyedge/cannyedge.html
*
* Note: T1 and T2 are lower and upper thresholds.
*/

void canny_edge_detection(pixel_t *in, pixel_t *out, int nx, int ny, int tmin, int tmax, float sigma) {

   int i, j, c, Gmax;
   int counter = 0;
   float Gx[]={
       -1, 0, 1,
       -2, 0, 2,
       -1, 0, 1};
   float Gy[]={
        1, 2, 1,
        0, 0, 0,
       -1,-2,-1};
   pixel_t *G    = calloc(nx*ny*sizeof(pixel_t), 1),
       *after_Gx = calloc(nx*ny*sizeof(pixel_t), 1),
       *after_Gy = calloc(nx*ny*sizeof(pixel_t), 1),
       *nms      = calloc(nx*ny*sizeof(pixel_t), 1);
   pixel_t *edges;
   int nedges, k, t;
   int nbs[8]; // neighbours

   gaussian_filter(in, out, nx, ny, sigma);

   convolution(out, after_Gx, Gx, nx, ny, 3, 0);
   convolution(out, after_Gy, Gy, nx, ny, 3, 0);

   for(i=1; i<nx-1; i++)
       for(j=1; j<ny-1; j++) {
           c = i + nx*j;
           //G[c] = abs(after_Gx[c]) + abs(after_Gy[c]);
           G[c] = hypot(after_Gx[c], after_Gy[c]);
           if (G[c]>Gmax)
               Gmax = G[c];
       }

   // Non-maximum suppression, straightforward implementation.
   for(i=1; i<nx-1; i++)
       for(j=1; j<ny-1; j++) {
           float dir;
           int nn,ss,ww,ee,nw,ne,sw,se;

           c = i + nx*j;
           nn = c - nx;
           ss = c + nx;
           ww = c + 1;
           ee = c - 1;
           nw = nn + 1;
           ne = nn - 1;
           sw = ss + 1;
           se = ss - 1;

           dir = fmod(atan2(after_Gy[c], after_Gx[c])+M_PI, M_PI)/M_PI*8;

           if( ((dir<=1 || dir>7) && G[c]>G[ee] && G[c]>G[ww])   // 0 deg
               ||
               ((dir>1 && dir<=3) && G[c]>G[nw] && G[c]>G[se])   // 45 deg
               ||
               ((dir>3 && dir<=5) && G[c]>G[nn] && G[c]>G[ss])   // 90 deg
               ||
               ((dir>5 && dir<=7) && G[c]>G[ne] && G[c]>G[sw])   // 135 deg
               )
               nms[c] = G[c];
           else
               nms[c] = 0;
       }

   // Reuse array
   edges = after_Gy;   // used as a stack
   memset(out, 0, sizeof(pixel_t)*nx*ny);
   memset(edges, 0, sizeof(pixel_t)*nx*ny);

   counter = 0;
   // Tracing edges with hysteresis . Non-recursive implementation.
   for(c=1, j=1; j<ny-1; j++)
       for(i=1; i<nx-1; i++) {
           if( nms[c] >= tmax && out[c] == 0 ) {       // trace edges
               out[c] = MAX_BRIGHTNESS;
               nedges = 1;
               edges[0] = c;

               do{
                   nedges--;
                   t = edges[ nedges ];

                   nbs[0] = t - nx;                    // nn
                   nbs[1] = t + nx;                    // ss
                   nbs[2] = t + 1;                     // ww
                   nbs[3] = t - 1;                     // ee
                   nbs[4] = nbs[0] + 1;                // nw
                   nbs[5] = nbs[0] - 1;                // ne
                   nbs[6] = nbs[1] + 1;                // sw
                   nbs[7] = nbs[1] - 1;                // se

                   for(k=0; k<8; k++)
                       if( nms[ nbs[k] ] >= tmin && out[ nbs[k] ] == 0 ) {
                           out[ nbs[k] ] = MAX_BRIGHTNESS;
                           edges[nedges++] = nbs[k];
                       }
               } while( nedges > 0 );
           }
           c++;
       }
   free(after_Gx);
   free(after_Gy);
   free(G);
   free(nms);

}

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

   pixel_t *bitmap_data, *temp_image;

   if(argc < 2){
       printf("Usage: %s image.bmp\n", argv[0]);
       exit(1);
   }

   bitmap_data = load_bmp(argv[1], &ih);
   temp_image = (pixel_t *) malloc(ih.bmp_bytesz*sizeof(pixel_t));

   printf("Info: %d x %d x %d\n", ih.width, ih.height, ih.bitspp);

   canny_edge_detection(bitmap_data, temp_image, ih.width, ih.height, 45, 50, 1);

   save_bmp("out.bmp", &ih, temp_image);
   free(bitmap_data);
   free(temp_image);
   return 0;

}</lang>