Merge branch 'master' into reinvert_SPECT_axis

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

@@ -111,20 +111,30 @@ jobs:
           compiler: gcc
           # compiler_version: 11
           cuda_version: "0"
+          # openmp=OFF as ON needs special stuff for Apple gcc which KT doesn't feel like getting into...
           BUILD_FLAGS: "-DSTIR_OPENMP=OFF"
           parallelproj: "OFF"
           BUILD_TYPE: "Debug"
           ROOT: "OFF"
           ITK: "OFF"
         - os: macOS-latest
           compiler: clang
-          compiler_version: 16
+          compiler_version: 18
           cuda_version: "0"
           BUILD_FLAGS: "-DSTIR_OPENMP=OFF"
           parallelproj: "OFF"
           BUILD_TYPE: "Release"
           ROOT: "OFF"
           ITK: "OFF"
+        - os: macOS-latest
+          compiler: clang
+          compiler_version: 21
+          cuda_version: "0"
+          BUILD_FLAGS: "-DSTIR_OPENMP=ON"
+          parallelproj: "OFF"
+          BUILD_TYPE: "Release"
+          ROOT: "OFF"
+          ITK: "OFF"
 
       # let's run all of them, as opposed to aborting when one fails
       fail-fast: false
@@ -191,7 +201,8 @@ jobs:
                 HOMEBREW_PREFIX=/usr/local
               fi
               LLVMDIR=$HOMEBREW_PREFIX/opt/llvm@${{ matrix.compiler_version }}
-              LDFLAGS="-L${LLVMDIR}/lib/c++ -Wl,-rpath,${LLVMDIR}/lib/c++"
+              # https://github.com/llvm/llvm-project/issues/94292#issuecomment-2149335628 says to remove this flag
+              # LDFLAGS="-L${LLVMDIR}/lib/c++ -Wl,-rpath,${LLVMDIR}/lib/c++"
               export SDKROOT=$(xcrun --show-sdk-path)
               echo SDKROOT="$SDKROOT"
               # make available to jobs below
@@ -230,7 +241,7 @@ jobs:
                   sudo apt install $CXX
               fi
               # other dependencies
-              sudo apt install libboost-dev libhdf5-serial-dev swig python3-dev nlohmann-json3-dev
+              sudo apt install libboost-dev libhdf5-serial-dev swig python3-dev nlohmann-json3-dev ninja-build
               if test "${{matrix.ITK}}XX" == "ONXX"; then
                 sudo apt install libinsighttoolkit5-dev
               fi
@@ -239,7 +250,16 @@ jobs:
               PYTHON_EXECUTABLE=$(which python3)
               ;;
             (macOS*)
-              brew install boost swig python nlohmann-json
+              if ! command -v ninja > /dev/null; then
+                  brew install ninja
+              fi
+              if ! command -v python3 > /dev/null; then
+                  brew install python
+              fi
+              if ! command -v swig > /dev/null; then
+                  brew install swig
+              fi
+              brew install boost nlohmann-json
               PYTHON_EXECUTABLE=$(which python3)
               ;;
             (windows*)
@@ -350,7 +370,7 @@ jobs:
           echo "cmake flags $BUILD_FLAGS $EXTRA_BUILD_FLAGS"          
           mkdir build
           cd build
-          cmake -S .. ${BUILD_FLAGS} ${EXTRA_BUILD_FLAGS}
+          cmake -G Ninja -S .. ${BUILD_FLAGS} ${EXTRA_BUILD_FLAGS}
           # report disk space
           case ${{matrix.os}} in
             (ubuntu* | macOS*)
@@ -375,7 +395,7 @@ jobs:
       run: |
           cd ${GITHUB_WORKSPACE}/build;
           source ${GITHUB_WORKSPACE}/my-env/bin/activate
-          cmake --build . -j 2 --config ${BUILD_TYPE}} --target install
+          cmake --build . --config ${BUILD_TYPE}} --target install
 
     - name: ctest
       shell: bash
@@ -519,5 +539,7 @@ jobs:
               ;;
             (*)
               python -m pytest .
+              cd ../examples/python
+              python recon_demo_gradient.py
               ;;
           esac

documentation/release_6.3.htm

@@ -19,6 +19,7 @@ <h2> Summary for end users (also to be read by developers)</h2>
 
 
   <h3>New functionality</h3>
+  <h4>General</h4>
   <ul>
     <li>
       The analytic Spline Reconstruction Technique (SRT) algorithm has been added in 2 different versions: for PET
@@ -51,9 +52,6 @@ <h3>New functionality</h3>
       However, projection data is currently still always returned as non-TOF (but list-mode data is read as TOF).<br>
       <a href=https://github.com/UCL/STIR/pull/1503>PR #1503</a>
     </li>
-    <li>
-      <tt>stir_timings</tt> has now an extra option to parse a par-file for a projector-pair.
-    </li>
     <li>
       Added the ability to set a forward projector for mask projection in the <code>ScatterEstimation</code> class.<br>
       <a href=https://github.com/UCL/STIR/pull/1530>PR #1530</a>
@@ -64,18 +62,42 @@ <h3>New functionality</h3>
       <a href=https://github.com/UCL/STIR/pull/1519>PR #1519</a>
     </li>
     <li>
-      <code>LmToProjData</code>/<tt>lm_to_projdata</tt> now no longer requires a template projection data. If none is
+      <code>LmToProjData</code> class/<tt>lm_to_projdata</tt> utility now no longer requires a template projection data. If none is
       specified, it will use the proj_data_info from the input list-mode.<br>
       <b>Warning</b> for some scanners with TOF capabilities, this will result in very large projection data (possibly
       larger than the default from the vendor).<br>
       <a href=https://github.com/UCL/STIR/pull/1315>PR #1315</a>
     </li>
+  </ul>
+  <h4>Python</h4>
+  <ul>
+    <li>
+      "Container" classes such as <code>FloatVoxelsOnCartesianGrid</code>, <code>ProjDataInMemory</code>
+      and array-classes now have numerical operations properly defined, for instance
+      <code>a = b + c</code> and <code>a -= 3</code>. Note that <code>a = 1 + b</code> is not
+      yet available.<br>
+      <a href=https://github.com/UCL/STIR/pull/1630>PR #1630</a>
+    </li>
+    <li>
+      The above "container" classes now have an extra member `as_array()` which returns a numpy <code>ndarray</code>. This
+      is equivalent to `stirextra.to_numpy()` which will become deprecated later. In addition, the
+      <code>fill()</code> method now directly accepts an <code>ndarray</code>, avoiding the need to go via an iterator.
+      These additions also make it easier to prt SIRF python code to STIR.<br>
+      <a href=https://github.com/UCL/STIR/pull/1632>PR #1632</a>
+    </li>
     <li>
-      Added a Python script to convert e7tools generated Siemens Biograph Vision 600 sinograms to STIR compatible format.
+      Added a Python script to convert e7tools generated Siemens Biograph Vision 600 sinograms to STIR compatible format.<br>
       <a href=https://github.com/UCL/STIR/pull/1593>PR #1593</a>
     </li>
   </ul>
 
+  <h4>Utilities</h4>
+  <ul>
+    <li>
+      <tt>stir_timings</tt> has now an extra option to parse a par-file for a projector-pair.
+    </li>
+  </ul>
+
   <h3>Changed functionality</h3>
   <ul>
     <li>
@@ -87,6 +109,14 @@ <h3>Changed functionality</h3>
     <li>
       Default ECAT scanner configurations updated to use a negative intrinsic tilt.
     </li>
+    <li>
+      When computing the sensitivity images in <code>PoissonLogLikelihoodWithLinearModelForMeanAndProjData</code>,
+      we now avoid creating an extra <code>ProjDataInMemory</code> in most cases (there are still some corner
+      cases for TOF data when using non-TOF projector for the sensitivity, but then the memory overhead is small).
+      This enables LAFOV PET reconstructions with the "ray tracing" matrix (we use less memory, and we avoid running into a current
+      <code>ProjDataInMemory</code> <a href=https://github.com/UCL/STIR/issues/1505> limitation on the number of bins</a>).<br>
+      <a href=https://github.com/UCL/STIR/pull/1617>PR #1716</a>
+      </li>
       <li>
         Boost format was replaced by `std::format` for formatting strings. If C++20 or newer is
         not yet used, a work-around is in place by using the [{fmt} library](https://github.com/fmtlib/fmt/)
@@ -154,6 +184,12 @@ <h3>New functionality</h3>
 
   <h3>Changed functionality</h3>
   <ul>
+    <li>
+      In C++, "Container" classes such as <code>VoxelsOnCartesianGrid</code>, <code>ProjDataInMemory</code>
+      and array-classes now have numerical operations that return objects of the correct type.
+      (Previously, it was a base-class such as <code>NumericVectorWithOffset<Array<2, float>, float></code>).<br>
+      <a href=https://github.com/UCL/STIR/pull/1630>PR #1630</a>
+    </li>
     <li>
       Made 2 (deprecated) members of <code>ProjDataInfoBlocksOnCylindricalNoArcCorr</code> <code>private</code>
       and do some clean-up of which files to include.
@@ -181,12 +217,12 @@ <h3>Changed functionality</h3>
 
   <h3>Bug fixes</h3>
   <ul>
-    <li>Fixed minor incompatibility with gcc-14 and clang-18 buy adding an extra include file<br>
+    <li>Fixed minor incompatibility with gcc-14 and clang-18 by adding an extra include file.<br>
       <a href=https://github.com/UCL/STIR/pull/1552>PR #1552</a>
     </li>
   </ul>
 
-  <h3>Deprecations</h3>
+  <h3>New deprecations</h3>
   <ul>
     <li>
       <code>truncate_end_planes</code> will be removed in v7.0
@@ -197,7 +233,13 @@ <h3>Other code changes</h3>
 
 
   <h3>Test changes</h3>
-
+  <h4>Python tests</h4>
+  <ul>
+    <li>
+      Tests for numerical operations for the "container" classes were added in the same PR as the code.<br>
+      <a href=https://github.com/UCL/STIR/pull/1630>PR #1630</a>
+    </li>
+  </ul>
 
   <h4>C++ tests</h4>
 

examples/python/recon_demo_gradient.py

@@ -1,18 +1,11 @@
 # %%
 import stir
 import stirextra
-import matplotlib.pyplot as pylab
+import matplotlib.pyplot as plt
 import os
 
 # go to directory with input files
 os.chdir('../recon_demo')
-
-# define a function due to current limitation on Swig. " multiply_images(a,b)"
-def multiply_images(a, b):
-    a.xapyb(a,b,a.get_empty_copy(),a)
-    return a
-
-
 # %%
 # initialise reconstruction object
 # we will do this here via a .par file
@@ -22,7 +15,7 @@ def multiply_images(a, b):
 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_subsensitivity_filenames('sens_subset{}.hv')
 poissonobj.set_recompute_sensitivity(True)
 
 
@@ -40,10 +33,6 @@ def multiply_images(a, b):
 s = recon.set_up(target)
 
 # %%
-
-# compute gradient of objective function 
-# put some data in the image
-target.fill(1);
 # create a copy to store the gradient
 gradient=target.get_empty_copy();
 # compute gradient
@@ -52,22 +41,24 @@ def multiply_images(a, b):
 
 # extract to python for plotting
 npimage = stirextra.to_numpy(gradient)
-pylab.plot(npimage[10, 30, :])
-pylab.show()
+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)
-pylab.plot(npimage[10, 30, :])
-pylab.show()
 
-# The followin is just the first iteration you need to create a for loop to run all the iterations and subsets
-EMupdate = multiply_images(target,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)
-pylab.plot(npimage[10, 30, :])
-pylab.show()
+plt.plot(npimage[10, 30, :])
+plt.show()
 
 
 

external_helpers/fmt

@@ -137,7 +137,7 @@ if(STIR_OPENMP)
   find_package(OpenMP REQUIRED)  
   add_definitions(${OpenMP_CXX_FLAGS})
 
-  if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+  if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" AND NOT APPLE)
     # work around https://gitlab.kitware.com/cmake/cmake/-/issues/26037
     set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_CXX -latomic)
     message(STATUS "OpenMP Linker flags for Clang: ${OpenMP_EXE_LINKER_FLAGS}")

src/CMakeLists.txt

@@ -19,7 +19,7 @@
     See STIR/LICENSE.txt for details
 */
 #include "stir/WienerArrayFilter2D.h"
-
+#include <vector>
 #include <algorithm>
 
 using std::nth_element;