/*      Author:         Jonathan Hardy
        file:           canny.cc
 
        This file reads in a PGM image file and processes it with the Canny
        edge detector algorithm.  The processed image is then written out to
        a new PGM image file with "_cnny" appended to the original filename.
        The arguments are as follows:
        canny <filename> <upperthresh> <lowerthresh> <filterwidth> <sigma>
*/

#include <stdlib.h> 
#include <stdio.h> 
#include <string.h> 
#include <math.h>

FILE *ptr;				// Image File
FILE *cnnyptr;				// Output File

void main(int argc, char**argv) 
{
char fname[13];				// Image filename
char cnnyname[18];			// Output filename
char width[4], height[4], numcol[4];	// Image information
char Type[3];				// Image type (P5)
int n, r, c, data;			// Counters and image data variable
char temp[100];				// Variable for input comment lines
char comment[] = "# New pgm file.\n";

const int ulimit = atoi(argv[2]);	// Get upper threshold
const int llimit = atoi(argv[3]);	// Get lower threshold
int filter_width = atoi(argv[4]);	// Get filter width
float sigma = atoi(argv[5]);		// Get sigma

/*      The following section reads the image file      */

strcpy(fname, argv[1]);			// Get filename from command line
strcat(fname, ".pgm");			// Append filename with ".pgm"
ptr = fopen(fname, "rb");		// Open image file
fread(&Type, 3, 1, ptr);		// Get image type
fread(&width, 4, 1, ptr);		// If no comment lines, get width
while (width[0] == '#'){		//  Otherwise, ignore comment lines
        fgets(temp, 100, ptr);		//  and get width
        fread(&width, 4, 1, ptr);
}
fread(&height, 4, 1, ptr);		// Get height
const int WIDTH = atoi(width);
const int HEIGHT = atoi(height);
int ImgArray[HEIGHT][WIDTH];		// Create image array
fread(&numcol, 4, 1, ptr);		// Get number of colors (255)
for (r = 0; r < HEIGHT; r++){		// Get image data and threshold it
	for (c = 0; c < WIDTH; c++){
		data = fgetc(ptr);
		if ((data < llimit) || (data > ulimit))
			ImgArray[r][c] = 0;
		else
			ImgArray[r][c] = data;
						
	}
}
fclose(ptr);				// Close image file

/*      The following iteration scans the image array and performs the 
	Canny algorithm on the data.  It calculates the magnitude
        response by taking the square root of the sum of the squares of 
	the dx and dy responses. 
*/

	float SmoothX[HEIGHT][WIDTH];	// Array for x-smoothed image
	float SmoothY[HEIGHT][WIDTH];	// Array for y-smoothed image
	float DiffX[HEIGHT][WIDTH];	// Array for x-differentiated image
	double t1, t2, ux, uy;		// Temp variables 
	double atmp, btmp, ctmp, dtmp, g0;	
	double gc, gn, gs, gw, ge, gne, gse, gsw, gnw;
/*	Arrays for the Gaussian and its derivatives	*/
	double G[filter_width], dG[filter_width], D2G[filter_width];

/*	This section calculates the values of the Gaussian and its derivatives
	based on the filter width and sigma provided in the command line
*/
	for (n=0; n<filter_width; n++){
		atmp = exp((-n*n)/(2*sigma*sigma));
		if (atmp>0.005 || n<2){
			btmp = exp((-((double)n-0.5)*((double)n-0.5)) 
				/ (2*sigma*sigma));
			ctmp = exp((-((double)n+0.5)*((double)n+0.5))
				/ (2*sigma*sigma));
			dtmp = exp((-n*n)/(0.5*sigma*sigma));
			G[n] = (atmp+btmp+ctmp)/3/(6.283185*sigma*sigma);
			dG[n] = ctmp - btmp;
			D2G[n] = 1.6*dtmp - atmp;
		}
		else break;
	}			
	filter_width = n;

/*	Set limits to avoid the image edges	*/

	int climit = WIDTH - (filter_width-1);
	int rlimit = HEIGHT - (filter_width-1);

/*	This section performs the Gaussian smoothing	*/

	for (c=filter_width-1; c<climit; c++){
		for (r=filter_width-1; r<rlimit; r++){
			t1 = ImgArray[r][c] * G[0];
			t2 = t1;
			for (n=1; n<filter_width; n++){
				t1 += G[n]*(ImgArray[r-n][c] +
					ImgArray[r+n][c]);
				t2 += G[n]*(ImgArray[r][c-n] +
					ImgArray[r][c+n]);
			}
			SmoothX[r][c] = t1;
			SmoothY[r][c] = t2;	// line 100
		}
	}

/*	This section performs the convolution of the smoothed image with 
	the derivative of the Gaussian
*/
	for (c=filter_width-1; c<climit; c++){
		for (r=filter_width-1; r<rlimit; r++){
			t1 = 0;
			for (n=1; n<filter_width; n++){
				t1 += dG[n]*(SmoothX[r][c-n] -
					SmoothX[r][c+n]);
			}
			DiffX[r][c] = t1;
		}
	}

	float DiffY[HEIGHT][WIDTH];
	climit = WIDTH - filter_width;
	for (c=filter_width; c<climit; c++){
		for (r=filter_width-1; r<rlimit; r++){
			t2 = 0;
			for (n=1; n<filter_width; n++){
				t2 += dG[n]*(SmoothY[r-n][c] -
					SmoothY[r+n][c]);
			}
			DiffY[r][c] = t2;
		}
	}

	int CnnyArray[HEIGHT][WIDTH];
	int DiffCnnyArray[HEIGHT][WIDTH];

	float nvlx, svlx, wvlx, evlx, nwvlx, nevlx, sevlx, swvlx;
	float nvly, svly, wvly, evly, nwvly, nevly, sevly, swvly;
	rlimit = HEIGHT - filter_width;	
	for (c=filter_width; c<climit; c++){
		for (r=filter_width; r<rlimit; r++){
			ux = DiffX[r][c];
			uy = DiffY[r][c];
			nvlx = DiffX[r-1][c];
			svlx = DiffX[r+1][c];
			wvlx = DiffX[r][c-1];
			evlx = DiffX[r][c+1];
			nwvlx = DiffX[r-1][c-1];
			nevlx = DiffX[r-1][c+1];
			sevlx = DiffX[r+1][c+1];
			swvlx = DiffX[r+1][c-1];	
			nvly = DiffY[r-1][c];
			svly = DiffY[r+1][c];
			wvly = DiffY[r][c-1];
			evly = DiffY[r][c+1];
			nwvly = DiffY[r-1][c-1];
			nevly = DiffY[r-1][c+1];
			sevly = DiffY[r+1][c+1];
			swvly = DiffY[r+1][c-1];	
			gc = sqrt(ux*ux + uy*uy);
			float t = gc*20.0;
			DiffCnnyArray[r][c] =(int)(t<256 ? t : 255);

/*	These lines find the magnitude response for a 3x3 neighborhood of 
	pixels
*/
			gn = sqrt(nvlx*nvlx + nvly*nvly);
			gs = sqrt(svlx*svlx + svly*svly);
			gw = sqrt(wvlx*wvlx + wvly*wvly);
			ge = sqrt(evlx*evlx + evly*evly);
			gne = sqrt(nevlx*nevlx + nevly*nevly);
			gse = sqrt(sevlx*sevlx + sevly*sevly);
			gsw = sqrt(swvlx*swvlx + swvly*swvly);
			gnw = sqrt(nwvlx*nwvlx + nwvly*nwvly);

/*	These lines find the local maxima	*/

			if (ux*uy>0){
				if (fabs(ux)<fabs(uy)){
					if((g0=fabs(uy*gc)) < 
					fabs(ux*gse+(uy-ux)*gs) ||
					g0<=fabs(ux*gnw+(uy-ux)*gn))
						continue;
					}
				else{
					if((g0=fabs(ux*gc)) <
					fabs(uy*gse+(ux-uy)*ge) ||
					g0<=fabs(uy*gnw+(ux-uy)*gw))
						continue;
				}
			}
			else{
				if (fabs(ux)<fabs(uy)){
					if((g0=fabs(uy*gc)) <
					fabs(ux*gne-(uy+ux)*gn) ||
					g0<=fabs(uy*gsw-(ux+uy)*gw))
						continue;
				}
				else{
					if((g0=fabs(ux*gc)) <
					fabs(uy*gne-(ux+uy)*ge) ||
					g0<=fabs(uy*gsw-(ux+uy)*gw))
						continue;
				}
			}
			CnnyArray[r][c] = DiffCnnyArray[r][c];
		}
	}

/*      This section writes the output to a new PGM file        */

	strcpy(cnnyname, argv[1]);
	strcat(cnnyname, "_cnny.pgm");
	cnnyptr = fopen(cnnyname, "wb");
	if (cnnyptr == (FILE*)NULL) {
	  fprintf(stdout, "Error opening file %s\n", cnnyname);
	  exit(-1);
	}

	fwrite(&Type, 3, 1, cnnyptr);	
        fwrite(&comment, 16, 1, cnnyptr);
        fwrite(&width, 4, 1, cnnyptr);
        fwrite(&height, 4, 1, cnnyptr);
        fwrite(&numcol, 4, 1, cnnyptr);

	for (r = 0; r < HEIGHT; r++){
	  for (c = 0; c < WIDTH; c++){
	    fputc(CnnyArray[r][c], cnnyptr);
	  }
	}
	fclose(cnnyptr); 
	
}