Merge images

/////////////////////////////////////////////////////////////////////
// (c) Gatan Inc.
/////////////////////////////////////////////////////////////////////
// Merges two images of different intensisties replacing saturated pixels
// in the stronger image with scaled values from the weaker image.
// Especially useful for merging diffraction patterns.
//	
// Asks for a threshold value which should be less than the value of the saturated
// pixels.  Choose a value that excludes the saturated values as well as nearby 
// pixels that will contain some of the spilled over saturated charge.  Improper 
// selection will result in an incorrect scaling ratio that should be observable
// using a line plot in the merged image.
/////////////////////////////////////////////////////////////////////
//	last modified 08-July-2014 BS



// This littel method creates dummy images for testing.
image GetDummyExposure( number exp )
{
	image img:=RealImage("Test",4,256,256)
	number r1=20
	number r2=30
	img = 1 + sqrt((mod(icol,2*r2)-r2)**2 + (mod(irow,2*r2)-r2)**2)
	img = ( img<r2 ? exp(-((img-r1))) : 0  )
	img = log(1+img*img*abs(itheta*icol))*icol*abs(itheta)+sin((iwidth/2-iradius)+itheta*Pi())*50
	img -= min(img)
	img /= max(img)	
	img *= exp
	return (img>6000?6000:img)+random()*300
}


If ( TwoButtonDialog( "Do you want to create 2 dummy images?", "Yes", "No" ) )
{
	image lowExp = GetDummyExposure(10000 )
	lowExp.SetName( "low exposure (no saturation, too noisy in background" )
	image highExp = GetDummyExposure(100000 )
	highExp .SetName( "high exposure (saturation, good background" )
	DisplayAt( lowExp, 100, 100 )
	DisplayAt( highExp, 100, 400 )
}


// Set the default threshold value 
number threshold = 3500.												

// Prompt the user to select the two images for merging
image img1, img2														
if ( !GetTwoImages( "Select the two images for merging", img1, img2 ))	
	exit(0)

// Check both images are 2D...if not, throw an appropriate error
if (img1.ImageGetNumDimensions() != 2 || img2.ImageGetNumDimensions() != 2) 			
	Throw("Both input images must be 2-dimensional.")					


// Check the dimension sizes of the input images
// If they are incompatible, throw an error message and halt
number xs1, xs2, ys1, ys2												
img1.Get2dSize( xs1, ys1 )												
img2.Get2dSize( xs2, ys2 )												
if ( xs1 != xs2 || ys1 != ys2 )											
	Throw("Input images are not the same size.")

// Prompt the user to specify  a cutoff threshold; exit if cancel is hit
if ( !GetNumber( "Threshold (discard pixels above this value): ", threshold, threshold ) )	
	exit(0)																

// Determine which is the saturated (highest intensity) and low intensity images 
image hi_img =  ( sum( img1 ) <  sum( img2 ) ) ? img2 : img1			
image lo_img =  ( sum( img1 ) <= sum( img2 ) ) ? img1 : img2			

// Calculate the ratio of intensities in regions below the threshold intensity
number ratio = sum(hi_img*(hi_img < threshold)) / sum(lo_img*(hi_img < threshold))			

// Merge the images; if the intensity in the saturated image is below the threshold, 
// use this pixel divided by the ratio.  Else, use the low intensity pixel.
image merged_img = ( hi_img <= threshold ) && (ratio > 0) ? hi_img/ratio  : lo_img   

// Set the name of the merged image...and display it  
merged_img.SetName("Merged image")										
ShowImage( merged_img )
													
// Finally, post the computed ratio
Result( "Scaling ratio = " + ratio + "\n" )