/*      Author:         Jonathan Hardy
        file:           hough.cc
 
        This file reads in a PGM image file of an image processed by an
        edge detector algorithm.  The Hough line detector algorithm is 
	performed on this data, and the detected lines are printed to the
	screen.  An image of the Hough space accumulator is written out to
        a new PGM image file with "_hough" appended to the original filename.
        The arguments are as follows:
        roberts <filename> <inputthreshold> <outputthreshold> 
*/
#include <stdlib.h> 
#include <stdio.h>
#include <string.h>
#include <math.h>

FILE *ptr;				// Image File
FILE *houghptr;				// Output File
void main(int argc, char**argv) 	
{
char fname[13];				// Image filename
char houghname[18];			// Output filename
char width[4], height[4], numcol[4];	// Image information
char Type[3];				// Image type (P5)
int r, c, data, accum;			// Counters and image data variables
char temp[100];				// Variable for input comment lines
char comment[17];			// Variable for output comment
char height_out[5];			// Variable for output height

/*      The following section reads the image file      */
strcpy(fname, argv[1]);			
strcat(fname, ".pgm");
int inthresh = atoi(argv[2]);
int outthresh = atoi(argv[3]);
ptr = fopen(fname, "rb");
fread(&Type, 3, 1, ptr);
fread(&width, 4, 1, ptr);
while (width[0] == '#'){
	fgets(temp, 100, ptr);
	fread(&width, 4, 1, ptr);
}						
fread(&height, 4, 1, ptr);
const int WIDTH = atoi(width);
const int HEIGHT = atoi(height);
int ImgArray[HEIGHT][WIDTH];		// Create image array
fread(&numcol, 4, 1, ptr);
for (r = 0; r < HEIGHT; r++){		// Read in data and threshold
	for (c = 0; c < WIDTH; c++){
		data = fgetc(ptr);	
		if (data > inthresh)
			ImgArray[r][c] = 255;
		else
			ImgArray[r][c] = 0;
	}
}
fclose(ptr);				// Close image file

/*  The following section performs the Hough algorithm on the image data  */

const float PI = 3.1415925;
const int rmax = WIDTH;
int HoughArray[181][rmax];		// Create Hough space
int rad, theta;
accum = 0;
double conv = 3.1415926535/180.0;	
for (theta=0; theta<181; theta++){	// Initialize Hough space to zero
	for (rad=0; rad<rmax; rad++){
		HoughArray[theta][rad] = 0;
	}				
}

/*	This section scans the image array for any pixel that is non-zero.
	For each pixel found, the radius is found for each line of angle 
	theta from 0 to 180 for lines of form r = ysin(theta) + xcos(theta)
*/
 
for (r=0; r<HEIGHT; r++){
        for (c=0; c<WIDTH; c++){
		if (ImgArray[r][c] != 0){
		for (theta=0; theta<181; theta++){
			rad = (int)(r*cos((double)(theta*conv)) +
				c*sin((double)(theta*conv)));
			if (rad != 0){
			if (HoughArray[theta][(int)((rad+HEIGHT)/2)] < 255) 
			HoughArray[theta][(int)((rad+HEIGHT)/2)] += 1;
			if (HoughArray[theta][(int)((rad+HEIGHT)/2)]>accum)
			accum=HoughArray[theta][(int)((rad+HEIGHT)/2)];
		}}}
        }
}

/*	This section finds local maxima in the Hough space and calculates 
	the endpoints for the lines represented.  These enpoints are 
	printed to the screen.
*/

int line = 0;
int endx1, endy1;
for (theta=0; theta<180; theta++){
for (rad=0; rad<rmax; rad++){
        if (HoughArray[theta][rad] > outthresh){
                line = 1;
                for (int r1=theta; r1<theta+7; r1++){
                for (int c1=rad; c1<rad+7; c1++){
                        if (HoughArray[theta][rad] < HoughArray[r1-3][c1-3]){
                                line = 0;
                                break;
                        }
                }}
        }
        if (line){
		printf("Line Found!\n");
		printf("Angle: %d\n", theta);
		printf("Endpoints: \n");
		if (theta != 90)
			endy1 = (int)((2*rad-HEIGHT)/cos(theta*conv));
			if (endy1 > 0 && endy1 < HEIGHT-1)
				printf("0,%d\n", endy1);
			endy1 = (int)((2*rad-HEIGHT-(WIDTH-1)*sin(theta*conv))
					/ cos(theta*conv));
			if (endy1 > 0 && endy1 < HEIGHT-1)
				printf("%d,%d\n",WIDTH-1,endy1);
		if ((theta != 0) && (theta != 180))
			endx1 = (int)((2*rad-HEIGHT)/sin(theta*conv));
			if (endx1 > 0 && endx1 < WIDTH-1)
				printf("%d,0\n", endx1);
			endx1 = (int)((2*rad-HEIGHT-(HEIGHT-1)*cos(theta*conv)) 
					/ sin(theta*conv));
			if (endx1 > 0 && endx1 < WIDTH-1)
				printf("%d,%d\n",endx1,HEIGHT-1);
	printf("\n");
	}
        line = 0;
}}

/*      This section writes the Hough space image to a new PGM file        */

	strcpy(houghname, argv[1]);
	strcat(houghname, "_hough.pgm");
        houghptr = fopen(houghname, "wb");
        fwrite(&Type, 3, 1, houghptr);
        comment = "# New pgm file.\n";
	height_out = "180\n";
        fwrite(&comment, 16, 1, houghptr);
        fwrite(&width, 4, 1, houghptr);
        fwrite(&height_out, 4, 1, houghptr);
        fwrite(&numcol, 4, 1, houghptr);
        for (r = 0; r < 180; r++){
                for (c = 0; c < rmax; c++){
                        fputc(HoughArray[r][c]*255/accum, houghptr);
                }
        }
        fclose(houghptr);
}