//
// Canny operator working on image sequence or camera picture
//
// Provides interactive control of three parameters:
// gaussian standard deviation for scale; 
// high, or starting, thresh - threshold on the filtered gradient value
// to start an edge line; 
// low/high - an edge line will be continued (threshold hysteresis) if 
// the gradient is above a secondary continuation threshold. This
// parameter sets the low threshold as a proportion of the high threshold.
//
// John Robinson 2004.
//
#include "picture.h"

double lowthresh;
int threshmaxima(int& outpel, const int& inpel) {
// Callback. On input, this picture is non-zero only at gradient maxima, 
// Other picture has gradient values.
// On output this picture has gradient values at maxima above lowthresh.
	if (outpel && (inpel > lowthresh))
		outpel = inpel;
	else
		outpel = 0;
	return 0;
	}

void followpoint(picture_of_int& out, picture_of_int& maxima,
		picture_of_int& mag, int lowthresh, int i, int j) {
// Recursive function to extend an edge from the current point so long
// as its value is above lowthresh.
	out[i][j] = 0;
	maxima[i][j] = 0;
	if (maxima[i-1][j-1])
		followpoint(out,maxima,mag,lowthresh,i-1,j-1);
	if (maxima[i-1][j])
		followpoint(out,maxima,mag,lowthresh,i-1,j);
	if (maxima[i-1][j+1])
		followpoint(out,maxima,mag,lowthresh,i-1,j+1);
	if (maxima[i][j-1])
		followpoint(out,maxima,mag,lowthresh,i,j-1);
	if (maxima[i][j+1])
		followpoint(out,maxima,mag,lowthresh,i,j+1);
	if (maxima[i+1][j-1])
		followpoint(out,maxima,mag,lowthresh,i+1,j-1);
	if (maxima[i+1][j])
		followpoint(out,maxima,mag,lowthresh,i+1,j);
	if (maxima[i+1][j+1])
		followpoint(out,maxima,mag,lowthresh,i+1,j+1);
	}


// Gradient callbacks
int gradv(int& pel, const int **buf)
	{
	pel = buf[0][1] - buf[0][-1];
	return 0;
	}
int gradh(int& pel, const int **buf) {
	pel = buf[1][0] - buf[-1][0];
	return 0;
	}


int main(int argc, char *argv[])
	{
	colour_picture incolour(argc, argv,480,640);	// Sequence processing
	double stddev = 1.4;			// Start variance of about 2
	parampic stddevparam(0.5, 4, 1.4, "Gaussian s.d.");
	double highthresh = 200.0;
	parampic highthreshparam(0, 1000, 200, "High (start) threshold");
	lowthresh = 100.0;
	double lowprop = 0.5;
	parampic lowpropparam(0, 1, 0.5, "Lowthresh/Highthresh");
	// The user doesn't control lowthresh directly, but as a proportion
	// of high thresh. This is so that only the highthresh slider
	// needs to be moved in most cases.
	int nr = incolour.nrows();
	int nc = incolour.ncols();
	picture_of_int in(nr,nc);
	picture_of_int filt(nr,nc);
	picture_of_int mag(nr,nc);
	picture_of_int maxima(nr,nc);
	picture_of_int h(nr,nc);
	picture_of_int v(nr,nc);
	picture_of_int out(nr,nc);
	incolour.show("input");
	out.show("output");
	while (!incolour.closerequested() && !out.closerequested()) {
	    incolour.next();
	    in = incolour;	// Convert to monochrome
	    in *= 10;	// So don't lose amplitude resolution in gaussian
	    incolour.reshow();
	    // Now get latest versions of parameters
	    stddevparam.X(stddev);
	    highthreshparam.X(highthresh);
	    lowpropparam.X(lowprop);
	    lowthresh = lowprop*highthresh;
	    // Make the gaussian
	    filt.gaussfilter(in,stddev);
	    // And the gradient components
	    h.neigh(filt,gradh);
	    v.neigh(filt,gradv);
	    int i;
	    maxima = 0;
	    // Now find the gradient magnitudes and test for maxima
	    for (i = 0; i < in.nrows(); i++)
	    for (int j = 0; j < in.ncols(); j++) {
		int hor = h[i][j];
		int ver = v[i][j];
		mag[i][j]=(int)(sqrt((double)(hor*hor+ver*ver)));
		// mag[i][j] = abs(hor) + abs(ver);
		if (ver == 0) {
			if ((hor>0)&&(hor>=h[i-1][j])&&(hor>h[i+1][j]))
				maxima[i][j] = 1;
			else if ((hor<0)&&(hor<=h[i-1][j])&&(hor<h[i+1][j]))
				maxima[i][j] = 1;
			}
		else if (hor == 0) {
			if ((ver>0)&&(ver>=v[i][j-1])&&(ver>v[i][j+1]))
				maxima[i][j] = 1;
			else if ((ver<0)&&(ver<=v[i][j-1])&&(ver<v[i][j+1]))
				maxima[i][j] = 1;
			}
		else {
			double offset;
			if (abs(hor) > abs(ver)) {
				offset = (double)ver/hor;
				if ((hor>0)&&(hor>=h[(double)i-1][j-offset])
				    &&(hor>h[(double)i+1][j+offset]))
					maxima[i][j] = 1;
				else if ((hor<0)&&
				    (hor<=h[(double)i-1][j-offset])
				    &&(hor<h[(double)i+1][j+offset]))
					maxima[i][j] = 1;
				}
			else {
				offset = (double)hor/ver;
				if ((ver>0)&&(ver>=v[i-offset][j-1])
				    &&(ver>v[i+offset][j+1]))
					maxima[i][j] = 1;
				else if ((ver<0)&&(ver<=v[i-offset][j-1])
				    &&(ver<v[i+offset][j+1]))
					maxima[i][j] = 1;
				}
			}
		}
	    // Remove all points below lowthresh from consideration
	    maxima.point(mag,threshmaxima);
	    out = 255;
	    // Start tracing lines from highthresh points, out so long as
	    // gradient value is above lowthresh.
	    for (i = 0; i < mag.nrows(); i++) {
		for (int j = 0; j < mag.ncols(); j++) {
			if (maxima[i][j] > highthresh)
				followpoint(out,maxima,mag,(int)lowthresh,i,j);
			}
		}
	    out.reshow();
	    }
	return 1;
	}