Merge pull request #1277 from KrisThielemans/CUDARDP

add STIR 6.2's new CudaRelativeDifferencePrior

CHANGES.md

@@ -7,6 +7,7 @@
   - fixed STIR linking problems by adding STIR_REGISTRIES to Reg executables
   - added `compute_attenuation_factors` method to AcquisitionSensitivityModel
   - added means for setting maximal and minimal value for scale factor in stir::ScatterEstimation
+  - added support for the STIR 6.2 (somewhat experimental) Cuda Relative Difference Prior
 
 * MR
   - fixed density_weight shape issue caused by upgrading numpy to version 2.0

cmake/config.py.in

@@ -4,4 +4,6 @@ SIRF_HAS_NiftyPET = @NiftyPET_BOOL_STR@
 
 SIRF_HAS_Parallelproj = @Parallelproj_BOOL_STR@
 
-SIRF_HAS_SPM = @SPM_BOOL_STR@
+STIR_WITH_CUDA = @STIR_WITH_CUDA_BOOL_STR@
+
+SIRF_HAS_SPM = @SPM_BOOL_STR@

src/CMakeLists.txt

@@ -54,6 +54,7 @@ if (DISABLE_STIR)
   message(STATUS "STIR support disabled.")
   set(NiftyPET_BOOL_STR "0" PARENT_SCOPE)
   set(Parallelproj_BOOL_STR "0" PARENT_SCOPE)
+  set(STIR_WITH_CUDA_BOOL_STR "0" PARENT_SCOPE)
 else()
   find_package(STIR REQUIRED)
   message(STATUS "STIR version found: ${STIR_VERSION}")
@@ -80,13 +81,20 @@ else()
     MESSAGE(STATUS "STIR not built with Parallelproj.")
     set(Parallelproj_BOOL_STR "0")
   endif()
+  if (STIR_WITH_CUDA)
+    set(STIR_WITH_CUDA_BOOL_STR "1")
+        message(STATUS "STIR was built with CUDA, corresponding functionality will be enabled.")
+  else()
+    MESSAGE(STATUS "STIR not built with CUDA.")
+    set(STIR_WITH_CUDA_BOOL_STR "0")
+  endif()
   ADD_SUBDIRECTORY(xSTIR)
   set(SIRF_BUILT_WITH_STIR TRUE PARENT_SCOPE)
   set(STIR_VERSION ${STIR_VERSION} PARENT_SCOPE)
 endif()
 set(NiftyPET_BOOL_STR ${NiftyPET_BOOL_STR} PARENT_SCOPE)
 set(Parallelproj_BOOL_STR ${Parallelproj_BOOL_STR} PARENT_SCOPE)
-
+set (STIR_WITH_CUDA_BOOL_STR ${STIR_WITH_CUDA_BOOL_STR} PARENT_SCOPE)
 
 
 ##########################################################################

src/xSTIR/cSTIR/cstir.cpp

@@ -186,6 +186,10 @@ void* cSTIR_newObject(const char* name)
 			return NEW_OBJECT_HANDLE(LogPrior3DF);
 		if (sirf::iequals(name, "RelativeDifferencePrior"))
 			return NEW_OBJECT_HANDLE(RDPrior3DF);
+#ifdef STIR_WITH_CUDA
+		if (sirf::iequals(name, "CudaRelativeDifferencePrior"))
+			return NEW_OBJECT_HANDLE(CudaRDPrior3DF);
+#endif
 		if (sirf::iequals(name, "PLSPrior"))
 			return NEW_OBJECT_HANDLE(PLSPrior3DF);
 		if (sirf::iequals(name, "TruncateToCylindricalFOVImageProcessor"))

src/xSTIR/cSTIR/cstir_p.cpp

@@ -703,9 +703,11 @@ void*
 sirf::cSTIR_setRelativeDifferencePriorParameter
 (DataHandle* hp, const char* name, const DataHandle* hv)
 {
-	auto& prior = objectFromHandle<xSTIR_RelativeDifferencePrior3DF>(hp);
-	if (sirf::iequals(name, "only_2D"))
-		prior.only2D(dataFromHandle<int>((void*)hv));
+        auto& prior = objectFromHandle<stir::RelativeDifferencePrior<float>>(hp);
+	if (sirf::iequals(name, "only_2D")) {
+		auto& xrdp = objectFromHandle<xSTIR_RelativeDifferencePrior3DF>(hp);
+		xrdp.only2D(dataFromHandle<int>((void*)hv));
+        }
 	else if (sirf::iequals(name, "kappa")) {
 		auto& id = objectFromHandle<STIRImageData>(hv);
 		prior.set_kappa_sptr(id.data_sptr());
@@ -723,7 +725,7 @@ void*
 sirf::cSTIR_RelativeDifferencePriorParameter
 (DataHandle* hp, const char* name)
 {
-	auto& prior = objectFromHandle<xSTIR_RelativeDifferencePrior3DF >(hp);
+	auto& prior = objectFromHandle<stir::RelativeDifferencePrior<float>>(hp);
 	if (sirf::iequals(name, "kappa")) {
 		auto sptr_im = std::make_shared<STIRImageData>(*prior.get_kappa_sptr());
 		return newObjectHandle(sptr_im);

src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h

@@ -1,7 +1,7 @@
 /*
 SyneRBI Synergistic Image Reconstruction Framework (SIRF)
-Copyright 2015 - 2019 Rutherford Appleton Laboratory STFC
-Copyright 2019 - 2020 University College London
+Copyright 2015 - 2024 Rutherford Appleton Laboratory STFC
+Copyright 2019 - 2024 University College London
 
 This is software developed for the Collaborative Computational
 Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
@@ -64,6 +64,9 @@ limitations under the License.
 #include "stir/recon_buildblock/QuadraticPrior.h"
 #include "stir/recon_buildblock/LogcoshPrior.h"
 #include "stir/recon_buildblock/RelativeDifferencePrior.h"
+#ifdef STIR_WITH_CUDA
+#include "stir/recon_buildblock/CUDA/CudaRelativeDifferencePrior.h"
+#endif
 #include "stir/SegmentBySinogram.h"
 #include "stir/Shape/Box3D.h"
 #include "stir/Shape/Ellipsoid.h"
@@ -116,7 +119,10 @@ namespace sirf {
 	typedef stir::QuadraticPrior<float> QuadPrior3DF;
 	typedef stir::LogcoshPrior<float> LogPrior3DF;
 	typedef stir::RelativeDifferencePrior<float> RDPrior3DF;
-	typedef stir::PLSPrior<float> PLSPrior3DF;
+#ifdef STIR_WITH_CUDA
+	typedef stir::CudaRelativeDifferencePrior<float> CudaRDPrior3DF;
+#endif
+        typedef stir::PLSPrior<float> PLSPrior3DF;
 	typedef stir::DataProcessor<Image3DF> DataProcessor3DF;
 	typedef stir::TruncateToCylindricalFOVImageProcessor<float> CylindricFilter3DF;
 

src/xSTIR/pSTIR/STIR.py

@@ -44,6 +44,7 @@
 import sirf.STIR_params as parms
 from sirf.config import SIRF_HAS_NiftyPET
 from sirf.config import SIRF_HAS_Parallelproj
+from sirf.config import STIR_WITH_CUDA
 
 if sys.version_info[0] >= 3 and sys.version_info[1] >= 4:
     ABC = abc.ABC
@@ -2529,6 +2530,24 @@ def get_kappa(self):
         check_status(image.handle)
         return image
 
+if STIR_WITH_CUDA:
+    class CudaRelativeDifferencePrior(RelativeDifferencePrior):
+        r"""Class for Relative Difference Prior using CUDA computations
+
+        Identical to RelativeDifferencePrior, but using STIR's CUDA implementation.
+        """
+
+        def __init__(self):
+            """init."""
+            self.name = 'CudaRelativeDifferencePrior'
+            self.handle = pystir.cSTIR_newObject(self.name)
+            check_status(self.handle)
+
+        def __del__(self):
+            """del."""
+            if self.handle is not None:
+                pyiutil.deleteDataHandle(self.handle)
+
 
 class PLSPrior(Prior):
     r"""Class for Parallel Level Sets prior.

src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py

@@ -15,10 +15,11 @@
 {licence}
 """
 import sirf.STIR
+import sirf.config
 from sirf.Utilities import runner, RE_PYEXT, __license__, examples_data_path, pTest
-__version__ = "2.0.0"
+__version__ = "2.1.0"
 __author__ = "Imraj Singh, Evgueni Ovtchinnikov, Kris Thielemans"
-  
+
 
 def Hessian_test(test, prior, x, eps=1e-3):
     """Checks that grad(x + dx) - grad(x) is close to H(x)*dx
@@ -52,18 +53,21 @@ def test_main(rec=False, verb=False, throw=True):
 
     _ = sirf.STIR.MessageRedirector(warn=None)
 
-    im_thorax = sirf.STIR.ImageData(examples_data_path('PET') + '/thorax_single_slice/emission.hv')
-    im_1 = im_thorax.get_uniform_copy(1)
-    im_2 = im_thorax.get_uniform_copy(2)
+    im_0 = sirf.STIR.ImageData(examples_data_path('PET') + '/brain/emission.hv')
+    im_1 = im_0.get_uniform_copy(1)
+    im_2 = im_0.get_uniform_copy(2)
 
-    for im in [im_thorax, im_1, im_2]:
+    for im in [im_0, im_1, im_2]:
       for penalisation_factor in [0,1,4]:
         for kappa in [True, False]:
-          for prior in [sirf.STIR.QuadraticPrior(), sirf.STIR.LogcoshPrior(), sirf.STIR.RelativeDifferencePrior()]:
+          priors = [sirf.STIR.QuadraticPrior(), sirf.STIR.LogcoshPrior(), sirf.STIR.RelativeDifferencePrior()]
+          if sirf.config.STIR_WITH_CUDA:
+              priors.append(sirf.STIR.CudaRelativeDifferencePrior())
+          for prior in priors:
             if kappa:
-              prior.set_kappa(im_thorax)
+              prior.set_kappa(im_0)
               # Check if the kappa is stored/returned correctly
-              test.check_if_equal_within_tolerance(im_thorax.norm(),prior.get_kappa().norm())
+              test.check_if_equal_within_tolerance(im_0.norm(),prior.get_kappa().norm())
 
             prior.set_penalisation_factor(penalisation_factor)
             prior.set_up(im)
@@ -78,7 +82,7 @@ def test_main(rec=False, verb=False, throw=True):
               if kappa:
                     # check if multiplying kappa and dividing penalisation factor gives same result
                    prior.set_penalisation_factor(penalisation_factor/4)
-                   prior.set_kappa(im_thorax*2)
+                   prior.set_kappa(im_0*2)
                    prior.set_up(im)
                    test.check_if_equal_within_tolerance(prior.get_gradient(im).norm(), grad_norm)