diff --git a/docs/rigid-iterative-ct-benchmark.md b/docs/rigid-iterative-ct-benchmark.md new file mode 100644 index 0000000..3c91499 --- /dev/null +++ b/docs/rigid-iterative-ct-benchmark.md @@ -0,0 +1,126 @@ +# Rigid Iterative CT Benchmark + +Date: 2026-06-16 + +Input: + +- `/workspace/localizer/Testing/CT.nii.gz` + +Build: + +```bash +docker compose -f docker-compose.core.yml run --rm dccc-core /workspace/scripts/docker-build-core.sh +``` + +Benchmark command shape: + +```bash +docker compose -f docker-compose.core.yml run --rm dccc-core \ + python3 -c '... subprocess.run(["./install/bin/DCCCcore", "rigid", "--input", "/workspace/localizer/Testing/CT.nii.gz", "--output", "/tmp/dccc_ct_rigid_*.nii", "--iterative"]) ...' +``` + +`/usr/bin/time` is not installed in the development container, so the benchmark used Python `time.perf_counter()` for wall time and `resource.getrusage(resource.RUSAGE_CHILDREN).ru_maxrss` for peak child-process RSS. + +## Baseline + +Baseline command: + +```bash +./install/bin/DCCCcore rigid \ + --input /workspace/localizer/Testing/CT.nii.gz \ + --output /tmp/dccc_ct_rigid_baseline.nii \ + --iterative +``` + +Result: + +- Return code: 0 +- Wall time: 32.233 seconds +- Max RSS: 6,810,932 KB + +The slow path was not caused by repeatedly reading the input file from disk. The input is loaded once. The main issues were: + +- `rigidOnly` requests still computed VoxelMorph spatial normalization, then discarded it. +- Rigid preprocessing computed intensity percentiles by sorting all high-resolution voxels as `double`. +- The iterative loop wrote `rigid_iter.nii` every pass, but that file was never read. + +## Changes Tested + +Implemented changes: + +- Added rigid-only normalization entry points so the `rigid` command stops before VoxelMorph. +- Split rigid alignment and VoxelMorph warping into separate lazy-loaded components, so `rigid` does not load VoxelMorph and `--manual-fov` does not load the rigid model. +- Removed the old combined `RegistrationPipeline` from the build to avoid reintroducing a runtime path that initializes both engines together. +- Removed the unused per-iteration `rigid_iter.nii` temporary write. +- Replaced full percentile sorting with `std::nth_element` over a `float` voxel buffer. + +An attempted low-resolution iterative working-image cache was rejected. It improved speed, but changed the final affine too much on `CT.nii.gz`; the old/new sform difference reached 19.5 mm in translation. The final implementation therefore preserves the original high-resolution iterative affine-estimation semantics. + +## Results + +Rejected low-resolution-cache experiment: + +- Return code: 0 +- Wall time: 16.848 seconds +- Max old/new sform absolute difference: 19.525314 mm + +Final result before component split: + +- Return code: 0 +- Wall time: 16.991 seconds +- Max old/new sform absolute difference: 0.0 +- Max old/new qoffset absolute difference: 0.0 +- Max old/new quaternion absolute difference: 0.0 +- Max RSS: 6,241,604 KB + +Final result after Rigid/VoxelMorph component split: + +- Return code: 0 +- Wall time: 6.270 seconds +- Max old/new sform absolute difference: 0.0 +- Max old/new qoffset absolute difference: 0.0 +- Max old/new quaternion absolute difference: 0.0 +- Max RSS: 2,009,480 KB + +An earlier final timing run before the affine correction measured 16.848 seconds and max RSS 6,603,216 KB, but that variant is not valid because the affine changed. + +After rigid-only branching, removal of the low-resolution-cache experiment, and component-level lazy loading, the accepted output preserves the old affine exactly for this CT benchmark while avoiding discarded VoxelMorph work and VoxelMorph model loading. + +Earlier intermediate result after rigid-only branching and low-resolution iterative working image, before the percentile change: + +- Return code: 0 +- Wall time: 23.528 seconds +- Max RSS: 6,543,312 KB + +This intermediate low-resolution result is included only for auditability and is not the accepted implementation. + +Compared with baseline: + +- Wall time improved by about 80.5%. +- Peak RSS improved because the `rigid` path no longer loads the VoxelMorph model or padded VoxelMorph processing buffers. + +## Verification + +Commands run: + +```bash +docker compose -f docker-compose.core.yml run --rm dccc-core /workspace/scripts/docker-build-core.sh +docker compose -f docker-compose.core.yml run --rm dccc-core pytest tests/test_normalize_cli.py -q +docker compose -f docker-compose.core.yml run --rm dccc-core pytest tests/test_acc_centiloid_centaurz_cli.py -q +``` + +Test result: + +- `tests/test_normalize_cli.py`: 4 passed in 32.11 seconds +- `tests/test_acc_centiloid_centaurz_cli.py`: 30 passed in 445.83 seconds +- `DCCCcore --version`: `4.2.3-alpha` + +Component split behavior check: + +- With a config whose `models.rigid` path points to a missing file, `normalize --manual-fov` returned 0. +- With the same config, `rigid` returned 1 while loading the missing rigid model. +- This confirms `--manual-fov` no longer initializes the rigid engine, while rigid paths still validate the rigid model. + +## Notes + +Further speed or memory reductions would require changing the iterative model input or the first rigid preprocessing pass so that very high-resolution images are downsampled before percentile clipping and smoothing. The low-resolution-cache experiment shows that this is a behavioral change and must be validated against affine accuracy before acceptance. diff --git a/localizer/src/CMakeLists.txt b/localizer/src/CMakeLists.txt index 6567911..24f1d31 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.2) +project(DCCCcore VERSION 4.2.3) set(DCCCCORE_VERSION_SUFFIX "-alpha") set(DCCCCORE_SOFTWARE_VERSION "${PROJECT_VERSION}${DCCCCORE_VERSION_SUFFIX}") @@ -53,8 +53,10 @@ target_sources(DCCCcore PRIVATE target_sources(DCCCcore PRIVATE core/normalizers/RigidVoxelMorphNormalizer.h core/normalizers/RigidVoxelMorphNormalizer.cpp - core/normalizers/RegistrationPipeline.h - core/normalizers/RegistrationPipeline.cpp + core/normalizers/RigidAlignmentNormalizer.h + core/normalizers/RigidAlignmentNormalizer.cpp + core/normalizers/VoxelMorphNormalizer.h + core/normalizers/VoxelMorphNormalizer.cpp core/normalizers/RigidRegistrationEngine.h core/normalizers/RigidRegistrationEngine.cpp core/normalizers/NonlinearRegistrationEngine.h diff --git a/localizer/src/conanfile.py b/localizer/src/conanfile.py index dbd9124..987db5a 100644 --- a/localizer/src/conanfile.py +++ b/localizer/src/conanfile.py @@ -4,7 +4,7 @@ class AppConan(ConanFile): name = "DCCCcore" - version = "4.2.2-alpha" + version = "4.2.3-alpha" settings = "os", "arch", "compiler", "build_type" generators = "CMakeDeps", "CMakeToolchain" diff --git a/localizer/src/core/config/Version.h b/localizer/src/core/config/Version.h index 5339b53..63a828d 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.2-alpha"; +const std::string SOFTWARE_VERSION = "4.2.3-alpha"; diff --git a/localizer/src/core/normalizers/RegistrationPipeline.cpp b/localizer/src/core/normalizers/RegistrationPipeline.cpp deleted file mode 100644 index 251e21b..0000000 --- a/localizer/src/core/normalizers/RegistrationPipeline.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include "RegistrationPipeline.h" - -RegistrationPipeline::RegistrationPipeline(const std::string& rigidModelPath, - const std::string& nonlinearModelPath) { - rigidEngine_ = std::make_unique(rigidModelPath); - nonlinearEngine_ = std::make_unique(nonlinearModelPath); -} - -ImageType::Pointer RegistrationPipeline::preprocess(ImageType::Pointer image) { - return ImagePreprocessor::preprocessForRigid(image); -} - -ImageType::Pointer RegistrationPipeline::preprocessVoxelMorph(ImageType::Pointer image) { - return ImagePreprocessor::preprocessForVoxelMorph(image); -} - -std::unordered_map> RegistrationPipeline::predict( - std::vector inputTensor, const std::vector inputShape) { - return rigidEngine_->predict(inputTensor, inputShape); -} - -std::unordered_map> RegistrationPipeline::predictVoxelMorph( - std::vector originalImg, std::vector movingImg, - std::vector templateImg) { - return nonlinearEngine_->predict(originalImg, movingImg, templateImg); -} - -std::tuple -RegistrationPipeline::getNewOriginAndDirection( - ImageType::Pointer preprocessedImage, ImageType::Pointer originalImage, - std::vector ac, std::vector pa, std::vector is) { - return rigidEngine_->getNewOriginAndDirection(preprocessedImage, originalImage, ac, pa, is); -} diff --git a/localizer/src/core/normalizers/RegistrationPipeline.h b/localizer/src/core/normalizers/RegistrationPipeline.h deleted file mode 100644 index b11018e..0000000 --- a/localizer/src/core/normalizers/RegistrationPipeline.h +++ /dev/null @@ -1,37 +0,0 @@ -#pragma once -#include "../common/ImageTypes.h" -#include "../preprocessing/ImagePreprocessor.h" -#include "RigidRegistrationEngine.h" -#include "NonlinearRegistrationEngine.h" -#include -#include - -/** - * @brief Complete registration pipeline combining rigid and nonlinear registration - * Replaces the old "Rigid" class with better separation of concerns - */ -class RegistrationPipeline { -public: - RegistrationPipeline(const std::string& rigidModelPath, const std::string& nonlinearModelPath); - ~RegistrationPipeline() = default; - - // Main interface methods (compatible with old Rigid class interface) - ImageType::Pointer preprocess(ImageType::Pointer image); - ImageType::Pointer preprocessVoxelMorph(ImageType::Pointer image); - - std::unordered_map> predict( - std::vector inputTensor, const std::vector inputShape); - - std::unordered_map> predictVoxelMorph( - std::vector originalImg, std::vector movingImg, - std::vector templateImg); - - std::tuple - getNewOriginAndDirection( - ImageType::Pointer preprocessedImage, ImageType::Pointer originalImage, - std::vector ac, std::vector pa, std::vector is); - -private: - std::unique_ptr rigidEngine_; - std::unique_ptr nonlinearEngine_; -}; diff --git a/localizer/src/core/normalizers/RigidAlignmentNormalizer.cpp b/localizer/src/core/normalizers/RigidAlignmentNormalizer.cpp new file mode 100644 index 0000000..b4c3cd2 --- /dev/null +++ b/localizer/src/core/normalizers/RigidAlignmentNormalizer.cpp @@ -0,0 +1,105 @@ +#include "RigidAlignmentNormalizer.h" + +#include "../common/Common.h" +#include "../preprocessing/ImagePreprocessor.h" + +#include +#include + +RigidAlignmentNormalizer::RigidAlignmentNormalizer(ConfigurationPtr config) + : config_(config) { + if (!config_) { + throw std::invalid_argument("RigidAlignmentNormalizer requires configuration"); + } + + rigidEngine_ = + std::make_unique(config_->getModelPath("rigid")); + paddedTemplate_ = Common::nifti::loadImage(config_->getTemplatePath("padded")); +} + +ImageType::Pointer RigidAlignmentNormalizer::align(ImageType::Pointer inputImage) { + ImageType::Pointer rigidImage = performAlignment(inputImage); + saveDebugImage(rigidImage, "rigid"); + return rigidImage; +} + +ImageType::Pointer RigidAlignmentNormalizer::alignIterative( + ImageType::Pointer inputImage, int maxIter, float threshold) { + ImageType::Pointer currentImage = performAlignment(inputImage, false); + saveDebugImage(currentImage, "rigid0"); + ImageType::PointType lastOrigin = currentImage->GetOrigin(); + + for (int i = 0; i < maxIter; ++i) { + currentImage = performAlignment(currentImage, true); + saveDebugImage(currentImage, "rigid" + std::to_string(i + 1)); + + float originShift = 0; + for (int j = 0; j < 3; ++j) { + originShift += + std::pow(currentImage->GetOrigin()[j] - lastOrigin[j], 2); + } + originShift = std::sqrt(originShift); + + if (originShift < threshold) { + break; + } + lastOrigin = currentImage->GetOrigin(); + } + + return currentImage; +} + +RigidAlignmentNormalizer::AlignmentEstimate RigidAlignmentNormalizer::estimate( + ImageType::Pointer inputImage, bool resampleFirst) { + ImageType::Pointer processedImage = inputImage; + + if (resampleFirst) { + processedImage = Common::image::resampleToMatch(paddedTemplate_, processedImage); + } + + processedImage = ImagePreprocessor::preprocessForRigid(processedImage); + saveDebugImage(processedImage, "rigid_preprocessed"); + + std::vector imageData; + Common::image::extractImageData(processedImage, imageData); + + auto orientation = rigidEngine_->predict(imageData, {1, 1, 64, 64, 64}); + return AlignmentEstimate{processedImage, orientation}; +} + +void RigidAlignmentNormalizer::apply(ImageType::Pointer targetImage, + const AlignmentEstimate& estimate) { + if (!targetImage) { + return; + } + + auto newOriginAndDirection = rigidEngine_->getNewOriginAndDirection( + estimate.preprocessedImage, + targetImage, + estimate.orientation.at("ac"), + estimate.orientation.at("pa"), + estimate.orientation.at("is")); + + targetImage->SetOrigin(std::get<0>(newOriginAndDirection)); + targetImage->SetDirection(std::get<1>(newOriginAndDirection)); +} + +ImageType::Pointer RigidAlignmentNormalizer::performAlignment( + ImageType::Pointer inputImage, bool resampleFirst) { + auto alignmentEstimate = estimate(inputImage, resampleFirst); + apply(inputImage, alignmentEstimate); + return inputImage; +} + +void RigidAlignmentNormalizer::setDebugMode(bool enable, const std::string& basePath) { + debugMode_ = enable; + debugBasePath_ = basePath; +} + +void RigidAlignmentNormalizer::saveDebugImage(ImageType::Pointer image, + const std::string& suffix) { + if (!debugMode_ || debugBasePath_.empty()) { + return; + } + Common::nifti::saveImage(image, debugBasePath_ + "_" + suffix + ".nii"); +} diff --git a/localizer/src/core/normalizers/RigidAlignmentNormalizer.h b/localizer/src/core/normalizers/RigidAlignmentNormalizer.h new file mode 100644 index 0000000..41ec1cd --- /dev/null +++ b/localizer/src/core/normalizers/RigidAlignmentNormalizer.h @@ -0,0 +1,39 @@ +#pragma once + +#include "../common/ImageTypes.h" +#include "../interfaces/IConfiguration.h" +#include "RigidRegistrationEngine.h" + +#include +#include +#include +#include + +class RigidAlignmentNormalizer { +public: + explicit RigidAlignmentNormalizer(ConfigurationPtr config); + + ImageType::Pointer align(ImageType::Pointer inputImage); + ImageType::Pointer alignIterative(ImageType::Pointer inputImage, + int maxIter = 5, + float threshold = 2.0f); + + void setDebugMode(bool enable, const std::string& basePath = ""); + +private: + struct AlignmentEstimate { + ImageType::Pointer preprocessedImage; + std::unordered_map> orientation; + }; + + ConfigurationPtr config_; + std::unique_ptr rigidEngine_; + ImageType::Pointer paddedTemplate_; + bool debugMode_ = false; + std::string debugBasePath_; + + AlignmentEstimate estimate(ImageType::Pointer inputImage, bool resampleFirst = false); + void apply(ImageType::Pointer targetImage, const AlignmentEstimate& estimate); + ImageType::Pointer performAlignment(ImageType::Pointer inputImage, bool resampleFirst = false); + void saveDebugImage(ImageType::Pointer image, const std::string& suffix); +}; diff --git a/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.cpp b/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.cpp index c26316e..fd7408e 100644 --- a/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.cpp +++ b/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.cpp @@ -1,218 +1,93 @@ #include "RigidVoxelMorphNormalizer.h" -#include "../common/Common.h" -#include -RigidVoxelMorphNormalizer::RigidVoxelMorphNormalizer(ConfigurationPtr config) +#include + +RigidVoxelMorphNormalizer::RigidVoxelMorphNormalizer(ConfigurationPtr config) : config_(config) { - initializeModel(); + if (!config_) { + throw std::invalid_argument("RigidVoxelMorphNormalizer requires configuration"); + } +} + +ImageType::Pointer RigidVoxelMorphNormalizer::normalize(ImageType::Pointer inputImage) { + ImageType::Pointer rigidImage = rigidNormalizer().align(inputImage); + return voxelMorphNormalizer().normalize(rigidImage); } -void RigidVoxelMorphNormalizer::initializeModel() { - std::string rigidModelPath = config_->getModelPath("rigid"); - std::string voxelMorphPath = config_->getModelPath("affine_voxelmorph"); - std::string templatePath = config_->getTemplatePath("padded"); - - registrationPipeline_ = std::make_unique(rigidModelPath, voxelMorphPath); - paddedTemplate_ = Common::nifti::loadImage(templatePath); +ImageType::Pointer RigidVoxelMorphNormalizer::normalizeRigidOnly( + ImageType::Pointer inputImage) { + return rigidNormalizer().align(inputImage); } -ImageType::Pointer RigidVoxelMorphNormalizer::normalize(ImageType::Pointer inputImage) { - // Standard spatial normalization pipeline - ImageType::Pointer rigidImage = performRigidAlignment(inputImage); - saveDebugImage(rigidImage, "rigid"); - return performVoxelMorphWarping(rigidImage); +ImageType::Pointer RigidVoxelMorphNormalizer::normalizeIterativeRigidOnly( + ImageType::Pointer inputImage, int maxIter, float threshold) { + return rigidNormalizer().alignIterative(inputImage, maxIter, threshold); } ImageType::Pointer RigidVoxelMorphNormalizer::normalizeIterative( ImageType::Pointer inputImage, int maxIter, float threshold) { - - ImageType::Pointer currentImage = performRigidAlignment(inputImage, false); - saveDebugImage(currentImage, "rigid0"); - ImageType::PointType lastOrigin = currentImage->GetOrigin(); - - for (int i = 0; i < maxIter; ++i) { - // Save temporary file - std::string tempPath = config_->getTempDirPath() + "/rigid_iter.nii"; - Common::nifti::saveImage(currentImage, tempPath); - - // Execute next rigid registration - currentImage = performRigidAlignment(currentImage, true); - saveDebugImage(currentImage, "rigid" + std::to_string(i + 1)); - - // Check convergence - float originShift = 0; - for (int j = 0; j < 3; ++j) { - originShift += std::pow(currentImage->GetOrigin()[j] - lastOrigin[j], 2); - } - originShift = std::sqrt(originShift); - - if (originShift < threshold) { - break; - } - lastOrigin = currentImage->GetOrigin(); - } - - return performVoxelMorphWarping(currentImage); + ImageType::Pointer rigidImage = + rigidNormalizer().alignIterative(inputImage, maxIter, threshold); + return voxelMorphNormalizer().normalize(rigidImage); } -ImageType::Pointer RigidVoxelMorphNormalizer::normalizeManualFOV(ImageType::Pointer inputImage) { - // Skip rigid registration, perform nonlinear registration directly - return performVoxelMorphWarping(inputImage); +ImageType::Pointer RigidVoxelMorphNormalizer::normalizeManualFOV( + ImageType::Pointer inputImage) { + return voxelMorphNormalizer().normalize(inputImage); } -RigidVoxelMorphNormalizer::NormalizationResult RigidVoxelMorphNormalizer::normalizeWithIntermediateResults(ImageType::Pointer inputImage) { +RigidVoxelMorphNormalizer::NormalizationResult +RigidVoxelMorphNormalizer::normalizeWithIntermediateResults(ImageType::Pointer inputImage) { NormalizationResult result; - - // Perform rigid alignment - result.rigidAlignedImage = performRigidAlignment(inputImage); - saveDebugImage(result.rigidAlignedImage, "rigid"); - - // Perform VoxelMorph warping - result.spatiallyNormalizedImage = performVoxelMorphWarping(result.rigidAlignedImage); - + result.rigidAlignedImage = rigidNormalizer().align(inputImage); + result.spatiallyNormalizedImage = + voxelMorphNormalizer().normalize(result.rigidAlignedImage); return result; } -RigidVoxelMorphNormalizer::NormalizationResult RigidVoxelMorphNormalizer::normalizeIterativeWithIntermediateResults(ImageType::Pointer inputImage, int maxIter, float threshold) { +RigidVoxelMorphNormalizer::NormalizationResult +RigidVoxelMorphNormalizer::normalizeIterativeWithIntermediateResults( + ImageType::Pointer inputImage, int maxIter, float threshold) { NormalizationResult result; - - ImageType::Pointer currentImage = performRigidAlignment(inputImage, false); - saveDebugImage(currentImage, "rigid0"); - ImageType::PointType lastOrigin = currentImage->GetOrigin(); - - for (int i = 0; i < maxIter; ++i) { - // Save temporary file - std::string tempPath = config_->getTempDirPath() + "/rigid_iter.nii"; - Common::nifti::saveImage(currentImage, tempPath); - - // Execute next rigid registration - currentImage = performRigidAlignment(currentImage, true); - saveDebugImage(currentImage, "rigid" + std::to_string(i + 1)); - - // Check convergence - float originShift = 0; - for (int j = 0; j < 3; ++j) { - originShift += std::pow(currentImage->GetOrigin()[j] - lastOrigin[j], 2); - } - originShift = std::sqrt(originShift); - - if (originShift < threshold) { - break; - } - lastOrigin = currentImage->GetOrigin(); - } - - result.rigidAlignedImage = currentImage; - result.spatiallyNormalizedImage = performVoxelMorphWarping(currentImage); - + result.rigidAlignedImage = + rigidNormalizer().alignIterative(inputImage, maxIter, threshold); + result.spatiallyNormalizedImage = + voxelMorphNormalizer().normalize(result.rigidAlignedImage); return result; } -ImageType::Pointer RigidVoxelMorphNormalizer::performRigidAlignment(ImageType::Pointer inputImage, bool resampleFirst) { - ImageType::Pointer processedImage = inputImage; - - if (resampleFirst) { - processedImage = Common::image::resampleToMatch(paddedTemplate_, processedImage); - } - - processedImage = registrationPipeline_->preprocess(processedImage); - saveDebugImage(processedImage, "rigid_preprocessed"); - - // Extract image data - std::vector imageData; - Common::image::extractImageData(processedImage, imageData); - - // Execute prediction - auto orientation = registrationPipeline_->predict(imageData, {1, 1, 64, 64, 64}); - - // Get new origin and direction - ImageType::Pointer originalImage = inputImage; // Use original image to get new origin and direction - auto newOriginAndDirection = registrationPipeline_->getNewOriginAndDirection( - processedImage, originalImage, - orientation["ac"], orientation["pa"], orientation["is"]); - - ImageType::PointType newOrigin = std::get<0>(newOriginAndDirection); - ImageType::DirectionType newDirection = std::get<1>(newOriginAndDirection); - - originalImage->SetDirection(newDirection); - originalImage->SetOrigin(newOrigin); - - return originalImage; -} - -ImageType::Pointer RigidVoxelMorphNormalizer::performVoxelMorphWarping(ImageType::Pointer rigidImage) { - // Resample to template space - ImageType::Pointer paddedImage = Common::image::resampleToMatch(paddedTemplate_, rigidImage); - - // Preprocessing - paddedImage = registrationPipeline_->preprocessVoxelMorph(paddedImage); - saveDebugImage(paddedImage, "elastic_preprocessed"); - - // Prepare data - std::vector paddedImageData, paddedTemplateData, paddedOriginalData; - Common::image::extractImageData(paddedImage, paddedImageData); - Common::image::extractImageData(paddedTemplate_, paddedTemplateData); - - // Load original resampled image - ImageType::Pointer paddedOriginalImage = Common::image::resampleToMatch(paddedTemplate_, rigidImage); - Common::image::extractImageData(paddedOriginalImage, paddedOriginalData); - - // Execute VoxelMorph prediction - auto warpedImageData = registrationPipeline_->predictVoxelMorph( - paddedOriginalData, paddedImageData, paddedTemplateData); - - // Create output image - ImageType::Pointer warpedImage = Common::image::createImageFromVector( - warpedImageData["warped"], paddedImage->GetLargestPossibleRegion().GetSize()); - - warpedImage->SetDirection(paddedTemplate_->GetDirection()); - warpedImage->SetOrigin(paddedTemplate_->GetOrigin()); - warpedImage->SetSpacing(paddedTemplate_->GetSpacing()); - - // Crop to MNI space - return cropMNI(warpedImage); -} - -ImageType::Pointer RigidVoxelMorphNormalizer::cropMNI(ImageType::Pointer image) { - ImageType::RegionType cropRegion; - ImageType::RegionType::IndexType start; - start[0] = config_->getInt("processing.crop_mni.start_x", 8); - start[1] = config_->getInt("processing.crop_mni.start_y", 16); - start[2] = config_->getInt("processing.crop_mni.start_z", 8); - - ImageType::RegionType::SizeType size; - size[0] = config_->getInt("processing.crop_mni.size_x", 79); - size[1] = config_->getInt("processing.crop_mni.size_y", 95); - size[2] = config_->getInt("processing.crop_mni.size_z", 79); - - cropRegion.SetSize(size); - cropRegion.SetIndex(start); - - using FilterType = itk::RegionOfInterestImageFilter; - FilterType::Pointer filter = FilterType::New(); - filter->SetRegionOfInterest(cropRegion); - filter->SetInput(image); - filter->Update(); - - return filter->GetOutput(); -} - std::string RigidVoxelMorphNormalizer::getName() const { return "RigidVoxelMorph"; } -bool RigidVoxelMorphNormalizer::isSupported(const std::string& modality) const { - // Support all modalities +bool RigidVoxelMorphNormalizer::isSupported(const std::string&) const { return true; } -void RigidVoxelMorphNormalizer::setDebugMode(bool enable, const std::string& basePath) { +void RigidVoxelMorphNormalizer::setDebugMode(bool enable, + const std::string& basePath) { debugMode_ = enable; debugBasePath_ = basePath; + if (rigidNormalizer_) { + rigidNormalizer_->setDebugMode(enable, basePath); + } + if (voxelMorphNormalizer_) { + voxelMorphNormalizer_->setDebugMode(enable, basePath); + } } -void RigidVoxelMorphNormalizer::saveDebugImage(ImageType::Pointer image, const std::string& suffix) { - if (!debugMode_ || debugBasePath_.empty()) return; - std::string filename = debugBasePath_ + "_" + suffix + ".nii"; - Common::nifti::saveImage(image, filename); +RigidAlignmentNormalizer& RigidVoxelMorphNormalizer::rigidNormalizer() { + if (!rigidNormalizer_) { + rigidNormalizer_ = std::make_unique(config_); + rigidNormalizer_->setDebugMode(debugMode_, debugBasePath_); + } + return *rigidNormalizer_; +} + +VoxelMorphNormalizer& RigidVoxelMorphNormalizer::voxelMorphNormalizer() { + if (!voxelMorphNormalizer_) { + voxelMorphNormalizer_ = std::make_unique(config_); + voxelMorphNormalizer_->setDebugMode(debugMode_, debugBasePath_); + } + return *voxelMorphNormalizer_; } diff --git a/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.h b/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.h index a48e703..f0d9001 100644 --- a/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.h +++ b/localizer/src/core/normalizers/RigidVoxelMorphNormalizer.h @@ -1,7 +1,9 @@ #pragma once #include "../interfaces/ISpatialNormalizer.h" #include "../interfaces/IConfiguration.h" -#include "RegistrationPipeline.h" // Use new registration pipeline +#include "RigidAlignmentNormalizer.h" +#include "VoxelMorphNormalizer.h" +#include /** * @brief Spatial normalizer based on rigid registration and VoxelMorph @@ -16,6 +18,8 @@ class RigidVoxelMorphNormalizer : public ISpatialNormalizer { bool isSupported(const std::string& modality) const override; // Extended functionality + ImageType::Pointer normalizeRigidOnly(ImageType::Pointer inputImage); + ImageType::Pointer normalizeIterativeRigidOnly(ImageType::Pointer inputImage, int maxIter = 5, float threshold = 2.0f); ImageType::Pointer normalizeIterative(ImageType::Pointer inputImage, int maxIter = 5, float threshold = 2.0f); ImageType::Pointer normalizeManualFOV(ImageType::Pointer inputImage); @@ -34,17 +38,13 @@ class RigidVoxelMorphNormalizer : public ISpatialNormalizer { private: ConfigurationPtr config_; - std::unique_ptr registrationPipeline_; - ImageType::Pointer paddedTemplate_; + std::unique_ptr rigidNormalizer_; + std::unique_ptr voxelMorphNormalizer_; // Debug parameters bool debugMode_ = false; std::string debugBasePath_ = ""; - - void initializeModel(); - ImageType::Pointer cropMNI(ImageType::Pointer image); - ImageType::Pointer performRigidAlignment(ImageType::Pointer inputImage, bool resampleFirst = false); - ImageType::Pointer performVoxelMorphWarping(ImageType::Pointer rigidImage); - void saveDebugImage(ImageType::Pointer image, const std::string& suffix); -}; + RigidAlignmentNormalizer& rigidNormalizer(); + VoxelMorphNormalizer& voxelMorphNormalizer(); +}; diff --git a/localizer/src/core/normalizers/VoxelMorphNormalizer.cpp b/localizer/src/core/normalizers/VoxelMorphNormalizer.cpp new file mode 100644 index 0000000..3641fff --- /dev/null +++ b/localizer/src/core/normalizers/VoxelMorphNormalizer.cpp @@ -0,0 +1,84 @@ +#include "VoxelMorphNormalizer.h" + +#include "../common/Common.h" +#include "../preprocessing/ImagePreprocessor.h" + +#include +#include +#include + +VoxelMorphNormalizer::VoxelMorphNormalizer(ConfigurationPtr config) + : config_(config) { + if (!config_) { + throw std::invalid_argument("VoxelMorphNormalizer requires configuration"); + } + + nonlinearEngine_ = std::make_unique( + config_->getModelPath("affine_voxelmorph")); + paddedTemplate_ = Common::nifti::loadImage(config_->getTemplatePath("padded")); +} + +ImageType::Pointer VoxelMorphNormalizer::normalize(ImageType::Pointer rigidImage) { + ImageType::Pointer paddedImage = + Common::image::resampleToMatch(paddedTemplate_, rigidImage); + + paddedImage = ImagePreprocessor::preprocessForVoxelMorph(paddedImage); + saveDebugImage(paddedImage, "elastic_preprocessed"); + + std::vector paddedImageData, paddedTemplateData, paddedOriginalData; + Common::image::extractImageData(paddedImage, paddedImageData); + Common::image::extractImageData(paddedTemplate_, paddedTemplateData); + + ImageType::Pointer paddedOriginalImage = + Common::image::resampleToMatch(paddedTemplate_, rigidImage); + Common::image::extractImageData(paddedOriginalImage, paddedOriginalData); + + auto warpedImageData = nonlinearEngine_->predict( + paddedOriginalData, paddedImageData, paddedTemplateData); + + ImageType::Pointer warpedImage = Common::image::createImageFromVector( + warpedImageData["warped"], paddedImage->GetLargestPossibleRegion().GetSize()); + + warpedImage->SetDirection(paddedTemplate_->GetDirection()); + warpedImage->SetOrigin(paddedTemplate_->GetOrigin()); + warpedImage->SetSpacing(paddedTemplate_->GetSpacing()); + + return cropMNI(warpedImage); +} + +ImageType::Pointer VoxelMorphNormalizer::cropMNI(ImageType::Pointer image) { + ImageType::RegionType cropRegion; + ImageType::RegionType::IndexType start; + start[0] = config_->getInt("processing.crop_mni.start_x", 8); + start[1] = config_->getInt("processing.crop_mni.start_y", 16); + start[2] = config_->getInt("processing.crop_mni.start_z", 8); + + ImageType::RegionType::SizeType size; + size[0] = config_->getInt("processing.crop_mni.size_x", 79); + size[1] = config_->getInt("processing.crop_mni.size_y", 95); + size[2] = config_->getInt("processing.crop_mni.size_z", 79); + + cropRegion.SetSize(size); + cropRegion.SetIndex(start); + + using FilterType = itk::RegionOfInterestImageFilter; + FilterType::Pointer filter = FilterType::New(); + filter->SetRegionOfInterest(cropRegion); + filter->SetInput(image); + filter->Update(); + + return filter->GetOutput(); +} + +void VoxelMorphNormalizer::setDebugMode(bool enable, const std::string& basePath) { + debugMode_ = enable; + debugBasePath_ = basePath; +} + +void VoxelMorphNormalizer::saveDebugImage(ImageType::Pointer image, + const std::string& suffix) { + if (!debugMode_ || debugBasePath_.empty()) { + return; + } + Common::nifti::saveImage(image, debugBasePath_ + "_" + suffix + ".nii"); +} diff --git a/localizer/src/core/normalizers/VoxelMorphNormalizer.h b/localizer/src/core/normalizers/VoxelMorphNormalizer.h new file mode 100644 index 0000000..a737795 --- /dev/null +++ b/localizer/src/core/normalizers/VoxelMorphNormalizer.h @@ -0,0 +1,26 @@ +#pragma once + +#include "../common/ImageTypes.h" +#include "../interfaces/IConfiguration.h" +#include "NonlinearRegistrationEngine.h" + +#include +#include + +class VoxelMorphNormalizer { +public: + explicit VoxelMorphNormalizer(ConfigurationPtr config); + + ImageType::Pointer normalize(ImageType::Pointer rigidImage); + void setDebugMode(bool enable, const std::string& basePath = ""); + +private: + ConfigurationPtr config_; + std::unique_ptr nonlinearEngine_; + ImageType::Pointer paddedTemplate_; + bool debugMode_ = false; + std::string debugBasePath_; + + ImageType::Pointer cropMNI(ImageType::Pointer image); + void saveDebugImage(ImageType::Pointer image, const std::string& suffix); +}; diff --git a/localizer/src/core/preprocessing/ImagePreprocessor.cpp b/localizer/src/core/preprocessing/ImagePreprocessor.cpp index a8c8345..bf105b6 100644 --- a/localizer/src/core/preprocessing/ImagePreprocessor.cpp +++ b/localizer/src/core/preprocessing/ImagePreprocessor.cpp @@ -24,10 +24,10 @@ ImageType::Pointer ImagePreprocessor::preprocessForVoxelMorph(ImageType::Pointer ImageType::Pointer ImagePreprocessor::clipIntensityPercentiles(ImageType::Pointer image, double lowerPercentile, double upperPercentile) { - auto sortedVoxelValue = getSortedPixelValues(image); + auto voxelValues = getPixelValues(image); - double lowerValue = getPercentileValue(sortedVoxelValue, lowerPercentile); - double upperValue = getPercentileValue(sortedVoxelValue, upperPercentile); + double lowerValue = getPercentileValue(voxelValues, lowerPercentile); + double upperValue = getPercentileValue(voxelValues, upperPercentile); using IntensityWindowingImageFilterType = itk::IntensityWindowingImageFilter; @@ -173,22 +173,27 @@ ImageType::Pointer ImagePreprocessor::resizeImage(ImageType::Pointer image, return resampler->GetOutput(); } -std::vector ImagePreprocessor::getSortedPixelValues(ImageType::Pointer image) { - itk::ImageRegionIterator it(image, image->GetRequestedRegion()); - std::vector pixelValues; +std::vector ImagePreprocessor::getPixelValues(ImageType::Pointer image) { + ImageType::RegionType region = image->GetRequestedRegion(); + itk::ImageRegionIterator it(image, region); + std::vector pixelValues; + pixelValues.reserve(region.GetNumberOfPixels()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { pixelValues.push_back(it.Get()); } - std::sort(pixelValues.begin(), pixelValues.end()); return pixelValues; } -double ImagePreprocessor::getPercentileValue(const std::vector& sortedValues, - double percentile) { - size_t index = static_cast(percentile * (sortedValues.size() - 1)); - return sortedValues[index]; -} +double ImagePreprocessor::getPercentileValue(std::vector& values, + double percentile) { + if (values.empty()) { + return 0.0; + } + size_t index = static_cast(percentile * (values.size() - 1)); + std::nth_element(values.begin(), values.begin() + index, values.end()); + return values[index]; +} diff --git a/localizer/src/core/preprocessing/ImagePreprocessor.h b/localizer/src/core/preprocessing/ImagePreprocessor.h index d1b19c5..7a42199 100644 --- a/localizer/src/core/preprocessing/ImagePreprocessor.h +++ b/localizer/src/core/preprocessing/ImagePreprocessor.h @@ -40,8 +40,7 @@ class ImagePreprocessor { const ImageType::SizeType& newSize); // Helper functions - static std::vector getSortedPixelValues(ImageType::Pointer image); - static double getPercentileValue(const std::vector& sortedValues, double percentile); + static std::vector getPixelValues(ImageType::Pointer image); + static double getPercentileValue(std::vector& values, double percentile); }; - diff --git a/localizer/src/core/services/SpatialNormalizationService.cpp b/localizer/src/core/services/SpatialNormalizationService.cpp index d40b913..12758ad 100644 --- a/localizer/src/core/services/SpatialNormalizationService.cpp +++ b/localizer/src/core/services/SpatialNormalizationService.cpp @@ -28,7 +28,15 @@ SpatialNormalizationOutput SpatialNormalizationService::normalize(const SpatialN normalizer->setDebugMode(true, request.options.debugOutputBasePath); } - if (request.options.useManualFOV) { + if (request.options.rigidOnly) { + output.rigidAlignedImage = request.options.useIterativeRigid + ? normalizer->normalizeIterativeRigidOnly( + inputImage, + request.options.maxIterations, + request.options.convergenceThreshold) + : normalizer->normalizeRigidOnly(inputImage); + output.spatiallyNormalizedImage = output.rigidAlignedImage; + } else if (request.options.useManualFOV) { output.rigidAlignedImage = inputImage; output.spatiallyNormalizedImage = normalizer->normalizeManualFOV(inputImage); } else if (request.options.useIterativeRigid) { @@ -37,15 +45,11 @@ SpatialNormalizationOutput SpatialNormalizationService::normalize(const SpatialN request.options.maxIterations, request.options.convergenceThreshold); output.rigidAlignedImage = normResult.rigidAlignedImage; - output.spatiallyNormalizedImage = request.options.rigidOnly - ? normResult.rigidAlignedImage - : normResult.spatiallyNormalizedImage; + output.spatiallyNormalizedImage = normResult.spatiallyNormalizedImage; } else { auto normResult = normalizer->normalizeWithIntermediateResults(inputImage); output.rigidAlignedImage = normResult.rigidAlignedImage; - output.spatiallyNormalizedImage = request.options.rigidOnly - ? normResult.rigidAlignedImage - : normResult.spatiallyNormalizedImage; + output.spatiallyNormalizedImage = normResult.spatiallyNormalizedImage; } }