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");
}