Viewing contents of file '../idllib/user_contrib/hamill/fbak_2d_multi.pro'
FUNCTION FBAK_2D_MULTI, SINO, X_COR=X_COR, N_HALF_ROT=N_HALF_ROT, $
X0=X0, PHI0=PHI0, QUIET=QUIET
;+
; 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
;
;-
;-------------------------------------------------------------------------------
; Check parameters.
;-------------------------------------------------------------------------------
on_error, 1
sizs = size(sino)
if sizs(0) ne 3 then stop,"Sinogram should be 3d: x-y-angle."
n = sizs(1)
n_slices = sizs(2)
n_phi = sizs(3)
case n_elements(x_cor) of
0: x_cor = replicate((n-1.0)/2.0,n_slices) ; DEFAULT
n_slices: begin ; OK; next check it's a vector
sizc = size(x_cor)
if sizc(0) ne 1 then message, "Please use a one-dimensional COR vector."
end
ELSE: message, "COR should agree with sinogram"
endcase
if keyword_set(n_half_rot) eq 0 then n_half_rot = 1
if keyword_set(x0) eq 0 then x0 = (n-1.0)/2.0 ; center of image
y0 = x0
if keyword_set(phi0) eq 0 then start_angle=0.0 else start_angle=phi0
if keyword_set(quiet) then display=0 else begin
; is there room to display all images or just some?
display = 1
n_im_x = !d.x_size/n
n_im_y = !d.y_size/n_slices
n_im_to_show = n_im_x*n_im_y
if n_im_to_show ge n_slices then begin
s1 = 0
s2 = n_slices - 1
endif else begin
s1 = (n_slices/2-n_im_to_show/2) > 0
s2 = s1 + n_im_to_show - 1
endelse
endelse
;-------------------------------------------------------------------------------
; Array of space-fixed x and y values, to be rotated slice by slice in the
; loop. Also, prepare a padded vector to store the projections during
; backprojection, removing all possibility of addressing outside of allocated
; memory.
;-------------------------------------------------------------------------------
x_array=fltarr(n,n) & for x=0,n-1 do x_array(x,*)=x-x0
y_array=fltarr(n,n) & for y=0,n-1 do y_array(*,y)=y-y0
offset = n ; pad by n on each end
projf = replicate ( 0.0, n + 2*offset ) ; vector to hold filtered projection
;-------------------------------------------------------------------------------
; Ramp filter in the space domain, ramp filter cut off at 1 x Nyquist. Filter
; values:
;
; at (i=0) = pi/4 ,
; at (odd i) = -1/(pi*i^2) ,
; at (even i) = 0.
;
;-------------------------------------------------------------------------------
still_busy
zero_value = !pi/4
vector = fltarr(n)
odd_ones = 2*lindgen(n/2)
vector(odd_ones) = -1.0/(!pi*(1+odd_ones)^2)
big_vector = [ reverse(vector), zero_value, vector ]
ramp_matrix = fltarr(n,n)
for y=0,n-1 do ramp_matrix(0,y) = big_vector((n-y):(n-y)+(n-1))
still_busy, "Filter matrix is generated ..."
;-------------------------------------------------------------------------------
; Loop over projections and slices. At each angle, we interpolate bilinearly.
;-------------------------------------------------------------------------------
recon = fltarr(n,n,n_slices)
for i_phi=0,n_phi-1 do begin
phi = start_angle + n_half_rot*!pi*i_phi/float(n_phi)
still_busy,"Backprojecting at phi="+ strtrim(phi*!radeg,2)+" degrees, . . ."
sin_phi = sin(phi)
cos_phi = cos(phi)
for slice=0,n_slices-1 do begin
xp = offset + x_cor(slice) + x_array*cos_phi + y_array*sin_phi
; rotated to space-fixed coordinates
ixp = fix(xp) ; lowest integer function "floor" for rotated x
ixp1 = ixp + 1 ; "ceil"
f_upper = xp - ixp ; coefficients above for linear interpolation
f_lower = 1.0 - f_upper ; coefficients below ...
projf(offset) = ramp_matrix#sino(*,slice,i_phi)
; filtration
recon(0,0,slice) = recon(*,*,slice) + $
f_lower*projf(ixp) + f_upper*projf(ixp1)
; add into the reconstruction array
endfor
if display then tvstack, recon(*,*,s1:s2)
endfor
;-------------------------------------------------------------------------------
; Remove the area outside the inscribed circle.
;-------------------------------------------------------------------------------
included = float(shift(dist(n),n/2,n/2) lt n/2)
for slice=0,n_slices-1 do recon(0,0,slice) = included*recon(*,*,slice)
return, recon
end