diff --git a/docker/Dockerfile.core b/docker/Dockerfile.core index 01aaf4a..cdaf371 100644 --- a/docker/Dockerfile.core +++ b/docker/Dockerfile.core @@ -27,6 +27,7 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ RUN python3 -m pip install --no-cache-dir --default-timeout=120 --retries=10 \ "conan>=2.0,<3" \ + numpy \ pandas \ pytest diff --git a/localizer/src/CMakeLists.txt b/localizer/src/CMakeLists.txt index 38745b5..afb7cea 100644 --- a/localizer/src/CMakeLists.txt +++ b/localizer/src/CMakeLists.txt @@ -1,7 +1,7 @@ # CMakeList.txt: DCCCcore Refactored Version cmake_minimum_required(VERSION 3.21) -project(DCCCcore VERSION 4.2.6) +project(DCCCcore VERSION 4.3.0) set(DCCCCORE_VERSION_SUFFIX "") set(DCCCCORE_SOFTWARE_VERSION "${PROJECT_VERSION}${DCCCCORE_VERSION_SUFFIX}") @@ -20,6 +20,7 @@ find_package(onnxruntime CONFIG REQUIRED) find_package(tomlplusplus CONFIG REQUIRED) find_package(Eigen3 CONFIG REQUIRED) find_package(rapidcsv CONFIG REQUIRED) +find_package(ZLIB REQUIRED) # Include directories # include_directories(${CMAKE_CURRENT_SOURCE_DIR}) @@ -67,6 +68,8 @@ target_sources(DCCCcore PRIVATE target_sources(DCCCcore PRIVATE core/preprocessing/ImagePreprocessor.h core/preprocessing/ImagePreprocessor.cpp + core/preprocessing/PetMotionCorrector.h + core/preprocessing/PetMotionCorrector.cpp ) # Add utility sources @@ -186,6 +189,8 @@ target_sources(DCCCcore PRIVATE spatialNormalizations/adni/AdniPetCoreCLI.cpp spatialNormalizations/rigid/RigidCLI.h spatialNormalizations/rigid/RigidCLI.cpp + spatialNormalizations/pet/PetMotionCorrectionCLI.h + spatialNormalizations/pet/PetMotionCorrectionCLI.cpp ) # Link libraries @@ -196,6 +201,7 @@ target_link_libraries(DCCCcore PRIVATE tomlplusplus::tomlplusplus Eigen3::Eigen rapidcsv::rapidcsv + ZLIB::ZLIB ) # Compile options diff --git a/localizer/src/README.md b/localizer/src/README.md index 27af390..da75eb7 100644 --- a/localizer/src/README.md +++ b/localizer/src/README.md @@ -110,12 +110,27 @@ Perform spatial standardization without metric calculation: # ADNI-style processing ./DCCCcore adni-pet-core --input pet.nii --output normalized.nii +# Multi-frame PET motion correction and arithmetic mean +./DCCCcore pet-motion-correct --input dynamic_pet.nii.gz --output averaged_pet.nii.gz + +# Optionally retain corrected frames and per-frame rigid motion parameters +./DCCCcore pet-motion-correct --input dynamic_pet.nii.gz --output averaged_pet.nii.gz \ + --save-corrected-dynamic corrected.nii.gz --motion-output motion.tsv + # Iterative rigid registration ./DCCCcore normalize --input pet.nii --output normalized.nii --iterative ./DCCCcore adni-pet-core --input pet.nii --output normalized.nii --iterative ./DCCCcore rigid --input pet.nii --output rigid.nii --iterative ``` +`pet-motion-correct` accepts general 4D `X × Y × Z × N` PET data. It uses frame 0 +as the fixed reference, independently registers every later frame to it with a +six-degree-of-freedom rigid transform, and writes the arithmetic mean as a 3D +float NIfTI. A 4D input passed to `adni-pet-core` is automatically motion-corrected +and averaged before the existing ADNI PET Core pipeline; 3D input follows the +existing pipeline unchanged. Motion TSV translations are in millimetres and +Euler rotations are in radians. + #### ADAD Analysis Run the ADAD decoupling-based metric: diff --git a/localizer/src/conanfile.py b/localizer/src/conanfile.py index 59f22f3..0b80149 100644 --- a/localizer/src/conanfile.py +++ b/localizer/src/conanfile.py @@ -4,7 +4,7 @@ class AppConan(ConanFile): name = "DCCCcore" - version = "4.2.6" + version = "4.3.0" settings = "os", "arch", "compiler", "build_type" generators = "CMakeDeps", "CMakeToolchain" @@ -15,6 +15,7 @@ def requirements(self): self.requires("onnxruntime/1.18.1") self.requires("tomlplusplus/3.4.0") self.requires("rapidcsv/8.84") + self.requires("zlib/1.3.2") # ---- conflict resolution: choose one Eigen for the whole graph ---- # If the graph shows ORT wants 3.4.0, prefer: self.requires("eigen/3.4.0") diff --git a/localizer/src/core/config/Version.h b/localizer/src/core/config/Version.h index 3ff495e..55ec70c 100644 --- a/localizer/src/core/config/Version.h +++ b/localizer/src/core/config/Version.h @@ -1,4 +1,4 @@ #pragma once #include -const std::string SOFTWARE_VERSION = "4.2.6"; +const std::string SOFTWARE_VERSION = "4.3.0"; diff --git a/localizer/src/core/preprocessing/PetMotionCorrector.cpp b/localizer/src/core/preprocessing/PetMotionCorrector.cpp new file mode 100644 index 0000000..b669b2d --- /dev/null +++ b/localizer/src/core/preprocessing/PetMotionCorrector.cpp @@ -0,0 +1,316 @@ +#include "PetMotionCorrector.h" + +#include "../common/NiftiIO.h" +#include "../common/PathUtils.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Pipeline::Preprocessing { + +namespace { + +using TransformType = itk::Euler3DTransform; + +std::uint16_t byteSwap16(std::uint16_t value) { + return static_cast((value >> 8U) | (value << 8U)); +} + +std::uint32_t byteSwap32(std::uint32_t value) { + return ((value & 0x000000FFU) << 24U) | + ((value & 0x0000FF00U) << 8U) | + ((value & 0x00FF0000U) >> 8U) | + ((value & 0xFF000000U) >> 24U); +} + +std::uint64_t byteSwap64(std::uint64_t value) { + return (static_cast(byteSwap32(static_cast(value))) << 32U) | + byteSwap32(static_cast(value >> 32U)); +} + +unsigned int declaredNiftiDimension(const std::string& inputPath) { + const auto fileName = Common::path::legacyFileName(Common::path::fromUtf8(inputPath)); + gzFile input = gzopen(fileName.c_str(), "rb"); + if (!input) { + throw std::runtime_error("Unable to open input image: " + inputPath); + } + std::array header{}; + const int bytesRead = gzread(input, header.data(), static_cast(header.size())); + gzclose(input); + if (bytesRead < 48) { + throw std::runtime_error("Input is too short to contain a valid NIfTI header: " + inputPath); + } + + std::uint32_t headerSize = 0; + std::memcpy(&headerSize, header.data(), sizeof(headerSize)); + bool swapBytes = false; + if (headerSize != 348U && headerSize != 540U) { + headerSize = byteSwap32(headerSize); + swapBytes = true; + } + + std::uint64_t dimension = 0; + if (headerSize == 348U) { + std::uint16_t nifti1Dimension = 0; + std::memcpy(&nifti1Dimension, header.data() + 40, sizeof(nifti1Dimension)); + dimension = swapBytes ? byteSwap16(nifti1Dimension) : nifti1Dimension; + } else if (headerSize == 540U) { + std::memcpy(&dimension, header.data() + 16, sizeof(dimension)); + if (swapBytes) { + dimension = byteSwap64(dimension); + } + } else { + throw std::runtime_error("Input does not contain a valid NIfTI-1 or NIfTI-2 header: " + + inputPath); + } + if (dimension < 1U || dimension > 7U) { + throw std::runtime_error("Input NIfTI declares an invalid image dimension: " + inputPath); + } + return static_cast(dimension); +} + +ImageType::Pointer extractFrame(DynamicImageType::Pointer dynamicImage, unsigned int frame) { + using ExtractFilterType = itk::ExtractImageFilter; + auto region = dynamicImage->GetLargestPossibleRegion(); + auto size = region.GetSize(); + auto index = region.GetIndex(); + size[3] = 0; + index[3] += static_cast(frame); + + auto extractor = ExtractFilterType::New(); + extractor->SetInput(dynamicImage); + extractor->SetExtractionRegion({index, size}); + extractor->SetDirectionCollapseToSubmatrix(); + extractor->Update(); + return extractor->GetOutput(); +} + +TransformType::Pointer registerFrame(ImageType::Pointer fixed, ImageType::Pointer moving) { + using MetricType = itk::MattesMutualInformationImageToImageMetricv4; + using OptimizerType = itk::RegularStepGradientDescentOptimizerv4; + using RegistrationType = itk::ImageRegistrationMethodv4; + using InitializerType = itk::CenteredTransformInitializer; + using ScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift; + + auto transform = TransformType::New(); + auto initializer = InitializerType::New(); + initializer->SetTransform(transform); + initializer->SetFixedImage(fixed); + initializer->SetMovingImage(moving); + initializer->MomentsOn(); + initializer->InitializeTransform(); + + auto metric = MetricType::New(); + metric->SetNumberOfHistogramBins(50); + + auto optimizer = OptimizerType::New(); + optimizer->SetLearningRate(2.0); + optimizer->SetMinimumStepLength(0.0005); + optimizer->SetNumberOfIterations(250); + optimizer->SetRelaxationFactor(0.5); + optimizer->SetReturnBestParametersAndValue(true); + + auto scalesEstimator = ScalesEstimatorType::New(); + scalesEstimator->SetMetric(metric); + scalesEstimator->SetTransformForward(true); + optimizer->SetScalesEstimator(scalesEstimator); + + auto registration = RegistrationType::New(); + registration->SetFixedImage(fixed); + registration->SetMovingImage(moving); + registration->SetMetric(metric); + registration->SetOptimizer(optimizer); + registration->SetInitialTransform(transform); + registration->InPlaceOn(); + registration->SetMetricSamplingStrategy(RegistrationType::MetricSamplingStrategyEnum::REGULAR); + registration->SetMetricSamplingPercentage(0.25); + + RegistrationType::ShrinkFactorsArrayType shrinkFactors; + shrinkFactors.SetSize(3); + shrinkFactors[0] = 4; + shrinkFactors[1] = 2; + shrinkFactors[2] = 1; + RegistrationType::SmoothingSigmasArrayType smoothingSigmas; + smoothingSigmas.SetSize(3); + smoothingSigmas[0] = 2; + smoothingSigmas[1] = 1; + smoothingSigmas[2] = 0; + registration->SetNumberOfLevels(3); + registration->SetShrinkFactorsPerLevel(shrinkFactors); + registration->SetSmoothingSigmasPerLevel(smoothingSigmas); + registration->SmoothingSigmasAreSpecifiedInPhysicalUnitsOn(); + registration->Update(); + return transform; +} + +ImageType::Pointer resampleFrame(ImageType::Pointer moving, + ImageType::Pointer fixed, + TransformType::Pointer transform) { + using ResampleFilterType = itk::ResampleImageFilter; + using InterpolatorType = itk::LinearInterpolateImageFunction; + auto resampler = ResampleFilterType::New(); + resampler->SetInput(moving); + resampler->SetTransform(transform); + resampler->SetInterpolator(InterpolatorType::New()); + resampler->SetUseReferenceImage(true); + resampler->SetReferenceImage(fixed); + resampler->SetDefaultPixelValue(0.0f); + resampler->Update(); + return resampler->GetOutput(); +} + +DynamicImageType::Pointer allocateCorrectedDynamic(DynamicImageType::Pointer source) { + auto output = DynamicImageType::New(); + output->SetRegions(source->GetLargestPossibleRegion()); + output->CopyInformation(source); + output->Allocate(); + output->FillBuffer(0.0f); + return output; +} + +void copyFrameIntoDynamic(ImageType::Pointer frame, + DynamicImageType::Pointer dynamic, + unsigned int frameNumber) { + itk::ImageRegionConstIterator sourceIt(frame, frame->GetLargestPossibleRegion()); + for (sourceIt.GoToBegin(); !sourceIt.IsAtEnd(); ++sourceIt) { + const auto sourceIndex = sourceIt.GetIndex(); + DynamicImageType::IndexType targetIndex; + targetIndex[0] = sourceIndex[0]; + targetIndex[1] = sourceIndex[1]; + targetIndex[2] = sourceIndex[2]; + targetIndex[3] = dynamic->GetLargestPossibleRegion().GetIndex()[3] + + static_cast(frameNumber); + dynamic->SetPixel(targetIndex, sourceIt.Get()); + } +} + +void addFrameToAverage(ImageType::Pointer frame, ImageType::Pointer average) { + itk::ImageRegionConstIterator sourceIt(frame, frame->GetLargestPossibleRegion()); + itk::ImageRegionIterator targetIt(average, average->GetLargestPossibleRegion()); + for (sourceIt.GoToBegin(), targetIt.GoToBegin(); !sourceIt.IsAtEnd(); ++sourceIt, ++targetIt) { + targetIt.Set(targetIt.Get() + sourceIt.Get()); + } +} + +} // namespace + +unsigned int PetMotionCorrector::inspectImageDimension(const std::string& inputPath) { + // ITK intentionally collapses trailing singleton axes. Read the declared NIfTI + // dimension first so X×Y×Z×1 remains distinguishable from a true 3D image. + return declaredNiftiDimension(inputPath); +} + +PetMotionCorrectionResult PetMotionCorrector::correct(const std::string& inputPath, + bool retainCorrectedDynamic) const { + const unsigned int dimension = inspectImageDimension(inputPath); + if (dimension != 4) { + if (dimension == 3) { + throw std::invalid_argument( + "pet-motion-correct requires a multi-frame/4D PET input; received a 3D image."); + } + throw std::invalid_argument( + "pet-motion-correct requires a multi-frame/4D PET input; received a " + + std::to_string(dimension) + "D image."); + } + + using ReaderType = itk::ImageFileReader; + auto reader = ReaderType::New(); + reader->SetFileName(Common::path::legacyFileName(Common::path::fromUtf8(inputPath))); + reader->Update(); + auto dynamicImage = reader->GetOutput(); + + const auto numberOfFrames = dynamicImage->GetLargestPossibleRegion().GetSize()[3]; + if (numberOfFrames == 0) { + throw std::invalid_argument("The 4D PET input contains no frames."); + } + + PetMotionCorrectionResult result; + if (retainCorrectedDynamic) { + result.correctedDynamicImage = allocateCorrectedDynamic(dynamicImage); + } + auto fixed = extractFrame(dynamicImage, 0); + result.averagedImage = ImageType::New(); + result.averagedImage->SetRegions(fixed->GetLargestPossibleRegion()); + result.averagedImage->CopyInformation(fixed); + result.averagedImage->Allocate(); + result.averagedImage->FillBuffer(0.0f); + + if (result.correctedDynamicImage) { + copyFrameIntoDynamic(fixed, result.correctedDynamicImage, 0); + } + addFrameToAverage(fixed, result.averagedImage); + result.motion.push_back({0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); + + for (unsigned int frame = 1; frame < numberOfFrames; ++frame) { + auto moving = extractFrame(dynamicImage, frame); + auto transform = registerFrame(fixed, moving); + auto corrected = resampleFrame(moving, fixed, transform); + if (result.correctedDynamicImage) { + copyFrameIntoDynamic(corrected, result.correctedDynamicImage, frame); + } + addFrameToAverage(corrected, result.averagedImage); + + const auto parameters = transform->GetParameters(); + result.motion.push_back({frame, + parameters[3], parameters[4], parameters[5], + parameters[0], parameters[1], parameters[2]}); + } + + const float scale = 1.0f / static_cast(numberOfFrames); + itk::ImageRegionIterator averageIt( + result.averagedImage, result.averagedImage->GetLargestPossibleRegion()); + for (averageIt.GoToBegin(); !averageIt.IsAtEnd(); ++averageIt) { + averageIt.Set(averageIt.Get() * scale); + } + return result; +} + +void PetMotionCorrector::saveAveragedImage(const PetMotionCorrectionResult& result, + const std::string& outputPath) { + Common::nifti::saveImage(result.averagedImage, outputPath); +} + +void PetMotionCorrector::saveCorrectedDynamicImage(const PetMotionCorrectionResult& result, + const std::string& outputPath) { + if (!result.correctedDynamicImage) { + throw std::invalid_argument( + "Corrected dynamic frames were not retained for this motion-correction run."); + } + using WriterType = itk::ImageFileWriter; + auto writer = WriterType::New(); + writer->SetFileName(Common::path::legacyFileName(Common::path::fromUtf8(outputPath))); + writer->SetInput(result.correctedDynamicImage); + writer->Update(); +} + +void PetMotionCorrector::saveMotionParameters(const PetMotionCorrectionResult& result, + const std::string& outputPath) { + std::ofstream output(Common::path::fromUtf8(outputPath)); + if (!output) { + throw std::runtime_error("Unable to open motion output: " + outputPath); + } + output << "frame\ttx\tty\ttz\trx\try\trz\n" << std::setprecision(12); + for (const auto& motion : result.motion) { + output << motion.frame << '\t' + << motion.tx << '\t' << motion.ty << '\t' << motion.tz << '\t' + << motion.rx << '\t' << motion.ry << '\t' << motion.rz << '\n'; + } +} + +} // namespace Pipeline::Preprocessing diff --git a/localizer/src/core/preprocessing/PetMotionCorrector.h b/localizer/src/core/preprocessing/PetMotionCorrector.h new file mode 100644 index 0000000..a01653f --- /dev/null +++ b/localizer/src/core/preprocessing/PetMotionCorrector.h @@ -0,0 +1,44 @@ +#pragma once + +#include "../common/ImageTypes.h" +#include +#include +#include +#include + +namespace Pipeline::Preprocessing { + +using DynamicImageType = itk::Image; + +struct PetMotionParameters { + unsigned int frame = 0; + double tx = 0.0; + double ty = 0.0; + double tz = 0.0; + double rx = 0.0; + double ry = 0.0; + double rz = 0.0; +}; + +struct PetMotionCorrectionResult { + ImageType::Pointer averagedImage; + DynamicImageType::Pointer correctedDynamicImage; + std::vector motion; +}; + +class PetMotionCorrector { +public: + static unsigned int inspectImageDimension(const std::string& inputPath); + + PetMotionCorrectionResult correct(const std::string& inputPath, + bool retainCorrectedDynamic = false) const; + + static void saveAveragedImage(const PetMotionCorrectionResult& result, + const std::string& outputPath); + static void saveCorrectedDynamicImage(const PetMotionCorrectionResult& result, + const std::string& outputPath); + static void saveMotionParameters(const PetMotionCorrectionResult& result, + const std::string& outputPath); +}; + +} // namespace Pipeline::Preprocessing diff --git a/localizer/src/spatialNormalizations/ModuleCatalog.cpp b/localizer/src/spatialNormalizations/ModuleCatalog.cpp index 4bae63a..8c31ec8 100644 --- a/localizer/src/spatialNormalizations/ModuleCatalog.cpp +++ b/localizer/src/spatialNormalizations/ModuleCatalog.cpp @@ -1,6 +1,7 @@ #include "ModuleCatalog.h" #include "standard/NormalizeCLI.h" #include "adni/AdniPetCoreCLI.h" +#include "pet/PetMotionCorrectionCLI.h" #include "rigid/RigidCLI.h" namespace Pipeline::SpatialNormalization { @@ -10,9 +11,8 @@ std::vector buildCLIModules() { modules.push_back(Standard::createCLI()); modules.push_back(Adni::createCLI()); modules.push_back(Rigid::createCLI()); + modules.push_back(Pet::createMotionCorrectionCLI()); return modules; } } // namespace Pipeline::SpatialNormalization - - diff --git a/localizer/src/spatialNormalizations/adni/AdniPetCoreCLI.cpp b/localizer/src/spatialNormalizations/adni/AdniPetCoreCLI.cpp index f6fa6d1..099f768 100644 --- a/localizer/src/spatialNormalizations/adni/AdniPetCoreCLI.cpp +++ b/localizer/src/spatialNormalizations/adni/AdniPetCoreCLI.cpp @@ -5,9 +5,12 @@ #include "../../core/common/PathUtils.h" #include "../../core/config/Version.h" #include "../../core/di/Bootstrap.h" +#include "../../core/preprocessing/PetMotionCorrector.h" #include "../../core/services/IFileService.h" #include "../../core/services/ISpatialNormalizationService.h" #include "../../metrics/shared/BatchLogging.h" +#include +#include #include #include #include @@ -26,6 +29,32 @@ struct RunConfig { constexpr const char* kBatchOutputSuffix = "_ADNI_style.nii"; +class TemporaryNiftiFile { +public: + TemporaryNiftiFile() = default; + TemporaryNiftiFile(const TemporaryNiftiFile&) = delete; + TemporaryNiftiFile& operator=(const TemporaryNiftiFile&) = delete; + + ~TemporaryNiftiFile() { + if (!path_.empty()) { + std::error_code ec; + std::filesystem::remove(path_, ec); + } + } + + std::string create() { + static std::atomic sequence{0}; + const auto timestamp = std::chrono::steady_clock::now().time_since_epoch().count(); + path_ = std::filesystem::temp_directory_path() / + ("dccc_adni_pet_motion_" + std::to_string(timestamp) + "_" + + std::to_string(sequence.fetch_add(1)) + ".nii.gz"); + return Common::path::toUtf8(path_); + } + +private: + std::filesystem::path path_; +}; + std::string resolveDebugBasePath(const NormalizeCommandOptions& options, const std::string& outputPath) { if (!options.enableDebugOutput || outputPath.empty()) { return {}; @@ -61,6 +90,24 @@ int processSingleImage(const NormalizeCommandOptions& options, auto fileService = container->resolve(); try { + TemporaryNiftiFile averagedDynamicInput; + const unsigned int inputDimension = + Pipeline::Preprocessing::PetMotionCorrector::inspectImageDimension(inputPath); + if (inputDimension == 4) { + std::cout << "[" << config.logTag + << "] 4D PET detected; motion-correcting each frame to frame 0 before ADNI processing." + << std::endl; + Pipeline::Preprocessing::PetMotionCorrector corrector; + auto motionResult = corrector.correct(inputPath); + request.inputPath = averagedDynamicInput.create(); + Pipeline::Preprocessing::PetMotionCorrector::saveAveragedImage( + motionResult, request.inputPath); + } else if (inputDimension != 3) { + throw std::invalid_argument( + "adni-pet-core requires a 3D or 4D PET input; received a " + + std::to_string(inputDimension) + "D image."); + } + auto output = spatialService->normalize(request); fileService->saveNormalizedImage({output.spatiallyNormalizedImage, outputPath}); if (logCompletion) { diff --git a/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.cpp b/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.cpp new file mode 100644 index 0000000..6ec0885 --- /dev/null +++ b/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.cpp @@ -0,0 +1,77 @@ +#include "PetMotionCorrectionCLI.h" + +#include "../../core/common/Filesystem.h" +#include "../../core/preprocessing/PetMotionCorrector.h" +#include +#include + +namespace Pipeline::SpatialNormalization::Pet { + +namespace { + +class PetMotionCorrectionCLI final : public ISpatialNormalizationCLI { +public: + std::string getSubcommandName() const override { + return "pet-motion-correct"; + } + + std::string getDescription() const override { + return "Rigidly align 4D PET frames to frame 0 and calculate their mean"; + } + + void configureArguments(argparse::ArgumentParser& parser) override { + parser.add_argument("--input") + .help("Input multi-frame/4D PET NIfTI path") + .required(); + parser.add_argument("--output") + .help("Output averaged 3D float NIfTI path") + .required(); + parser.add_argument("--save-corrected-dynamic") + .help("Optional output path for the corrected 4D PET") + .default_value(std::string{}); + parser.add_argument("--motion-output") + .help("Optional TSV output path for per-frame rigid motion parameters") + .default_value(std::string{}); + } + + int execute(const argparse::ArgumentParser& parser, const std::string&) override { + const auto inputPath = parser.get("--input"); + const auto outputPath = parser.get("--output"); + const auto correctedPath = parser.get("--save-corrected-dynamic"); + const auto motionPath = parser.get("--motion-output"); + + if (!Common::fs::ensureParentDirectory(outputPath) || + (!correctedPath.empty() && !Common::fs::ensureParentDirectory(correctedPath)) || + (!motionPath.empty() && !Common::fs::ensureParentDirectory(motionPath))) { + std::cerr << "[pet-motion-correct] Failed to prepare an output directory." << std::endl; + return EXIT_FAILURE; + } + + try { + Pipeline::Preprocessing::PetMotionCorrector corrector; + auto result = corrector.correct(inputPath, !correctedPath.empty()); + Pipeline::Preprocessing::PetMotionCorrector::saveAveragedImage(result, outputPath); + if (!correctedPath.empty()) { + Pipeline::Preprocessing::PetMotionCorrector::saveCorrectedDynamicImage( + result, correctedPath); + } + if (!motionPath.empty()) { + Pipeline::Preprocessing::PetMotionCorrector::saveMotionParameters(result, motionPath); + } + std::cout << "[pet-motion-correct] Corrected " << result.motion.size() + << " frame(s); averaged 3D PET saved to " << outputPath << std::endl; + return EXIT_SUCCESS; + } catch (const std::exception& ex) { + std::cerr << "[pet-motion-correct] Processing failed: " << ex.what() << std::endl; + return EXIT_FAILURE; + } + } +}; + +} // namespace + +SpatialNormalizationCLIPtr createMotionCorrectionCLI() { + return std::make_shared(); +} + +} // namespace Pipeline::SpatialNormalization::Pet diff --git a/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.h b/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.h new file mode 100644 index 0000000..dfe63b2 --- /dev/null +++ b/localizer/src/spatialNormalizations/pet/PetMotionCorrectionCLI.h @@ -0,0 +1,9 @@ +#pragma once + +#include "../../core/interfaces/ISpatialNormalizationCLI.h" + +namespace Pipeline::SpatialNormalization::Pet { + +SpatialNormalizationCLIPtr createMotionCorrectionCLI(); + +} // namespace Pipeline::SpatialNormalization::Pet diff --git a/localizer/src/tests/test_adad_cli.py b/localizer/src/tests/test_adad_cli.py index 4c7bf56..1e5c54b 100644 --- a/localizer/src/tests/test_adad_cli.py +++ b/localizer/src/tests/test_adad_cli.py @@ -78,6 +78,7 @@ def test_basic(self, run_subprocess, tmp_path, test_files): f"STDOUT:\n{result.stdout}\nSTDERR:\n{result.stderr}" ) assert output_path.exists(), "adni-pet-core command did not create the expected output file." + assert "4D PET detected" not in result.stdout def test_manual_fov_iterative(self, run_subprocess, tmp_path, test_files): output_path = tmp_path / "adni_pet_core_iter_manual.nii" diff --git a/localizer/src/tests/test_pet_motion_correct_cli.py b/localizer/src/tests/test_pet_motion_correct_cli.py new file mode 100644 index 0000000..3c2a49a --- /dev/null +++ b/localizer/src/tests/test_pet_motion_correct_cli.py @@ -0,0 +1,205 @@ +import gzip +import struct + +import numpy as np + + +_NIFTI_DTYPES = { + 2: np.uint8, + 4: np.int16, + 8: np.int32, + 16: np.float32, + 64: np.float64, + 512: np.uint16, +} + + +def _open_nifti(path, mode): + return gzip.open(path, mode) if str(path).endswith(".gz") else open(path, mode) + + +def _write_float_nifti(path, data): + """Write a small identity-geometry NIfTI without adding a test-only imaging dependency.""" + data = np.asarray(data, dtype="" + ndim = struct.unpack_from(f"{endian}h", payload, 40)[0] + dimensions = struct.unpack_from(f"{endian}{ndim}h", payload, 42) + datatype = struct.unpack_from(f"{endian}h", payload, 70)[0] + offset = int(struct.unpack_from(f"{endian}f", payload, 108)[0]) + dtype = np.dtype(_NIFTI_DTYPES[datatype]).newbyteorder(endian) + count = int(np.prod(dimensions)) + data = np.frombuffer(payload, dtype=dtype, count=count, offset=offset).copy() + slope = struct.unpack_from(f"{endian}f", payload, 112)[0] + intercept = struct.unpack_from(f"{endian}f", payload, 116)[0] + if slope != 0.0: + data = data.astype(np.float32) * slope + intercept + return data.reshape(tuple(reversed(dimensions))), ndim + + +def _nifti_datatype(path): + with _open_nifti(path, "rb") as source: + header = source.read(74) + endian = "<" if struct.unpack_from("" + return struct.unpack_from(f"{endian}h", header, 70)[0] + + +def _synthetic_pet(shape=(32, 32, 32)): + z, y, x = np.indices(shape, dtype=np.float32) + return ( + 1.1 * np.exp(-((x - 10.0) ** 2 / 20.0 + (y - 13.0) ** 2 / 35.0 + (z - 17.0) ** 2 / 26.0)) + + 0.8 * np.exp(-((x - 22.0) ** 2 / 12.0 + (y - 20.0) ** 2 / 20.0 + (z - 10.0) ** 2 / 18.0)) + + 0.45 * np.exp(-((x - 18.0) ** 2 / 8.0 + (y - 8.0) ** 2 / 9.0 + (z - 23.0) ** 2 / 12.0)) + ).astype(np.float32) + + +def _transformed_pet(angle_degrees, translation, shape=(32, 32, 32)): + z, y, x = np.indices(shape, dtype=np.float32) + center = (np.asarray(shape[::-1], dtype=np.float32) - 1.0) / 2.0 + points = np.stack((x, y, z), axis=-1) + angle = np.deg2rad(angle_degrees) + rotation = np.array( + [[np.cos(angle), -np.sin(angle), 0.0], + [np.sin(angle), np.cos(angle), 0.0], + [0.0, 0.0, 1.0]], + dtype=np.float32, + ) + source = (points - center - np.asarray(translation, dtype=np.float32)) @ rotation + center + sx, sy, sz = source[..., 0], source[..., 1], source[..., 2] + return ( + 1.1 * np.exp(-((sx - 10.0) ** 2 / 20.0 + (sy - 13.0) ** 2 / 35.0 + (sz - 17.0) ** 2 / 26.0)) + + 0.8 * np.exp(-((sx - 22.0) ** 2 / 12.0 + (sy - 20.0) ** 2 / 20.0 + (sz - 10.0) ** 2 / 18.0)) + + 0.45 * np.exp(-((sx - 18.0) ** 2 / 8.0 + (sy - 8.0) ** 2 / 9.0 + (sz - 23.0) ** 2 / 12.0)) + ).astype(np.float32) + + +class TestPetMotionCorrectCLI: + def test_rejects_3d_input(self, run_subprocess, tmp_path): + input_path = tmp_path / "single_pet.nii.gz" + output_path = tmp_path / "average.nii.gz" + _write_float_nifti(input_path, _synthetic_pet()) + + result = run_subprocess([ + "pet-motion-correct", "--input", str(input_path), "--output", str(output_path) + ]) + + assert result.returncode != 0 + assert "multi-frame/4D PET" in result.stderr + assert not output_path.exists() + + def test_single_frame_4d_is_preserved(self, run_subprocess, tmp_path): + frame = _synthetic_pet((20, 21, 22)) + input_path = tmp_path / "one_frame.nii.gz" + output_path = tmp_path / "average.nii.gz" + corrected_path = tmp_path / "corrected.nii.gz" + motion_path = tmp_path / "motion.tsv" + _write_float_nifti(input_path, frame[np.newaxis, ...]) + + result = run_subprocess([ + "pet-motion-correct", + "--input", str(input_path), + "--output", str(output_path), + "--save-corrected-dynamic", str(corrected_path), + "--motion-output", str(motion_path), + ]) + + assert result.returncode == 0, result.stderr + averaged, averaged_dimension = _read_nifti(output_path) + corrected, corrected_dimension = _read_nifti(corrected_path) + assert averaged_dimension == 3 + assert _nifti_datatype(output_path) == 16 + assert corrected_dimension == 4 + np.testing.assert_allclose(averaged, frame, rtol=0.0, atol=1e-6) + np.testing.assert_allclose(corrected[0], frame, rtol=0.0, atol=1e-6) + assert motion_path.read_text().splitlines() == [ + "frame\ttx\tty\ttz\trx\try\trz", + "0\t0\t0\t0\t0\t0\t0", + ] + + def test_known_rigid_motion_improves_alignment(self, run_subprocess, tmp_path): + fixed = _synthetic_pet() + frames = np.stack([ + fixed, + _transformed_pet(6.0, (2.0, -1.5, 1.0)), + _transformed_pet(-5.0, (-1.5, 2.0, -1.0)), + ]) + input_path = tmp_path / "moving_frames.nii.gz" + output_path = tmp_path / "average.nii.gz" + corrected_path = tmp_path / "corrected.nii.gz" + motion_path = tmp_path / "motion.tsv" + _write_float_nifti(input_path, frames) + + result = run_subprocess([ + "pet-motion-correct", + "--input", str(input_path), + "--output", str(output_path), + "--save-corrected-dynamic", str(corrected_path), + "--motion-output", str(motion_path), + ]) + + assert result.returncode == 0, result.stderr + corrected, dimension = _read_nifti(corrected_path) + averaged, averaged_dimension = _read_nifti(output_path) + assert dimension == 4 + assert averaged_dimension == 3 + assert _nifti_datatype(output_path) == 16 + before_mse = np.mean((frames[1:] - fixed) ** 2) + after_mse = np.mean((corrected[1:] - fixed) ** 2) + assert after_mse < before_mse * 0.35 + np.testing.assert_allclose(averaged, corrected.mean(axis=0), rtol=1e-5, atol=1e-6) + rows = motion_path.read_text().splitlines() + assert len(rows) == 4 + assert rows[1] == "0\t0\t0\t0\t0\t0\t0" + + +class TestAdniPetCoreDynamicInput: + def test_single_frame_4d_runs_automatic_preprocessing( + self, run_subprocess, tmp_path, test_files + ): + original = test_files["input"].read_bytes() + dynamic_header = bytearray(original) + endian = "<" if struct.unpack_from("" + original_dimension = struct.unpack_from(f"{endian}h", dynamic_header, 40)[0] + assert original_dimension == 3 + struct.pack_into(f"{endian}h", dynamic_header, 40, 4) + struct.pack_into(f"{endian}h", dynamic_header, 48, 1) + dynamic_path = tmp_path / "single_frame_dynamic.nii" + dynamic_path.write_bytes(dynamic_header) + output_path = tmp_path / "adni_dynamic_output.nii" + + result = run_subprocess([ + "adni-pet-core", "--input", str(dynamic_path), "--output", str(output_path) + ]) + + assert result.returncode == 0, ( + f"STDOUT:\n{result.stdout}\nSTDERR:\n{result.stderr}" + ) + assert "4D PET detected" in result.stdout + _, output_dimension = _read_nifti(output_path) + assert output_dimension == 3