Merge branch 'UCL:master' into SRT2D

CMakeLists.txt

@@ -54,9 +54,9 @@ endif()
 
 ####### Set Version number etc
 set(VERSION_MAJOR   6)
-set(VERSION_MINOR   0)
+set(VERSION_MINOR   1)
 set(VERSION_PATCH   0)
-set(VERSION   060000) # only used in STIRConfig.h.in and swig/CMakeLists.txt
+set(VERSION   060100) # only used in STIRConfig.h.in and swig/CMakeLists.txt
 
 set(STIR_VERSION
   ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -121,7 +121,13 @@ if(NOT DISABLE_CERN_ROOT)
   find_package(CERN_ROOT)
   if (CERN_ROOT_FOUND)
     message(STATUS "ROOT Version: ${CERN_ROOT_VERSION}")
-    if (${CERN_ROOT_VERSION} VERSION_GREATER 6.23.99)
+    if (${CERN_ROOT_VERSION} VERSION_GREATER 6.29.99)
+        message(STATUS "ROOT Version is >= 6.30. Setting the minimum CXX version to 20.")
+        UseCXX(20)
+    elseif (${CERN_ROOT_VERSION} VERSION_GREATER 6.27.99)
+        message(STATUS "ROOT Version is >= 6.28. Setting the minimum CXX version to 17.")
+        UseCXX(17)
+    elseif (${CERN_ROOT_VERSION} VERSION_GREATER 6.23.99)
         message(STATUS "ROOT Version is >= 6.24. Setting the minimum CXX version to 14.")
         UseCXX(14)
     endif()

VERSION.txt

@@ -1 +1 @@
-6.0.0
+6.1.0

documentation/STIR-UsersGuide.tex

@@ -44,7 +44,7 @@
 \\[3cm]
 
 \textbf{{\huge User's Guide\\
- Version 6.0}}
+ Version 6.1}}
 \end{center}
 
 \end{spacing}
@@ -444,7 +444,8 @@ \subsubsection{
 
 { \subsubsubsection{Cygwin on Windows}
 }
-\textit{Cygwin support has not been tested since about 2018 but likely still works. Instructions below might be out-of-date though.}
+\textit{Cygwin support has not been tested since about 2018 but likely still works. Instructions below are likely out-of-date though.
+We highly recommend to use WSL instead.}
 
 If you are using Windows but would like 
 to have nearly everything that Linux/Unix has to offer, Cygwin could help. Check out http://cygwin.com.

documentation/STIR-developers-overview.tex

@@ -60,8 +60,8 @@ \section{
 doxygen on the source files that come in the STIR distribution). 
 
 \subsection{Language support}
-STIR is written in C++ and currently requires C++-11, but it is compatible with newer versions of the C++ standard.
-We will enforce C++-14 from STIR 6.0.
+STIR is written in C++ and currently requires C++-14, but it is compatible with newer versions of the C++ standard.
+We will enforce C++-16 from STIR 6.2.
 
 Python and MATLAB support is provided via \R2Lurl{http://www.swig.org/}{SWIG}. This means that Python/MATLAB interfaces follow
 the C++ classes closely, although some differences are required as these languages do not support templates for instance.
@@ -170,8 +170,8 @@ \subsection{
 view.
 \end{itemize} 
 
-Note that in STIR version 6.0, Time-of-Flight (TOF) will be supported. This introduces another
-index. However, \texttt{Sinogram} and \texttt{Viewgram} will remain 2D objects, and \texttt{Segment*} 3D.
+Note that since STIR version 6.0, Time-of-Flight (TOF) is supported. This introduces another
+index. However, \texttt{Sinogram} and \texttt{Viewgram} remain 2D objects, and \texttt{Segment*} 3D.
 This will also be the case once we have layers and energy windows.
 In STIR 5.2, we have therefore introduced new classes
 \texttt{SinogramIndices}, \texttt{ViewgramIndices}\footnote{Replacing \texttt{ViewSegmentNumbers} in previous versions of STIR.}
@@ -897,7 +897,7 @@ \subsubsection{
 \end{itemize}
 
 
-At the moment, we have three derived classes:
+At the moment, we have the following derived classes:
 \begin{itemize}
 \item 
 \texttt{ProjMatrixByBinUsingRayTracing} uses essentially the same 
@@ -916,6 +916,10 @@ \subsubsection{
 part of the projection matrix needs to be stored. Our current 
 implementation does not yet provide a very compact format for 
 storing the elements (although they are of course stored sparsely).
+\item
+\texttt{ProjMatrixByBinSPECTUB} is for parallel hole SPECT
+\item
+\texttt{ProjMatrixByBinPinholeSPECTUB} is for multi-pinhole SPECT
 \end{itemize}
 
 \begin{figure}[htbp]
@@ -956,8 +960,10 @@ \subsection{Objective functions}
 that is computed is generally not the gradient of the
 log-likelihood that corresponds to the forward projector.
 However, one hopes that it still points towards the optimum.
-The corresponding objective function is implemented in the class\\
-\texttt{PoissonLogLikelihoodWithLinearModelForMeanAndProjData}.
+The corresponding objective function is implemented in the classes\\
+\texttt{PoissonLogLikelihoodWithLinearModelForMeanAndProjData} for
+projection data and \texttt{PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin}
+for list-mode data. There are classes for dynamic and gated data as well.
 
 There can be different objective function that use common operations.
 For instance, the objective function could implement a least squares

documentation/release_6.1.htm

@@ -7,27 +7,30 @@
   <body>
     <h1>Summary of changes in STIR release 6.1</h1>
 
-<p>This version is 95% backwards compatible with STIR 6.0 for the user (see below).
-Developers might need to make code changes as 
-detailed below.
-</p>
-<h2>Overall summary</h2>
-<p>
-</p>
+    <p>
+      This version is 100% backwards compatible with STIR 6.0 for the user, except for the bug-fix in the RDP (see below).
+    </p>
 
+<h2>Overall summary</h2>
+    <p>
+      This version adds capability of using Parallelproj (CPU and GPU versions) for TOF data. In addition,
+      the list-mode objective function has several improvements, including speed-up by using multi-threading
+      if caching was not enabled.
+    </p>
 
-<p>Of course, there is also the usual code-cleanup and 
-improvements to the documentation.
-</p>
+    <p>
+      Of course, there is also the usual code-cleanup and improvements to the documentation.
+    </p>
 
-<p>This release contains mainly code written by Kris Thielemans (UCL).
-</p>
+    <p>
+      This release contains mainly code written by Nicole Jurjew (UCL) and Kris Thielemans (UCL).
+    </p>
 
 <h2>Patch release info</h2>
 <ul>
   <li>
-    6.1.0 released ?/?/2020<br>
-    <a href="https://github.com/UCL/STIR/milestone/4">GitHub Milestone 6.0</a>
+    6.1.0 released 15/05/2024<br>
+    <a href="https://github.com/UCL/STIR/milestone/10">GitHub Milestone 6.1</a>
   </li>
   <!--
       <li> 4.0.1 released 28/04/2020
@@ -39,22 +42,14 @@ <h2>Patch release info</h2>
 
 <h2> Summary for end users (also to be read by developers)</h2>
 
-<h3>Changes breaking backwards compatibility from a user-perspective</h3>
-<h4>General</h4>
-<h4>Python (and MATLAB)</h4>
-
-<ul>
-  <li> </li>
-</ul>
-
 
 <h3>New functionality</h3>
 <ul>
 <li>
   Add TOF capability of the parallelproj projector (see <a href=https://github.com/UCL/STIR/pull/1356>PR #1356</a>)
 </li>
 <li>
-  Read TOF bin order from interfile header (see <a href=https://github.com/UCL/STIR/pull/1389> PR #1389</a>)
+  It is now possible to read TOF bin order from the interfile header (see <a href=https://github.com/UCL/STIR/pull/1389> PR #1389</a>)
 </li>
 <li>
   <code>PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin</code> can now compute the value
@@ -88,7 +83,7 @@ <h3>Changed functionality</h3>
 
 <h3>Bug fixes</h3>
 <ul>
-  <li>The Relative Difference Prior gave incorrect results, probably since switching to C++-14 in version 6.0, although we are not sure.
+  <li>The <strong>Relative Difference Prior gave incorrect results</strong>, probably since switching to C++-14 in version 6.0, although we are not sure.
     See <a href=https://github.com/UCL/STIR/pull/1410>PR #1410</a> and associated <a href=https://github.com/UCL/STIR/pull/1409>issue #1409</a>.
   </li>
   <li>
@@ -102,13 +97,6 @@ <h3>Known problems</h3>
     <p>See <a href=https://github.com/UCL/STIR/labels/bug>our issue tracker</a>.</p>
 
 
-<h3>Documentation changes</h3>
-<ul>
-<li>Added documentation on new features</li> 
-<li>Also check the wiki in addition to the provided PDFs.
-</li>
- </ul>
-
 <H2>What's new for developers (aside from what should be obvious
 from the above):</H2>
 
@@ -129,28 +117,25 @@ <h3>New functionality</h3>
 
 <h3>Other code changes</h3>
 <ul>
-<li>
-</li>
+  <li>
+    Fixes an incompatibility with C++20.
+  </li>
 </ul>
 
 <h3>Build system</h3>
-<ul>
-<li>
-</li>
-</ul>
-
-<h3>Test changes</h3>
-<h4>recon_test_pack changes</h4>
 <ul>
   <li>
-    updated version number and added some clarification to the README.txt
+    Force C++ version according to CERN ROOT versions: ROOT 6.28.10 needs C++17 and 6.30.2 needs C++20.
+    Also some fixes when relying on <code>root-config</code>.
   </li>
 </ul>
 
+<h3>Test changes</h3>
+
 <h4>C++ tests</h4>
 <ul>
   <li>
-    objective function and priors now have a numerical test for <code>accumulate_Hessian_times_input</code>
+    Objective functions (both projection-data and list-mode) and priors now have a numerical test for <code>accumulate_Hessian_times_input</code>
     <br>
     <a href=https://github.com/UCL/STIR/pull/1418> PR #1418</a>
   </li>

examples/python/.gitignore

@@ -0,0 +1,2 @@
+*.png
+

recon_test_pack/README.txt

@@ -1,5 +1,5 @@
 
-README file for STIR recon_test_pack version 3.0 (and later versions?)
+README file for STIR recon_test_pack version 6.1 (and later versions?)
 ----------------------------------------------------------------------
 
 This test pack runs some simple tests to check if various STIR reconstruction

recon_test_pack/run_test_listmode_recon.sh

@@ -20,7 +20,7 @@ if [ -n "$TRAVIS" -o -n "$GITHUB_WORKSPACE" ]; then
     set -e
 fi
 
-echo This script should work with STIR version ">=" 6.0. If you have
+echo This script should work with STIR version ">=" 6.1. If you have
 echo a later version, you might have to update your test pack.
 echo Please check the web site.
 echo

src/buildblock/ExamInfo.cxx

@@ -55,14 +55,14 @@ ExamInfo::parameter_info() const
 bool
 ExamInfo::operator==(const ExamInfo& p1) const
 {
-  return abs(this->up_energy_thres - p1.up_energy_thres) <= 1 &&   // keV
-         abs(this->low_energy_thres - p1.low_energy_thres) <= 1 && // keV
+  return std::abs(this->up_energy_thres - p1.up_energy_thres) <= 1 &&   // keV
+         std::abs(this->low_energy_thres - p1.low_energy_thres) <= 1 && // keV
          this->radionuclide == p1.radionuclide && this->time_frame_definitions == p1.time_frame_definitions &&
          //              this->branching_ratio==p1.branching_ratio &&
          ((this->calibration_factor <= 0 && p1.calibration_factor <= 0)
-          || abs(this->calibration_factor / p1.calibration_factor - 1.) <= 1E-3)
+          || std::abs(this->calibration_factor / p1.calibration_factor - 1.) <= 1E-3)
          && this->imaging_modality == p1.imaging_modality && this->patient_position == p1.patient_position
-         && abs(this->start_time_in_secs_since_1970 - p1.start_time_in_secs_since_1970) <= .5; // sec
+         && std::abs(this->start_time_in_secs_since_1970 - p1.start_time_in_secs_since_1970) <= .5; // sec
 }
 
 END_NAMESPACE_STIR

src/buildblock/ProjDataInfo.cxx

@@ -92,24 +92,24 @@ ProjDataInfo::get_sampling_in_k(const Bin& bin) const
 float
 ProjDataInfo::get_sampling_in_t(const Bin& bin) const
 {
-  return abs(get_t(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() + 1, bin.tangential_pos_num()))
-             - get_t(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() - 1, bin.tangential_pos_num())))
+  return std::abs(get_t(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() + 1, bin.tangential_pos_num()))
+                  - get_t(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() - 1, bin.tangential_pos_num())))
          / 2;
 }
 
 float
 ProjDataInfo::get_sampling_in_m(const Bin& bin) const
 {
-  return abs(get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() + 1, bin.tangential_pos_num()))
-             - get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() - 1, bin.tangential_pos_num())))
+  return std::abs(get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() + 1, bin.tangential_pos_num()))
+                  - get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() - 1, bin.tangential_pos_num())))
          / 2;
 }
 
 float
 ProjDataInfo::get_sampling_in_s(const Bin& bin) const
 {
-  return abs(get_s(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num(), bin.tangential_pos_num() + 1))
-             - get_s(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num(), bin.tangential_pos_num() - 1)))
+  return std::abs(get_s(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num(), bin.tangential_pos_num() + 1))
+                  - get_s(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num(), bin.tangential_pos_num() - 1)))
          / 2;
 }
 

src/buildblock/Radionuclide.cxx

@@ -96,9 +96,9 @@ Radionuclide::get_modality(bool check) const
 bool
 Radionuclide::operator==(const Radionuclide& r) const
 {
-  return (abs(energy - r.energy) <= 1E-1 || (energy <= 0 && r.energy <= 0))
-         && (abs(branching_ratio - r.branching_ratio) <= 1E-1 || (branching_ratio <= 0 && r.branching_ratio <= 0))
-         && (abs(half_life - r.half_life) <= 1 || (r.half_life <= 0 && r.half_life <= 0)) && (modality == r.modality);
+  return (std::abs(energy - r.energy) <= 1E-1 || (energy <= 0 && r.energy <= 0))
+         && (std::abs(branching_ratio - r.branching_ratio) <= 1E-1 || (branching_ratio <= 0 && r.branching_ratio <= 0))
+         && (std::abs(half_life - r.half_life) <= 1 || (r.half_life <= 0 && r.half_life <= 0)) && (modality == r.modality);
 }
 
 END_NAMESPACE_STIR

src/buildblock/TimeFrameDefinitions.cxx

@@ -268,8 +268,8 @@ TimeFrameDefinitions::operator==(const TimeFrameDefinitions& t) const
   for (int frame = 0; frame < frame_times.size(); frame++)
     {
 
-      const bool is_identical = (abs(frame_times.at(frame).first - t.frame_times.at(frame).first) <= 10e-5)
-                                && (abs(frame_times.at(frame).second - t.frame_times.at(frame).second) <= 10e-5);
+      const bool is_identical = (std::abs(frame_times.at(frame).first - t.frame_times.at(frame).first) <= 10e-5)
+                                && (std::abs(frame_times.at(frame).second - t.frame_times.at(frame).second) <= 10e-5);
       if (!is_identical)
         return false;
     }

src/buildblock/extend_projdata.cxx

@@ -74,9 +74,9 @@ extend_segment(const SegmentBySinogram<float>& segment,
   const auto average_phi_sampling = phi_range / (segment.get_proj_data_info_sptr()->get_num_views() - 1);
   // check if 360 or 180 degrees
   // use a rather large tolerance to cope with non-uniform sampling in BlocksOnCylindrical
-  if (abs(phi_range - 2 * _PI) < 5 * average_phi_sampling)
+  if (std::abs(phi_range - 2 * _PI) < 5 * average_phi_sampling)
     flip_views = false; // if views cover 360°, we can simply wrap around
-  else if ((abs(phi_range - _PI) < 5 * average_phi_sampling) && (segment.get_segment_num() == 0))
+  else if ((std::abs(phi_range - _PI) < 5 * average_phi_sampling) && (segment.get_segment_num() == 0))
     flip_views = true; // if views cover 180°, the tangential positions need to be flipped
   else
     {
@@ -97,7 +97,7 @@ extend_segment(const SegmentBySinogram<float>& segment,
             }
           else if (flip_views)
             {
-              const int sym_dim = std::min(abs(min_dim[3]), max_dim[3]);
+              const int sym_dim = std::min(std::abs(min_dim[3]), max_dim[3]);
               for (int tang_pos = -sym_dim; tang_pos <= sym_dim; tang_pos++)
                 {
                   out[axial_pos][min_dim[2] + view_edge][tang_pos]

src/buildblock/find_fwhm_in_image.cxx

@@ -200,14 +200,14 @@ interpolate_line(const Array<3, elemT>& input_array,
 
     line_110 = extract_line(input_array, location_110, dimension);
     line_111 = extract_line(input_array, location_111, dimension);
-    line = line_000 * (1 - abs(z0)) * (1 - abs(y0)) * (1 - abs(x0));
-    line += line_001 * (1 - abs(z0)) * (1 - abs(y0)) * abs(x0);
-    line += line_010 * (1 - abs(z0)) * abs(y0) * (1 - abs(x0));
-    line += line_100 * abs(z0) * (1 - abs(y0)) * (1 - abs(x0));
-    line += line_011 * (1 - abs(z0)) * abs(y0) * abs(x0);
-    line += line_101 * abs(z0) * (1 - abs(y0)) * abs(x0);
-    line += line_110 * abs(z0) * abs(y0) * (1 - abs(x0));
-    line += line_111 * abs(z0) * abs(y0) * abs(x0);
+    line = line_000 * (1 - std::abs(z0)) * (1 - std::abs(y0)) * (1 - std::abs(x0));
+    line += line_001 * (1 - std::abs(z0)) * (1 - std::abs(y0)) * std::abs(x0);
+    line += line_010 * (1 - std::abs(z0)) * std::abs(y0) * (1 - std::abs(x0));
+    line += line_100 * std::abs(z0) * (1 - std::abs(y0)) * (1 - std::abs(x0));
+    line += line_011 * (1 - std::abs(z0)) * std::abs(y0) * std::abs(x0);
+    line += line_101 * std::abs(z0) * (1 - std::abs(y0)) * std::abs(x0);
+    line += line_110 * std::abs(z0) * std::abs(y0) * (1 - std::abs(x0));
+    line += line_111 * std::abs(z0) * std::abs(y0) * std::abs(x0);
   }
   return line;
 }

src/buildblock/inverse_SSRB.cxx

@@ -87,12 +87,12 @@ inverse_SSRB(ProjData& proj_data_4D, const ProjData& proj_data_3D)
                   // for the first slice there is no previous
                   const auto distance_to_previous = in_ax_pos_num == proj_data_3D.get_min_axial_pos_num(0)
                                                         ? std::numeric_limits<float>::max()
-                                                        : abs(out_m - in_m.at(in_ax_pos_num - 1));
-                  const auto distance_to_current = abs(out_m - in_m.at(in_ax_pos_num));
+                                                        : std::abs(out_m - in_m.at(in_ax_pos_num - 1));
+                  const auto distance_to_current = std::abs(out_m - in_m.at(in_ax_pos_num));
                   // for the last slice there is no next
                   const auto distance_to_next = in_ax_pos_num == proj_data_3D.get_max_axial_pos_num(0)
                                                     ? std::numeric_limits<float>::max()
-                                                    : abs(out_m - in_m.at(in_ax_pos_num + 1));
+                                                    : std::abs(out_m - in_m.at(in_ax_pos_num + 1));
                   if (distance_to_current <= distance_to_previous && distance_to_current <= distance_to_next)
                     {
                       if (distance_to_current <= 1E-4)

src/include/stir/numerics/FastErf.h

@@ -65,7 +65,7 @@ class FastErf
   inline void set_num_samples(int num_samples);
 
   //! Returns the maximum sample value
-  inline int get_maximum_sample_value() const;
+  inline double get_maximum_sample_value() const;
   //! Sets the maximum sample value
   inline void set_maximum_sample_value(double maximum_sample_value);