Canny edge detector

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={\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$ - current pixel, $p_{a},\;p_{b}$ - 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>
include <stdio.h>
include <stdlib.h>
include <float.h>
include <math.h>
include <string.h>

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

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 int i; unsigned char c; 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 };

FILE* fp = fopen(filename,"w"); if (fp == NULL) return 1; struct bmpfile_magic mag = Template:0x42, 0x4d;

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 rgb_t color = {0, 0, 0, 0}; int i; 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); unsigned char c; for(i=0; i<bmp_ih->bmp_bytesz; i++){ 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 dir; 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);

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]; }

int nn,ss,ww,ee,nw,ne,sw,se; // Non-maximum suppression, straightforward implementation. for(i=1; i<nx-1; i++) for(j=1; j<ny-1; j++) { 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 = atan2(after_Gy[c], after_Gx[c]) * 8.0/M_PI;

if( (( fabs(dir)<=1 || fabs(fabs(dir)-8)<=1 ) && G[c] > G[ee] && G[c] > G[ww]) /* 0 deg */ || (( fabs(dir-2)<=1 || fabs(dir+6)<=1 ) && G[c] > G[nw] && G[c] > G[se]) /* 45 deg */ || (( fabs(fabs(fabs(dir)-4))<=1 ) && G[c] > G[nn] && G[c] > G[ss]) /* 90 deg */ || (( fabs(dir-6)<=1 || fabs(dir+2)<=1 ) && G[c] > G[ne] && G[c] > G[sw]) /* 135 deg */ ) nms[c] = G[c]; else nms[c] = 0; } // Reuse array pixel_t *edges = after_Gy; // used as a stack memset(out, 0, sizeof(pixel_t)*nx*ny); memset(edges, 0, sizeof(pixel_t)*nx*ny); int nedges, k, t; int nbs[8]; // neighbours

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) { if(argc < 2){ printf("Usage: %s image.bmp\n", argv[0]); exit(1); } pixel_t *bitmap_data, *temp_image;

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>