/* Edge detect an image */

/* mark@topic.demon.co.uk */
/* mpaddock@cix.compulink.co.uk */

#include <proto/exec.h>
#include <proto/dos.h>
#include <proto/utility.h>

#include <dos/dos.h>

#include <utility/tagitem.h>
#include <libraries/MPImage.h>
#include <pragmas/MPImage_pragmas.h>
#include <clib/MPImage_protos.h>
#include <dos.h>
#include <stdlib.h>
#include <dos/rdargs.h>
#include <stdio.h>
#include <mffp.h>

#define TEMPLATE "FROM1/K/A,FROM2/K,TO/K,PERCENT1=PERC1/N/K/A,PERCENT2=PERC2/N/K,MINPOINTS=MINP/N/K,MAXPOINTS=MAXP/N/K,MINDIFF/N/K,X/N/K/A,Y/K/N/A,PUBSCREEN/K,GRID/S,HOOK/K/N"

#define OPT_FILE1			0
#define OPT_FILE2			1
#define OPT_TO				2
#define OPT_PERC1			3
#define OPT_PERC2			4
#define OPT_MIN			5
#define OPT_MAX			6
#define OPT_MINDIFF		7
#define OPT_X				8
#define OPT_Y				9
#define OPT_PUB			10
#define OPT_GRID			11
#define OPT_HOOK			12

#define OPT_COUNT			13

static LONG opts[OPT_COUNT] = {
	0
};

//extern long __oslibversion=39;

extern long __stack = 16000;
struct Library *MPImageBase;

struct Point {
	int x,y;
	int newx1,newy1;
	int newx2,newy2;
	int diff1;
	int diff2;
};

const char Version[]="$VER: EdgePoints 4.4 (23.2.97)";

void SetMessage(char *message) {
	if (opts[OPT_HOOK]) {
		CallHookPkt(*((struct Hook **)opts[OPT_HOOK]),MPIP_MESSAGE,(APTR)message);
	}
}

void SetMax(int max) {
	if (opts[OPT_HOOK]) {
		CallHookPkt(*((struct Hook **)opts[OPT_HOOK]),MPIP_MAX,(APTR)max);
	}
}

void SetCur(int cur) {
	if (opts[OPT_HOOK]) {
		CallHookPkt(*((struct Hook **)opts[OPT_HOOK]),MPIP_CURR,(APTR)cur);
	}
}

int
cmp(const struct Point *a,const struct Point *b) {
	if ((a->diff1+a->diff2) < (b->diff1+b->diff2)) {
		return -1;
	}
	if ((a->diff1+a->diff2) == (b->diff1+b->diff2)) {
		return 0;
	}
	return 1;
}

struct MPImage *
EdgeDetect(char *file,int perc,int x,int y) {
	struct MPImage *MPi;
	int i,j,k;
	UWORD *pa,*op,*pw1,*pw2,*pw3;
	UBYTE *p1,*p2,*p3,*opb;
	UWORD mymax;
	int bar[255*4+1] = {0};
	int top,newtop;
	SetMessage("Loading Image");
	if (MPi=LoadMPImage(file,NULL,
								MPIF_GREY)) {
		if (((x == MPi->Width) && (y == MPi->Height)) || RescaleMPImage(MPi,x,y)) {
			if ((pa = calloc(sizeof(UWORD),x * y))) {
//				Printf("1.1\n");
				op = pa+1+x;
				p1 = MPi->Red;
				p3 = p1+(2*x);
				SetMessage("Processing Red 1");
				SetMax(y - 1);
				for (i = 2; i < y; ++i) {
					SetCur(i);
					for (j = 2; j < x; ++j) {
						k = (*p3 + (2*p3[1]) + p3[2] - *p1 - (2*p1[1]) - p1[2]);
						if (k<0) {
							*op++ = -k;
						}
						else {
							*op++ = k;
						}
						++p1;
						++p3;
					}
					op+=2;
					p1+=2;
					p3+=2;
				}
				if (!MPi->GreyScale) {
//					Printf("1.2\n");
					op = pa+1+x;
					p1 = MPi->Green;
					p3 = p1+(2*x);
					SetMessage("Processing Green 1");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (*p3 + (2*p3[1]) + p3[2] - *p1 - (2*p1[1]) - p1[2]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p3;
						}
						op+=2;
						p1+=2;
						p3+=2;
					}
//					Printf("1.3\n");
					op = pa+1+x;
					p1 = MPi->Blue;
					p3 = p1+(2*x);
					SetMessage("Processing Blue 1");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (*p3 + (2*p3[1]) + p3[2] - *p1 - (2*p1[1]) - p1[2]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p3;
						}
						op+=2;
						p1+=2;
						p3+=2;
					}
				}
//				Printf("2.1\n");
				op = pa+1+x;
				p1 = MPi->Red;
				p2 = p1+x;
				p3 = p2+x;
				SetMessage("Processing Red 2");
				for (i = 2; i < y; ++i) {
					SetCur(i);
					for (j = 2; j < x; ++j) {
						k = (*p2 + (2* *p3) + p3[1] - p1[1] - (2*p1[2]) - p2[2]);
						if (k < 0) {
							if (-k < *op) {
								op++;
							}
							else {
								*op++ = -k;
							}
						}
						else {
							if (k > *op) {
								*op++ = k;
							}
							else {
								++op;
							}
						}
						++p1;
						++p2;
						++p3;
					}
					op+=2;
					p1+=2;
					p2+=2;
					p3+=2;
				}
				if (!MPi->GreyScale) {
//					Printf("2.2\n");
					op = pa+1+x;
					p1 = MPi->Green;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Green 2");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (*p2 + (2* *p3) + p3[1] - p1[1] - (2*p1[2]) - p2[2]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
//					Printf("2.3\n");
					op = pa+1+x;
					p1 = MPi->Blue;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Blue 2");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (*p2 + (2* *p3) + p3[1] - p1[1] - (2*p1[2]) - p2[2]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
				}
//				Printf("3.1\n");
				op = pa+1+x;
				p1 = MPi->Red;
				p2 = p1+x;
				p3 = p2+x;
				SetMessage("Processing Red 3");
				for (i = 2; i < y; ++i) {
					SetCur(i);
					for (j = 2; j < x; ++j) {
						k = (p1[2] + (2*p2[2]) + p3[2] - *p1 - (2* *p2) - *p3);
						if (k < 0) {
							if (-k < *op) {
								op++;
							}
							else {
								*op++ = -k;
							}
						}
						else {
							if (k > *op) {
								*op++ = k;
							}
							else {
								++op;
							}
						}
						++p1;
						++p2;
						++p3;
					}
					op+=2;
					p1+=2;
					p2+=2;
					p3+=2;
				}
				if (!MPi->GreyScale) {
//					Printf("3.2\n");
					op = pa+1+x;
					p1 = MPi->Green;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Green 3");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (p1[2] + (2*p2[2]) + p3[2] - *p1 - (2* *p2) - *p3);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
//					Printf("3.3\n");
					op = pa+1+x;
					p1 = MPi->Blue;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Blue 3");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (p1[2] + (2*p2[2]) + p3[2] - *p1 - (2* *p2) - *p3);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
				}
//				Printf("4.1\n");
				op = pa+1+x;
				p1 = MPi->Red;
				p2 = p1+x;
				p3 = p2+x;
				SetMessage("Processing Red 4");
				for (i = 2; i < y; ++i) {
					SetCur(i);
					for (j = 2; j < x; ++j) {
						k = (p2[2] + (2*p3[2]) + p3[1] - *p2 - (2* *p1) - p1[1]);
						if (k < 0) {
							if (-k < *op) {
								op++;
							}
							else {
								*op++ = -k;
							}
						}
						else {
							if (k > *op) {
								*op++ = k;
							}
							else {
								++op;
							}
						}
						++p1;
						++p2;
						++p3;
					}
					op+=2;
					p1+=2;
					p2+=2;
					p3+=2;
				}
				if (!MPi->GreyScale) {
//					Printf("4.2\n");
					op = pa+1+x;
					p1 = MPi->Green;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Green 4");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (p2[2] + (2*p3[2]) + p3[1] - *p2 - (2* *p1) - p1[1]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
//					Printf("4.3\n");
					op = pa+1+x;
					p1 = MPi->Blue;
					p2 = p1+x;
					p3 = p2+x;
					SetMessage("Processing Blue 4");
					for (i = 2; i < y; ++i) {
						SetCur(i);
						for (j = 2; j < x; ++j) {
							k = (p2[2] + (2*p3[2]) + p3[1] - *p2 - (2* *p1) - p1[1]);
							if (k < 0) {
								if (-k < *op) {
									op++;
								}
								else {
									*op++ = -k;
								}
							}
							else {
								if (k > *op) {
									*op++ = k;
								}
								else {
									++op;
								}
							}
							++p1;
							++p2;
							++p3;
						}
						op+=2;
						p1+=2;
						p2+=2;
						p3+=2;
					}
				}
				SetMessage("Computing Histogram");
				pw1 = pa;
				for (i=0; i<y; ++i) {
					for (j=0; j<x; ++j) {
						++bar[*pw1++];
					}
				}
				newtop = (x*y * perc)/100;
				top = 0;
				for (i = (255*4); (i > 0) && (top < newtop); --i) {
					top += bar[i];
				}
				if (i) {
					mymax = i-1;
				}
				else {
					mymax = 0;
				}
				op = pa;
				for (i=0; i<y; ++i) {
					for (j=0; j<x; ++j) {
						if (*op > mymax) {
							*op++ = 255;
						}
						else {
							*op++ = 0;
						}
					}
				}
				opb = MPi->Red;
				pw1 = pa;
				pw2 = pw1+x;
				pw3 = pw2+x;
				for (j = 0; j<x; ++j) {
					*opb++=255;
				}
				SetMessage("Computing Edges");
				for (i=2; i<y; ++i) {
					SetCur(i);
					*opb++=255;
					for (j=2; j<x; ++j) {
						if ((*pw1 && pw1[2] && *pw3 && pw3[2]) ||
							 (pw1[1] && *pw2 && pw2[2] && pw3[1])) {
							*opb++ = 0;
						}
						else {
							if (pw2[1]) {
								if ((*pw1 + pw1[1] + pw1[2] + *pw2 + pw2[2] + *pw3 + pw3[1] + pw3[2]) < 512) {
									*opb++ = 0;
								}
								else {
									*opb++ = 255;
								}
							}
							else {
								*opb++ = 0;
							}
						}
						++pw1;
						++pw2;
						++pw3;
					}
					*opb++=255;
					pw1+=2;
					pw2+=2;
					pw3+=2;
				}
				for (j = 0; j<x; ++j) {
					*opb++=255;
				}
				free(pa);
			}
			else {
				Printf("Error allocating work buffer");
				FreeMPImage(MPi);
				MPi = NULL;
			}
		}
		else {
			Printf(MPImageErrorMessage());
			Printf("\n");
			FreeMPImage(MPi);
			MPi = NULL;
		}
	}
	else {
		Printf(MPImageErrorMessage());
		Printf("\n");
	}
	return MPi;
}

int
main(int argc,char **argv) {
	struct RDArgs *rdargs;
	UBYTE *op;
	FILE *fp;
	int resx = RETURN_OK;
	int i, j, k, l, m;
	int x, y;
	int pixel=1;
	struct MPImage *MPi1,*MPi2;
	struct Point *Points;
	int maxp,minp;
	int randx,randy;
	BOOL skip;
	int numx=0,numy=0;

	if (!(rdargs = ReadArgs((char *)TEMPLATE, opts, NULL))) {
		PrintFault(IoErr(), NULL);
		return RETURN_ERROR;
	}
	if (MPImageBase = OpenLibrary("MPImage.library",6)) {
		SetMPImageScreen((char *)opts[OPT_PUB],opts[OPT_HOOK]?0:MPIF_PROGRESS);
		if (opts[OPT_HOOK]) {
			MPProgressHook(*((struct Hook **)opts[OPT_HOOK]));
		}
		x = *((ULONG *)opts[OPT_X]);
		y = *((ULONG *)opts[OPT_Y]);
		if (!opts[OPT_FILE2]) {
			SetMessage("Edge Detecting");
			if (MPi1 = EdgeDetect((char *)opts[OPT_FILE1],*((ULONG *)opts[OPT_PERC1]),x,y)) {
				SetMessage("Saving Image");
				if (!SaveMPImage((char *)opts[OPT_TO],MPi1->Red,MPi1->Red,MPi1->Red,MPi1->Width,MPi1->Height,
								MPIS_FORMAT,	MPI_COLOUR,
								MPIS_COLOURS,	2,
								MPIS_12BIT,		TRUE,
								TAG_END)) {
					resx = RETURN_FAIL;
					Printf(MPImageErrorMessage());
					Printf("\n");
				}					
				FreeMPImage(MPi1);
			}
			else {
				resx = RETURN_FAIL;
				Printf(MPImageErrorMessage());
				Printf("\n");
			}
			CloseLibrary(MPImageBase);
			return resx;
		}
		if (!opts[OPT_TO] || !opts[OPT_PERC2] || !opts[OPT_MIN] || !opts[OPT_MAX] ||
			 !opts[OPT_MINDIFF]) {
			Printf("If FROM2 supplied then all parameters must be supplied\n");
			CloseLibrary(MPImageBase);
			return RETURN_FAIL;
		}
		SetMessage("Edge Detecting 1");
		if (MPi1 = EdgeDetect((char *)opts[OPT_FILE1],*((ULONG *)opts[OPT_PERC1]),x,y)) {
			SetMessage("Edge Detecting 2");
			if (MPi2 = EdgeDetect((char *)opts[OPT_FILE2],*((ULONG *)opts[OPT_PERC2]),x,y)) {
				maxp = *((ULONG *)opts[OPT_MAX]);
				minp = *((ULONG *)opts[OPT_MIN]);
				if (opts[OPT_GRID]) {
					SetMessage("Computing Grid");
					pixel = (x + y)/sqrt((double)maxp);
					numx = x/pixel;
					numy = y/pixel;
					maxp = (numx + 1) * (numy + 1);
				}
				if (Points = calloc(maxp,sizeof(struct Point))) {
					Points[0].x = 0;
					Points[0].y = 0;
					Points[0].newx1 = 0;
					Points[0].newy1 = 0;
					Points[0].newx2 = 0;
					Points[0].newy2 = 0;
					Points[1].x = x-1;
					Points[1].y = 0;
					Points[1].newx1 = x-1;
					Points[1].newy1 = 0;
					Points[1].newx2 = x-1;
					Points[1].newy2 = 0;
					Points[2].x = 0;
					Points[2].y = y-1;
					Points[2].newx1 = 0;
					Points[2].newy1 = y-1;
					Points[2].newx2 = 0;
					Points[2].newy2 = y-1;
					Points[3].x = x-1;
					Points[3].y = y-1;
					Points[3].newx1 = x-1;
					Points[3].newy1 = y-1;
					Points[3].newx2 = x-1;
					Points[3].newy2 = y-1;
					if (opts[OPT_GRID]) {
						k = 4;
						for (j = 1; j<numy; ++j) {
							Points[k].x = 0;
							Points[k].y = j * pixel;
							Points[k].diff1 = 0x7fffffff;
							Points[k].diff2 = 0x7fffffff;
							++k;
							Points[k].x = x-1;
							Points[k].y = j * pixel;
							Points[k].diff1 = 0x7fffffff;
							Points[k].diff2 = 0x7fffffff;
							++k;
						}
						for (i = 1; i<numx; ++i) {
							for (j = 0; j<numy; ++j) {
								Points[k].x = i * pixel;
								Points[k].y = j * pixel;
								Points[k].diff1 = 0x7fffffff;
								Points[k].diff2 = 0x7fffffff;
								++k;
							}
							Points[k].x = i * pixel;
							Points[k].y = y-1;
							Points[k].diff1 = 0x7fffffff;
							Points[k].diff2 = 0x7fffffff;
							++k;
						}
					}
					else {
						SetMessage("Computing Random Points");
						randy = RAND_MAX/MPi1->Height;
						randx = RAND_MAX/MPi1->Width;
						for (k=4; k<maxp; ++k) {
							Points[k].x = rand()/randx;
							Points[k].y = rand()/randy;
							Points[k].diff1 = 0x7fffffff;
							Points[k].diff2 = 0x7fffffff;
						}
					}
					op = MPi1->Red;
					SetMessage("Computing Image 1");
					SetMax(y - 1);
					for (i=0; i<y; ++i) {
						SetCur(i);
//						Printf("%ld\n",i);
						for (j=0; j<x; ++j) {
							if (*op++) {
								for (k=4; k<maxp; ++k) {
									int diff1;
									diff1 = abs(i-Points[k].y);
									if (diff1 < Points[k].diff1) {
										int diff2;
										diff2 = abs(j-Points[k].x);
										if (diff2 < Points[k].diff1) {
											diff1 *= diff1;
											if (diff1 < Points[k].diff1) {
												diff2 *= diff2;
												if (diff2 < Points[k].diff1) {
													diff1 += diff2;
													if (diff1 < Points[k].diff1) {
														Points[k].diff1 = diff1;
														Points[k].newx1 = j;
														Points[k].newy1 = i;
													}
												}
											}
										}
									}
								}
							}
						}
					}
					op = MPi2->Red;
					SetMessage("Computing Image 2");
					for (i=0; i<y; ++i) {
						SetCur(i);
//						Printf("%ld\n",i);
						for (j=0; j<x; ++j) {
							if (*op++) {
								for (k=4; k<maxp; ++k) {
									int diff1;
									diff1 = abs(i-Points[k].y);
									if (diff1 < Points[k].diff2) {
										int diff2;
										diff2 = abs(j-Points[k].x);
										if (diff2 < Points[k].diff2) {
											diff1 *= diff1;
											if (diff1 < Points[k].diff2) {
												diff2 *= diff2;
												if (diff2 < Points[k].diff2) {
													diff1 += diff2;
													if (diff1 < Points[k].diff2) {
														Points[k].diff2 = diff1;
														Points[k].newx2 = j;
														Points[k].newy2 = i;
													}
												}
											}
										}
									}
								}
							}
						}
					}
					SetMessage("Sorting Points");
					qsort(&(Points[4]),maxp-4,sizeof(struct Point),cmp);
//					Printf("Sorted\n");
					SetMessage("Writing Points");
					SetMax(minp - 1);
					if (fp = fopen((char *)opts[OPT_TO],"w")) {
						int mindiff;
						mindiff = *((ULONG *)opts[OPT_MINDIFF]);
						fprintf(fp,"/* EdgePoints AddPoints file */\n");
						for (k = 0, m = 0; (k<maxp) && (m<minp); ++k) {
							SetCur(m);
							skip = FALSE;
							for (l = 0; (l < k) && !skip ; ++l) {
								if (((abs(Points[l].newx1 - Points[k].newx1) +
									   abs(Points[l].newy1 - Points[k].newy1)) < mindiff) ||
									 ((abs(Points[l].newx2 - Points[k].newx2) +
									   abs(Points[l].newy2 - Points[k].newy2)) < mindiff)) {
									skip = TRUE;
								}
							}
							if (!skip) {
								fprintf(fp,"'addpoint x1=%d y1=%d x2=%d y2=%d'\n",Points[k].newx1,Points[k].newy1,Points[k].newx2,Points[k].newy2);
								++m;
							}
						}
						SetCur(minp - 1);
						fclose(fp);
					}
					else {
						Printf("Failure opening '%s'\n",opts[OPT_TO]);
					}
				}
				else {
					Printf("Failure allocating Points");
				}
				FreeMPImage(MPi2);
			}
			else {
				resx = RETURN_FAIL;
			}
			FreeMPImage(MPi1);
		}
		else {
			resx = RETURN_FAIL;
		}
		CloseLibrary(MPImageBase);
	}
	else {
		Printf((char *)"Error opening MPImage.library\n");
		resx = RETURN_FAIL;
	}
	FreeArgs(rdargs);
	return resx;
}