/*
	File: mh.c

	Copyright 1999 by Ralph Pugmire

	Marr Hildreth Edge detection filter
*/

#include <SetUpA4.h>
#include "FilterInterface.h"
#include "Utilities.h"
#include "mh.h"
#include "math.h"
#include <stdlib.h>

#define Length(string) (*(unsigned char *)(string))

typedef struct TParameters{
	int		wsize;
	Boolean	sbsflag,mtflag;
	double 	sigma, zband, total;
} TParameters, *PParameters, **HParameters;

short 			gResult;
FilterRecordPtr gStuff;
int				mx,my,wsize,
			    midx,midy,	/* offset to middle of convolving mask */
			    xmax,ymax,
			    nxmax,nymax;

Boolean 		sbsflag,mtflag;
double 			sigma,
				zband, 
				total;
imgdr 			c_img, m_img;
imgd 			i_img,o_img;	/* pointer to my style image descriptor */

/*****************************************************************************/
/* My routines*/
double calctotal(Boolean mtflag){
int    
	x,y;
double
    zoffset,dtotal,
    dsigma,	    /* sigma as double */
    rs,		    /* distance from centre of mask*/
    pis4,	    /* 1/pi s^4 */
    ro2s,	    /* r^2 / 2sigma^2 */
    midx,midy;	    /* middle of mask */
    
    midx = (m_img.x - 1)/2.0;
    midy = (m_img.y - 1)/2.0;
    dsigma = sigma;
    pis4 = 1.0 / 3.14159 / pow(dsigma,4.0);
    
    for (y=dtotal=0; y<m_img.y; y++){
		for (x=0; x<m_img.x; x++){
		    rs = (midx-x)*(midx-x) + (midy-y)*(midy-y);
		    ro2s = rs / 2 / dsigma / dsigma;
		    m_img.e[y][x] = 10 * pis4 * (1.0 - ro2s) * exp(-ro2s);
		    dtotal += m_img.e[y][x] ;
		}
	}
	if(mtflag){
		zoffset = (dtotal) / wsize / wsize;
		for (y=dtotal=0; y<m_img.y; y++){
		    for (x=0; x<m_img.x; x++){
				m_img.e[y][x] -= zoffset;
				dtotal += m_img.e[y][x] ;
		    }
		}
	}         
	return dtotal;
}

/*****************************************************************************/

/* Calls the host's TestAbort function */

Boolean TestAbort (void){
	return CallPascalB (gStuff->abortProc);
	}

/*****************************************************************************/

/* Calls the host's UpdateProgress procedure */

void UpdateProgress (long done, long total){
	CallPascal (done, total, gStuff->progressProc);
	}

/*****************************************************************************/

/* Centers a dialog template 1/3 of the way down on the main screen. */

#define menuHeight 20
void CenterDialog (DialogTHndl dt){
	Rect r;
	short width;
	short height;

	width  = screenBits.bounds.right;
	height = screenBits.bounds.bottom;

	r = (**dt).boundsRect;
	OffsetRect (&r, -r.left, -r.top);
	OffsetRect (&r, (width - r.right) / 2,
					(height - r.bottom - menuHeight) / 3 + menuHeight);
	(**dt).boundsRect = r;
	}

#undef menuHeight

/*****************************************************************************/

/* Displays the about dialog box for the plug-in module. */

void DoAbout (void){
	short item;
	DialogPtr dp;
	DialogTHndl dt;

	dt = (DialogTHndl) GetResource ('DLOG', AboutdialogID);
	HNoPurge ((Handle) dt);
	CenterDialog (dt);

	dp = GetNewDialog (AboutdialogID, nil, (WindowPtr) -1);

	ModalDialog (nil, &item);

	DisposDialog (dp);
	HPurge ((Handle) dt);
	}

#undef dialogID

/*****************************************************************************/

/* UserItem to outline the OK button in a dialog box. */

pascal void OutlineOK (DialogPtr dp, short item){
	Rect r;
	Handle h;
	short itemType;
	
	SetUpA4 ();

	item = OK;

	GetDItem (dp, item, &itemType, &h, &r);

	PenNormal ();
	PenSize (3, 3);
	
	InsetRect (&r, -4, -4);
	FrameRoundRect (&r, 16, 16);

	PenNormal ();
	
	RestoreA4 ();
	}

/*****************************************************************************/

/* Asks the user for the plug-in filter module's parameters. Note that
   the image size information is not yet defined at this point. Also, do
   not assume that the calling program will call this routine every time the
   filter is run (it may save the data held by the parameters handle in
   a macro file). */


void DoParameters (void){
	Rect r;
	short j;
	Str255 s;
	unsigned char *p, *pp;
	Handle h;
	short item;
	DialogPtr dp;
	short itemType;
	DialogTHndl dt;
	int x,y;
	float prevtotal;
	Boolean afterpeak;
	
	if (!gStuff->parameters)
		{
		gStuff->parameters = NewHandle ((long) sizeof (TParameters));

		if (!gStuff->parameters)
			{
			gResult = memFullErr;
			return;
			}

		((PParameters) *gStuff->parameters)->sigma = 1.2;
		((PParameters) *gStuff->parameters)->zband = 0.01;
		((PParameters) *gStuff->parameters)->total = 0.0;
		((PParameters) *gStuff->parameters)->wsize = 11;
		((PParameters) *gStuff->parameters)->mtflag = false;
		((PParameters) *gStuff->parameters)->sbsflag = true;
		}

	dt = (DialogTHndl) GetResource ('DLOG', ParamsdialogID);
	HNoPurge ((Handle) dt);
	CenterDialog (dt);

	dp = GetNewDialog (ParamsdialogID, nil, (WindowPtr) -1);
	
	RememberA4 ();

	GetDItem (dp, hookItem, &itemType, &h				  , &r);
	SetDItem (dp, hookItem,  itemType, (Handle) &OutlineOK, &r);

	sigma = ((PParameters) *gStuff->parameters)->sigma;
	zband = ((PParameters) *gStuff->parameters)->zband;
	wsize = ((PParameters) *gStuff->parameters)->wsize;
	mtflag = ((PParameters) *gStuff->parameters)->mtflag;
	sbsflag = ((PParameters) *gStuff->parameters)->sbsflag;
	gResult = dimimgr(&m_img,wsize,wsize); /* create and init mask image */
	if (gResult != noErr) return;
	total = calctotal(mtflag);
	freeimgr(&m_img);
	
	SetDReal(dp,sigmaItem,sigma,2);
	SetDReal(dp,zbandItem,zband,3);
	SetDReal(dp,totalItem,total,6);
	SetDNum(dp,wsizeItem,wsize);
	SetDialogItem(dp,mtflagItem,mtflag);
	SetDialogItem(dp,sbsflagItem,sbsflag);
	SelIText (dp, sigmaItem, 0, 32767);
	do{

		ModalDialog (nil, &item);

		switch (item){
		case wsizeItem:
			/* Beep and remove anything other than 0-9 */
			GetDString(dp,wsizeItem,s);
			for (p=s+1; p<=s+s[0]; p++){
				if (*p>'9' || *p<'0'){
					SysBeep(1);
					s[0]--;
					for (pp=p; pp<=s+s[0]; pp++)
						*pp = *(pp+1);
				}
				SetDString(dp,wsizeItem,s);
			}
			break;
		case calcwsizeItem:
			//Find odd value for which total is less than 0.01 after reached a max
			sigma = GetDReal(dp,sigmaItem);
			afterpeak = false;
			for (wsize=3; wsize<25; wsize += 2){
				SetDNum(dp,wsizeItem,wsize);
				gResult = dimimgr(&m_img,wsize,wsize); /* create and init mask image */
				if (gResult != noErr) return;
				total = calctotal(false);
				SetDReal(dp,totalItem,total,6);
				freeimgr(&m_img);
				if (total<prevtotal)  afterpeak=true;
				if ((total<0.01) && afterpeak) break;
				prevtotal = total;
			}
			gResult = dimimgr(&m_img,wsize,wsize); /* create and init mask image */
			if (gResult != noErr) return;
			total = calctotal(mtflag);
			SetDReal(dp,totalItem,total,6);
			freeimgr(&m_img);
			break;
			
		case sbsflagItem:
			sbsflag = !sbsflag;
			SetDialogItem(dp,sbsflagItem,sbsflag);
			break;
			
		case mtflagItem:
			mtflag = !mtflag;
			SetDialogItem(dp,mtflagItem,mtflag);
			if (!mtflag)
				SelIText (dp, zbandItem, 0, 32767);
			else
				SelIText (dp, sigmaItem, 0, 32767);
			// Then fall into calcItem	

		case calcItem:
			sigma = GetDReal(dp,sigmaItem);
			wsize = GetDNum(dp,wsizeItem);
			gResult = dimimgr(&m_img,wsize,wsize); /* create and init mask image */
			if (gResult != noErr) return;
			total = calctotal(mtflag);
			SetDReal(dp,totalItem,total,6);
			freeimgr(&m_img);
			break;
		}
	}while (item != OK && item != Cancel);

	if (item == OK){
		sigma = GetDReal(dp,sigmaItem);
		wsize = GetDNum(dp,wsizeItem);
		zband = GetDReal(dp,zbandItem);
		((PParameters) *gStuff->parameters)->sigma = sigma;
		((PParameters) *gStuff->parameters)->zband = zband;
		((PParameters) *gStuff->parameters)->total = total;
		((PParameters) *gStuff->parameters)->wsize = wsize;
		((PParameters) *gStuff->parameters)->mtflag = mtflag;
		((PParameters) *gStuff->parameters)->sbsflag = sbsflag;
	} 

	DisposDialog (dp);
	HPurge ((Handle) dt);

	if (item == Cancel){
		gResult = 1;
		return;
	}

}
#undef dialogID
#undef hookItem

/*****************************************************************************/

/* Prepare to filter an image.	If the plug-in filter needs a large amount
   of buffer memory, this routine should set the bufferSpace field to the
   number of bytes required. */

void DoPrepare (void){
	typedef short double *psdouble;
	double zoffset;
	int i,x,y;
	unsigned char 
	 	*ptr;

	sigma = ((PParameters) *gStuff->parameters)->sigma;
	zband = ((PParameters) *gStuff->parameters)->zband;
	wsize = ((PParameters) *gStuff->parameters)->wsize;
	mtflag = ((PParameters) *gStuff->parameters)->mtflag;
	sbsflag = ((PParameters) *gStuff->parameters)->sbsflag;

	my = gStuff->filterRect.bottom - gStuff->filterRect.top + 1;
	mx = gStuff->filterRect.right - gStuff->filterRect.left + 1;
    xmax = mx - wsize + 1;
    ymax = my - wsize + 1;
    midx = (wsize - 1)/2.0;
    midy = (wsize - 1)/2.0;
    nxmax = mx - midx -1;
    nymax = my - midy -1;

/* tell them how much space we want */
	gStuff->bufferSpace = 	(long)wsize * wsize * sizeof(short double) 
							+ (long)wsize * sizeof(psdouble)
							+(long)mx * my * sizeof(short double) 
							+ (long)my * sizeof(psdouble);

}

/*****************************************************************************/

/* Requests pointer to the first part of the image to be filtered. */
/* Lets try asking for all the image in one go */
void DoStart (void){
	
    gResult = dimimgr(&c_img,my,mx); /* create and init real conv img*/
	if (gResult != noErr) return;

	gResult = dimimgr(&m_img,wsize,wsize); /* create and init kernel image */
	if (gResult != noErr){
		freeimgr(&c_img);
 		return;
 	}
	
	total = calctotal(mtflag);

	gStuff->inRect = gStuff->filterRect;

	gStuff->inLoPlane = 0;
	
	switch (gStuff->imageMode){
		case filterModeGrayScale:
			gStuff->inHiPlane = 0;
			break;
		default:
			gResult = filterBadMode;
	}

	gStuff->outRect = gStuff->inRect;
	gStuff->outLoPlane = gStuff->inLoPlane;
	gStuff->outHiPlane = gStuff->inHiPlane;
}

/*****************************************************************************/

/* Filters the area and requests the next area. */

void DoContinue (void){
/* convolve with input image */
#define max(a,b) ((a>b) ? (a) : (b))
#define min(a,b) ((a<b) ? (a) : (b))
    int 
		i,d,x,y,sj,j,jj;
    short double
		*re,smin,smax,cmin,cmax;
    short double
		c,mzband;
	unsigned char
		*tpe,*pe,*ptr;

	/*convert the input and output roi's to my style images*/
	i_img.y = my;
	i_img.x = mx;
	if (!(i_img.e = (ppc)NewPtr(sizeof(ptr) * my))){
		gResult = memFullErr; return;}
	ptr=(unsigned char *) gStuff->inData;
	for (i=0; i<my; i++){
		i_img.e[i] = ptr;
		ptr += gStuff->inRowBytes;
	}
	o_img.y = my;
	o_img.x = mx;
	if (!(o_img.e = (unsigned char **)NewPtr(sizeof(ptr) * my))){
		gResult = memFullErr; return;}
	ptr=(unsigned char *) gStuff->outData;
	for (i=0; i<my; i++){
		o_img.e[i] = ptr;
		ptr += gStuff->outRowBytes;
	}
	mzband = zband * 255.0;
    for (y=0;y < ymax; y++){
		UpdateProgress ((long) y,(long) ymax+nymax/10.0);
		if (TestAbort ()){
			gResult = 1;
			return;
		}
		for (x=0; x < xmax; x++){
		    sj = y + wsize;
		    for (jj=0,j=y,c=0; j < sj ; j++,jj++){
				re = m_img.e[jj];
				pe = &i_img.e[j][x];
				tpe = &i_img.e[j][x + wsize];
				for (; pe < tpe; c += *re++ * *pe++);
		    }	
		    if (c>0 && c<=mzband)
			    c=0;
		    if (c<0 && c>=mzband)
			    c=0;
		    c_img.e[y+midy][x+midx] = c;
		}
	}

/*Detect zero crossings */
    for (y=midy+1,smax=0; y < nymax-1; y++){
		UpdateProgress ((long) ymax+y/10,(long) ymax+nymax/10.0);
		for (x=midx+1; x < nxmax-1; x++){
		    short double a,b,c,s;
		    c = c_img.e[y][x];
		    if (c<0) c = -c;
		    s = 0;
		    a = c_img.e[y-1][x-1];
		    b = c_img.e[y+1][x+1];
		    if (a>0 && b<0) s = max(s,a-b);
		    else if (a<0 && b>0) s = max(s,b-a);
		    a = c_img.e[y-1][x];
		    b = c_img.e[y+1][x];
		    if (a>0 && b<0) s = max(s,a-b);
		    else if (a<0 && b>0) s = max(s,b-a);
		    a = c_img.e[y-1][x+1];
		    b = c_img.e[y+1][x-1];
		    if (a>0 && b<0) s = max(s,a-b);
		    else if (a<0 && b>0) s = max(s,b-a);
		    a = c_img.e[y][x-1];
		    b = c_img.e[y][x+1];
		    if (a>0 && b<0) s = max(s,a-b);
		    else if (a<0 && b>0) s = max(s,b-a);
	        if (s!=0){
	        	if (c>s) s=0;
	        	else if (sbsflag) s = s - c; else s = 1 - c/s;
			}
			smax = max(s,smax);
			c_img.e[y-2][x-2] = s;
		}
    }

/* put it in the ouput image and scale to 0-255 */
    /*printf("Max edge strength = %f\n",smax);*/
    for (y=midy+1; y < nymax-1; y++)
		for (x=midx+1; x < nxmax-1; x++)
		    o_img.e[y][x] = 255.0 - (c_img.e[y-2][x-2]*255.0/smax);
		
	SetRect (&gStuff->inRect, 0, 0, 0, 0); /* let em know we finished */
	gStuff->outRect = gStuff->inRect;
}

/*****************************************************************************/

/* This routine will always be called if DoStart does not return an error
   (even if DoContinue returns an error or the user aborts the operation).
   This allows the module to perform any needed cleanup.  None is required
   in this example. */

void DoFinish (void){
    DisposPtr((Ptr)i_img.e);   /*Free dynamic space used by my img desc array of pts*/
    DisposPtr((Ptr)o_img.e);   /*Free dynamic space used by my img desc array of pts*/
	freeimgr(&m_img);
	freeimgr(&c_img);

}

/*****************************************************************************/

/* Main dispatching routine.  Initializes and sets up the global variables,
   and performs the operation specified by the selector. */

pascal void main (short selector, FilterRecordPtr stuff, long *data, short *result){

	/* Allow access to global variables */

	RememberA0 ();
	SetUpA4 ();

	/* Copy the current quickdraw globals into the plug-in local copy */

	asm
		{
			MOVE.L	(A5),A0 		; Get address of real quickdraw globals
			SUB.L	#126,A0 		; Move to start
			LEA 	randSeed,A1 	; Get address of local copy
			MOVE.W	#64,D0			; Globals are 65 words long
		@1	MOVE.W	(A0)+,(A1)+ 	; Copy a word
			DBF 	D0,@1			; Move to next word
		}

	/* Perform the requested operation */

	gStuff	= stuff;
	gResult = noErr;

	switch (selector)
		{

		case filterSelectorAbout:
			DoAbout ();
			break;

		case filterSelectorParameters:
			DoParameters ();
			break;

		case filterSelectorPrepare:
			DoPrepare ();
			break;

		case filterSelectorStart:
			DoStart ();
			break;

		case filterSelectorContinue:
			DoContinue ();
			break;

		case filterSelectorFinish:
			DoFinish ();
			break;

		default:
			gResult = filterBadParameters;

		}

	*result = gResult;

	/* Restore the application's A4 register */

	RestoreA4 ();

	}