Appendix

A Novel Approach to Spectral MRI Tiffany A. Fetzner

Appendix A: IDL code for running experimental tests with Projection Repetition used

PRO Exp_Rie
;pro Exp_Rie is an IDL procedure for computing the inverse radon transform of experimental test images

TargetDegrees=180.0/256.0
NVIEWS=fix((180/TargetDegrees)+1)
ImageAngles=(180*asin(fix(256*sin(findgen(NVIEWS)*TargetDegrees*!PI/180)+ $
0.5+(findgen(NVIEWS) eq 0))/256.0)/!PI)* $
(1-(2*((findgen(NVIEWS)*TargetDegrees) gt 90.0)))+ $
(180.0*((findgen(NVIEWS)*TargetDegrees) gt 90.0))

N = 256L    ;Define number of columns in A.
M = 256L    ;Define number of rows in A.
;nviews = 255    ;Number of views.
K = CEIL(SQRT(N^2 + M^2))

;A=FLOAT(READ_SKW(PICKFILE(FILTER='*.skw')))
A_tiff=Findgen(3,256,256)
A_tiff=FLOAT(READ_TIFF(PICKFILE(FILTER='*.tif')))
;A=(findgen(256,256) eq 32639)*255
A=findgen(256,256)
FOR xxx=0, 255 DO BEGIN
FOR yyy=0, 255 DO BEGIN
A[xxx,yyy]=(A_tiff[0,xxx,yyy]+A_tiff[0,xxx,yyy]+A_tiff[0,xxx,yyy])/3
ENDFOR
ENDFOR
A=A-Min(A)
A=A*(255.5/MAX(A))
tv,FIX(A)
pp = FLTARR(K, nviews)    ;Make array for sinogram.
pp_th = FLTARR(K, nviews)    ;Make array for sinogram.
FOR I=0, NVIEWS-1 DO BEGIN
RIEMANN, pp, A, ((ImageAngles[I]*!PI)/180), ROW=i, cubic=-0.5
RIEMANN, pp_th, A, I * !PI/nviews, ROW=i, cubic=-0.5
;if ((I mod 8) eq 0) Then tvscl,pp
ENDFOR
      ;Compute each view.
temp_pp=pp-Min(pp)
temp_pp = temp_pp * (255.5/MAX(temp_pp))
TV, temp_pp
Print, 'Experimental Radon I5'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), pp, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), temp_pp
         ;Show sinogram.

temp_pp_th=pp_th-Min(pp_th)
temp_pp_th = temp_pp_th * (255.5/MAX(temp_pp_th))
TV, temp_pp_th
Print, 'Theoretical Radon I2'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), pp, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), temp_pp_th
;Show sinogram.

Y=(findgen(256)-128)
Ramlak=(abs(Y mod 2 ) * ((-4*1)/((!PI*Y^2)+((Y-1) mod 2)))) +(Y EQ 0)*(!PI)

wide_pp=findgen(fix(512*sqrt(2)),nviews)*0.0
wide_pp_th=findgen(fix(512*sqrt(2)),nviews)*0.0

   ;FOR xx=0,(fix(256*sqrt(2))-1),1 DO BEGIN
   ; ;FOR yy=0,(nviews-1),1 DO BEGIN
   ; ; wide_pp(xx+fix(127*sqrt(2)),yy)=pp(xx,yy)
   ; ;ENDFOR
   ; wide_pp(xx+fix(128*sqrt(2)),0:(nviews-1))=pp(xx,0:(nviews-1))
   ;ENDFOR
   wide_pp[181:542,0:(nviews-1)]=pp[0:361,0:(nviews-1)]
   wide_pp_th[181:542,0:(nviews-1)]=pp_th[0:361,0:(nviews-1)]
tvscl,wide_pp
tvscl,wide_pp_th
INV_FT_RamFilter=convol(wide_pp,Ramlak)
INV_FT_RamFilter_th=convol(wide_pp_th,Ramlak)
temp_RamFilter = INV_FT_RamFilter - MIN(INV_FT_RamFilter)
temp_RamFilter = temp_RamFilter * (255.5/MAX(temp_RamFilter))
tv,temp_RamFilter

;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), temp_RamFilter, CENTER=1
;WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), temp_RamFilter
;temp_RamFilter_th = INV_FT_RamFilter_th - MIN(INV_FT_RamFilter_th)
;temp_RamFilter_th = temp_RamFilter_th * (255.5/MAX(temp_RamFilter_th))
;tv,temp_RamFilter_th
;
;
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), temp_RamFilter, CENTER=1
;WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), temp_RamFilter

B = FLTARR(N,M) ;Initial reconstructed image.
B_th = FLTARR(N,M) ;Initial reconstructed image.
FOR I=0, nviews-1 DO BEGIN
RIEMANN, INV_FT_RamFilter, B, ((ImageAngles[I]*!PI)/180), /BACKPROJECT, ROW=i,cubic=-0.5
RIEMANN, INV_FT_RamFilter_th, B_th, I * !PI/nviews, /BACKPROJECT, ROW=i,cubic=-0.5
;if ((I mod 8) eq 0) THEN TVSCL,B
ENDFOR

;adjustments=findgen(17,3)
;checked=-2
;for b_shift=28000L,36000L, 1L do BEGIN
; for scl_pwr=120,120 do Begin
; print, ((scl_pwr+(160*(256.0+b_shift)))*10/8208), '% done'
; b_scale=(2.0)^((scl_pwr-80)/4)
B_new=((B-29625)/1024)
B_temp = B_new - MIN(B_new)
B_new_th=((b_th-29625)/1024)
B_temp_th = B_new_th - MIN(B_new_th)
;print, (max(B_temp)/255.5)
;mom_b=median(b[0:64,192:255])
;mom_a=median(a[0:64,192:255])
;print, mom_b, mom_a,mom_b/mom_a
B_temp = B_temp / (MAX(B_temp)/255.5)
;Do the backprojection for each view.
TV, B_temp ;Show reconstructed array.
Print, 'Experimental Reconstruction I6'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), B, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), B_temp

B_temp_th = B_temp_th / (MAX(B_temp_th)/255.5)
;Do the backprojection for each view.
TV, B_temp_th ;Show reconstructed array.
Print, 'Theoretical reconstruction I3'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), B, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), B_temp_th
c = FLTARR(N,M)
e = FLTARR(N,M)
g = FLTARR(N,M)
c=A-B_new
c_tmp=c-Min(c)
c_tmp=c_tmp*(255.5/Max(c_tmp))
tv, c_tmp
Print, 'Difference Image I7'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), c_tmp, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), c_tmp
d=(Total(c^2))
print, 'Standard deviation I7:',d
e=B_new-B_new_th
e_tmp=e-Min(e)
e_tmp=e_tmp*(255.5/Max(e_tmp))
tv, e_tmp
Print, 'Difference Image I8'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), e_tmp, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), e_tmp
f=(Total(e^2))
print, 'Standard deviation I8:',f

g=A-B_new_th
g_tmp=g-Min(g)
g_tmp=g_tmp*(255.5/Max(g_tmp))
tv, g_tmp
print, 'Difference Image I4:'
;WRITE_SKW, PICKFILE(FILTER='*.skw', /WRITE), g_tmp, CENTER=1
WRITE_TIFF, PICKFILE(FILTER='*.tif', /write), g_tmp
H=(Total(g^2))
print, 'Standard deviation:',H
;if (checked eq (0-2) ) then begin
; adjustments[0:15,0]=(findgen(16,1)*0.0)+d
; adjustments[0:15,1]=(findgen(16,1)*0.0)+b_scale
; adjustments[0:15,2]=(findgen(16,1)*0.0)+b_shift
; print, adjustments
; print, 'you should only see this once!'
; checked=-1
;endif Else Begin
; checked =-1
; while ((adjustments[(1+checked),0] lt d) and (checked lt 15)) DO checked = checked + 1
; for counter=16,min([16, (checked+2)]), -1 DO BEGIN
; adjustments[counter,0:2]=adjustments[(counter-1),0:2]
; ENDFOR
; adjustments[(checked+1),0]=d
; adjustments[(checked+1),1]=b_scale
; adjustments[(checked+1),2]=b_shift
; if (checked eq -1) then print, adjustments
; checked=-1
;endelse
;endfor
;endfor
;print, adjustments
END

Table of Contents | Appendix | Thesis