//	Butterworth
//
//	last modified 10-Sept-2002 RTH
//		
//		Interactively applies a Butterworth filter to an image. 
//		This filter works as a high or low pass filter.
//		Also works on selections and keeps filter data with image.  
//		
//		For DM 3.3.1 or higher! 
//		
//	The Butterworth filter is defined as: 
//		H = 1/[1 + C*(R/R0)**2*n]   => 
//		R = R0 ((1-H)/H)**1/(2*n) 
//			as in "The Image processing handbook", by J. Russ, p198 
//			note:  here C==1 is used (=> cutoff R0 is at 50% value) 
//		
//		
//		
	number size, sizeX, sizeY, posX, posY, key, key2, orderOfFilter, R0, sizeX2, sizeY2 
	number xcv, xcu, rad2, zoom, ii, mini, top, left, bottom, right, annot, annot1, annot2
	number displayzoom, calibrate, scale, scaleY, imageDataTypeIs, version, jj 
	number WindowSizeX, WindowSizeY 
	string nr, pixelSizeUnit, pixelSizeUnitF, filtertype, Neg_Val 
	string platform, ApplyToFullImage 
	number uparrow, downarrow, leftarrow, rightarrow, tops, lefts, bottoms, rights 
	number sizeXFull, sizeYFull, DMversion 
	ComplexImage frontF, frontF1, frontF2, frontC 
	Image front, display, filter, test, final, finalFull 
	CloseProgressWindow() 
		
//		
//		Solve version issues: 
//		
	GetPersistentNumberNote("Private:Configuration:ApplicationVersion", DMversion) 
	if( DMVersion <= 330 ) 
		Throw("Needs DigitalMicrograph version 3.3.1 or higher") 
	uparrow = 30;	downarrow = 31; 	leftarrow = 28;		rightarrow = 29;
		
//		
//		Put up info window 
//		
	OkDialog("Press space bar to end script at any time." + \
			 "\n\n(For additional info see \nResults- and Progress- windows)")
		
//		
//		Check on front Image:  only for real/integer images
//		
	front := GetFrontImage() 
	try
		mini = min(front) 
	catch
		Throw("Only for real or integer images.") 
		
//		
//		Check for selection as well as square and power of 2 feature: 
//		
	GetWindowSize(front, WindowSizeX, WindowSizeY) 
	GetSelection( front, tops, lefts, bottoms, rights) 
	GetSize(front, sizeXFull, sizeYFull) 
	if((rights-lefts) != bottoms - tops) 
		Throw( "Only for square images (selections) of the power of two!" ) 
	size = bottoms - tops 
	while( size / 2 >= 1 ) 
		size = size/2 
	if( size != 1 ) 
		Throw( "Only for square images (selections) of the power of two!" )
	sizeX = rights - lefts
	sizeY = bottoms - tops
		
	OpenResultsWindow() 
	Result("\n\nButterworth filter: " + datestamp() + "\n")
		
//		
//		Determine high / or low pass filter 
//		
	filtertype = "low" 
	if( !TwoButtonDialog( "\nHigh- or low- pass filter? ", "Low", "High" ) ) 
		filtertype = "high"
	Result( "  image: " + GetName( front ) + "\n") 
	Result( "  " + filtertype + "pass" + " filter\n") 
		
//		
//		Get image info: 
//		
	OpenAndSetProgressWindow("Running FFT", "...", "")
	GetWindowPosition( front, posX, posY)
	GetScale( front, scale, scaleY) 
	scale = Round(10000/(scale* sizeX/2))/10000 
	pixelsizeunit = GetUnitString( front ) 
	if(pixelsizeunit != "") 
	{	
		calibrate = 1 
		if(pixelSizeUnit == "nm") 
			pixelSizeUnitF = "1/nm" 
		if(pixelSizeUnit == "�m") 
			pixelSizeUnitF = "1/�m" 
		if(pixelSizeUnit == "�") 
			pixelSizeUnitF = "1/�" 
	}	
	zoom = GetZoom( front ) 
	displayzoom = zoom //sizeX / 256 
		
//		
//		test if image is o.k. for filtering 
//		
	mini = min(front) 
	Neg_Val = "" 
	if(mini < 0 ) 
	{	
		if(! TwoButtonDialog("Negatives values in image found!" + \
			 "\n\n  This may cause sever artefacts!" + \
			 "\n\n     Continue anyway? ", \
			 "YES", "NO") ) 
		{	
			CloseProgressWindow() 
			Result("exit: Negative values found in image\n")
			Throw(-128)
		}	
		Result("  (!) negative values found in image ..... continued \n")
		Neg_Val = "found" 
	}	
		
//		
//		Display FFT of front image and draw oval annotations 
//		
	frontF := RealFFT(front[tops, lefts, bottoms, rights]) 
	OpenAndSetProgressWindow("Getting ready", "for display", "(modulus(log10())) ... ")
	display = log10( modulus(frontF) ) 
	orderOfFilter = 3 
	R0 = sizeX / 8 
	sizeX2 = sizeX / 2 
	sizeY2 = sizeY / 2  
	top = sizeY / 2 - R0 
	left = sizeX / 2 - R0
	bottom = sizeY / 2 + R0
	right = sizeX / 2 + R0
	SetZoom( display, displayzoom ) 
	DisplayAt(display, posX + WindowSizeX * sizeX/SizeXFull + 24, posY + 17) 
	xcu = Round( R0 * 3**(-1/orderOfFilter) ) - R0  // 90% cutoff 
	xcv = Round( R0 * 3**( 1/orderOfFilter) ) - R0  // 10% cutoff 
	annot  = CreateOvalAnnotation(display, top+1, left+1, bottom, right ) 
	annot1 = CreateOvalAnnotation(display, top+1-xcu, left+1-xcu, bottom+xcu, right+xcu ) 
	annot2 = CreateOvalAnnotation(display, top+1-xcv, left+1-xcv, bottom+xcv, right+xcv ) 
	SetAnnotationBackground( display, annot, 4 ) 
	SetAnnotationBackground( display, annot1, 4 ) 
	SetAnnotationBackground( display, annot2, 4 ) 
	updateImage( display ) 
		
//		
//		Create and apply filter, draw filtered image to screen: 
//		
	filter := CreateFloatImage("filter", sizeX, sizeY)
	rad2 = R0 * R0 
	filter = (icol - sizeX/2) * (icol - sizeX/2) + (irow - sizeX/2) * (irow - sizeX/2) 
	filter = filter / rad2 
	if( orderOfFilter == 1 ) 
	{	
		filter = 1 + filter / rad2 
		filter = 1 / filter 
	}	
	else 
	{	
		test = filter 
		ii = 1 
		while( ii < orderOfFilter )
		{	
			test = test * filter 
			ii += 1 
		}	
		test = test + 1 
		filter = 1/test 
		deleteImage(test) 
	}	
	if(filtertype == "high") 
		filter = 1 - filter 
	SetPixel(filter, sizeX/2, sizeX/2, 1)
	frontF1 = frontF * filter 
	OpenAndSetProgressWindow("Running invers FFT", "...", "")
	final = modulus( IFFT(frontF1) ) 
	if(sizeXFull > sizeX || sizeYFull > sizeY) 
		SetName( final, GetName( front ) + ".sBW" ) 
	else 
		SetName( final, GetName( front ) + ".BW" ) 
	SetZoom( final, zoom ) 
	DisplayAt( final, posX + 10, posY + 17 )	
	updateImage( final ) 
	updateImage(display)
			
//			
	if(calibrate == 1 ) 
		OpenAndSetProgressWindow("R0 = " + R0*scale + "  " + pixelSizeUnitF, "order of filter: " + orderOfFilter, "spacebar to end") 
	else 
		OpenAndSetProgressWindow("R0 = " + R0 + " pixel", "order of filter: " + orderOfFilter, "spacebar to end")
		
//		
//			W H I L E    L O O P
//		
	Result("     up  / down -> change diameter\n");
	Result("     left/right -> change order\n");
	Result("     press space bar to continue\n")
	jj = 0;
	while( key != 32 && key2 != 32 && key != 13 && key2 != 13  && key != 3 && key2 != 3 ) 
	{		
		key = GetKey() 
		Delay(1) 
		key2 = GetKey() 
			
		if(key == 0 && key2 == 0) 
		{	
			jj += 1;
			if(jj == 2)
			{	
				rad2 = R0 * R0 
				filter = (icol - sizeX/2) * (icol - sizeX/2) + (irow - sizeX/2) * (irow - sizeX/2) 
				filter = filter / rad2 
				if( orderOfFilter == 1 ) 
				{	
					filter = 1 + filter / rad2 
					filter = 1 / filter 
				}	
				else 
				{	
					test = filter 
					ii = 1 
					while( ii < orderOfFilter )
					{	
						test = test * filter 
						ii += 1 
					}	
					test = test + 1 
					filter = 1/test 
					deleteImage(test) 
				}	
				if(filtertype == "high") 
					filter = 1 - filter 
				SetPixel(filter, sizeX/2, sizeX/2, 1)
				frontF1 = frontF * filter 
				final = modulus( IFFT(frontF1) ) 
				updateImage( final ) 
				updateImage(display) 
//				OpenAndSetProgressWindow("arrow keys: continue", "space bar: end", "")
			}	
			DoEvents()
		}	
		else
		if( key == rightarrow || key2 == rightarrow && orderOfFilter > 1 )
		{	
			jj = 0;
			if(orderOfFilter < 4.1) 
				orderOfFilter += -.25 
			else 
				orderOfFilter += -1 
			xcu = Round( R0 * 3**(-1/orderOfFilter) ) - R0  // 90% cutoff 
			xcv = Round( R0 * 3**( 1/orderOfFilter) ) - R0  // 10% cutoff 
			MoveAnnotation( display, annot1, top+1-xcu, left+1-xcu, bottom+xcu, right+xcu ) 
			MoveAnnotation( display, annot2, top+1-xcv, left+1-xcv, bottom+xcv, right+xcv ) 
			updateImage( display )
			if(calibrate == 1 ) 
				OpenAndSetProgressWindow("R0 = " + R0*scale + "  " + pixelSizeUnitF, "order of filter: " + orderOfFilter, "spacebar to end") 
			else 
				OpenAndSetProgressWindow("R0 = " + R0 + " pixel", "order of filter: " + orderOfFilter, "spacebar to end")
		}	
		else 
		if( key == leftarrow || key2 == leftarrow && orderOfFilter < 20 )
		{	
			jj = 0;
			if(orderOfFilter < 4) 
				orderOfFilter += .25 
			else
				orderOfFilter += 1 
			xcu = Round( R0 * 3**(-1/orderOfFilter) ) - R0  // 90% cutoff 
			xcv = Round( R0 * 3**( 1/orderOfFilter) ) - R0  // 10% cutoff 
			MoveAnnotation( display, annot1, top+1-xcu, left+1-xcu, bottom+xcu, right+xcu ) 
			MoveAnnotation( display, annot2, top+1-xcv, left+1-xcv, bottom+xcv, right+xcv ) 
			updateImage( display )
			if(calibrate == 1 ) 
				OpenAndSetProgressWindow("R0 = " + R0*scale + "  " + pixelSizeUnitF, "order of filter: " + orderOfFilter, "spacebar to end") 
			else 
				OpenAndSetProgressWindow("R0 = " + R0 + " pixel", "order of filter: " + orderOfFilter, "spacebar to end")
		}	
		else
		if( key == uparrow || key2 == uparrow && R0 <= sizeX2  )
		{	
			jj = 0;
			if( key == uparrow && key2 == uparrow && R0 <= sizeX2 - 4 ) 
				R0 += 4 
			else 
				R0 += 1 
			top = sizeY2 - R0 
			left = sizeX2 - R0
			bottom = sizeY2 + R0
			right = sizeX2 + R0
			xcu = Round( R0 * 3**(-1/orderOfFilter) ) - R0  // 90% cutoff 
			xcv = Round( R0 * 3**( 1/orderOfFilter) ) - R0  // 10% cutoff 
			MoveAnnotation( display, annot, top+1, left+1, bottom, right ) 
			MoveAnnotation( display, annot1, top+1-xcu, left+1-xcu, bottom+xcu, right+xcu ) 
			MoveAnnotation( display, annot2, top+1-xcv, left+1-xcv, bottom+xcv, right+xcv ) 
			updateImage( display )
			if(calibrate == 1 ) 
				OpenAndSetProgressWindow("R0 = " + R0*scale + "  " + pixelSizeUnitF, "order of filter: " + orderOfFilter, "spacebar to end") 
			else 
				OpenAndSetProgressWindow("R0 = " + R0 + " pixel", "order of filter: " + orderOfFilter, "spacebar to end")
		}
		else
		if( key == downarrow || key2 == downarrow && R0 > 1 )
		{	
			jj = 0;
			if( key == downarrow && key2 == downarrow && R0 > 4 ) 
				R0 += -4 
			else 
				R0 += -1 
			top = sizeY2 - R0 
			left = sizeX2 - R0
			bottom = sizeY2 + R0
			right = sizeX2 + R0
			xcu = Round( R0 * 3**(-1/orderOfFilter) ) - R0  // 90% cutoff 
			xcv = Round( R0 * 3**( 1/orderOfFilter) ) - R0  // 10% cutoff 
			MoveAnnotation( display, annot, top+1, left+1, bottom, right ) 
			MoveAnnotation( display, annot1, top+1-xcu, left+1-xcu, bottom+xcu, right+xcu ) 
			MoveAnnotation( display, annot2, top+1-xcv, left+1-xcv, bottom+xcv, right+xcv ) 
			updateImage( display )
			if(calibrate == 1 ) 
				OpenAndSetProgressWindow("R0 = " + R0*scale + "  " + pixelSizeUnitF, "order of filter: " + orderOfFilter, "spacebar to end") 
			else 
				OpenAndSetProgressWindow("R0 = " + R0 + " pixel", "order of filter: " + orderOfFilter, "spacebar to end")
		}
//		else
//		if(key == 13 || key2 == 13 || key == 3 || key2 == 3)
//		{	
//		}
	}	
	DeleteImage( display ) 
	nr = "th" 
	if(orderOfFilter == 1) 
		nr = "st" 
	if(orderOfFilter == 2) 
		nr = "nd" 
	if(orderOfFilter == 3) 
		nr = "rd" 
	if(orderOfFilter > 3 && orderOfFilter < 21) 
		nr = "th" 
	if(orderOfFilter > 20) 
	nr = "." 
				
//		
//		Finishing off: 
//		
	SetScale( final, scaleY, scaleY) 
	SetUnitString( final, pixelsizeunit) 
	CloseProgressWindow() 
	Result("  order: " + orderOfFilter + "\n") 
	if(calibrate == 1) 
		Result("  radius: " + R0*scale + "  " + pixelSizeUnitF + " \n") 
	else 
		Result("  radius: " + R0 + " pixel\n") 
	Result("output: " + GetName( front ) + ".BW\n\n") 
		
///////////////////////////////////
//		
//		If working on selection, should filter be applied to total image?  	
//		(only if full image is square and of the power of two)
//		
	if((sizeXFull != sizeX || sizeYFull != SizeY) && sizeXFull == SizeYFull) 
	{	
		size = sizeXFull
		while( size / 2 >= 1 ) 
			size = size/2 
		if( size == 1 ) 
		{	
			if(TwoButtonDialog("Apply filter to full image too?","YES","no")) 
				ApplyToFullImage = "yes"
			else 
				Throw(-128)
		}
	}	
	else 
		Throw(-128) 
		
/////////////////////////////////////////////////////////////////////
	OpenAndSetProgressWindow("applying filter", "    ...", "")
	R0 = R0 * sizeXFull/sizeX 
	rad2 = R0 * R0 
	filter := CreateFloatImage("filter", sizeXFull, sizeYFull)
	filter = (icol - sizeXFull/2) * (icol - sizeXFull/2) + (irow - sizeXFull/2) * (irow - sizeXFull/2) 
	filter = filter / rad2 
	if( orderOfFilter == 1 ) 
	{	
		filter = 1 + filter / rad2 
		filter = 1 / filter 
	}	
	else 
	{	
		test := filter 
		ii = 1 
		while( ii < orderOfFilter )
		{	
			test = test * filter 
			ii += 1 
		}	
		test = test + 1 
		filter = 1/test 
		deleteImage(test) 
	}	
	if(filtertype == "high") 
		filter = 1 - filter 
	SetPixel(filter, sizeXFull/2, sizeXFull/2, 1) 
	frontF2 := RealFFT(front[0,0,sizeYFull,sizeXFull]) * filter 
	OpenAndSetProgressWindow("Running invers FFT", "...", "")
	finalFull = modulus( IFFT(frontF2) ) 
	ShowImage(finalFull)
	updateImage( finalFull ) 
		
		
//		
//		Finishing off: 
//		
	SetScale( finalFull, scaleY, scaleY) 
	SetUnitString( finalFull, pixelsizeunit) 
	SetName( finalFull, GetName( front ) + ".BW" ) 
	CloseProgressWindow()