Hamill Contrib Library

This page is a listing of the entire contents of this library for IDL. This listing is the long version. Viewing the much more compact listing may be handier.

[Go Back to Main IDL Libraries Search Page]


Last modified: Thu Dec 21 21:32:51 2000.

List of Routines


Routine Descriptions

CHECK_TYPE_DIM

[Next Routine] [List of Routines]
 NAME:
	CHECK_TYPE_DIM

 PURPOSE:
	This function returns a TRUE or FALSE value (1 or 0) depending
	on the type of the variable, e.g. is it of type integer.

 CALLING SEQUENCE:
	RESULT = CHECK_DATA_TYPE ( VARIABLE )
	But the result only makes sense if one keyword or more are set.

 INPUTS:
	VARIABLE: an IDL variable to be checked.
	
 KEYWORD PARAMETERS:
	There are two classes of keywords.  First, there is ...
		DIMENSIONS	the number of dimensions; check that the
				input variable is as specified; this parameter
				can be a scalar (only one dimensionality is
				expected) or a vector (if several dimensionali-
				ties are allowed)
	Next, there are keywords that specify the data types expected for
	VARIABLE.  The choices are: UNDEFINED, BYTE, INTEGER, LONG, FLOAT,
	DOUBLE, COMPLEX, STRING, and STRUCTURE; and there is also the choice
	NUMERIC, which is the same as setting the BYTE, INTEGER, LONG, FLOAT,
	and DOUBLE keywords.

 OUTPUTS:
	The function returns 1 if the dimensions are correct and if the data
	type matches any set keyword, otherwise it returns 0.

 EXAMPLE:
	On entering a procedure, you can check that the input parameter P1 is
	a 1, 2, or 3-dimensional floating or double array, as follows:

	    if not check_data_type ( p1, dim=[1,2,3], /float, /double ) then $
		message, "Invalid parameter P1"

 MODIFICATION HISTORY:
 	Written by:	James Hamill
			Siemens Medical Systems
			2501 N. Barrington Rd.
			Hoffman Estates, IL  60195-7372
			(708)304-7760
			hamill@sgi.siemens.com
			October, 1993

(See /host/bluemoon/usr2/idllib/user_contrib/hamill/check_type_dim.pro)


FBAK_2D_MULTI

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	FBAK_2D_MULTI

 PURPOSE:
	This procedure will reconstruct a 3D volume from a set of projection
	images, using the filtered backprojection method and the formalism of
	nuclear medicine.  In nuclear medicine, one refers to a "stack" of 2D
	images that comprise the 3D one; hence the name.

 CATEGORY:
	Image reconstruction from projections.

 CALLING SEQUENCE:
	VOLUME = FBAK_2D_MULTI ( SINO )

 INPUTS:
	SINO: a stack of sinograms (N)x(n_slices)x(N_phi), where n_slices is
	the number of slices to reconstruct and N_phi is the number of view
	angles, which are assumed to be equally spaced.
	
 KEYWORD PARAMETERS:
	X_COR: A vector of center of rotation offsets, as required in nuclear
		medicine.  Default=ALL 0.0, meaning that each projection image
		is centered at x=(N/2).
	N_HALF_ROT: How many half rotations are represented by the input
		sinograms.  1 implies a 180 degree orbit, 2 implies a 360 degree
		orbit.  Default=1.
	X0: The image center in X and Y.  Default=(N-1)/2.0.
	PHI0: Starting view angle.  Default=0.0.
	QUIET: If this keyword is not set then by default we display the
		incomplete reconstructions as each view is processed.

 OUTPUT:
	The reconstructed volume.

 DEPENDENCIES:
	This procedure calls two user-contributed IDL procedures:
		STILL_BUSY
		and TVSTACK
	Their presence is only cosmetic, so you might want to delete the code
	that uses these procedures.

 SIDE EFFECTS:
	If the /QUIET keyword is not used, the image display will be used
	to display results, requiring the TVSTACK procedure.  The STILL_BUSY
	procedure is called to announce that this sometimes long-running program
	is still at work.

 PROCEDURE:
	The projection images are assumed to have been filtered in 2D before
	this routine is called, thereby eliminating the need for filters other
	than the ramp filter.  First, we compute the 1D space-domain ramp filter
	and we place it in a 2D array.  We loop over view angles (BY ANSATZ
	THESE ARE SPACED EQUALLY FROM PHI0 TO THE MAXIMUM ANGLE) and at each
	angle we: filter by multiplying each row of the projection image by this
	array, then we backproject into the image, using linear interpolation.

	This procedure is by no means optimized for speed, but the code is
	pretty streamlined.

 EXAMPLE:
	If the following commands are combined as a procedure and compiled and
	executed, we will generate projection images for a volume that contains
	three discs of radius 10, 20, and 30 pixels, and reconstruct them.  On a
	SUN IPX workstation, the example runs in ~ 2 minutes.

	n = 128L
	n_slices = 3L
	n_views = 90L
	x0 = [ 50, 60, 70 ]	; centers of the discs, x direction
	y0 = x0			; centers of the discs, y direction
	r0 = [ 10, 20, 30 ]	; disc radii
	dx = x0 - (n-1)/2.0	; centers relative to center of proj'n
	dy = y0 - (n-1)/2.0	; centers relative to center of proj'n
	
	;  Make sinogram by averaging over fine (expanded) samples.
	
	expander = 4L
	np = n*expander
	xp = findgen(np) - (np-1)/2.0	; rotated x, expanded coord.
	dxp = expander*dx		; relative centers, expanded
	dyp = expander*dy
	r0p = r0*expander		; radii in expanded coordinates

	sino = fltarr(n,n_slices,n_views)

	for view=0,n_views-1 do begin
	  angle = !pi*view/n_views
	  for slice=0,n_slices-1 do begin
	    dx_rotated_expanded = $
		dxp(slice)*cos(angle) + dyp(slice)*sin(angle)
	    expanded_projection = $
		sqrt(( r0p(slice)^2 - (xp-dx_rotated_expanded)^2 ) > 0 )
	    sino(0,slice,view) = rebin(expanded_projection,n)
	  endfor
	endfor

	;  Reconstruct.

	recon = fbak_2d_multi ( sino )

	end

 MODIFICATION HISTORY:
 	Written by:	James Hamill
			Siemens Medical Systems
			2501 N. Barrington Rd.
			Hoffman Estates, IL  60195-7372
			(708)304-7760
			hamill@sgi.siemens.com
			March, 1991

(See /host/bluemoon/usr2/idllib/user_contrib/hamill/fbak_2d_multi.pro)


IM_POISSON

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	IM_POISSON

 PURPOSE:
	Add Poisson noise to an array.

 CATEGORY:
	Image simulation.

 CALLING SEQUENCE:
	RESULT = POISSON_IMAGE ( IMAGE[, SEED] )

 INPUT:
	IMAGE: A numeric array (byte, integer, long, or float) of arbitrary
	dimensionality.  This is the array of values around which values in the
	result will be Poisson-distributed.

 OPTIONAL INPUTS:
	SEED: A longword seed for the random number generator.  If this is not
	supplied, the value -123456789L is used for generating the first random
	value.
	
 KEYWORD PARAMETER:
	OUTPUT_KIND: The data type of the output array, that is byte, integer,
	longword, or float.  The words "byte", "int", "integer","long",
	"longword", and "float", in upper or lower case, are accepted, as are
	the numeric IDL values 1,2,3,4 for byte, integer, longword, and float.

 OUTPUT:
	POISSON_IMAGE returns a copy of the input array, Poisson noise added.

 DEPENDENCIES:
	This procedure calls a user-contributed IDL procedure:
		STILL_BUSY

 RESTRICTIONS:
	Negative input values are mapped to a result of 0.  A call is made to
	the routine STILL_BUSY to provide messages once a minute during the
	processing of large arrays.  You can eliminate that line of code if you
	want to.

 PROCEDURE:
	Create an image with values Poisson-distributed around the mean values
	IMAGE, using Knuth's "Algorithm Q".  (Donald E. Knuth, The Art of
	Computer Programming, Volume 2, "Seminumerical Algorithms", Addison-
	Wesley (1969), page 117.  This routine IS NOT VECTORIZED, AND SHOULD RUN
	SLOWLY.  A deft IDL'er could probably vectorize the algorithm, and
	anyone who does so is entitled to a gold star.  Where the gaussian and
	Poisson distributions are essentially identical (mean value > 50) a
	normal, that is gaussian, distribution is used.

 EXAMPLE:
	Here is how you can create a 100x100 array of values Poisson-distributed
	around the mean value 5.0, and check the empirical probability against
	the Poisson distribution:

		n = 100
		a = replicate(5.0,n,n)
		b = poisson_image ( a, out="byte" )
		tvscl, b
		print, stdev(b,mean), mean, sqrt(5.0)
		h = histogram ( b )
		prob = float(h)/n_elements(b)
		probi = exp(-5.0)
		for i=0,10 do begin print, i, prob(i), probi & $
			probi=probi*5.0/(i+1)

 MODIFICATION HISTORY:
 	Written by:	James Hamill
			Siemens Medical Systems
			2501 N. Barrington Rd.
			Hoffman Estates, IL  60195-7372
			(708)304-7760
			hamill@sgi.siemens.com
			February, 1992

(See /host/bluemoon/usr2/idllib/user_contrib/hamill/poisson_image.pro)


STILL_BUSY

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	STILL_BUSY

 PURPOSE:
	This function can be embedded inside an IDL code that is expected to
	run for a long time, announcing with printed messages that the code is
	still running.  Unlike a PRINT statement in a loop, STILL_BUSY makes its
	announcements just once a minute.

 CATEGORY:
	Calming down impatient users like the author.

 CALLING SEQUENCE:
	STILL_BUSY[, ANNOUNCEMENT]

 OPTIONAL INPUT:
	ANNOUNCEMENT: An array of messages or values to be printed.  If no
	parameter is used, the message reads "STILL BUSY" plus the time and
	date.
	
 KEYWORD PARAMETERS:
	NONE.

 OUTPUT:
	NONE.

 COMMON BLOCKS:
	BIZZY

 SIDE EFFECTS:
	The overhead of calling this function many times might be onerous; so
	don't bury it deep inside nested loops.

 PROCEDURE:
	Each time the procedure is invoked, it decodes the !STIME system vari-
	able to see if the minute has changed since the last time.  If so, then
	the program prints the announcement.  If not, it returns and does
	nothing.

 EXAMPLE:
	The following code might be used in a coded procedure that runs for
	twenty minutes.  About twenty times out of the 1000 passages through the
	loop, a message appears on the screen.

		for iter = 1, 1000 do begin
		  number_cruncher, data, i
		  STILL_BUSY, "Processing iteration number " + strtrim(i,2)
		endfor

 MODIFICATION HISTORY:
 	Written by:	James Hamill
			Siemens Medical Systems
			2501 N. Barrington Rd.
			Hoffman Estates, IL  60195-7372
			(708)304-7760
			hamill@sgi.siemens.com
			October, 1989

(See /host/bluemoon/usr2/idllib/user_contrib/hamill/still_busy.pro)


TVSTACK

[Previous Routine] [List of Routines]
 NAME:
	TVSTACK

 PURPOSE:
	This procedure displays a multidimensional array (e.g. 3D) as a stack of
	images lain side by side, and by default scaled between the minimum and
	maximum values in the array.  This is a method of visualizing 3D data
	sets, for example reconstructions of tracer uptake in nuclear medicine.
	For an image measuring NX by NY by NZ, the effect is somewhat the same
	as typing

		for z=0,nz-1 do tvscl, image(*,*,z), z

	except that they aren't necessarily scaled to the frame maxima, and
	some options are available through the keywords.

 CATEGORY:
	IMAGE DISPLAY

 CALLING SEQUENCE:
	TVSTACK, IMAGES

 INPUTS:
	IMAGES: The multidimensional array to be displayed.  It may have any
	number of dimensions but is displayed as if it is 3D.
	
 KEYWORD PARAMETERS:
	W_BORDER: Width of the border to be drawn around each frame displayed.
		If framing is performed, each NX by NY image takes up
		(NX+2*W_BORDER) pixels in the X direction and (NY+2*W_BORDER)
		pixels in the Y direction.  Default = NO BORDER.
	C_BORDER: Color of the border, if the W_BORDER keyword is set.  Default=
		!d.n_colors/3.
	TVSCL: If this keyword is set, each 2D image is scaled to its own max
		and min.  By default, the max and min of the multidimensional
		array are used.
	POSITION: The screen position used by the first image frame in the
		sequence.  Thus, images are positioned as if the command

			for z=0,nz-1 do tvscl, image(*,*,z), position+z

		were used.  Default = 0 (start in upper left corner of th
		display window.)
	REVERSE_DIM: which dimension of the image is to be reversed (using the
		REVERSE function from the USERLIB); 1 means X, 2 means Y.
		Default = NO REVERSAL.
	ROTATE_PAR: what parameter to use to rotate the image frames, with the
		IDL ROTATE function (q.v.) in which 0=no rotation, 1=90 degrees
		CCW without transposing, and so on.  Defaul=NO ROTATION.
	SKIP: Use this keyword if you desire to display only some of the images
		for example every other one, which would be represented by
		SKIP=2.  Images are positioned as if the command

			for z=0,nz-1,SKIP do tvscl, image(*,*,z), z/SKIP

		were used.  Default = 1, that is, do not skip.  If you specify
		/SKIP you end up not skipping (!!!) because in this case you
		specify steps of 1.

 OUTPUTS:
	None; the procedure writes into the current display window.

 EXAMPLE:  A 3D gaussian can be displayed as follows:

	n = 50
	xysq = shift(dist(n)^2,n/2,n/2)	; x^2 + y^2
	rsq = fltarr(n,n,n)		; array to hold r^2
	for z=0,n-1 do rsq(0,0,z) = (z-n/2.0)^2 + xysq
					; r^2 = x^2 + y^2 + z^2
	fwhm=10.0 & sig=fwhm/(2*sqrt(2*alog(2))) & gauss_3D=exp(-rsq/(2*sig^2))
					; 3D gaussian of 10 channels FWHM
	tvstack, gauss_3D, /w_b		; display it

 MODIFICATION HISTORY:
 	Written by:	James Hamill
			Siemens Medical Systems
			2501 N. Barrington Rd.
			Hoffman Estates, IL  60195-7372
			(708)304-7760
			hamill@sgi.siemens.com
			April, 1992

(See /host/bluemoon/usr2/idllib/user_contrib/hamill/tvstack.pro)