Rotational average

/////////////////////////////////////////////////////////////////////
// (c) Gatan Inc.
/////////////////////////////////////////////////////////////////////
// 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. 
/////////////////////////////////////////////////////////////////////
//	last modified 08-July-2014 BS (comments only)


/*	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))