SPECTGPU

examples/python/recon_demo_SPECTGPU.py

@@ -0,0 +1,76 @@
+# Demo of how to use STIR from python to reconstruct some data
+# To run in "normal" Python, you would type the following in the command line
+#  execfile('recon_demo.py')
+# In ipython, you can use
+#  %run recon_demo.py
+
+# Copyright 2012-06-05 - 2013 Kris Thielemans
+# Copyright 2015 University College London
+
+# This file is part of STIR.
+#
+# SPDX-License-Identifier: Apache-2.0
+#
+# See STIR/LICENSE.txt for details
+
+import stir
+import matplotlib.pyplot as pylab
+import os
+
+# go to directory with input files
+os.chdir('../recon_demo')
+
+# initialise reconstruction object
+# we will do this here via a .par file
+recon = stir.OSMAPOSLReconstruction3DFloat('recon_demo_OSEM_SPECTGPU.par')
+# now modify a few settings from in Python for illustration
+recon.set_num_subsets(2)
+num_subiterations = 4
+# set filenames to save subset sensitivities (for illustration purposes)
+poissonobj = recon.get_objective_function()
+poissonobj.set_subsensitivity_filenames('sens_subset%d.hv')
+poissonobj.set_recompute_sensitivity(True)
+
+# get initial image
+target = recon.get_initial_data()#stir.FloatVoxelsOnCartesianGrid.read_from_file('init.hv')
+# we will just fill the whole array with 1 here
+target.fill(1)
+
+# run a few iterations and plot intermediate results
+
+# Switch 'interactive' mode on for pylab.
+# Without it, the python shell will wait after every pylab.show() for you
+# to close the window.
+try:
+    pylab.ion()
+except:
+    print("Enabling interactive-mode for plotting failed. Continuing.")
+
+s = recon.set_up(target)
+if (s.succeeded()):
+    pylab.figure()
+    for i in range(1,num_subiterations+1):
+        print('\n--------------------- Subiteration ', i)
+        recon.set_start_subiteration_num(i)
+        recon.set_num_subiterations(i)
+        s = recon.reconstruct(target)
+        # currently we need to explicitly prevent recomputing sensitivity
+        # when we call reconstruct() again in the next iteration
+        poissonobj.set_recompute_sensitivity(False)
+        # extract to python for plotting
+        npimage = target.as_array()
+        pylab.plot(npimage[10, 30, :])
+        pylab.show()
+
+    # Save the final image
+    recon_filename = f"recon_image_{num_subiterations}.hv"
+    print(f"Saving Image as: {recon_filename}")
+    target.write_to_file(recon_filename)
+
+    # plot slice of final image
+    pylab.figure()
+    pylab.imshow(npimage[10, :, :])
+    # Keep figures open until user closes them
+    pylab.show(block=True)
+else:
+    print('Error setting-up reconstruction object')

examples/python/recon_demo_gradient_SPECTGPU.py

@@ -0,0 +1,66 @@
+# %%
+import stir
+import matplotlib.pyplot as plt
+import os
+
+# go to directory with input files
+os.chdir('../recon_demo')
+
+# %%
+# initialise reconstruction object
+# we will do this here via a .par file
+recon = stir.OSMAPOSLReconstruction3DFloat('recon_demo_OSEM_SPECTGPU.par')
+# now modify a few settings from in Python for illustration
+recon.set_num_subsets(2)
+num_subiterations = 4
+# set filenames to save subset sensitivities (for illustration purposes)
+poissonobj = recon.get_objective_function()
+poissonobj.set_subsensitivity_filenames('sens_subset{}.hv')
+poissonobj.set_recompute_sensitivity(True)
+help(recon)
+
+
+# %%
+# construct image related to the data to reconstruct
+projdata=recon.get_input_data()#recon.tir.ProjData.read_from_file('smalllong.hs');
+# use smaller voxels than the default
+# zoom=2.216842;
+target=recon.get_initial_data() #stir.FloatVoxelsOnCartesianGrid(projdata.get_proj_data_info(), zoom);
+# get initial image
+# help(target)
+# target = stir.FloatVoxelsOnCartesianGrid.read_from_file('init.hv')
+# we will just fill the whole array with 1 here
+# target.fill(1)
+s = recon.set_up(target)
+
+# %%
+
+# compute gradient of objective function
+# create a copy to store the gradient
+gradient=target.get_empty_copy();
+# compute gradient
+subset_num=1;
+poissonobj.compute_sub_gradient(gradient,target,subset_num)
+
+# extract to python for plotting
+npimage = gradient.as_array()
+plt.plot(npimage[10, 30, :])
+plt.show()
+
+# this is useful to find the EM update (i.e. multiply with image)
+poissonobj.compute_sub_gradient_without_penalty_plus_sensitivity(gradient,target,subset_num)
+# extract to python for plotting
+npimage = gradient.as_array()
+plt.plot(npimage[10, 30, :])
+plt.show()
+
+# The followin is just the first iteration you need to create a for
+# loop to run all the iterations and subsets
+EMupdate = target*gradient
+# extract to python for plotting
+npimage = EMupdate.as_array()
+plt.plot(npimage[10, 30, :])
+plt.show()
+
+
+

examples/recon_demo/recon_demo_OSEM_SPECTGPU.par

@@ -0,0 +1,37 @@
+OSMAPOSLParameters :=
+
+objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+input file := ../../recon_test_pack/SPECT/SPECTUB/input.hs
+projector pair type := SPECTGPU
+              Projector Pair Using SPECTGPU Parameters :=
+                  ; width of PSF
+
+                  ;PSF type of correction {  3D  }
+                  ;
+                  ; These values are mostly dependent on your collimator.
+                  ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
+                  ; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+                  collimator slope := 0.0163
+                  collimator sigma 0(cm) := 0.1466
+                  number of PSF sigmas := 2 ; default 2 (95% of the PSF)
+
+
+                  attenuation image filename := ../../recon_test_pack/SPECT/SPECTUB/attMapRec.hv
+
+         End Projector Pair Using SPECTGPU Parameters :=
+
+end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+initial estimate := ../../recon_test_pack/SPECT/SPECTUB/init.hv
+output filename prefix := my_SPECTGPU_image
+
+number of subsets:= 4
+number of subiterations:= 16
+save estimates at subiteration intervals:= 4
+
+
+
+
+END :=

examples/samples/OSMAPOSL_osem_SPECTGPU.par

@@ -0,0 +1,46 @@
+OSMAPOSLParameters :=
+  ; sample .par file to use OSEM on SPECT data. 
+  ; Any of the algorithm parameters illustrated for PET (such as filtering, prior etc)
+  ; will work for SPECT as well, as would OSSPS.
+  ; The only thing different here is the projector.
+
+  objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+
+  PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+    input file := peak.hs
+
+    projector pair type := SPECTGPU
+      Projector Pair Using SPECTGPU Parameters :=
+            ; width of PSF
+
+            ;PSF type of correction {  3D  }
+            ; next 2 parameters define the PSF. They are ignored if psf_type is "Geometrical"
+            ; These values are mostly dependent on your collimator.
+            ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
+            ; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+            collimator slope := 0.0163
+            collimator sigma 0(cm) := 0.1466
+            number of PSF sigmas := 2 ; default 2 (95% of the PSF)
+
+            attenuation image filename := umap_zoomed.hv
+
+     End Projector Pair Using SPECTGPU Parameters :=
+  subset sensitivity filenames := sens_%d
+  recompute sensitivity := 1
+  z zoom := 0.5
+  xy output image size (in pixels) := 128
+  z output image size (in pixels) := 128
+  end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+  ; best to specifiy an initial image (e.g. filled with 1) for sizes
+  ; or see OSMAPOSL_osem_with_interfiltering.par and FBP2D_SPECT.par for some keywords 
+  ;initial estimate:= init.hv
+  output filename prefix := outgpu
+
+  ; needs to be a divisor of the number of views
+  number of subsets:= 8
+  number of subiterations:= 16
+  Save estimates at subiteration intervals:= 16
+
+END :=

examples/samples/OSMAPOSL_osem_SPECTUB.par

@@ -8,17 +8,17 @@ OSMAPOSLParameters :=
 
   PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
 
-    input file := input.hs
+    input file := peak.hs
 
     projector pair type := Matrix
       Projector Pair Using Matrix Parameters :=
         Matrix type := SPECT UB
         Projection Matrix By Bin SPECT UB Parameters:=
             ; width of PSF
-            maximum number of sigmas:= 2.0
+            maximum number of sigmas:= 1.3
 
             ;PSF type of correction { 2D // 3D // Geometrical }
-            psf type:= 2D
+            psf type:= Geometrical
             ; next 2 parameters define the PSF. They are ignored if psf_type is "Geometrical"
             ; These values are mostly dependent on your collimator.
             ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
@@ -29,25 +29,30 @@ OSMAPOSLParameters :=
             ;Attenuation correction { Simple // Full // No }
             attenuation type := Simple      
             ;Values in attenuation map in cm-1
-            attenuation map := attMapRec.hv
+            attenuation map := umap_zoomed.hv
 
             ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No}
-            mask type := Explicit Mask
-            mask file := mask.hv
+            ;mask type := Explicit Mask
+            ;mask file := mask.hv
 
            ; if next variable is set to 0, only a single view is kept in memory
            keep all views in cache:=1
 
         End Projection Matrix By Bin SPECT UB Parameters:=
 
      End Projector Pair Using Matrix Parameters :=
+     subset sensitivity filenames := UBsens_%d
+  recompute sensitivity := 1
+  z zoom := 0.5
+  xy output image size (in pixels) := 128
+  z output image size (in pixels) := 128
 
   end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
 
   ; best to specifiy an initial image (e.g. filled with 1) for sizes
   ; or see OSMAPOSL_osem_with_interfiltering.par and FBP2D_SPECT.par for some keywords 
-  initial estimate:= init.hv
-  output filename prefix := out
+  ;initial estimate:= init.hv
+  output filename prefix := outUB
 
   ; needs to be a divisor of the number of views
   number of subsets:= 8

examples/samples/backward_project_SPECTGPU.par

@@ -0,0 +1,9 @@
+Back Projector parameters:=
+   type := SPECTGPU
+   Back projector Using SPECTGPU Parameters :=
+      attenuation image filename := umap_inverted.hv
+      collimator slope := 0.01
+      collimator sigma 0(cm) := 0.07
+      number of PSF sigmas := 2 ; default 2 (95% of the PSF)
+      End Back Projector Using SPECTGPU Parameters :=
+End:=

examples/samples/forward_project_SPECTGPU.par

@@ -0,0 +1,9 @@
+Forward Projector parameters:=
+   type := SPECTGPU
+   Forward projector Using SPECTGPU Parameters :=
+      attenuation image filename := umap_inverted.hv
+      collimator slope := 0.01
+      collimator sigma 0(cm) := 0.07
+      number of PSF sigmas := 2 ; default 2 (95% of the PSF)
+      End Forward Projector Using SPECTGPU Parameters :=
+End:=

examples/samples/projector_pair_SPECTGPU.par

@@ -0,0 +1,19 @@
+ projector pair parameters:=
+   projector pair type := SPECTGPU
+     Projector Pair Using SPECTGPU Parameters :=
+                  ; width of PSF
+
+                  ;PSF type of correction {  3D  }
+                  ;
+                  ; These values are mostly dependent on your collimator.
+                  ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
+                  ; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+                  ;collimator slope := 0.0163
+                  ;collimator sigma 0(cm) := 0.1466
+                  ;number of PSF sigmas := 2 ; default 2 (95% of the PSF)
+
+
+                  ;attenuation image filename := ${CMAKE_SOURCE_DIR}/recon_test_pack/SPECT/SPECTUB/attMapRec.hv
+
+         End Projector Pair Using SPECTGPU Parameters :=
+End projector pair parameters:=   

examples/samples/projector_pair_SPECTUB.par

@@ -0,0 +1,32 @@
+ projector pair parameters:=
+    projector pair type:= Matrix
+       Projector Pair Using Matrix Parameters :=
+          Matrix type := SPECT UB
+            Projection Matrix By Bin SPECT UB Parameters:=
+              ; width of PSF
+              ;maximum number of sigmas:= 1.3
+
+              ;PSF type of correction { 2D // 3D // Geometrical }
+              psf type:= Geometrical
+              ; next 2 parameters define the PSF. They are ignored if psf_type is "Geometrical"
+              ; These values are mostly dependent on your collimator.
+              ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
+              ; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+              ;collimator slope := 0.0163
+              ;collimator sigma 0(cm) := 0.1466
+
+              ;Attenuation correction { Simple // Full // No }
+              ;attenuation type := Simple      
+              ;Values in attenuation map in cm-1
+              ;attenuation map := umap_zoomed.hv
+
+              ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No}
+              ;mask type := Explicit Mask
+              ;mask file := mask.hv
+
+           ; if next variable is set to 0, only a single view is kept in memory
+           keep all views in cache:=1
+
+        End Projection Matrix By Bin SPECT UB Parameters:=
+    End Projector Pair Using Matrix Parameters :=
+End projector pair parameters:=   

recon_test_pack/SPECT/SPECTGPU/FBP2D.par

@@ -0,0 +1,21 @@
+FBP2DParameters :=
+output filename prefix := out/FBP
+input file := ../SPECTUB/input.hs
+Z output image size (in pixels):=64
+z zoom:=.5
+xy output image size (in pixels) := 129
+Back Projector type:= matrix
+  Back Projector Using Matrix Parameters :=
+	Matrix type := SPECT UB
+	Projection Matrix By Bin SPECT UB Parameters:=
+		psf type:= Geometrical
+		attenuation type := No	
+		mask type := No
+                ;mask file := ../SPECTUB/mask.hv
+
+                ; no point in keeping things in memory as FBP runs the projector only once
+                keep all views in cache := 0
+
+	End Projection Matrix By Bin SPECT UB Parameters:=
+  End Back Projector Using Matrix Parameters :=
+end:=