/*****************************************************
Program name:  	coins_detect
Description:
	read an input image and find US coins in this
	image, count this coins
Usage:
    coins_detect
To compile (on a Sun workstation):
    g++ -Wall -o coins_detect coins_detect.c -ltiff
*****************************************************/

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

#define d1             169                         // diameter of the 1 cent coin in pixels
#define d5             185                         // diameter of the 5 cent coin in pixels
#define d10            153                         // diameter of the 10 cent coin in pixels
#define d25            205                         // diameter of the 25 cent coin in pixels
#define t1             80                          // threshold for the 1 cent coin
#define t5             90                          // threshold for the 5 cent coin
#define t10            80                          // threshold for the 10 cent coin
#define t25            90                          // threshold for the 25 cent coin
#define size_connect   5                           // 1/2 of size of squares for regions_connect
#define size_search    20                          // 1/2 of size of squares for search_maximum
#define sig			   2                           // sigma for the canny edge detector
#define th             100                         // higher threshold for canny edge detector
#define tl             40                          // lower threshold for canny edge detector
#define GV             256                         // number of gray values
#define pi             3.141592654                 // value of pi


struct edge                                        // struct to hold magnitude and direction
{                                                  // of edges
    double mag;
    int dr, dc;
};

void canny_edge_detect (unsigned char *in_image, double *out_image, float sigma, float tlow,
						float thigh, uint32 rows, uint32 columns);
void print_error (char *err_message);
void gaussian_filtering (unsigned char *in, double *out, float sigma, uint32 rows,
						 uint32 columns);
void roberts (double *in, struct edge *out, uint32 rows, uint32 columns);
void scale (struct edge *in, uint32 rows, uint32 columns);
void suppress_nonmaxima (struct edge *in, double *out, uint32 rows, uint32 columns);
void follow_edge (uint32 er, uint32 ec, double *in, double *out, uint32 rows, uint32 columns,
				  float tlow);
void hough_transform_circles (double *in, int *hough_space, uint32 rows, uint32 columns);
void hough_threshold (int *hough_space, uint32 rows, uint32 columns);
void regions_connect (int *hough_space, int *hough_connect, uint32 rows, uint32 columns);
void hough_label (int *hough_space, int *hough_connect, int label, int *sum, double *rc,
				  double *cc, int *n, uint32 r, uint32 c, int i, uint32 rows, uint32 columns);
void search_maximum (int *hough_space, uint32 rows, uint32 columns);
float count_coins (int *hough_space, uint32 rows, uint32 columns);
void generate_output_image(unsigned char *in_image, unsigned char *out_image, int *hough_space, uint32 rows, uint32 columns);


int main()
{
    uint32 r;        	            		   // row index
    uint32 c; 	       	            		   // column index
    uint32 rows;           		        	   // number of rows in image
    uint32 columns;		                       // number of columns in image
    uint16 BitsPerSample;              		   // normally 8 for grayscale image
    uint16 SamplesPerPixel;    	               // normally 1 for grayscale image
    uint16 PhotoMetric;            		       // normally 1 for grayscale image
    unsigned char *in_image;            	   // pointer for input image array
    unsigned char *out_image;            	   // pointer for output image array
    TIFF *in_filep;                			   // handle for input image file
    TIFF *out_filep;                           // handle for output image file
    int	check_io;           		           // status of I/O operation
    char input[20];                            // input string
    char in_name[24];                          // input file name
    char out_name[28];                         // output file name
	char mess[55];                             // error message
	double *tmp1_image;                        // holds temporary image 1
	int *hough_space;                          // parameter space for the Hough transform
	int *hough_connect;                        // holds the components of the Hough space
	                                           // belonging together
	int label;                                 // holds the next label to be applied
	                                           // to a connected component
	int sum;                                   // sum of the values in the Hough space
	                                           // belonging to one component
	int n;                                     // number of pixels belonging to one component
	int i;                                     // parameter for the diameter in the Hough space
	double rc, cc;                             // row and column index of the centroid
	uint32 r_centroid, c_centroid;             // row and column index of the centroid
	float coins;                               // value of the detected coins




    // input of the input file name
	printf("Input file name (without .tif):      ");
	scanf("%s", input);
    strcpy(in_name, input);
	strcat(in_name, ".tif");
    strcpy(out_name, input);
	strcat(out_name, "_out.tif");


	// Open input image file
    in_filep = TIFFOpen(in_name, "r");
    if (in_filep == NULL)
	{
	    strcpy(mess, "Could not open input file ");
		strcat(mess, in_name);
        print_error (mess);
	}

    // Open output image file
    out_filep = TIFFOpen (out_name, "w"); 
    if (out_filep == NULL)
	{
	    strcpy(mess, "Could not open output file ");
		strcat(mess, out_name);
        print_error (mess);
	}


    // Determine the size of the input image
    TIFFGetField(in_filep, TIFFTAG_IMAGELENGTH, &rows);
    TIFFGetField(in_filep, TIFFTAG_IMAGEWIDTH, &columns);
    TIFFGetField(in_filep, TIFFTAG_BITSPERSAMPLE, &BitsPerSample);
    TIFFGetField(in_filep, TIFFTAG_SAMPLESPERPIXEL, &SamplesPerPixel);
    TIFFGetField(in_filep, TIFFTAG_PHOTOMETRIC, &PhotoMetric);

    // Specify TIFF header fields for output image
    TIFFSetField(out_filep, TIFFTAG_IMAGELENGTH, rows);
    TIFFSetField(out_filep, TIFFTAG_IMAGEWIDTH, columns);
    TIFFSetField(out_filep, TIFFTAG_BITSPERSAMPLE, 8);
    TIFFSetField(out_filep, TIFFTAG_SAMPLESPERPIXEL, 1);
    TIFFSetField(out_filep, TIFFTAG_PHOTOMETRIC, 1);
    TIFFSetField(out_filep, TIFFTAG_PLANARCONFIG, 1);


    // Allocate memory to hold image arrays
    in_image = (unsigned char *) _TIFFmalloc(rows*columns);
    if (in_image == NULL) 
	    print_error ("Could not allocate memory!");
    out_image = (unsigned char *) _TIFFmalloc(rows*columns);
    if (out_image == NULL) 
	    print_error ("Could not allocate memory!");
	tmp1_image = (double *) calloc(rows*columns, sizeof(double));
    if(tmp1_image == NULL)
	    print_error ("Could not allocate memory!");

    // Allocate memory to hold Hough spaces
	hough_space = (int *) calloc(4*rows*columns, sizeof(int));
    if(hough_space == NULL)
	    print_error ("Could not allocate memory!");
	hough_connect = (int *) calloc(4*rows*columns, sizeof(int));
    if(hough_connect == NULL)
	    print_error ("Could not allocate memory!");

	// Read image pixel values from file, row by row 
    for (r = 0; r < rows; r++)
    {
        check_io = TIFFReadScanline(in_filep, &in_image[r*columns], r, 1);
        if (check_io != 1)
	        print_error ("Could not read image from file!");
    }


	///////////////////////
    // Process the image //
	///////////////////////

	// detect edges using a canny edge detector
    canny_edge_detect (in_image, tmp1_image, sig, tl, th, rows, columns);

	// applying the Hough transform for circles
    hough_transform_circles (tmp1_image, hough_space, rows, columns);

	// set values less than or equal to the threshold to 0
	hough_threshold (hough_space, rows, columns);

	// connect the regions in the Hough space
    regions_connect (hough_space, hough_connect, rows, columns);

	// label the regions and compute centroid and add the values
	// in the Hough space belonging to these region
	printf("labeling\n");
	label = 0;
	for (i = 0; i < 4; i++)
		for (r = 0; r < rows; r++)
		    for (c = 0; c < columns; c++)
				if (hough_connect[i*rows*columns+r*columns+c] == -1)
				{
					label++;
					rc = 0;
					cc = 0;
					sum = 0;
					n = 0;
					hough_label (hough_space, hough_connect, label, &sum, &rc, &cc, &n,
						         r, c, i, rows, columns);
					r_centroid = rc;
					c_centroid = cc;
					hough_space[i*rows*columns+r_centroid*columns+c_centroid] = sum;
				}

	// search for maximum
	search_maximum (hough_space, rows, columns);

	// count coins
	coins = count_coins (hough_space, rows, columns);

	printf("\n");
	printf("detected value:             %.2f $\n", coins);
	printf("\n");

	// generate output image
	generate_output_image(in_image, out_image, hough_space, rows, columns);


    // Write output image to file, row by row 
    for (r = 0; r < rows; r++)
    {
        check_io = TIFFWriteScanline(out_filep, &out_image[r*columns], r, 1);
        if (check_io != 1)
	        print_error ("Could not write image to file!");
    }


    // Deallocate the image memory, close file streams, and exit
    _TIFFfree(in_image);
    _TIFFfree(out_image);
    free(tmp1_image);
	free(hough_space);
	free(hough_connect);
    TIFFClose(in_filep);
    TIFFClose(out_filep);

	printf("ready\n");

    exit (0);
	return 0;
}



/*****************************************************
Function name:   canny_edge_detect
Description:
    detects edges in an input image using a Canny edge
	detector and stores the result in an output image
Input parameters:
    in_image      --  pointer to input image
    out_image     --  pointer to output image
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void canny_edge_detect (unsigned char *in_image, double *out_image, float sigma, float tlow,
						float thigh, uint32 rows, uint32 columns)
{
    uint32 r;       		           		   // row index
    uint32 c; 	        	           		   // column index
    double *tmp1_image;                        // holds temporary image 1
    struct edge *tmp2_image;                   // holds temporary image 2


    // allocate memory to hold temporary images
	tmp1_image = (double *) calloc(rows*columns, sizeof(double));
    if(tmp1_image == NULL)
	    print_error ("Could not allocate memory!");
    tmp2_image = (struct edge *) calloc(rows*columns, sizeof(struct edge));
    if(tmp2_image == NULL)
	    print_error ("Could not allocate memory!");

    // initialize output image
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			out_image[r*columns+c] = 0;

    // applying of the Gaussian mask to the image
	gaussian_filtering (in_image, tmp1_image, sigma, rows, columns);

    // computing of magnitude and direction of edges using the Roberts operator
    roberts (tmp1_image, tmp2_image, rows, columns);

	// scale the temporary image
    scale (tmp2_image, rows, columns);

	// suppress nonmaxima
    suppress_nonmaxima (tmp2_image, tmp1_image, rows, columns);

	// edge detect
	printf("edge detect\n");
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			if(tmp1_image[r*columns+c] >= thigh)     // if pixel value over the higher threshold,
				follow_edge(r, c, tmp1_image, out_image, rows, columns, tlow); // follow the edge

    // deallocate the memory
    free(tmp1_image);
    free(tmp2_image);
}


/*****************************************************
Function name:   print_error
Description:
    Print an error message, and exit the program.
Input parameters:
    err_message -- pointer to string to be printed
Returned value:
    none
*****************************************************/
void print_error (char *err_message)
{
    fprintf(stderr, "Error:  %s\n", err_message);
    exit (1);
}


/*****************************************************
Function name:   gaussian_filtering
Description:
    applies a Gaussian filter to an image
Input parameters:
    in            --  pointer to input image
    out           --  pointer to output image
	sigma         --  sigma of the gaussian mask
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void gaussian_filtering (unsigned char *in, double *out, float sigma, uint32 rows,
						 uint32 columns)
{
    uint32 r;       	            		   // row index
    uint32 c; 	       	            		   // column index
    int hr, hc;                                // holds the indexes used for the computation
	float sum;                                 // result of the Gaussian filter at one pixel
    int s;                                     // size of the Gaussian filter
	int x;                                     // coordinate for Gaussian mask
	int gc;                                    // column index of the Gaussian mask
    double *g_filter;                          // holds the mask for the Gaussian filter
    double *tmpg_image;                        // holds temporary image


	// calculate the size of the for the Gaussian filter
	s = 8*sqrt(2)*sigma;
	if (div(s,2).rem == 0)
		s = s+1;

    if (((s+1)/2 > 2*rows) || ((s+1)/2 > 2*columns))
        print_error ("Gaussian mask too large!");

	printf("size of the Gaussian filter:         %d X %d\n", s, s);

    // allocate memory
    tmpg_image = (double *) calloc(rows*columns, sizeof(double));
    if(tmpg_image == NULL)
	    print_error ("Could not allocate memory!");
    g_filter = (double *) calloc(s, sizeof(double));
    if(g_filter == NULL)
	    print_error ("Could not allocate memory!");

	printf("computing the Gaussian mask\n");

    // compute the mask for the Gaussian filter
	for (gc = 0; gc < s; gc++)
		{
			x = gc - (s-1)/2;

            g_filter[gc]= 1/(sqrt(2*pi) * sigma) * exp(-pow(x,2) / (2*pow(sigma,2)));
		}

	printf("applying the Gaussian mask\n");

    // applying of the Gaussian mask to the image
	// horizontal
    for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
		{
			sum = 0;

            for (gc = 0; gc < s; gc++)
				{
					hc = c + gc - (s-1)/2;             // index of the input image

					if (hc < 0)                        // if the index is not in the input
						hc = -hc -1;                   // image, reflect the image (virtually)
					if (hc > (columns - 1))
						hc = columns -(hc-columns) -1;

                    sum += g_filter[gc] * in[r*columns+hc];     // result of the Gaussian
				}                                               // mask at one pixel

			tmpg_image[r*columns+c] = sum;           // storing the result to an temporary image
		}

	// vertical	
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
		{
			sum = 0;

            for (gc = 0; gc < s; gc++)
				{
					hr = r + gc - (s-1)/2;

					if (hr < 0)                        // if the index is not in the input
						hr = -hr -1;                   // image, reflect the image (virtually)
					if (hr > (rows - 1))
						hr = rows -(hr-rows) -1;

                    sum += g_filter[gc] * tmpg_image[hr*columns+c];  // result of the Gaussian
				}                                                    // mask at one pixel

			out[r*columns+c] = sum;                  // storing the result to the output image
		}

    // deallocate memory
    free(tmpg_image);
    free(g_filter);
}


/*****************************************************
Function name:   roberts
Description:
    applies Roberts operator to an image and stores
	the magnitude and dirction in an output image
Input parameters:
    in            --  pointer to input image
    out           --  pointer to output image
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void roberts (double *in, struct edge *out, uint32 rows, uint32 columns)
{
    uint32 r;      		            		   // row index
    uint32 c; 	        	           		   // column index
    int hr, hc;                                // holds the indexes used for the computation
	float sum1, sum2;                          // result of the Roberts operator at one pixel
	double alpha;                              // holds  the angle
	int rr, rc;                                // row and column index of the Roberts operator
    const int roberts1[2][2]                   // holds the 1. Roberts operator
		= {{1,0},{0,-1}};
    const int roberts2[2][2]                   // holds the 2. Roberts operator
		= {{0,-1},{1,0}};
	
	
	printf("computing magnitude and direction\n");

    // computing of magnitude and direction of edges using the Roberts operator
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
		{
			sum1 = 0;
			sum2 = 0;

			for (rr = 0; rr < 2; rr++)                   // indexes of the Roberts operator
				for (rc = 0; rc < 2; rc++)
				{
					hc = c + rc;                         // indexes of the temporary image
					hr = r + rr;

					if (hc > (columns -1))
						hc = columns - 1;
					if (hr > (rows -1))
						hr = rows - 1;

					sum1 += roberts1[rr][rc] * in[hr*columns+hc]; // results of applying
					sum2 += roberts2[rr][rc] * in[hr*columns+hc]; // the Roberts operators
				}

			out[r*columns+c].mag = sqrt(pow(sum1,2)+pow(sum2,2));  // storing the magnitude
                                                                   // to temporary image 2
			if (sum2 != 0)                       // calculation of the direction
			    alpha = atan(sum1/sum2);
			else
				alpha = pi/2;

			if ((alpha < -3*pi/8) || (alpha >= 3*pi/8))   // storing the direction to temporary
			{                                             // image 2 in terms of +1, 0, -1
			    out[r*columns+c].dr = 1;           // for rows and columns
			    out[r*columns+c].dc = 1;           // => 4 directions between -pi/2 and +pi/2
			}
			if ((alpha >= -3*pi/8) && (alpha < -pi/8))
			{
			    out[r*columns+c].dr = 1;
			    out[r*columns+c].dc = 0;
			}
			if ((alpha >= -pi/8) && (alpha < pi/8))
			{
			    out[r*columns+c].dr = 1;
			    out[r*columns+c].dc = -1;
			}
			if ((alpha >= pi/8) && (alpha < 3*pi/8))
			{
			    out[r*columns+c].dr = 0;
			    out[r*columns+c].dc = -1;
			}
		}
}


/*****************************************************
Function name:   scale
Description:
    scales an image
Input parameters:
    in            --  pointer to input image
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void scale (struct edge *in, uint32 rows, uint32 columns)
{
    uint32 r;      		            		   // row index
    uint32 c;         	            		   // column index
	double scale_factor;                       // scaling factor
    double max_val;                            // hold the maximal value of the temporary images

	printf("scaling\n");

	// scale the temporary image
	max_val = 0;                     // initialization of maximum value

	for (r = 0; r < rows; r++)       // calcultion of the maximum value
        for (c = 0; c < columns; c++)
			if (in[r*columns+c].mag > max_val)
				max_val = in[r*columns+c].mag;

	scale_factor = (GV - 1)/max_val;

	for (r = 0; r < rows; r++)                       // scaling of the temporary image 2
        for (c = 0; c < columns; c++)
			in[r*columns+c].mag = in[r*columns+c].mag * scale_factor;
}


/*****************************************************
Function name:   suppress_nonmaxima
Description:
    suppresses the non maximum values
Input parameters:
    in            --  pointer to input image
    out           --  pointer to output image
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void suppress_nonmaxima (struct edge *in, double *out, uint32 rows, uint32 columns)
{
    uint32 r;       		          		   // row index
    uint32 c; 	        	           		   // column index
    int hr, hc;                                // holds the indexes used for the computation


	printf("suppress nonmaxima\n");

	// suppress nonmaxima
	for (r = 0; r < rows; r++)                      // indexes of the temporary images
        for (c = 0; c < columns; c++)
		{
			hc = c + in[r*columns+c].dc;    // indexes of the pixel in the direction
			hr = r + in[r*columns+c].dr;    // of the gradient

			if (hc < 0)
				hc = 0;
			if (hc > (columns -1))
				hc = columns - 1;
			if (hr < 0)
				hr = 0;
			if (hr > (rows -1))
				hr = rows - 1;

			if (in[r*columns+c].mag <= in[hr*columns+hc].mag)  // suppress nonmaxima
				out[r*columns+c] = 0;
			else
				out[r*columns+c] = in[r*columns+c].mag;

			hc = c - in[r*columns+c].dc;    // indexes of the pixel in the directon
			hr = r - in[r*columns+c].dr;    // of the gradient
		
			if (hc < 0)
				hc = 0;
			if (hc > (columns -1))
				hc = columns - 1;
			if (hr < 0)
				hr = 0;
			if (hr > (rows -1))
				hr = rows - 1;

			if (in[r*columns+c].mag <= in[hr*columns+hc].mag)  // suppress nonmaxima
				out[r*columns+c] = 0;
		}
}


/*****************************************************
Function name:   follow_edge
Description:
    follows an edge in an input image and prints
	this edge to an output image (has to be initialized)
Input parameters:
    er            --  row
	ec            --  column
	in            --  pointer to input image
	out           --  pointer to output image
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void follow_edge (uint32 er, uint32 ec, double *in, double *out, uint32 rows, uint32 columns,
				  float tlow)
{
	int hr, hc;								   // row and column index


	out[er*columns+ec] = GV-1;                 // Pixel gets highest value

	hr = er + 1;
	if (hr < rows)          // test all neigbors, if their values are over the threshold,
	{						// then follow this edge
		hc = ec + 1;
		if (hc < columns)
			if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
				follow_edge (hr, hc, in, out, rows, columns, tlow);

		hc = ec;
		if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
			follow_edge (hr, hc, in, out, rows, columns, tlow);

		hc = ec - 1;
		if (hc >= 0)
			if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
				follow_edge (hr, hc, in, out, rows, columns, tlow);
	}

	hr = er - 1;
	if (hr >= 0)
	{
		hc = ec + 1;
		if (hc < columns)
			if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
				follow_edge (hr, hc, in, out, rows, columns, tlow);

		hc = ec;
		if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
			follow_edge (hr, hc, in, out, rows, columns, tlow);

		hc = ec - 1;
		if (hc >= 0)
			if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
				follow_edge (hr, hc, in, out, rows, columns, tlow);
	}

	hr = er;
	hc = ec + 1;
	if (hc < columns)
		if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
			follow_edge (hr, hc, in, out, rows, columns, tlow);

	hc = ec - 1;
	if (hc >= 0)
		if ((in[hr*columns+hc] >= tlow) && (out[hr*columns+hc] != GV-1))
			follow_edge (hr, hc, in, out, rows, columns, tlow);
}


/*****************************************************
Function name:   hough_transform_circles
Description:
    applies the Hough transform for circles on an
	image and stores the result in a 3D array
	holding the	Hough space
Input parameters:
    in            --  pointer to input image
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void hough_transform_circles (double *in, int *hough_space, uint32 rows, uint32 columns)
{
    int r;       		            	   	   // row index
    int c; 	        	            	       // column index
	int rh;                                    // row index of the Hough space
	int ch;                                    // column index of the Hough space
	int i;                                     // diameter index of the Hough space
	int chmin, chmax;                          // minimum and maximum values of the column value
										       // in the Hough space
	const int d[4] = {d10, d1, d5, d25};       // array holding the diameters


	printf("performing the Hough transform\n");

    // initialize the Hough space
	for (rh = 0; rh < rows; rh++)
        for (ch = 0; ch < columns; ch++)
			for (i = 0; i < 4; i++)
				hough_space[i*rows*columns+rh*columns+ch] = 0;

    // perform the Hough transform
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			if (in[r*columns+c] == (GV-1))      // compute a circle for every pixel in the
				for (i = 0; i < 4; i++)         // edge image
				{
					chmin = c - d[i]/2;         // minimum column value of the circle
					chmax = c + d[i]/2;         // maximum column value of the circle

					if (chmin < 0)              // column values of the circle must
						chmin = 0;              // lie inside the Hough space
					if (chmax >= columns)
						chmax = columns - 1;

					for (ch = chmin; ch <= chmax; ch++)    // compute the circle
					{
						if (((ch-c) <= d[i]/2) && ((c-ch) <= d[i]/2))
						{
							rh = r + sqrt(pow(d[i]/2, 2) - pow(ch-c, 2));
							if (rh < rows)
								hough_space[i*rows*columns+rh*columns+ch]
								    = hough_space[i*rows*columns+rh*columns+ch] +1;

							rh = r - sqrt(pow(d[i]/2, 2) - pow(ch-c, 2));
							if (rh >= 0)
								hough_space[i*rows*columns+rh*columns+ch]
								    = hough_space[i*rows*columns+rh*columns+ch] +1;
						}
					}
				}
}


/*****************************************************
Function name:   hough_threshold
Describtion:
    sets the values less than or equal to
	the threshold to 0
Input parameters:
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void hough_threshold (int *hough_space, uint32 rows, uint32 columns)
{
	uint32 r;								// row index
	uint32 c;								// column index
	int i;									// index for the diameters of the circles
	const int t[4] = {t10, t1, t5, t25};    // holds the different thresholds for the
											// different diameters


	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			for (i = 0; i < 4; i++)
				if (hough_space[i*rows*columns+r*columns+c] <= t[i])
					hough_space[i*rows*columns+r*columns+c] = 0;
}


/*****************************************************
Function name:   regions_connect
Description:
    connect the regions in the Hough space by placing
	a square labeled -1 over every pixel in the
	Hough space being larger than 0, the result is
	stored in hough_connect
Input parameters:
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
    hough_connect --  pointer to the 3D array holding
	                  the conncted Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void regions_connect (int *hough_space, int *hough_connect, uint32 rows, uint32 columns)
{
	uint32 r;								// row index
	uint32 c;								// column index
	int i;									// index for the diameters of the circles
	uint32 rs, cs;							// row and column index of the squares
	uint32 rs_min, rs_max;					// minimum and maximum index (row)
	uint32 cs_min, cs_max;					// minimum and maximum index (column)
	const int size = size_connect;			// 1/2 of size of the squares


	printf("connect regions\n");

	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			for (i = 0; i < 4; i++)
				if (hough_space[i*rows*columns+r*columns+c] != 0)
				{
					rs_min = r - size;		// the square must lie inside the Hough space
					if (rs_min < 0)
						rs_min = 0;
					rs_max = r + size;
					if (rs_max >= rows)
						rs_max = rows - 1;
					cs_min = c - size;
					if (cs_min < 0)
						cs_min = 0;
					cs_max = c + size;
					if (cs_max >= columns)
						cs_max = columns - 1;

					for (rs = rs_min; rs <= rs_max; rs++)
						for (cs = cs_min; cs <= cs_max; cs++)
							hough_connect[i*rows*columns+rs*columns+cs] = -1;
				}
}


/*****************************************************
Function name:   hough_label
Description:
    label the connected regions in hough_connect,
	add the values in the Hough space belonging to
	these regions, compute the centroid of each region
	and write the value to the centroids of these
	regions
Input parameters:
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
    hough_connect --  pointer to the 3D array holding
	                  the conncted Hough space
	label         --  label to apply to the region
	r             --  row index of the current pixel
	c             --  column index of the current pixel
	i			  --  diameter index of the current
					  pixel
	rows          --  number of rows
	columns       --  number of columns
Returned values:
    *sum          --  sum of the values in the
	                  Hough space belonging to the
					  region
	*rc           --  row index of the centroid
	*cc           --  column index of teh centroid
	*n            --  number of pixels belonging to
	                  region
*****************************************************/
void hough_label (int *hough_space, int *hough_connect, int label, int *sum, double *rc,
				  double *cc, int *n, uint32 r, uint32 c, int i, uint32 rows, uint32 columns)
{
	uint32 rh, ch;							// row and column index


	hough_connect[i*rows*columns+r*columns+c] = label;		// apply the label
	*sum = *sum + hough_space[i*rows*columns+r*columns+c];	// compute the sum
	*rc = (*rc**n+r)/(*n+1);			// compute the row index of the centroid
	*cc = (*cc**n+c)/(*n+1);			// compute the column index of the centroid
	*n = *n+1;							// compute the number of pixels belonging to that region
    hough_space[i*rows*columns+r*columns+c] = 0; // set used pixel in the Hough space to 0
	
	// look at the 4-neighbors
	// if neighbor belongs to the region, then call function again
	rh = r;
	ch = c+1;
	if (ch >= columns)
		ch = columns-1;
	if (hough_connect[i*rows*columns+rh*columns+ch] == -1)
		hough_label (hough_space, hough_connect, label, &*sum, &*rc, &*cc, &*n, rh, ch, i,
					 rows, columns);

	ch = c-1;
	if (ch < 0)
		ch = 0;
	if (hough_connect[i*rows*columns+rh*columns+ch] == -1)
		hough_label (hough_space, hough_connect, label, &*sum, &*rc, &*cc, &*n, rh, ch, i,
					 rows, columns);

	ch = c;
	rh = rh+1;
	if (rh >= rows)
		rows = rows-1;
	if (hough_connect[i*rows*columns+rh*columns+ch] == -1)
		hough_label (hough_space, hough_connect, label, &*sum, &*rc, &*cc, &*n, rh, ch, i,
					 rows, columns);

	rh = rh-1;
	if (rh < 0)
		rows = 0;
	if (hough_connect[i*rows*columns+rh*columns+ch] == -1)
		hough_label (hough_space, hough_connect, label, &*sum, &*rc, &*cc, &*n, rh, ch, i,
					 rows, columns);
}


/*****************************************************
Function name:   search_maximum
Description:
    compares the value at one layer in the Hough space 
	with the values nearby in the other layers by
	comparing the value with the values in a square
	region around this pixel in the other layers,
	non maximum values are deleted
Input parameters:
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void search_maximum (int *hough_space, uint32 rows, uint32 columns)
{
	uint32 r;								// row index
	uint32 c;								// column index
	int i;									// index for the diameter (layer)
	int j;									// index for the diameters
	uint32 rs, cs;							// row and column indexes of the square region
	uint32 rs_min, rs_max;					// border of the square region
	uint32 cs_min, cs_max;					// border of the square region
	const int size = size_search;			// 1/2 of size of the square region


	printf("searching for maximum\n");

	for (i = 0; i < 4; i++)
		for (r = 0; r < rows; r++)
		    for (c = 0; c < columns; c++)
				if (hough_space[i*rows*columns+r*columns+c] != 0)
				{
					rs_min = r - size;		// square region must lie inside the Hough space
					if (rs_min < 0)
						rs_min = 0;
					rs_max = r + size;
					if (rs_max >= rows)
						rs_max = rows - 1;
					cs_min = c - size;
					if (cs_min < 0)
						cs_min = 0;
					cs_max = c + size;
					if (cs_max >= columns)
						cs_max = columns - 1;

					for (rs = rs_min; rs <= rs_max; rs++)         // search in the region
						for (cs = cs_min; cs <= cs_max; cs++)
							for (j = i+1; j < 4; j++)
								if (hough_space[i*rows*columns+r*columns+c]
									< hough_space[j*rows*columns+rs*columns+cs])
									hough_space[i*rows*columns+r*columns+c] = 0;
								else
									hough_space[j*rows*columns+rs*columns+cs] = 0;
				}
}


/*****************************************************
Function name:   count_coins
Description:
    counts the coins detected in the Hough space
Input parameters:
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    value of the detected coins
*****************************************************/
float count_coins (int *hough_space, uint32 rows, uint32 columns)
{
	uint32 r;										  // row index
	uint32 c;										  // column index
	int i;											  // index of the diameters
	int n;											  // number of coins
	float count;									  // detected value
	const float value[4] = {0.10, 0.01, 0.05, 0.25};  // holds the values of the different coins


	printf("counting coins\n");
	printf("\n");

	count = 0;

	// search in Hough space
	// if coin detected, add value to count
	for (i = 0; i < 4; i++)
	{
		n = 0;
		for (r = 0; r < rows; r++)
		    for (c = 0; c < columns; c++)
				if (hough_space[i*rows*columns+r*columns+c] != 0)
					n++;

		count += n * value[i];
		printf("number of coins %.2f $:     %d\n", value[i], n);
	}

	printf("\n");

	return count;
}


/*****************************************************
Function name:   generate_output_image
Description:
    counts the coins detected in the Hough space
Input parameters:
    in_image      --  pointer to the input image
	out_image     --  pointer to the output image
    hough_space   --  pointer to the 3D array holding
	                  the Hough space
	rows          --  number of rows
	columns       --  number of columns
Returned value:
    none
*****************************************************/
void generate_output_image(unsigned char *in_image, unsigned char *out_image, int *hough_space, uint32 rows, uint32 columns)
{
	uint32 r;									// row index
	uint32 c;									// column index
	int i;										// index of the diameters
	uint32 hr, hc;								// indexes for the template
	uint32 hr_min, hc_min;						// maximal and minimum values for the template
	const int number[4][24][26]					// holds the templates for the numbers
		= {{{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,0,0,0,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1},
			{1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1},
			{1,1,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}},

           {{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}},

           {{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}},

           {{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1},
			{1,1,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1},
			{1,1,1,1,0,0,0,1,1,0,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,0,0,1,1,1,1,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,0,0,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,0,0,0,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1},
			{1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1},
			{1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1},
			{1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1},
			{1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1},
			{1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,1,1,1,1},
			{1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,1,1},
			{1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
			{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
            {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}}};


	// write input image in output image		
	for (r = 0; r < rows; r++)
        for (c = 0; c < columns; c++)
			out_image[r*columns+c] = in_image[r*columns+c];

	// write value numbers belonging to detected coins in output image
	for (i = 0; i < 4; i++)
		for (r = 0; r < rows; r++)
		    for (c = 0; c < columns; c++)
			{
				if (hough_space[i*rows*columns+r*columns+c] != 0)
				{
					hr_min = r - 12;
					hc_min = c - 13;

					if (hr_min >= rows)
						hr_min = rows - 25;
					if (hr_min < 0)
						hr_min = 0;

					if (hc_min >= columns)
						hr_min = columns - 27;
					if (hc_min < 0)
						hc_min = 0;

					for (hr = hr_min; hr < hr_min+24; hr++)
						for (hc = hc_min; hc < hc_min+26; hc++)
							if((hr >= 0) && (hr < rows) && (hc >= 0) && (hc < columns))
								out_image[hr*columns+hc] = (GV-1) * number[i][hr-r+12][hc-c+13];
				}
			}
}