#include <iostream>
#include <string>
#include <iomanip>
#include <fstream>
#include <map>
#include <cmath>

using namespace std;

//struct for holding rgb color values
typedef struct {
int r,g,b;
} rgb;

const int NUMANGLES = 240;
const int NUMDIST = 240;
const double RELATIVE_LINE_THRESHOLD = 0.2; //relative percentage of the image's diagonal

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

string magicnum, comment;
comment = "#created by Andrew Stebbins";
int height, width, range;
int edge_threshold = 120;
int hysteresis_threshold = 90;
int line_threshhold = 0; //this will change as soon we get the width/height of the image
ifstream in;

cout << "opening file: " << argv[1] << endl;

in.open(argv[1], ios::in);

getline(in, magicnum);

if(!(magicnum == "P3")){
        cout << "error: image is not of type PPM!" << endl;
        return 0;
}

getline(in, comment);

in >> width;
in >> height;
in >> range;
cout << "width: " << width << " height: " << height << " range: " << range << endl; 

double MAXDIST = width*cos(3.14159/4.0) + height*sin(3.14159/4.0);
line_threshhold = (int)(sqrt((double)height*(double)height+(double)width*(double)width)*RELATIVE_LINE_THRESHOLD);
cout << "threshold: " << line_threshhold << endl;

cout << "initializing arrays... " << endl;

cout << "ppm" << endl;
rgb** ppm = new rgb*[height];
for(int k = 0; k < height; k++){
	ppm[k] = new rgb[width];
	for(int c = 0; c < width; c++){
		ppm[k][c].r = 0;
		ppm[k][c].g = 0;
		ppm[k][c].b = 0;
	}
}

cout << "pgm" << endl;
int** pgm = new int*[height];
for(int k = 0; k < height; k++){
        pgm[k] = new int[width];
	for(int c = 0; c < width; c++){
		pgm[k][c] = 0;
	}
}

cout << "smoothed" << endl;
int** smoothed = new int*[height];
for(int k = 0; k < height; k++){
        smoothed[k] = new int[width];
	for(int c = 0; c < width; c++){
		smoothed[k][c] = 0;
	}
}

cout << "edges" << endl;
int** edges = new int*[height];
for(int k = 0; k < height; k++){
        edges[k] = new int[width];
	for(int c = 0; c < width; c++){
		edges[k][c] = 0;
	}
}

int** gx = new int*[height];
for(int k = 0; k < height; k++)
        gx[k] = new int[width];
int** gy = new int*[height];
for(int k = 0; k < height; k++)
        gy[k] = new int[width];
int** edge_vals = new int*[height];
for(int k = 0; k < height; k++)
        edge_vals[k] = new int[width];

cout << "vertices" << endl;
int** vertices = new int*[height];
for(int k = 0; k < height; k++){
        vertices[k] = new int[width];
        for(int c = 0; c < width; c++){
                vertices[k][c] = 0;
        }
}

cout << "lines" << endl;
int** lines = new int*[height];
for(int k = 0; k < height; k++){
        lines[k] = new int[width];
	for(int c = 0; c < width; c++){
		lines[k][c] = 0;
	}
}

cout << "accum" << endl;
int** accumulator = new int*[NUMANGLES];
for(int k = 0; k < NUMANGLES; k++){
        accumulator[k] = new int[NUMDIST];
	for(int c = 0; c < NUMDIST; c++){
		accumulator[k][c] = 0;
	}
}

cout << "done." << endl;

cout << "reading file... " << endl;
//read the input file into an array of rgb data structures
for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
		in >> ppm[y][x].r;
		in >> ppm[y][x].g;
		in >> ppm[y][x].b;
	}
}	
in.close();

cout << "done.\ngenerating grayscale image... " << endl; 
//generate a greyscale representation of the input file
for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
		pgm[y][x] = (ppm[y][x].r + ppm[y][x].g + ppm[y][x].b)/3;
	}
}

cout << "done.\napplying gaussian filtering... " << endl;
//generate a flattened grayscale image using gaussian filtering
for(int y = 1; y < height-1; y++){
        for(int x=1; x < width-1; x++){
		smoothed[y][x]=(pgm[y][x]*4+pgm[y-1][x]*2+pgm[y+1][x]*2+pgm[y][x-1]*2+pgm[y][x+1]*2+pgm[y-1][x+1]*1+pgm[y-1][x-1]*1+pgm[y+1][x+1]*1+pgm[y+1][x-1]*1)/16;
	}
	}
cout << "done.\ndetecting edges... " << endl;
//detect edges
for(int y = 2; y < height-2; y++){
        for(int x=2; x < width-2; x++){
		if(y>0){
			gx[y][x] = smoothed[y][x-1]*(-2)+smoothed[y-1][x-1]*(-1)+smoothed[y+1][x-1]*(-1)+smoothed[y+1][x]*(2)+smoothed[y+1][x-1]*(1)+smoothed[y+1][x+1]*(1);
                        gy[y][x] = smoothed[y-1][x]*(-2)+smoothed[y-1][x-1]*(-1)+smoothed[y-1][x+1]*(-1)+smoothed[y+1][x]*(2)+smoothed[y+1][x-1]*(1)+smoothed[y+1][x+1]*(1);
			if(abs(gx[y][x])+abs(gy[y][x])>edge_threshold)
	                        edges[y][x] = abs(gx[y][x])+abs(gy[y][x]);
		}
	}
}

//apply non-maximum suppression
for(int y = 2; y < height-2; y++){
        for(int x = 2; x < width-2; x++){
                if(abs(gx[y][x])>abs(gy[y][x])){
                        if(gx[y][x]<0){
                                if(gx[y][x-1]<gx[y][x] || gx[y][x+1]<gx[y][x])
                                        edges[y][x] = 0;
                        }
                        else {
                                if(gx[y][x+1]>gx[y][x] || gx[y][x-1]>gx[y][x])
                                        edges[y][x] = 0;
                        }
                }
                else if(abs(gx[y][x])<abs(gy[y][x])){
                        if(gy[y][x]<0){
                                if(gy[y-1][x]<gy[y][x] || gy[y+1][x]<gy[y][x])
                                        edges[y][x] = 0;
                        }
                        else {
                                if(gy[y+1][x]>gy[y][x] || gy[y-1][x]>gy[y][x])
                                        edges[y][x] = 0;
                        }
                }
        }
}

//hysteresis

cout << "done.\nfilling hough space... " << endl;

//fill hough space
for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
		if(edges[y][x] > edge_threshold){
			for(int a = 0; a < NUMANGLES; a++){
				double ad = 360.0*(double)a/(double)NUMANGLES;
				double ar = (double)(3.14159)*((double)(ad/180.0));
				int d = (int)((double)NUMDIST * (((double)x*cos(ar)+(double)y*sin(ar))/(double)MAXDIST));
				accumulator[a][abs(d)]++;
			}	
		}
	}
}

cout << "done.\ninitializing array to store line data... " << endl;

//initialize array to store detected lines
for(int y = 0; y < height; y++){
	for(int x=0; x < width; x++){
		lines[y][x] = 0;	
	}
}

cout << "done.\ndetecting lines... " << endl;
//use hough space accumulator to detect lines
for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
			for(int a = 0; a < NUMANGLES; a++){
				double ad = 360.0*(double)a/(double)NUMANGLES;
				double ar = (double)(3.14159)*((double)(ad/180.0));
				int d = (int)((double)NUMDIST * (((double)x*cos(ar)+(double)y*sin(ar))/(double)MAXDIST));
				if(accumulator[a][abs(d)] > line_threshhold){
					lines[y][x] = 1;
				}
			}	
	}
}

cout << "done.\ndetecting vertices... " << endl;

for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
                        for(int a = 0; a < NUMANGLES; a++){
                                double ad = 360.0*(double)a/(double)NUMANGLES;
                                double ar = (double)(3.14159)*((double)(ad/180.0));
                                int d = (int)((double)NUMDIST * (((double)x*cos(ar)+(double)y*sin(ar))/(double)MAXDIST));
				if(accumulator[a][abs(d)] > line_threshhold){
					for(int a2 = 0; a2 < NUMANGLES; a2++){
                        		        double ad2 = 360.0*(double)a2/(double)NUMANGLES;
                	                	double ar2 = (double)(3.14159)*((double)(ad2/180.0));
		                                int d2 = (int)((double)NUMDIST * (((double)x*cos(ar2)+(double)y*sin(ar2))/(double)MAXDIST));
						if(accumulator[a2][abs(d2)] > line_threshhold){
							if(abs(ad - ad2) >= 75 && abs(ad - ad2) <= 105)
								vertices[y][x] = 1;
						}
					}
				}
                        }
        }
}

cout << "done.\nbeginning output... " << endl;
//output image of Hough space
fstream outh;
outh.open("houghspace.pgm", fstream::out);
outh << "P2" << endl;
outh << comment << endl;
outh << NUMANGLES-1 << " " << NUMDIST-1 << endl;
outh << range << endl;

for(int y = 0; y < NUMANGLES; y++){
        for(int x=0; x < NUMDIST; x++){
		outh << accumulator[y][x] << endl;
        }
}
outh.close();

//output the smoothed grayscale image with no edge superposition
fstream out1;
out1.open("smoothed.pgm", fstream::out);
out1 << "P2" << endl;
out1 << comment << endl;
out1 << width << " " << height << endl;
out1 << range << endl;

for(int y = 0; y < height; y++){
        for(int x = 0; x < width; x++){
		out1 << smoothed[y][x] << endl;	
	}
}
out1.close();

//output the detected edges
fstream out2;
out2.open("edges.ppm", fstream::out);

out2 << magicnum << endl;
out2 << comment << endl;
out2 << width << " " << height << endl;
out2 << range << endl;

for(int y = 0; y < height; y++){
	for(int x = 0; x < width; x++){
			out2 << edges[y][x] << endl;
                        out2 << edges[y][x] << endl;
                        out2 << edges[y][x] << endl;
	}
}	
out2.close();

//output the detected lines
fstream outl;
outl.open("lines.ppm", fstream::out);

outl << magicnum << endl;
outl << comment << endl;
outl << width << " " << height << endl;
outl << range << endl;

for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
                //if pixel is on a detected line, make red
                if(lines[y][x] == 1){
                        outl << 255 << endl;
                        outl << 0 << endl;
                        outl << 0 << endl;
                }
                else{
                      	outl << ppm[y][x].r << endl;
                        outl << ppm[y][x].g << endl;
                        outl << ppm[y][x].b << endl;
                }
        }
}
outl.close();

//output the detected lines
fstream outv;
outv.open("vertices.ppm", fstream::out);

outv << magicnum << endl;
outv << comment << endl;
outv << width << " " << height << endl;
outv << range << endl;

for(int y = 0; y < height; y++){
        for(int x=0; x < width; x++){
                //if pixel is on a detected line, make red
                if(vertices[y][x] == 1){
                        outv << 255 << endl;
                        outv << 0 << endl;
                        outv << 0 << endl;
                }
                else{
                        outv << ppm[y][x].r << endl;
                        outv << ppm[y][x].g << endl;
                        outv << ppm[y][x].b << endl;
                }
        }
}
outv.close();

cout << "done.\n";
return 0;
}