diff --git a/CHANGES.md b/CHANGES.md index 78260c398..d63c644d0 100644 --- a/CHANGES.md +++ b/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 diff --git a/cmake/config.py.in b/cmake/config.py.in index b4e70f9c4..8b25da9b5 100644 --- a/cmake/config.py.in +++ b/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@ \ No newline at end of file +STIR_WITH_CUDA = @STIR_WITH_CUDA_BOOL_STR@ + +SIRF_HAS_SPM = @SPM_BOOL_STR@ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1a5e2f850..c41c2bd98 100644 --- a/src/CMakeLists.txt +++ b/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) ########################################################################## diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index cae22fa01..7757d665e 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/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")) diff --git a/src/xSTIR/cSTIR/cstir_p.cpp b/src/xSTIR/cSTIR/cstir_p.cpp index 35b5cf2bb..7f3d29781 100644 --- a/src/xSTIR/cSTIR/cstir_p.cpp +++ b/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(hp); - if (sirf::iequals(name, "only_2D")) - prior.only2D(dataFromHandle((void*)hv)); + auto& prior = objectFromHandle>(hp); + if (sirf::iequals(name, "only_2D")) { + auto& xrdp = objectFromHandle(hp); + xrdp.only2D(dataFromHandle((void*)hv)); + } else if (sirf::iequals(name, "kappa")) { auto& id = objectFromHandle(hv); prior.set_kappa_sptr(id.data_sptr()); @@ -723,7 +725,7 @@ void* sirf::cSTIR_RelativeDifferencePriorParameter (DataHandle* hp, const char* name) { - auto& prior = objectFromHandle(hp); + auto& prior = objectFromHandle>(hp); if (sirf::iequals(name, "kappa")) { auto sptr_im = std::make_shared(*prior.get_kappa_sptr()); return newObjectHandle(sptr_im); diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h index 7d88ee17d..a1e8b5025 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h +++ b/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 QuadPrior3DF; typedef stir::LogcoshPrior LogPrior3DF; typedef stir::RelativeDifferencePrior RDPrior3DF; - typedef stir::PLSPrior PLSPrior3DF; +#ifdef STIR_WITH_CUDA + typedef stir::CudaRelativeDifferencePrior CudaRDPrior3DF; +#endif + typedef stir::PLSPrior PLSPrior3DF; typedef stir::DataProcessor DataProcessor3DF; typedef stir::TruncateToCylindricalFOVImageProcessor CylindricFilter3DF; diff --git a/src/xSTIR/pSTIR/STIR.py b/src/xSTIR/pSTIR/STIR.py index 5fddaa42d..efbc63c00 100644 --- a/src/xSTIR/pSTIR/STIR.py +++ b/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. diff --git a/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py b/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py index c54eab676..ef7741bba 100644 --- a/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py +++ b/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)