Trouble with peak finding in canny algorithm

I'm writing a canny edge detector for my intro to robotic vision class. There are output samples to check if we have it right, and my magnitude image is coming out exact but my peak finder is not complete. The peaks I have found appear to be in the correct output but I am missing a significant amount. I have looked through the peaks part of the code many times and still cannot find where the problem is. I was hoping that maybe someone with more experience with the canny algorithm could help me out. I've looked around the web for examples of canny algorithms but can't find the difference that is causing the problem. Any help or info would be appreciated. Thank you.

Code:

#include <stdio.h> #include <math.h> #define PICSIZE 256 #define MAXMASK 100 int pic[PICSIZE][PICSIZE]; double outpicx[PICSIZE][PICSIZE]; double outpicy[PICSIZE][PICSIZE]; double ival[PICSIZE][PICSIZE]; double maskx[MAXMASK][MAXMASK]; double masky[MAXMASK][MAXMASK]; double peaks[PICSIZE][PICSIZE]; main() { int i,j,p,q,mr,centx,centy, sig; double maskval,sum1, sum2, maxival, dir; FILE *fo1, *fo2,*fp1, *fo3; char str [80]; fp1=fopen("garb34.pgm","rb"); fscanf(fp1, "%s", str); fscanf(fp1, "%s", str); fscanf(fp1, "%s", str); //the 3 output files fo1=fopen("garb34mag.pgm","wb"); fo2=fopen("garb34peaks.pgm","wb"); fo3=fopen("garb34dblthrsh.pgm","wb"); printf("Enter the value for sigma.\n"); scanf("%d", &sig); mr = (int)(sig * 3); centx = (MAXMASK / 2); centy = (MAXMASK / 2); for (i=0;i<256;i++) { for (j=0;j<256;j++) { pic[i][j] = getc (fp1); } } for (p=-mr;p<=mr;p++) { for (q=-mr;q<=mr;q++) { maskval = ((-q)*(exp(-1*(((p*p)+(q*q))/(2*(sig*sig)))))); (maskx[p+centy][q+centx]) = maskval; maskval = ((-p)*(exp(-1*(((p*p)+(q*q))/(2*(sig*sig)))))); (masky[p+centy][q+centx]) = maskval; } } for (i=mr;i<=255-mr;i++) { for (j=mr;j<=255-mr;j++) { sum1 = 0; sum2 = 0; for (p=-mr;p<=mr;p++) { for (q=-mr;q<=mr;q++) { sum1 += pic[i+p][j+q] * maskx[p+centy][q+centx]; sum2 += pic[i+p][j+q] * masky[p+centy][q+centx]; } } outpicx[i][j] = sum1; outpicy[i][j] = sum2; } } maxival = 0; for (i=mr;i<256-mr;i++) { for (j=mr;j<256-mr;j++) { ival[i][j]=sqrt((double)((outpicx[i][j]*outpicx[i][j]) + (outpicy[i][j]*outpicy[i][j]))); if (ival[i][j] > maxival) maxival = ival[i][j]; } } //print the headers for all of the files. fprintf(fo1,"P5\n"); fprintf(fo1,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo1,"255\n"); fprintf(fo2,"P5\n"); fprintf(fo2,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo2,"255\n"); fprintf(fo3,"P5\n"); fprintf(fo3,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo3,"255\n"); //loop through each pixel and output its according value for each //of the 3 images. for (i=0;i<256;i++) { for (j=0;j<256;j++) { //gradient magnitude. ival[i][j] = (ival[i][j] / maxival) * 255; fprintf(fo1,"%c",(char)((int)(ival[i][j]))); if (outpicx[i][j] == 0) { //printf("here\n"); outpicx[i][j] = 10^-15; } dir = (atan2(outpicy[i][j],outpicx[i][j]))/3.14159 * 180.0; // Calculate actual direction of edge peaks[i][j] = 0; /* Convert actual edge direction to approximate value */ if (((dir < 22.5) && (dir > -22.5) ) || (dir > 157.5) || (dir < -157.5) ) { if((ival[i][j] > ival[i][j-1]) && (ival[i][j] > ival[i][j+1])) peaks[i][j] = 255; } else if(((dir > 22.5) && (dir < 67.5)) || ((dir < -112.5) && (dir > -157.5))) { if((ival[i][j] > ival[i+1][j+1]) && (ival[i][j] > ival[i-1][j-1])) peaks[i][j] = 255; } else if (((dir > 112.5) && (dir < 157.5)) || ((dir < -22.5) && (dir > -67.5))) { if((ival[i][j] > ival[i-1][j+1]) && (ival[i][j] > ival[i+1][j-1])) peaks[i][j] = 255; } else { if((ival[i][j] > ival[i+1][j]) && (ival[i][j] > ival[i-1][j])) peaks[i][j] = 255; } //low threshold of 40 /*if(ival[i][j] > 40) fprintf(fo2,"%c",(char)((int)(255))); else fprintf(fo2,"%c",(char)((int)(0))); */ //high threshold of 100 if(ival[i][j] > 100) fprintf(fo3,"%c",(char)((int)(255))); else fprintf(fo3,"%c",(char)((int)(0))); } } for (i=0;i<PICSIZE;i++) { for (j=0;j<PICSIZE;j++) { if (peaks[i][j] == 255) fprintf(fo2,"%c",(char)((int)(255))); else fprintf(fo2,"%c",(char)((int)(0))); } } //close all of the files. fclose(fp1); fclose(fo1); fclose(fo2); fclose(fo3); system("PAUSE"); }

Can't attach the pgm files but

here the input file can be found:
http://www.cs.ucf.edu/courses/cap4453/inputpics/ under garb34.pgm

here the correct output can be found:
http://www.cs.ucf.edu/courses/cap4453/outputpics/ under cannypeaks.pgm

Run my code to see what i am getting which is a partial output of the correct one it seems. (i.e. what outputs is correct its just missing parts)

What works better than random chunks of code is to actually describe the process in words, and then to turn those words into pseudocode, include BOTH of those with your code, so people understand what it is you want to have done.

Quzah.

Sorry couldn't find the edit button so i guess ill just repost. I've added some comments to clear it up as best i can but a lot of this code is just pieced together from the professors code of sobel and marr algorithms. Basically i want to take an image and output 3 images from it. The first is to be a magnitude image which my code does correctly. The 2nd is an image of all peaks (with no thresholding of any kind at this point) which is where I am having the problem. My code is finding peaks and those peaks are correct but it is not finding all peaks within the image.

A peak is where the two neighboring pixels are of lesser value than it. The direction in which to check is based on the inverse tangent of the y convolution (outpicy) divided by the x convolution (outpicx) It is supposed to be 100% just canny algorithm not doing anything fancy with it or modified in any way. So i guess the problem would be that im not correctly searching for all of the peaks. I hope this clears it up.

Code:

#include <stdio.h> #include <math.h> #define PICSIZE 256 #define MAXMASK 100 int pic[PICSIZE][PICSIZE]; double outpicx[PICSIZE][PICSIZE]; double outpicy[PICSIZE][PICSIZE]; double ival[PICSIZE][PICSIZE]; double maskx[MAXMASK][MAXMASK]; double masky[MAXMASK][MAXMASK]; double peaks[PICSIZE][PICSIZE]; main() { int i,j,p,q,mr,centx,centy, sig; double maskval,sum1, sum2, maxival, dir; FILE *fo1, *fo2,*fp1, *fo3; char str [80]; fp1=fopen("garb34.pgm","rb"); fscanf(fp1, "%s", str); fscanf(fp1, "%s", str); fscanf(fp1, "%s", str); //the 3 output files fo1=fopen("garb34mag.pgm","wb"); fo2=fopen("garb34peaks.pgm","wb"); fo3=fopen("garb34dblthrsh.pgm","wb"); printf("Enter the value for sigma.\n"); scanf("%d", &sig); //create the mask raidus based on the sigma value entered. mr = (int)(sig * 3); //create offsets based on the defined maximum mask. centx = (MAXMASK / 2); centy = (MAXMASK / 2); //get each pixel of the picture and place it into a 256x256 array. for (i=0;i<256;i++) { for (j=0;j<256;j++) { pic[i][j] = getc (fp1); } } //loop through from -mask radius to mask radius and create the x and y //masks using a gaussian curve. for (p=-mr;p<=mr;p++) { for (q=-mr;q<=mr;q++) { maskval = ((-q)*(exp(-1*(((p*p)+(q*q))/(2*(sig*sig)))))); (maskx[p+centy][q+centx]) = maskval; maskval = ((-p)*(exp(-1*(((p*p)+(q*q))/(2*(sig*sig)))))); (masky[p+centy][q+centx]) = maskval; } } //perform the convolution and store the output into outpicx and outpicy. for (i=mr;i<=255-mr;i++) { for (j=mr;j<=255-mr;j++) { sum1 = 0; sum2 = 0; for (p=-mr;p<=mr;p++) { for (q=-mr;q<=mr;q++) { sum1 += pic[i+p][j+q] * maskx[p+centy][q+centx]; sum2 += pic[i+p][j+q] * masky[p+centy][q+centx]; } } outpicx[i][j] = sum1; outpicy[i][j] = sum2; } } //loop through again and get the magnitude of each pixel. maxival = 0; for (i=mr;i<256-mr;i++) { for (j=mr;j<256-mr;j++) { ival[i][j]=sqrt((double)((outpicx[i][j]*outpicx[i][j]) + (outpicy[i][j]*outpicy[i][j]))); if (ival[i][j] > maxival) maxival = ival[i][j]; } } //print the headers for all of the files. fprintf(fo1,"P5\n"); fprintf(fo1,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo1,"255\n"); fprintf(fo2,"P5\n"); fprintf(fo2,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo2,"255\n"); fprintf(fo3,"P5\n"); fprintf(fo3,"%d %d\n", PICSIZE, PICSIZE); fprintf(fo3,"255\n"); //loop through each pixel and output its according value for each //of the 3 images. for (i=0;i<256;i++) { for (j=0;j<256;j++) { //gradient magnitude. ival[i][j] = (ival[i][j] / maxival) * 255; fprintf(fo1,"%c",(char)((int)(ival[i][j]))); if (outpicx[i][j] == 0) { //printf("here\n"); outpicx[i][j] = 10^-15; } /*The following code is where I would imagine the problem to be * * * * */ dir = (atan2(outpicy[i][j],outpicx[i][j]))/3.14159 * 180.0; // Calculate actual direction of edge //initialize the peaks of the image to zero. peaks[i][j] = 0; //check which sector the value falls under and see if its a peak. //if so turn the peak on (set to 255). if (((dir < 22.5) && (dir > -22.5) ) || (dir > 157.5) || (dir < -157.5) ) { //if its bigger than its two neighboring pixels if((ival[i][j] > ival[i][j-1]) && (ival[i][j] > ival[i][j+1])) peaks[i][j] = 255; } else if(((dir > 22.5) && (dir < 67.5)) || ((dir < -112.5) && (dir > -157.5))) { if((ival[i][j] > ival[i+1][j+1]) && (ival[i][j] > ival[i-1][j-1])) peaks[i][j] = 255; } else if (((dir > 112.5) && (dir < 157.5)) || ((dir < -22.5) && (dir > -67.5))) { if((ival[i][j] > ival[i-1][j+1]) && (ival[i][j] > ival[i+1][j-1])) peaks[i][j] = 255; } else { if((ival[i][j] > ival[i+1][j]) && (ival[i][j] > ival[i-1][j])) peaks[i][j] = 255; } //low threshold of 40 /*if(ival[i][j] > 40) fprintf(fo2,"%c",(char)((int)(255))); else fprintf(fo2,"%c",(char)((int)(0))); */ //high threshold of 100 //unimportant for my question just a placeholder for the 3rd output. if(ival[i][j] > 100) fprintf(fo3,"%c",(char)((int)(255))); else fprintf(fo3,"%c",(char)((int)(0))); } } //output the peaks image. for (i=0;i<PICSIZE;i++) { for (j=0;j<PICSIZE;j++) { if (peaks[i][j] == 255) fprintf(fo2,"%c",(char)((int)(255))); else fprintf(fo2,"%c",(char)((int)(0))); } } //close all of the files. fclose(fp1); fclose(fo1); fclose(fo2); fclose(fo3); system("PAUSE"); }

You should really split this all up into functions and test them one at a time. It makes it much easier. I would also just make a sample input that was easy to work with, and see if it finds what you want. Just take your picture array, and randomly populate it with some 1s and the rest 0s. Does it find all the 1s?

Quzah.