CM5 Program Listing
C *******************************************************************
C PROGRAM FOR TOMOGRAPHIC RECONSTRUCTON OF PARALLEL BEAM DATA
C ON THE CONNECTION MACHINE CM5.
C
C PROGRAMMER: RAMAN RAO
C DATE OF LAST REVISION: 1/2/95
C VIRGINIA POLYTECHNIC INSTITUTE AND STATE UNIVERSITY
C BLACKSBURG, VA 24061
C********************************************************************
C This program carries out CT reconstruction on the CM5
C for parallel beam projection data.
C
C odim=object dimension
C rdim=reconstruction dimension
C object=object (2D) to be forward projected
C lb1 = stores low bins for the forward projections
C hb1=stores high bins for the forward-projections
C mb1=stores the mid bin values for the forward projections
C frac1=stores the fractional part of the mid bin for forward projections
C step = the width of each block in the checker board
C angles=number of angles at which the forward projections
C are performed
C mat = projection matrix (proj. data is stored here)
C lb2 = stores low bins for the back-projections
C hb2=stores high bins for the back-projections
C mb2=stores the mid bin values for the back-projections
C frac2=stores the fractional part of the mid bin for back projections
C f = final density matrix for the reconstructed object
C Ham=hamming window function is generated and stored here
C H2=temporary row vector for storing a row
C lbmat=stores the filtered proj. values corresponding to the low bin
C hbmat=stores the filtred proj. values corresponding to the high bin
C psum=stores the interpolated values corresponding to the mid-bins
C xcen1, ycen1= center of the object dimension
C xcen2, ycen2=center of the reconstructed dimension
C xmax,ymax = maximum size of the z and y dimensions
C mp1= mid point of the rays
C tw = adjustment for keeping the mb inside the object(adjust until
C the program gives no error
C H1 stores the ramp function (response function)
C evalh = subroutine which generates the response function
C theta= angles which are used in the projection
C The CMF$ commands enable the respective arrays to be executed
C in parallel.
C The data of dimension odimxodim is reconstructed to rdimxrdim.
C In this process, zero padding is done to fill the frequency values
C for dimensions greater than odim (p. 75 kak) to reduce the dishing
C and dc artifacts(fig. 3.17,p.76 kak)
C
program filtered_bak
implicit none
C include the cmssl routines
include '/usr/include/cm/cmssl-cmf.h'
C include the timer routines
include '/usr/include/cm/timer-fort.h'
integer odim,rdim,step
parameter (odim=16,rdim=64,step=4,angles=128)
integer rows,cols,i,j,k,setup_id,ier,n1,rays,blank,i1,icen
integer hb1(1:odim,1:odim,1:128),lb1(1:odim,1:odim,1:128),j1
CMF$ LAYOUT hb1(:news,:news,:news),lb1(:news,:news,:news)
integer hb2(1:rdim,1:rdim,1:128),lb2(1:rdim,1:rdim,1:128)
CMF$ LAYOUT hb2(:news,:news,:news),lb2(:news,:news,:news)
complex mat(1:128,1:rdim),H1(1:rdim),H2(1:rdim)
CMF$ LAYOUT mat(:news,:news),H1(:news),H2(:news)
real tau1,tw,xcen1,xmax,ycen1,ymax,f(1:rdim,1:rdim),Ham(1:rdim)
CMF$ LAYOUT f(:news,:news),Ham(1:news)
real theta(1:128),frac1(1:odim,1:odim,1:128),
$ mb1(1:odim,1:odim,1:128)
CMF$ LAYOUT theta(:news),frac1(:news,:news,:news)
CMF$ LAYOUT mb1(:news,:news,:news)
real frac2(1:rdim,1:rdim,1:128),mb2(1:rdim,1:rdim,1:128)
CMF$ LAYOUT frac2(:news,:news,:news),mb2(:news,:news,:news)
real PI,psum(1:rdim,1:rdim,1:128),p1(1:odim,1:odim,1:128)
CMF$ LAYOUT psum(:news,:news,:news),p1(:news,:news,:news)
real p2(1:rdim,1:rdim,1:128)
CMF$ LAYOUT p2(:news,:news,:news)
real hbmat(1:rdim,1:rdim,1:128),lbmat(1:rdim,1:rdim,1:128),const
CMF$ LAYOUT hbmat(:news,:news,:news), lbmat(:news,:news,:news)
real object(1:odim,1:odim),den,frac
CMF$ LAYOUT object(:news,:news)
real xcen2,ycen2,mp1,mp2
integer lb,hb,timer1
timer1=0
C clear the timer
call CM_timer_clear(timer1)
C start the timer
call CM_timer_start(timer1)
xcen1=real(odim)/2.
ycen1=real(odim)/2.
xcen2 = real(rdim)/2.
ycen2 = real(rdim)/2.
tw = 1.0/(1.15 * sqrt(2.))
rays = odim
mp1=real(rays)/2.
mp2 = real(rdim)/2.
n1=6
C rows = # of angles
rows=128
cols=rdim
PI = acos(-1.0)
const = PI/real(rows)
C initialize the mat, H1 (ramp function), and object
forall(i=1:rows,j=1:cols) mat(i,j) = (0.0,0.0)
forall(j=1:cols) H1(j) = (0.0,0.0)
forall(i=1:odim,j=1:odim) object(i,j) =0.0
print *, "mat and H1 init completed"
call CM_timer_stop(timer1)
C this routine projects the data . The data is 128 128 x odim rows,
C and the data is supposed to look like a checker board.
C making the checker board
blank =0
do i =1,odim,step
if (blank .eq. 0) then
blank = 255
else
blank = 0
endif
do j=1,odim,step
if (blank .eq. 0) then
blank = 255
else
blank = 0
endif
forall(i1=i:i+step,j1=j:j+step)
$ object(i1,j1) = real(blank)
enddo
enddo
call CM_timer_start(timer1)
C Begins the forward projecton process. Projections are found by adding
C the densities along parallel rays (since this is parallel
C beam tomography)
C finds theta for (in radians) for all the steps
theta = real([1:128])*PI/128.
print *, 'theta evaluated'
C stop timer
call CM_timer_stop(timer1)
C finds the value of (xcos(theta)+ysin(theta)) for all the
C angles and all x,y. The values below are normalized to
C make the forward projection space to lie in the space (-1,1)
forall(i=1:odim,j=1:odim,k=1:128)
$ p1(i,j,k) = ((real(i)-xcen1)/xcen1)*cos(theta(k))
$ + ((real(j)-ycen1)/ycen1)*sin(theta(k))
print *, 'p1 evaluated'
C find the mid-bin of the forward projections
forall(i=1:odim,j=1:odim,k=1:128)
$ mb1(i,j,k) = (p1(i,j,k)*mp1*tw) + mp1
C find the low bins
lb1 = int(mb1)
print *, 'lb1 evaluated'
C find the fractional values
frac1 = mb1-real(lb1)
print *, 'frac1 evaluated'
C find the high bins
forall(i=1:odim,j=1:odim,k=1:128) hb1(i,j,k) = lb1(i,j,k) + 1
print *, 'hb1 evaluated'
C following is executed serially: not timed as these carry out the forward
C projections
do i=1,odim
do j=1,odim
do k=1,128
lb = lb1(i,j,k)
hb = hb1(i,j,k)
frac = frac1(i,j,k)
den = object(i,j)
mat(k,lb) = mat(k,lb)+cmplx((1.0-frac)*den)
mat(k,hb) = mat(k,hb)+cmplx(frac*den)
enddo
enddo
enddo
call CM_timer_start(timer1)
C evaluate the response function(ref. p72 kak)
call evalh(H1,tau1,n1)
print *, 'H matrix evaluation completed'
C Find the FFT of the projection matrix row by row.
C The parallel FFT on the CM5 is used to find the fourier
C transforms row by row.
do i=1,rows
forall(j=1:cols) H2(j) = mat(i,j)
C setup the FFT variables inside the CM5
setup_id = fft_setup(H2,'CTOC',ier)
call fft(H2,'CTOC',CMSSL_f_xform,setup_id,ier)
if (ier .ne. 0) then
write(*,*) 'error in fft','i=',i,'ier=',ier
endif
C deallocate the FFT setup variables
call deallocate_fft_setup( setup_id)
forall(j=1:cols) mat(i,j) = H2(j)
enddo
print *, 'forward transform of MAT completed'
C point to point multiply of the two fft's-the fft of the
C proj. matrix and the response function.
C Multiplies each row in parallel one row at a time
C Also multiplied point to point by the Hamming window function to
C reduce the gibbs phenomenon (p.75 kak)
forall(i = 1:cols)
$ Ham(i)=(0.54-0.46*cos((2.0*(i-1)*PI)/real(cols)))
do i=1,128
mat(i,:) = mat(i,:) * h1 * cmplx(Ham)
enddo
print *, 'point to point multiply completed'
C find the inverse FFT
C Inverse FFT is found by finding the complex conj of the
C matrrix row by row, finding the forward FFT , again finding
C the complex conj. and then scaling(ref. oppenheim and schafer
C : digital signal processing, or any other dsp book)
do i =1,rows
forall(j=1:cols) H2(j) = mat(i,j)
setup_id = fft_setup(H2,'CTOC',ier)
call fft(H2,'CTOC',CMSSL_i_xform,setup_id,ier)
if (ier .ne. 0) then
write(*,*) 'error in ifft','i=',i,'ier=',ier
endif
call deallocate_fft_setup( setup_id)
forall(j=1:cols) mat(i,j) = H2(j)
enddo
print *, 'inverse fft completed'
call CM_timer_stop(timer1)
write(*,*) 'time for calculating filtered projections'
call CM_timer_print(timer1)
call CM_timer_start(timer1)
C Interpolate to find the actual values. The interpolated
C values are stored in the variable 'psum'
C The code below actually does the backprojection, but some
C of the code which has been executed earlier is also
C timed as it is required for back-projection.
C Finds the value of (xcos(theta)+ysin(theta)) with x and y
C normalized to lie in the space (-1,1)
forall(i=1:rdim,j=1:rdim,k=1:128)
$ p2(i,j,k) = ((real(i)-xcen2)/xcen2)*cos(theta(k))
$ + ((real(j)-ycen2)/ycen2)*sin(theta(k))
print *, 'p2 evaluated'
C finds the mid bins in parallel
forall(i=1:rdim,j=1:rdim,k=1:128)
$ mb2(i,j,k) = (p2(i,j,k)*mp1*tw) + mp1
C finds the low bins in parallel
lb2 = int(mb2)
print *, 'lb evaluated'
C finds the fractional values of the mid bins
frac2 = mb2-real(lb2)
print *, 'frac2 evaluated'
C finds the high bins from the low bins
forall(i=1:rdim,j=1:rdim,k=1:128) hb2(i,j,k) = lb2(i,j,k) + 1
print *, 'hb2 evaluated'
C finds the filtered values corresponding to the low bins
forall(i=1:rdim,j=1:rdim,k=1:128)
$ lbmat(i,j,k) = real(mat(k,lb2(i,j,k)))
print *, 'lbmat2 evaluated'
C finds the filtered values corresponding to the high bins
forall(i=1:rdim,j=1:rdim,k=1:128)
$ hbmat(i,j,k) = real(mat(k,hb2(i,j,k)))
print *, 'hbmat2 evaluated'
C find the interpolated values corresponding to the fractional
C values just evaluated
psum = ((1.-frac2)*lbmat) + (frac2*hbmat)
print *, 'psum evaluated'
C finds the total summation over all the angles by summing over
C all the angles
f = SUM(psum,dim=3)
call CM_timer_stop(timer1)
call CM_timer_print(timer1)
print *,'f matrix evaluated'
C write the 'f' matrix to a file
open(unit=11,file='cm5_f.fil',status='old',form='formatted')
write(11,*), ((f(i,j),j=1,rdim),i=1,rdim)
close(11)
end
C finds the response function directly in the frequency
C domain. (Just a ramp function over the reconstruction space)
subroutine evalh(H,tau,N)
real tau, PI
integer N,j,i,k,ii,rdim
parameter (rdim=64)
complex H(1:rdim),a1(32)
CMF$ LAYOUT h(:news),a1(:news)
PI = acos(-1.0)
forall(j=1:rdim) H(j) = (0.,0.)
forall (i=1:32)
$ a1(i) = real(i)
do i =1,rdim
if (i .le. 32) then
H(i) = real(i)
else
H(i) = 65.-real(i)
endif
enddo
C can multipy with a Hamming function if necessary
c do i = 1,rdim
c H(i) = H(i)*(.54 + .46*cos(2*pi*real(i)/real(rdim)))
c enddo
return
end