#!/bin/csh -e
# script,recon/sim2d
#
# This is an example script that generates simulated PET transmission
# and emission scans, and performs 2D transmission and emission reconstruction
# using penalized-likelihood methods.
# It uses many ASPIRE functions, so should be useful to new users
# as a template and for testing new installations of ASPIRE.
#
# The final output images will be put on the web site, for comparison.
#
# The only software needed is the three ASPIRE binaries: "op, wt, i"
# and the "xv" program (or something equivalent for image display),
# and of course the unix "csh" shell.
#
# Jeff Fessler, fessler@umich.edu, June 1998, University of Michigan


#set data = /users/fessler/data/sim,tran
set data = .	;		# where should output data go?
set exit = exit	;		# halt execution after each step?
alias j 'op -chat 0 disp \!*'	# jiffy display

set arg = ""
if ($#argv < 1) then
	echo "usage: $0 [arg]"
else
	set arg = $1
endif

# transmission scan parameters
set	tcount	= 3e6
set	trandpercent	= 5
set	seed	= $trandpercent
set	nrep	= 1
set	trun	= ,$tcount,r$trandpercent,s$seed,$nrep

set	ps	= 4.2;	# pix_size
set	rs	= 3.4;	# ray_spacing
set	sw	= $rs;	# strip_width
set	nx	= 128
set	ny	= 64
set	nb	= 160	# sinogram radial bins
set	na	= 192	# projection angles
set	tdsc	= $data/t,$nx,$ny,$nb,$na.dsc
set	twtf	= ${tdsc:r}.wtf
set	twtr	= ${tdsc:r},row.wtf
set	tpsf	= ${tdsc:r},psf.fld
set	tmap	= $data/tmap,$nx,$ny.fld	# attenuation map
set	tproj	= $data/tproj,$nb,$na.fld	# li
set	ttrues	= $data/ttrues,$tcount.fld	# bi exp(-li)
set	tblank	= $data/tblank,$tcount.fld	# bi
set	tprompt	= $data/tprompt$trun.fld	# yi
set	tline	= -;	# $data/tline$trun.fld	# not needed now
#set	tinfo	= $data/tinfo$trun.fld
#set	tmask	= $data/tmask.fld
set	tmask	= -;	# not needed

set	t = "`echo $tcount | sed 's/e/ * 10^/'`"
set	trandmean = `echo "$t * $trandpercent / 100 / $nb / $na" | bc -l`
echo trandmean=$trandmean tcount=$tcount

if $arg == reset then
	rm $data/*.wtf $data/*.dsc $data/*.fld
endif

# clean work files
touch tmppop.pgm; rm -f t.fld t.mat t.dat tmppop*.pgm


#
# build transmission system matrix file
#
if !(-e $tdsc) then
	echo "build system matrix"
	wt -chat 0 dsc 2 nx $nx ny $ny nb $nb na $na \
		pixel_size $ps ray_spacing $rs strip_width $sw \
		scale 0 xrad 62 yrad 30 >! $tdsc
endif
if !(-e $twtf) then
	wt gen $tdsc
	wt col2row $twtr $twtf;	# convert to row grouped for OS methods
$exit
endif

#
# synthesize attenuation map using ellipses
# with lung=0.002 bone=0.016 water=0.0096
#
if (!(-e $tmap) || $arg == tmap) then
	echo "build attenuation map"
	op ellipse $tmap $nx $ny \
		-25 0	15 18	0 -0.0076 2 \
		25 0	15 18	0 -0.0076 2 \
		0 -15	5 5	0 0.0064 2 \
		0 0	57 27	0 0.0096 2
	j $tmap
$exit
endif

#
# transmission noiseless projection
#
if (!(-e $tproj) || $arg == tproj) then
	echo "project attenuation map"
	i proj2 $tproj $tmap $twtf
	j $tproj
$exit
endif

#
# blank scan (noiseless), with log-normal detector efficiencies
#
if (!(-e $tblank) || $arg == tblank) then
	echo "make blank scan"
	op sim blank $ttrues $tblank $tproj 0.3 $tcount $seed
	j $tblank $ttrues
$exit
endif

#
# noisy transmission scan
#
if (!(-e $tprompt) || $arg == tprompt) then
	echo "make transmission scan"
	op sim pet $tprompt - $ttrues $nrep $seed - $trandpercent 1
	j -red $tprompt $ttrues
$exit
endif

#
# Noisy FBP transmission reconstruction (to initialize iteration)
#
	set tfbpwin = "cls3sinc,13,1"
	set tfbp = $data/tfbp,$tfbpwin$trun.fld
if (!(-e $tfbp) || $arg == tfbp) then
	echo "transmission FBP"
	i fbp2t dsc $tfbp - $tprompt $tblank 1 - $trandmean $tdsc $tfbpwin -

#	echo y | op nonlin max $tfbp $tfbp 0 0
	j -red $tfbp $tmap
$exit
endif

#
# make support mask from transmission FBP image
# not needed, not tested.
#
if ("$tmask" != "-" && (!(-e "$tmask") || $arg == tmask)) then
	echo y | op pre attn mask $tmask $tfbp 0.001 1 "f e1-+3,3 d2-+3,3"
	# make sure within .wtf support
	i support t.fld $twtf
	echo y | op mul $tmask $tmask t.fld
	j -s -m2 t.fld $tmask $tmap
$exit
endif

#
# \hat{l_i} and weights for PWLS, note no smoothing!
# not needed, not tested.
#
if ("$tline" != "-" && (!(-e "$tline") || $arg == tinfo)) then
	op pre trans blank $tline $tinfo $tprompt - $tblank - $trandmean 0 1 1 0 0 0 0
	j $tinfo 0 0
	j --red -a $tline $tproj
$exit
endif

#
# psf of A'A
# not essential
#
if (0 && !(-e $tpsf) || $arg == tpsf) then
	i proj2 $tpsf - $tdsc -p
	echo y | op psfsym $tpsf $tpsf
	j $tpsf
#	op psfpls - - $tpsf $tq_l2b
$exit
endif


#
# transmission penalized-likelihood reconstruction
#

# transmission algorithm (OSTR, Erdogan, 1999, PMB)
set	nsub	= 8;	# 8 subsets
set	talg	= osc,$nsub,1.0

# transmission penalty (edge preserving)
#set	l2b	= 18;	set proot = quad,1,-
#set	ureg	= $data/ureg.wtf
#set	ureg	= $data/ureg,all-support-only.wtf
#set	l2b	= 21;	set proot = lange3,1,ureg:${ureg}:,0.0003,4
#set	l2b	= 18;	set proot = quad,0,ureg:${ureg}:
#set	l2b	= 21;	set proot = lange3,1,-,0.0003,4
#set	l2b	= 13;	set proot = quad,1,center
#set	tq_l2b	= 13
set	tn_l2b	= 21;	set proot = huber,2,-,0.0003,ih,4; # play with these...
set	tpenal = $tn_l2b,$proot;	unset proot 

set init	= $tfbp
set niter	= 6
set saver	= -
set method	= @$niter@$talg@$tpenal
set tpl		= $data/tpl,$talg,$tpenal$trun.fld
set slice	= -
if (!(-e $tpl) || $arg == tpl) then
	set tmask = -
	i trpl2 $tpl $init $tprompt $tblank 1 - $trandmean $twtr $tmask \
		$method $saver 1 1 1e9 0 $slice
	j -red $tpl $tfbp $tmap; # look how nice iterative is vs FBP!
$exit
endif


#
##################################################################
# Emission Simulation
##################################################################
#

set	ecount	= 1e6
set	erandpercent	= 15
set	eseed	= $erandpercent
set	erun	= ,$ecount,r$erandpercent,s$eseed,$nrep

set	edsc	= $tdsc
set	ewtf	= $twtf
set	ewtr	= $twtr
set	emap	= $data/emap,$nx,$ny.fld;	# true emission map
set	eproj	= $data/eproj,$nb,$na.fld
#set	line	= $data/line,$nb,$na.fld
set	etrues	= $data/etrues,$ecount.fld
set	ecalib	= $data/ecalib,$ecount.fld;	# ray-dependent factors (c_i)
set	ecalib_noisy	= $data/ecalib,noisy,$ecount.fld
set	eprompt	= $data/eprompt$erun.fld;	# raw noisy prompt coincidences
#set	info	= - # $data/info$erun.fld

set	t = "`echo $ecount | sed 's/e/ * 10^/'`"
set	erandmean = `echo "$t * $erandpercent / 100 / $nb / $na" | bc -l`
echo erandmean=$erandmean ecount=$ecount

# emission system can just match transmission model for simplicity
#if !(-e $edsc) then
#	wt -chat 0 dsc 2 nx $nx ny $ny nb $nb na $na \
#		pixel_size $ps ray_spacing $rs strip_width $sw \
#		xrad 62 yrad 30 >! $edsc
#endif
#if !(-e $ewtf) then
#	wt gen $edsc
#endif

#
# synthesize emission object - silly heart image
#
if (!(-e $emap) || $arg == emap) then
	op ellipse t.fld $nx $ny \
		0 10	8 10	-20 4 2 \
		0 10	5 7	-20 -4 2
	echo y | op rect $emap $nx $ny 1  64 17 0  22 16 1 1
	echo y | op mul t.fld t.fld $emap
	echo y | op ellipse $emap $nx $ny \
		-25 0	15 18	0 -1 2 \
		25 0	15 18	0 -1 2 \
		0 -15	5 5	0 -2 2 \
		0 0	57 27	0 2 2
	echo y | op add $emap $emap t.fld
	j -m 2 $emap
$exit
endif

#
# emission noiseless sinogram
#
if (!(-e $eproj) || $arg == eproj) then
	i proj2 $eproj $emap $ewtf
	j $eproj
$exit
endif

#
# 'calibration' sinogram (c_i), based on true survival probabilities
# including log-normal detector efficiencies
# c_i = T d_i exp(-l_i)	where d_i is det. effic.
# and l_i is line integral through attenuation map and T is scan time
#
if (!(-e $ecalib) || $arg == ecalib) then
	op nonlin exp t.fld $tproj -1 1		# exp(-li) survival prob.
	op range t.fld
	op sim calib $etrues $ecalib $eproj t.fld 0.3 $ecount $eseed
	op range $etrues
	j $ecalib; # that is not "noise" but is random detector efficiencies.
$exit
endif

#
# noisy sinograms
#
if (!(-e $eprompt) || $arg == eprompt) then
	op sim pet $eprompt - $etrues $nrep $eseed - $erandpercent 1
	j -red $eprompt; # show zero pixels as 'blue'
$exit
endif

#
# Use noisy ACFs from estimated attenuation map.
# set to "if 0" to use ideal ACFs
#
set	acf = acf0	# noiseless ACFs
if 0 then
	if (!(-e $ecalib_noisy) || $arg == calib_noisy) then
		i proj2 t.fld $tpl $twtf
		echo y | op sub t.fld t.fld $tproj	# li_hat - li_true
		echo y | op nonlin exp t.fld t.fld -1 1	# exp(-"")
		op mul $ecalib_noisy $ecalib t.fld
		j $ecalib $ecalib_noisy
	$exit
	endif
	set ecalib = $ecalib_noisy
	set acf = acf1	# noisy ACFs
endif

#
# Noiseless emission FBP
#
set efbpwin = cls3sinc,12,1,1
set efbp0 = $data/efbp0,$efbpwin.fld
if (!(-e $efbp0) || $arg == efbp0) then
	i fbp2e dsc $efbp0 - $eproj - 1 - 0 $edsc $efbpwin -
	j -red $efbp0 $emap
$exit
endif

#
# Noisy emission FBP
# note, this is incorrect since emission smoothing should be applied
# before attenuation correction to avoid resolution mismatch.
#
set efbp = $data/efbp,$acf,$efbpwin$erun.fld
if (!(-e $efbp) || $arg == efbp) then
	i fbp2e dsc $efbp - $eprompt $ecalib 1 - $erandmean $edsc $efbpwin -
	j -red $efbp $emap $efbp0
$exit
endif


#
# Everything here and below is various iterative reconstruction methods
#


#
# OSEM (unregularized)
#
set method = @2@osema,16@-	# 2 iterations, 16 subsets
set osem = $data/osem,$acf,$method$erun.fld
if (!(-e $osem) || $arg == osem) then
	set saver = -
	set init = -
	i empl2 $osem $init $eprompt $ecalib - $erandmean $ewtr - \
		$method $saver 1 1e9 0 -
	j -red $osem $emap $efbp
$exit
endif


# which beta
set eq_l2b = 13

# emission algorithm
#set	ealg	= sage,3,raster1	# recommended 2d algorithm, for now :-)
set	ealg	= osdpc,8;		# depierro regularized osem
set	proot	= quad,1,b2info		# recommended penalty, for now :-)
set	epenal	= $eq_l2b,$proot; unset proot
#set	proot	= huber,2,center,0.5,ih,4	# try this penalty?
#set	epenal	= $en_l2b,$proot; unset proot

#
# emission PL recon
#
set niter	= 10
set method	= @$niter@$ealg@$epenal
set init	= $efbp
set epl		= $data/epl,$acf,$ealg,$epenal,$niter$erun.fld
set slice	= -
if $arg == epl then
	if (-e $epl) then
		echo $epl exists
	else
		set saver = disp,1
		set saver = -
		i empl2 $epl $init $eprompt $ecalib - $erandmean $ewtr - \
			$method $saver 1 1e9 0 $slice
		j -red $epl $emap $efbp
	endif
#	j -a $epl $emap
$exit
endif

#
# make final picture in GIF format for web site / verification
#
if $arg == pic then
	j -0 -nrow 2 $epl $osem $efbp $emap
	op conv sim2d.gif tmp,op,disp.pgm
endif

exit ############# nothing below here has been tested recently ###########
exit ############# nothing below here has been tested recently ###########
exit ############# nothing below here has been tested recently ###########

#
# Noiseless EM
#
	set em0	= t,em0.fld
if $arg == em0 then
	set method = @30@em,1@-
	set saver = stack,1
	if 0 then
		echo y | op div t.fld $trues $calib
		set inputs = "t.fld - - 0e-5"
	else
		echo y | op add t.fld $trues - $randmean
		set inputs = "t.fld $calib - $randmean"
	endif
	i empl2 $em0 - $inputs $ewtf - $method $saver 1 1e9 0 -
	j --red -a $emap $fbp0 $em0
exit
endif

#
# Noiseless OSEM
#
if $arg == osem0 then
	set method = @30@osema,1@-
	set method = @30@osemb,1@-
	set method = @20@osema,15@-
	set method = @120@osemb,12@-
	set osem0 = t,$method.fld
	set saver = stack,5
	set saver = stack,1
	if 0 then
		echo y | op div t.fld $trues $calib
		set inputs = "t.fld - - 0e-5"
	else
		echo y | op add t.fld $trues - $randmean
		set inputs = "t.fld $calib - $randmean"
	endif
	set init = $emap
	set init = $fbp0
	set init = -
	i empl2 $osem0 $init $inputs $ewtr - $method $saver 1 1e9 0 -
	j --red -a $emap $fbp0 $osem0
exit
endif


#
# Noiseless psca
#
	set alg		= psca,1,raster1
	set niter	= 10
	set method	= @$niter@$alg@$penal
	set psca0	= $data/$alg,$penal,$niter,nonoise.fld
if ($arg == psca0) then
	# add randoms to trues
	set prompt = t.fld
	echo y | op add t.fld $trues - $randmean
	set init = -
	set init = $fbp0
	set saver = disp,5
	echo y | \
	i empl2 $psca0 $init $prompt $calib - $randmean $ewtf - \
		$method $saver 1 1e9 0 - \
		| grep dlike
	j -a $sage0 $fbp0 $emap
$exit
endif

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

#
# precorrect data for fbp (or pwls)
#
if (!(-e $meas) || $arg == pre) then
	op pre emis calib $meas $info $prompt - $calib - $randmean 7
	echo y | op mean t.fld $meas
	j -a t.fld $proj
$exit
endif

if $arg == stat then
	set dir=$data
	foreach file ($data/sage*$run.fld $data/fbp*$run.fld)
		set root = $dir/`basename $file .fld`
		set mean = $root,mean.fld
		set std = $root,std.fld
		if !(-e $mean) then
			op stat $mean $std - - $file
			j -a $mean $emap $sage0 $fbp0
		endif
	end
endif

#
# psf of A'A
#
if (!(-e $epsf) || $arg == epsf) then
	i proj2 $epsf - $edsc -p
	echo y | op psfsym $epsf $epsf
	j $epsf
endif

#
# predict variance of one pixel
#
	# diagonal of fisher info: calib^2 / ymean
	set diag = $data/diag,$count,$randpercent.fld
if $arg == var then
	if !(-e $diag) then
		echo y | op add $diag $trues - $randmean
		echo y | op div $diag $calib $diag
		echo y | op mul $diag $calib $diag
		j $diag
	exit
	endif

	echo $penal
	set btype = `echo $penal | sed s/$l2b,quad,'[12],//'`
#	set btype = `echo $btype | sed 's/(.*//'`
	echo btype = $btype

#	var,emis -ata $epsf -btype b2info
	var,roi -ata $epsf -btype $btype \
	-ix 22 -iy 37 -chat 1 \
	-diag $diag -hood 1 -l2b $l2b -show 0 -wtf $ewtf # -root Var
endif

#
# predict variance of roi
#
	set roi = $data/roi1.fld
if $arg == roi then
	echo y | op ellipse $roi $nx $ny -10 -5 10 6 0 1 3
	echo y | op conv $roi $roi byte		# full pixels only
#	echo y | op add t.fld $emap t.fld -9
#	j -s -m3 t.fld $emap

	var,roi -ata $epsf -btype b2info \
	-roi $roi -chat 1 \
	-diag $diag -hood 1 -l2b $l2b -show 0 -wtf $ewtf # -root Var
endif

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