use different min/max thresholds for scatter scale for mMR example [c…

…i skip]

Fixes UCL#1163

documentation/release_5.2.htm

@@ -127,7 +127,10 @@ <h3>Changed functionality</h3>
   <br /><a href="https://github.com/UCL/STIR/pull/1243/">PR #1243</a>.
   </li>
   <li>The <code>Succeeded</code> class has a new method <code>bool succeeded()</code> enabling more concise code (avoiding the need for comparing with <code>Succeeded::yes</code> which is especially verbose in Python).
-    </li>
+  </li>
+  <li>The example files for the Siemens mMR now use lower min/max thresholds for the (single) scatter scale. This gives better results, see <a href="https://github.com/UCL/STIR/issues/1163/">Issue #1163</a>.
+  <br /><a href="https://github.com/UCL/STIR/pull/1279/">PR #1279</a>.
+  </li>
 </ul>
 
 <h3>Deprecated functionality</h3>

examples/Siemens-mMR/scatter_and_recon.sh

@@ -75,7 +75,7 @@ echo "Estimating scatter (be patient). Log saved in output/scatter.log"
 # filename-prefix for additive sino (i.e. "precorrected" sum of scatter and randoms)
 total_additive_prefix=output/total_additive
 num_scat_iters=3
-scatter_pardir=${pardir}/../samples/scatter_estimation_par_files
+scatter_pardir=${pardir}/scatter_estimation_par_files
 # you might have to change this for a different scanner than the mMR
 scatter_recon_num_subiterations=21
 scatter_recon_num_subsets=21

examples/Siemens-mMR/scatter_estimation_par_files/README.md

@@ -0,0 +1,19 @@
+# Example files for running scatter estimation for the Siemens mMR
+
+Files made by Nikos Efthimou and fine-tuned by Kris Thielemans.<br>
+Copyright University of Hull 2018-2019<br>
+copyright University College London 2016, 2020<br>
+Distributed under the Apache 2.0 License
+
+These files are almost identical to those in
+[examples/samples/scatter_estimation_par_files/](../../samples/scatter_estimation_par_files/README.md),
+see there for some more information.
+
+Currently the only difference are the lower values for
+```
+maximum scatter scaling factor := .5
+minimum scatter scaling factor := 0.1
+```
+
+These have been shown to work better for mMR data, see e.g.
+[STIR issue #1163](https://github.com/UCL/STIR/issues/1163).

examples/Siemens-mMR/scatter_estimation_par_files/postfilter_Gaussian_for_mask.par

@@ -0,0 +1,12 @@
+PostFilteringParameters :=
+  Postfilter type :=  Separable Gaussian
+Separable Gaussian Filter Parameters :=
+x-dir filter FWHM (in mm):= 20
+y-dir filter FWHM (in mm):= 20
+z-dir filter FWHM (in mm):= 15
+; optionally restrict kernel sizes
+; x-dir maximum kernel size := 129
+; y-dir maximum kernel size := 129
+; z-dir maximum kernel size := 31
+END Separable Gaussian Filter Parameters :=
+End PostFiltering Parameters:=

examples/Siemens-mMR/scatter_estimation_par_files/run_reconstruction.par

@@ -0,0 +1,44 @@
+Reconstruction Parameters :=
+reconstruction type := OSMAPOSL
+OSMAPOSLParameters :=
+
+objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+maximum absolute segment number to process := -1
+
+projector pair type := Matrix
+  Projector Pair Using Matrix Parameters :=
+  Matrix type := Ray Tracing
+  Ray tracing matrix parameters :=
+   number of rays in tangential direction to trace for each bin:= 5
+  End Ray tracing matrix parameters :=
+  End Projector Pair Using Matrix Parameters :=
+
+;recompute sensitivity := 0
+;subset sensitivity filenames := scatter_subset_sens_%d.hv
+
+; reconstruct at large voxel size to save time
+zoom := 0.2
+
+end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+; initial estimate := 
+enforce initial positivity condition:=1
+
+number of subsets:= ${scatter_recon_num_subsets}
+number of subiterations:=${scatter_recon_num_subiterations}
+;save estimates at subiteration intervals:= ${scatter_recon_num_subiterations}
+
+; smooth a bit as we use a down-sampled scanner (during the scatter estimation resolution can be low)
+post-filter type :=  Separable Gaussian
+Separable Gaussian Filter Parameters :=
+  x-dir filter FWHM (in mm):= 15
+  y-dir filter FWHM (in mm):= 15
+  z-dir filter FWHM (in mm):= 15
+END Separable Gaussian Filter Parameters :=
+;
+; Disable output
+;
+disable output := 1
+End OSMAPOSLParameters:=
+End reconstruction Parameters:=

examples/Siemens-mMR/scatter_estimation_par_files/scatter_estimation.par

@@ -0,0 +1,86 @@
+Scatter Estimation Parameters :=
+
+;Run in debug mode
+;; A new folder called extras will be created, in which many 
+;; extra files will be stored 
+run in debug mode := 1
+
+; Measured data
+input file := ${sino_input}
+
+; Attenuation Image
+attenuation image filename := ${atnimg}
+
+; Normalisation coefficients & attenuation data
+Normalisation type := from ProjData
+   Bin Normalisation From ProjData :=
+            normalisation projdata filename:= ${NORM}
+   End Bin Normalisation From ProjData:=
+
+attenuation correction factors filename := ${acf3d}
+
+;; Background data (not  normalised).
+; Should be set to the randoms estimate (unless you precorrected, but we haven't tested that)
+background projdata filename := ${randoms3d}
+
+; Mask for tail-fitting
+; It will be computed by masking the attenuation image, and forward projecting that.
+; If !recompute mask projdata then the filename must be set.
+recompute mask projdata := 1
+mask projdata filename := ${mask_projdata_filename}
+
+; Input or output filename - depends on recompute  
+recompute mask image := 1
+mask image filename := ${mask_image}
+; threshold to be applied after filtering (in cm^-1). Default value is below
+mask attenuation image min threshold  := 0.003
+; optional filename to specify a filter before thresholding the attenuation image
+; By default a Gaussian filter with FWHM (15,20,20) will be used. Here we use an explicit file as an example.
+mask attenuation image filter filename := ${scatter_pardir}/postfilter_Gaussian_for_mask.par
+;End of Mask
+
+;Parameter file for the tail fitting of the scatter data (within the mask)
+tail fitting parameter filename := ${scatter_pardir}/tail_fitting.par
+
+; Run simulation and reconstruction in 2D and export SSRB sinograms (currently required)
+run in 2d projdata := 1
+
+; ScatterSimulation parameters 
+; could read from a file, but instead we have them below
+; scatter simulation parameter filename := ${scatter_pardir}/scatter_simulation.par
+Scatter Simulation type := PET Single Scatter Simulation
+  PET Single Scatter Simulation Parameters :=
+    ; could change some parameters here if you need to (not recommended)
+  End PET Single Scatter Simulation Parameters:=
+
+; next option is the default
+use scanner downsampling in scatter simulation := 1
+
+; could add parameters below, but reading it from file
+; reconstruction type := ...
+reconstruction parameter filename := ${scatter_pardir}/run_reconstruction.par 
+
+;
+; This is the number of times which the Scatter Estimation will
+; iterate. Default is 5
+
+number of scatter iterations := ${num_scat_iters}
+
+; Average the first two activity images
+do average at 2 := 1
+
+; Export scatter estimates of each iteration 
+export scatter estimates of each iteration := 1
+
+output scatter estimate name prefix := ${scatter_prefix}
+output additive estimate name prefix:= ${total_additive_prefix}
+
+maximum scatter scaling factor := 0.4
+minimum scatter scaling factor := 0.1
+
+;Upsample and fit 
+; defaults to 3.
+upsampling half filter width := 3
+remove interleaving before upsampling := 1
+
+End Scatter Estimation Parameters :=

examples/Siemens-mMR/scatter_estimation_par_files/tail_fitting.par

@@ -0,0 +1,4 @@
+CreateTailMaskFromACFs :=
+        ACF-threshold := 1.1
+        safety-margin := 4
+END CreateTailMaskFromACFs :=