#!/bin/csh -e
# spect,sim,sphere
# simulations of 3D SPECT using ASPIRE.
# examine count recovery for spheres in ellipse background
# for FBP/Chang vs OS-EM with full 3D detector response
# Copyright 2002-04-20, Jeff Fessler, The University of Michigan

which sed > /dev/null	# avoids 'sed not found' message
alias prom 'echo -n " "; echo -n "hit space: "; set dummy = $<'
alias prom 'exit'

set arg = ""
if ($#argv >= 1) then
	set arg = $1 
endif

#set	nx	= 128
set	nx	= 64
set	na	= 60	# projection views
set	nz	= `echo "$nx / 2" | bc`
if 1 then
	set scatperc	= 20	# 20% scatter
	set	vol	= 100	# sphere volume in cc = cm^3 ALL UNITS IN CM!
	set	backg	= 20	# elliptical cylinder background intensity
	set	foreg	= 100	# sphere intensity
	set	dsamp	= 1	# subsampling for partial volume
	set	case	= d$dsamp,v$vol,b$backg
endif

set	ny	= $nx	# required by spect3
set	nu	= $nx	# required by spect3
set	nv	= $nz	# required by spect3
set	orbit	= 360
set	ostart	= 0
set	dx	= `echo "0.7196 * 64 / $nx" | bc -l`	# pixel size in cm
set	dx	= `/usr/bin/printf %g $dx`		# neater
#set	du	= $dx	# ray spacing, must equal dx for 3s
set	count	= 1e6	# irrelevant since noise free

# show images
alias	ji	'echo DISPLAY \!*; j --red -b -1 -nrow 8 \!*'
alias	ji	'echo DISPLAY \!*; op disp -red -b -1 -ncol 8 \!*'
# display sinograms or projection views
#alias	js	'op -chat 0 transpose t9 \!* 1,2; j --red -b -1 -nrow 10'
alias	js	'echo DISPLAY \!*; j --red -b -1 -nrow 10 \!*' # show views
if 0 then	# for slow remote logins, or for users without 'j' script
	alias	j	'op range'
	alias	js	j
	alias	ji	j
endif
set	i	= "i -nthread 2" # because i have dual processors

set	data	= ./out,n$nx,$case
if !(-d $data) mkdir $data

#set	psf	= ./dat/psf,uhe,$nx.fld
set	gpsf	= gauss,0.5,1.2 # fwhm in cm at detector and isocenter

set	voi	= $data/voi.fld
set	xtrue	= $data/xtrue.fld
set	mask	= $data/mask.fld
set	mumap	= $data/mumap.fld
set	proj	= $data/proj,$gpsf.fld
set	ci	= $data/ci.fld
set	ybi	= $data/ybi.fld
set	chang	= $data/chang.fld

set	sfilter = 1
set	spre	= 3s@$dx,$dx,$dx,$orbit,$ostart,$sfilter
#set	stype	= $spre,conv1sym@$mumap@$psf@-$nx,$nz,$na
set	stype	= $spre,$gpsf@$mumap@-@-$nx,$nz,$na
set	stype0	= $spre,none@$mumap@-@-$nx,$nz,$na	# no psf for chang!

#
# voi
#
if ($arg == voi || !(-e $voi)) then
	# sphere radius in pixels from volume and pixel size:
	set r = `echo "e((1/3) * l(3/4 * $vol)) / $dx" | bc -l`
	set x = `echo "-20 * $nx" / 128 | bc -l`
	echo "sphere radius = $r"
	op ellipsoid $voi $nx $ny $nz \
		$x 0 0	$r $r $r	0 0 1 $dsamp
prom
endif


#
# object: background + voi
#
set ell_x = `echo "48 * $nx / 128" | bc`
set ell_y = `echo "32 * $nx / 128" | bc`
if ($arg == xtrue || !(-e $xtrue)) then

	# background ellipses
	set n0 = `echo "$nx / 16" | bc`
	set n1 = `echo "$nx / 16 * 6" | bc`
	op unif t0 $nx $ny $n0 0
	op ellipse t1 $nx $ny 0 0 $ell_x $ell_y 0 $backg $dsamp
	op rep t1 t1 $n1
	op stack t1 float t0 t1 t0	# 8 48 8 -> 64

	set lam = `echo "$foreg - $backg" | bc -l`	# incremental intensity

	op mul t0 $voi - $lam
	op add $xtrue t0 t1

	ji $xtrue
prom
endif


#
# mask
#
if ($arg == mask || (("$mask" != "-") && !(-e "$mask"))) then
	set rx = `echo "$ell_x + 4 * $nx / 128" | bc`
	set ry = `echo "$ell_y + 4 * $nx / 128" | bc`
	op ellipse t0 $nx $ny 0 0 $rx $ry 0 1 1
	op conv $mask t0 byte
	op add t0 $xtrue $mask
	ji t0
prom
endif


#
# attenuation map (mumap): a uniform ellipse
#
set mu = 0.11	# 0.11/cm for iodine 360KeV in water
if ($arg == mumap || (("$mumap" != "-") && !(-e "$mumap"))) then
	op ellipse $mumap $nx $ny 0 0 $ell_x $ell_y 0 $mu 3
	j --red $mumap
prom
endif


#
# PSFs
# here is where we would make customized psfs if needed.
# instead, use built-in gaussian psf model.
#
if 0 then #&& ($arg == psfs || (("$psf" != "-") && !(-e "$psf"))) then
	mat << END
	make_psf_uhe_22cm_128x64
END
	j --red $psf
	op range $psf
prom
endif


#
# noise-free projections
#
if (0 || $arg == proj || !(-e $proj)) then
	time $i proj3 $proj $xtrue $stype $mask
	js $proj
prom
endif


#
# Chang attenuation correction factors (for FBP)
#
if ($arg == chang || !(-e $chang)) then
	$i back3 t0 $nx $ny $nz - - $stype0 -
	op div0 $chang - t0 $na
	op range $chang
	ji $chang
	# compare to 2d approach - agrees since uniform ellipse all slices
	if 0 then
		set rx = `echo "$ell_x * $dx" |  bc`	# radius in mm
		set ry = `echo "$ell_y * $sx" |  bc`	# radius in mm
		op chang t0 $nx $ny $mu 0 0 $rx $ry $na $sx $ostart $orbit
		op range t0 $chang
	endif
prom
endif


#
# FBP/Chang with ramp filter
# trick: use pixel size = 1 for consistency with 3s projector 
#
set fbpwin = boxcar,1,1
set fbp0 = $data/fbp0,$fbpwin.fld
if ($arg == fbp0 || !(-e $fbp0)) then
	op transpose t0 $proj 1,2	# make sinograms nu,na,nv
	op fbp t0 t0 user,t $nx $ny $fbpwin 1 1 0  0 0 $orbit $ostart
	op mul $fbp0 t0 $chang
	op range $xtrue $fbp0
	ji $fbp0
prom
endif


#
# look at sensitivity pattern
# fix: the end slices are about 1% more sensitive - WHY?
#
if ($arg == sens) then
	op unif t0 $nu $nv $na
	$i back3 t0 $nx $ny $nz t0 - $stype $mask
	op mul t1 t0 $chang	# apply chang for fun
	op sub t2 t1 - $na	# subtract # of angles to make near 0
	j --red t2
prom
endif


#
# noise-free iterative reconstructions
# initialize with uniform image?????????????????????????????
#
if ($backg == 0) then
	set	init	= -$nx,$ny,$nz,1
else
	set	init	= -$nx,$ny,$nz,$backg
endif


#
# emulate effects of scatter
#
if 1 then
	set ri = $data/ri,$scatperc.fld
	if !(-e $ri) then
		op blur t0 $proj 10,20,zer 10,20,zer 1,1,zer
		op mul t0 t0 - 0.01
		op mul $ri t0 - $scatperc

		op add t99 $proj $ri
	endif
	set yi = t99
else
	set yi = $proj
endif


#
# EM and OSEM noise-free
# note: EM converges *much* faster for the point source
# with shift=tiny rather than shift=1 !
# shift=0 causes numerical precision problems (overflow?)
# set nsubset=1 for EM
#
set	niter	= 20
set	niter	= 2
set	penal	= -
set	nsubset	= 10
set	nsubset	= 5
set	otype	= osemc	# usually use this!
set	alg	= $otype,fast,$nsubset,$na,1.0
set	em0	= $data/${otype},s$scatperc,n$nsubset,i$niter.fld
set	method	= @$niter@$alg@$penal
set	saver	= stack,1
if ($arg == em0 || !(-e $em0)) then
	i -chat 15 -nthread 2 empl3 $em0 $init $yi - 1 $ri 1 $stype $mask \
		$method $saver 0 1 1e9 0 -
	op range $em0 $xtrue
	ji $em0
#prom
endif

exit ##################################################################