adding a Python demo for EM step with gradient (#1620)

* adding a demo for EM step with gradient

* adding reconddemo to Actions workflow

* disable appveyor test as it seems to hang (maybe because of plotting)

---------

Co-authored-by: Kris Thielemans <[email protected]>

Author: danieldeidda

.appveyor.yml

@@ -68,6 +68,10 @@ test_script:
   - cd ..\src
   - "set PYTHONPATH=C:\\projects\\stir\\install\\python"
   - python -m pytest .
+  - echo "Running Python examples"
+  - cd ..\examples\python
+  - rem Disabled as fails, probably because of plotting
+  - rem python recon_demo_gradient.py
   - conda deactivate
 
 on_finish:

.github/workflows/build-test.yml

@@ -539,5 +539,7 @@ jobs:
               ;;
             (*)
               python -m pytest .
+              cd ../examples/python
+              python recon_demo_gradient.py
               ;;
           esac

examples/python/recon_demo_gradient.py

@@ -0,0 +1,66 @@
+# %%
+import stir
+import stirextra
+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.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)
+
+
+# %%
+# construct image related to the data to reconstruct
+projdata=stir.ProjData.read_from_file('smalllong.hs');
+# use smaller voxels than the default
+zoom=2.216842;
+target=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 = stirextra.to_numpy(gradient)
+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 = stirextra.to_numpy(gradient)
+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 = stirextra.to_numpy(EMupdate)
+plt.plot(npimage[10, 30, :])
+plt.show()
+
+
+