/*----------------------------------------------------------------------------------------------

Hit ctrl-enter to execute.

Description:

	This script performs rotational average of an image. An image is passed in. If a diffractogram, it 
	is packed complex as a result of FFT and data type is changed to real by means of modulus extraction. 
	Then the image dimension (hight and width) is found from the image. This information is used 
	to declare several intermediate images such as "dst", "line_projection", and "rotational_average"

	The image is displayed in cardinal coordinates, and needs to be transformed into polar coordinates. 
	The passed variable "sample" defines the number of segment of the 360 degree angular range. 
	The transformation of coordinates ("img" in cardinal, "dst" in polar) is realized by bilinear inter-
	polation (the use of warp function). Then line projection is easily calculated by adding up all
	the columns in "dst" (polar coordinates). The averaged line intensity is also normalized by the
	number of segment ("sample"). Finally this (1D) line profile is used to generate a (2D) image --
	"rotational_average" by using again the bilinear interpolation technique (the warp function). The
	last line is to display the rotational averaged image using the iradius intrinsic variable. 

	By Ming Pan,  Gatan, Inc.
	Modified and annotated by Paul Thomas, March 2001
	BugFix by RTH, June 2002

-----------------------------------------------------------------------------------------------*/


/*	The Rotational_Average method performs the rotational averaging.
	It is passed a source image and the sampling density 
	i.e. the number of segments to split the full 360 degree
	angular range into.  It returns the rotationally averaged image */ 

image Rotational_Average(image img, number samples)
{
	// Define neccessary parameters and constants
	number pi = 3.1416
	number xscale, yscale, xsize, ysize
	number centerx, centery, halfMinor
	number scale = img.ImageGetDimensionScale(0)
	string unit = img.ImageGetDimensionUnitString(0)

	// Likewise, declare intermediate images
	image rotational_average, dst, line_projection
	
	// If the source image is complex, take the modulus						
	if ( img.ImageIsDataTypeComplex( ))
		img := modulus(img)

	// Get the dimension sizes, and determine half the smallest dimension 
	img.Get2dSize( xsize, ysize )
	halfMinor = min( xsize, ysize )/2

	// Find the centre of the image
	centerx = xsize / 2
	centery = ysize / 2

	// Convert the image to polar co-ordinates...
	dst := RealImage( "dst", 4, halfMinor, samples )
	dst = warp( img, icol*sin(irow*2*pi/samples) + \
			centerx, icol*cos(irow*2*pi/samples) + centery )

	// and create a line projection using the icol intrinsic variable, 
	// normalising with the sampling density
	line_projection := RealImage( "line projection", 4, halfMinor, 1 )
	line_projection = 0
	line_projection[icol,0] += dst
	line_projection /= samples

	// Create rotationally averaged image using the line-projection and the iradius intrinsic variable
	rotational_average := RealImage( "Rotational Average", 4, halfMinor * 2, halfMinor * 2 )
	rotational_average = warp( line_projection, iradius, 0 )
	rotational_average.ImageSetDimensionScale( 0, scale )
	rotational_average.ImageSetDimensionScale( 1, scale )
	rotational_average.ImageSetDimensionUnitString( 0, unit )
	rotational_average.ImageSetDimensionUnitString( 1, unit )

	// Finally, return the rotationally averaged image
	return rotational_average
}

//  Call the method using the front image and a sampling density of 256
//  Display the result  
ShowImage( Rotational_Average(GetFrontImage(), 256))

