From 725bb4d37280aac8b141c1abfe35052b6b1030f4 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Sat, 3 Oct 2026 10:22:17 +0200 Subject: [PATCH] STYLE: Mark the sources for the ITK style checks and format them The clang-format linter action now fails when no tracked file carries hooks.style=clangformat; without a .gitattributes it had been passing without checking a file. Mark the C++ sources as ITK marks its own and apply clang-format 19.1.7 with the module's .clang-format: whitespace and line breaks only. --- .gitattributes | 6 + examples/ImpactMetricExample.cxx | 18 +- include/ImpactLoss.h | 6 +- include/itkImageToFeaturesMap.hxx | 5 +- include/itkImageToTensorFilter.h | 24 +- include/itkImpactCoarseRegistration.hxx | 14 +- include/itkImpactFineRegistration.h | 18 +- include/itkImpactFineRegistration.hxx | 67 ++-- include/itkImpactImageToImageMetricv4.h | 14 +- include/itkImpactImageToImageMetricv4.hxx | 105 +++--- ...ageMetricv4GetValueAndDerivativeThreader.h | 23 +- ...eMetricv4GetValueAndDerivativeThreader.hxx | 88 +++-- include/itkImpactModelConfiguration.h | 16 +- include/itkImpactModelConfigurationDetail.h | 2 +- include/itkImpactOnlineInference.h | 66 ++-- include/itkImpactPatchTiling.h | 56 +-- include/itkImpactTorchRegistrationHelpers.h | 35 +- include/itkInterpolateVectorImageFunction.h | 19 +- include/itkInterpolateVectorImageFunction.hxx | 8 +- include/itkTensorToImageFilter.h | 11 +- include/itkTensorToImageFilter.hxx | 23 +- src/itkImpactModelConfiguration.cxx | 18 +- test/itkImageToFeaturesMapTest.cxx | 4 +- test/itkImageToTensorFilterTest.cxx | 18 +- test/itkImpactBackendGTest.cxx | 318 +++++++++--------- .../itkInterpolateVectorImageFunctionTest.cxx | 4 +- test/itkTensorToImageFilterTest.cxx | 4 +- 27 files changed, 500 insertions(+), 490 deletions(-) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..46c7b7b --- /dev/null +++ b/.gitattributes @@ -0,0 +1,6 @@ +# The sources the ITK style checks read (KWStyle, clang-format), as ITK marks its own. +[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat + +*.h our-c-style +*.cxx our-c-style +*.hxx our-c-style diff --git a/examples/ImpactMetricExample.cxx b/examples/ImpactMetricExample.cxx index 7524051..cd0c819 100644 --- a/examples/ImpactMetricExample.cxx +++ b/examples/ImpactMetricExample.cxx @@ -49,11 +49,11 @@ main(int argc, char * argv[]) std::cerr << " device: \"cpu\" (default), \"cuda\", \"cuda:0\", ...\n"; return EXIT_FAILURE; } - const char * const modelPath = argv[1]; - const char * const fixedPath = argv[2]; - const char * const movingPath = argv[3]; - const char * const outputPath = argv[4]; - const std::string device = (argc > 5) ? argv[5] : "cpu"; + const char * const modelPath = argv[1]; + const char * const fixedPath = argv[2]; + const char * const movingPath = argv[3]; + const char * const outputPath = argv[4]; + const std::string device = (argc > 5) ? argv[5] : "cpu"; constexpr unsigned int Dimension = 3; using PixelType = float; @@ -84,14 +84,14 @@ main(int argc, char * argv[]) // --- 2. metric: register moving onto fixed by comparing anatomical features ------- using MetricType = itk::ImpactImageToImageMetricv4; - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector models{ config }; metric->SetModelsConfiguration(models); - metric->SetDistance({ "L2" }); // per-layer loss: L1, L2, NCC, Cosine, Dice, ... + metric->SetDistance({ "L2" }); // per-layer loss: L1, L2, NCC, Cosine, Dice, ... metric->SetLayersWeight({ 1.0f }); - metric->SetSubsetFeatures({ 4 }); // random channel subset for speed (0 = all) + metric->SetSubsetFeatures({ 4 }); // random channel subset for speed (0 = all) metric->SetPCA({ 0 }); - metric->SetMode("Static"); // "Static" (precomputed features) or "Jacobian" + metric->SetMode("Static"); // "Static" (precomputed features) or "Jacobian" metric->SetDevice(device); using TransformType = itk::TranslationTransform; diff --git a/include/ImpactLoss.h b/include/ImpactLoss.h index 0489490..817b684 100644 --- a/include/ImpactLoss.h +++ b/include/ImpactLoss.h @@ -537,8 +537,7 @@ class L1Cosine : public Loss } }; -inline RegisterLoss L1Cosine_reg( - "L1Cosine"); +inline RegisterLoss L1Cosine_reg("L1Cosine"); /** * \class Cosine @@ -763,8 +762,7 @@ class NCC : public Loss if (N <= 0) return 0.0; torch::Tensor u = this->m_sfm - (this->m_sf * this->m_sm / N); - torch::Tensor v = - Denominator(Variance(this->m_sff, this->m_sf, N), Variance(this->m_smm, this->m_sm, N)); + torch::Tensor v = Denominator(Variance(this->m_sff, this->m_sf, N), Variance(this->m_smm, this->m_sm, N)); return 1.0 - (u / v).mean().item(); } diff --git a/include/itkImageToFeaturesMap.hxx b/include/itkImageToFeaturesMap.hxx index 721624d..19171c6 100644 --- a/include/itkImageToFeaturesMap.hxx +++ b/include/itkImageToFeaturesMap.hxx @@ -53,14 +53,13 @@ ImageToFeaturesMap::ImageToFeaturesMap() template void -ImageToFeaturesMap -::PrintSelf(std::ostream & os, Indent indent) const +ImageToFeaturesMap::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); } template -void +void ImageToFeaturesMap::AddInput(const TInputImage * input) { if (!m_Interpolator) diff --git a/include/itkImageToTensorFilter.h b/include/itkImageToTensorFilter.h index eb333ef..d431187 100644 --- a/include/itkImageToTensorFilter.h +++ b/include/itkImageToTensorFilter.h @@ -54,10 +54,10 @@ class ITK_TEMPLATE_EXPORT ImageToTensorFilter : public ProcessObject using InputImagePointer = typename InputImageType::Pointer; using InputImageConstPointer = typename InputImageType::ConstPointer; - using TensorType = torch::Tensor; + using TensorType = torch::Tensor; using TensorHandleType = std::shared_ptr; using TensorDataObject = itk::SimpleDataObjectDecorator; - + using InterpolatorPointer = typename TInterpolator::Pointer; using InputImagePixelType = typename InputImageType::PixelType; @@ -89,9 +89,14 @@ class ITK_TEMPLATE_EXPORT ImageToTensorFilter : public ProcessObject itkSetMacro(OutputSpacing, InputSpacingType); itkGetConstReferenceMacro(OutputSpacing, InputSpacingType); itkSetVectorMacro(OutputSpacing, const float, ImageDimension); - - void SetInterpolator(typename TInterpolator::Pointer interp) { m_Interpolator = interp; } - void SetTransform(std::function fct) + + void + SetInterpolator(typename TInterpolator::Pointer interp) + { + m_Interpolator = interp; + } + void + SetTransform(std::function fct) { m_Transform = fct; } @@ -141,15 +146,16 @@ class ITK_TEMPLATE_EXPORT ImageToTensorFilter : public ProcessObject void PrintSelf(std::ostream & os, Indent indent) const override; - + void VerifyPreconditions() const override; - void GenerateData() override; + void + GenerateData() override; private: - InterpolatorPointer m_Interpolator; - InputSpacingType m_OutputSpacing{ MakeFilled(1.0) }; + InterpolatorPointer m_Interpolator; + InputSpacingType m_OutputSpacing{ MakeFilled(1.0) }; std::function m_Transform; }; diff --git a/include/itkImpactCoarseRegistration.hxx b/include/itkImpactCoarseRegistration.hxx index a7fe085..04d784c 100644 --- a/include/itkImpactCoarseRegistration.hxx +++ b/include/itkImpactCoarseRegistration.hxx @@ -276,8 +276,7 @@ ImpactCoarseRegistration::GetDisplacementField() -> D template auto -ImpactCoarseRegistration::GetDisplacementFieldTransform() - -> DisplacementFieldTransformType * +ImpactCoarseRegistration::GetDisplacementFieldTransform() -> DisplacementFieldTransformType * { return m_DisplacementFieldTransform.GetPointer(); } @@ -421,7 +420,9 @@ ImpactCoarseRegistration::GenerateData() if (kept > 0 && kept < fixedLayers[l].size(1)) { const torch::Tensor channels = - torch::randperm(fixedLayers[l].size(1), torch::TensorOptions().dtype(torch::kLong)).narrow(0, 0, kept).to(device); + torch::randperm(fixedLayers[l].size(1), torch::TensorOptions().dtype(torch::kLong)) + .narrow(0, 0, kept) + .to(device); fixedLayers[l] = fixedLayers[l].index_select(1, channels).contiguous(); movingLayers[l] = movingLayers[l].index_select(1, channels).contiguous(); } @@ -741,8 +742,7 @@ ImpactCoarseRegistration::GenerateData() torch::Tensor f2 = (dispBack / scaleT).flip(1); // Identity sampling grid at coarse resolution, channel-first {1, Dim, coarse...}, x,y,z. - torch::Tensor idAffine = - torch::eye(ImageDimension, torch::TensorOptions().dtype(torch::kFloat32).device(device)); + torch::Tensor idAffine = torch::eye(ImageDimension, torch::TensorOptions().dtype(torch::kFloat32).device(device)); idAffine = torch::cat({ idAffine, torch::zeros({ static_cast(ImageDimension), 1 }, idAffine.options()) }, 1) .unsqueeze(0); @@ -753,7 +753,7 @@ ImpactCoarseRegistration::GenerateData() { gridSize.push_back(c); } - torch::Tensor idGrid = torch::affine_grid_generator(idAffine, gridSize, /*align_corners=*/true); + torch::Tensor idGrid = torch::affine_grid_generator(idAffine, gridSize, /*align_corners=*/true); std::vector toChannelFirst; toChannelFirst.push_back(0); toChannelFirst.push_back(idGrid.dim() - 1); @@ -803,7 +803,7 @@ ImpactCoarseRegistration::GenerateData() std::vector componentShape(ImageDimension + 2, 1); // {1, Dim, 1, ...}: one factor per component componentShape[1] = static_cast(ImageDimension); const torch::Tensor cellsToVoxels = torch::tensor(cellVoxels, torch::kLong).to(disp.options()).view(componentShape); - torch::Tensor dispFull; + torch::Tensor dispFull; if constexpr (ImageDimension == 3) dispFull = F::interpolate( disp * cellsToVoxels, F::InterpolateFuncOptions().size(exactSize).mode(torch::kTrilinear).align_corners(false)); diff --git a/include/itkImpactFineRegistration.h b/include/itkImpactFineRegistration.h index ab45929..5c9e01e 100644 --- a/include/itkImpactFineRegistration.h +++ b/include/itkImpactFineRegistration.h @@ -287,18 +287,18 @@ class ITK_TEMPLATE_EXPORT ImpactFineRegistration GenerateData() override; private: - typename FixedImageType::ConstPointer m_FixedImage{ nullptr }; - typename MovingImageType::ConstPointer m_MovingImage{ nullptr }; - typename MaskImageType::ConstPointer m_FixedMask{ nullptr }; - typename MaskImageType::ConstPointer m_MovingMask{ nullptr }; - typename DisplacementFieldType::Pointer m_InitialDisplacementField{ nullptr }; + typename FixedImageType::ConstPointer m_FixedImage{ nullptr }; + typename MovingImageType::ConstPointer m_MovingImage{ nullptr }; + typename MaskImageType::ConstPointer m_FixedMask{ nullptr }; + typename MaskImageType::ConstPointer m_MovingMask{ nullptr }; + typename DisplacementFieldType::Pointer m_InitialDisplacementField{ nullptr }; std::vector m_FixedModelsConfiguration; std::vector m_MovingModelsConfiguration; - std::vector m_Distance; - std::vector m_LayersWeight; - std::vector m_SubsetFeatures; - std::vector m_PCA; + std::vector m_Distance; + std::vector m_LayersWeight; + std::vector m_SubsetFeatures; + std::vector m_PCA; std::string m_Device{ "cpu" }; unsigned int m_Seed{ 0 }; diff --git a/include/itkImpactFineRegistration.hxx b/include/itkImpactFineRegistration.hxx index a2c4d08..175ed85 100644 --- a/include/itkImpactFineRegistration.hxx +++ b/include/itkImpactFineRegistration.hxx @@ -181,8 +181,7 @@ ImpactFineRegistration::GetDisplacementField() -> Dis template auto -ImpactFineRegistration::GetDisplacementFieldTransform() - -> DisplacementFieldTransformType * +ImpactFineRegistration::GetDisplacementFieldTransform() -> DisplacementFieldTransformType * { return m_DisplacementFieldTransform.GetPointer(); } @@ -311,25 +310,27 @@ ImpactFineRegistration::GenerateData() if (coarseSpatial != spatial) { if constexpr (ImageDimension == 3) - initField = torch::nn::functional::interpolate( - initField, - torch::nn::functional::InterpolateFuncOptions().size(coarseSpatial).mode(torch::kTrilinear).align_corners(true)); + initField = torch::nn::functional::interpolate(initField, + torch::nn::functional::InterpolateFuncOptions() + .size(coarseSpatial) + .mode(torch::kTrilinear) + .align_corners(true)); else - initField = torch::nn::functional::interpolate( - initField, - torch::nn::functional::InterpolateFuncOptions().size(coarseSpatial).mode(torch::kBilinear).align_corners(true)); + initField = torch::nn::functional::interpolate(initField, + torch::nn::functional::InterpolateFuncOptions() + .size(coarseSpatial) + .mode(torch::kBilinear) + .align_corners(true)); } theta = initField.set_requires_grad(true); } else { - theta = - torch::zeros(fieldShape, torch::TensorOptions().dtype(torch::kFloat32).device(device).requires_grad(true)); + theta = torch::zeros(fieldShape, torch::TensorOptions().dtype(torch::kFloat32).device(device).requires_grad(true)); } // ---- 3. Base identity sampling grid (normalized [-1,1], last-dim order x,y,z), align_corners=true. ---- - torch::Tensor idAffine = - torch::eye(ImageDimension, torch::TensorOptions().dtype(torch::kFloat32).device(device)); + torch::Tensor idAffine = torch::eye(ImageDimension, torch::TensorOptions().dtype(torch::kFloat32).device(device)); idAffine = torch::cat({ idAffine, torch::zeros({ static_cast(ImageDimension), 1 }, idAffine.options()) }, 1) .unsqueeze(0); // {1, N, N+1} std::vector gridSize; @@ -339,8 +340,7 @@ ImpactFineRegistration::GenerateData() { gridSize.push_back(s); } - torch::Tensor grid0 = - torch::affine_grid_generator(idAffine, gridSize, /*align_corners=*/true); // {1, z,y,x, N} + torch::Tensor grid0 = torch::affine_grid_generator(idAffine, gridSize, /*align_corners=*/true); // {1, z,y,x, N} // Per-component normalization in (z, y, x) order: the field units (s_min) per grid_sample unit, (size-1)/2 voxels // (exact for align_corners=true) of s_a / s_min units each. @@ -394,7 +394,8 @@ ImpactFineRegistration::GenerateData() return field; } if constexpr (ImageDimension == 3) - return F::interpolate(field, F::InterpolateFuncOptions().size(target).mode(torch::kTrilinear).align_corners(true)); + return F::interpolate(field, + F::InterpolateFuncOptions().size(target).mode(torch::kTrilinear).align_corners(true)); else return F::interpolate(field, F::InterpolateFuncOptions().size(target).mode(torch::kBilinear).align_corners(true)); }; @@ -510,15 +511,15 @@ ImpactFineRegistration::GenerateData() // Feature-mode setup: extract the fixed/moving feature layers (constants, not differentiated // through), optionally PCA-reduce them (fit on fixed), and build one loss per kept layer. // Intensity mode skips this and compares raw voxels. - std::vector fixedLayers; - std::vector movingLayers; - std::vector fixedLayersOnline; // "Jacobian" mode: fixed features, pooled - std::vector pcaBasis; // per kept layer; undefined entry = no PCA - std::vector> losses; - std::vector layerWeights; + std::vector fixedLayers; + std::vector movingLayers; + std::vector fixedLayersOnline; // "Jacobian" mode: fixed features, pooled + std::vector pcaBasis; // per kept layer; undefined entry = no PCA + std::vector> losses; + std::vector layerWeights; // SubsetFeatures: per kept layer, that many of its channels drawn at random at every iteration (0 = all). std::vector subsets; // per layer; an undefined entry keeps every channel - auto drawSubsets = [&](const std::vector & layers) { + auto drawSubsets = [&](const std::vector & layers) { subsets.assign(layers.size(), torch::Tensor()); for (size_t l = 0; l < layers.size(); ++l) { @@ -731,8 +732,7 @@ ImpactFineRegistration::GenerateData() } for (size_t l = 0; l < layerCount; ++l) { - const std::string name = - m_Distance.empty() ? std::string("L2") : m_Distance[std::min(l, m_Distance.size() - 1)]; + const std::string name = m_Distance.empty() ? std::string("L2") : m_Distance[std::min(l, m_Distance.size() - 1)]; losses.push_back(Impact::LossFactory::Instance().Create(name)); layerWeights.push_back(l < m_LayersWeight.size() ? m_LayersWeight[l] : 1.0f); if (sampled && losses.back()->IsSpatial()) @@ -951,10 +951,8 @@ ImpactFineRegistration::GenerateData() } // ---- 4. Adam loop (entirely on device; no host copies). ---- - torch::optim::Adam optimizer({ theta }, - torch::optim::AdamOptions(m_LearningRate) - .betas(std::make_tuple(m_Beta1, m_Beta2)) - .eps(m_Epsilon)); + torch::optim::Adam optimizer( + { theta }, torch::optim::AdamOptions(m_LearningRate).betas(std::make_tuple(m_Beta1, m_Beta2)).eps(m_Epsilon)); m_MetricValuesPerIteration.clear(); m_MetricValuesPerIteration.reserve(m_NumberOfIterations); @@ -962,7 +960,7 @@ ImpactFineRegistration::GenerateData() // Every layer's loss divided by its value at the stage's first iteration, so each starts at 1 and // LayersWeight weighs comparable quantities (see Impact::LossNormalization). Impact::LossNormalization normalization; - auto normalized = [&](size_t l, const torch::Tensor & value) -> torch::Tensor { + auto normalized = [&](size_t l, const torch::Tensor & value) -> torch::Tensor { if (!m_NormalizeLosses) { return value; @@ -1032,7 +1030,7 @@ ImpactFineRegistration::GenerateData() const torch::Tensor smoothedControl = smoothControl(theta); // graph to theta torch::Tensor smoothed = smoothedControl.detach().requires_grad_(true); smoothed.mutable_grad() = torch::zeros_like(smoothed); - int64_t voxels = 1; + int64_t voxels = 1; for (const int64_t extent : spatial) { voxels *= extent; @@ -1439,9 +1437,11 @@ ImpactFineRegistration::GenerateData() if (cur != tgt) { if constexpr (ImageDimension == 3) - mll = F::interpolate(mll, F::InterpolateFuncOptions().size(tgt).mode(torch::kTrilinear).align_corners(true)); + mll = + F::interpolate(mll, F::InterpolateFuncOptions().size(tgt).mode(torch::kTrilinear).align_corners(true)); else - mll = F::interpolate(mll, F::InterpolateFuncOptions().size(tgt).mode(torch::kBilinear).align_corners(true)); + mll = + F::interpolate(mll, F::InterpolateFuncOptions().size(tgt).mode(torch::kBilinear).align_corners(true)); } torch::Tensor counted; if (masked) @@ -1575,7 +1575,8 @@ ImpactFineRegistration::GenerateData() m_WarpedMovingImage->SetDirection(m_FixedImage->GetDirection()); m_WarpedMovingImage->Allocate(); - ImageRegionIteratorWithIndex wit(m_WarpedMovingImage, m_WarpedMovingImage->GetLargestPossibleRegion()); + ImageRegionIteratorWithIndex wit(m_WarpedMovingImage, + m_WarpedMovingImage->GetLargestPossibleRegion()); for (wit.GoToBegin(); !wit.IsAtEnd(); ++wit) { const auto idx = wit.GetIndex(); diff --git a/include/itkImpactImageToImageMetricv4.h b/include/itkImpactImageToImageMetricv4.h index 2519fd8..022bf6d 100644 --- a/include/itkImpactImageToImageMetricv4.h +++ b/include/itkImpactImageToImageMetricv4.h @@ -303,18 +303,18 @@ class ITK_TEMPLATE_EXPORT ImpactImageToImageMetricv4 int m_FeaturesMapUpdateInterval{ 0 }; /** Derivative evaluations since Initialize(), the metric's stand-in for an iteration count. * Mutable because GetValueAndDerivative() is const, as the base metric declares it. */ - mutable unsigned long m_CurrentIteration{ 0 }; - std::string m_Mode; - std::string m_FeatureMapsPath; - bool m_NormalizeLosses{ true }; + mutable unsigned long m_CurrentIteration{ 0 }; + std::string m_Mode; + std::string m_FeatureMapsPath; + bool m_NormalizeLosses{ true }; /** Latched at the first evaluation after Initialize(); mutable because the evaluation is const. */ mutable Impact::LossNormalization m_LossNormalization; - std::string m_Device = "cpu"; + std::string m_Device = "cpu"; // Zero means "seed from the clock", which is what the per-work-unit generator does with it. // It must have a value even when the user never calls SetSeed(): it is read on every // evaluation, and an indeterminate one would make the metric unreproducible at random. - unsigned int m_Seed{ 0 }; - unsigned int m_BatchSize{ 64 }; + unsigned int m_Seed{ 0 }; + unsigned int m_BatchSize{ 64 }; std::vector> m_features_indexes; }; diff --git a/include/itkImpactImageToImageMetricv4.hxx b/include/itkImpactImageToImageMetricv4.hxx index 143f807..650b34d 100644 --- a/include/itkImpactImageToImageMetricv4.hxx +++ b/include/itkImpactImageToImageMetricv4.hxx @@ -25,7 +25,7 @@ #include "itkImageToFeaturesMapInternals.h" #include "itkImpactTorchRegistrationHelpers.h" // Impact::PythonGilReleaseGuard #include "itkImpactModelConfigurationDetail.h" // ComputeImageMetadata -#include "itkImpactBatchBudget.h" // ConfigureBatchSize +#include "itkImpactBatchBudget.h" // ConfigureBatchSize #include namespace itk @@ -46,7 +46,8 @@ struct ImpactImageToImageMetricv4>; + itk::InterpolateVectorImageFunction>; /** A feature map image together with its interpolator. */ struct FeaturesMap @@ -80,20 +81,16 @@ template -ImpactImageToImageMetricv4::ImpactImageToImageMetricv4() +ImpactImageToImageMetricv4:: + ImpactImageToImageMetricv4() : m_Internals(std::make_shared()) { // We have our own GetValueAndDerivativeThreader's that we want // ImageToImageMetricv4 to use. - using ImpactDenseGetValueAndDerivativeThreaderType = - ImpactImageToImageMetricv4GetValueAndDerivativeThreader< - ThreadedImageRegionPartitioner, - Superclass, - Self>; + using ImpactDenseGetValueAndDerivativeThreaderType = ImpactImageToImageMetricv4GetValueAndDerivativeThreader< + ThreadedImageRegionPartitioner, + Superclass, + Self>; using ImpactSparseGetValueAndDerivativeThreaderType = ImpactImageToImageMetricv4GetValueAndDerivativeThreader; this->m_DenseGetValueAndDerivativeThreader = ImpactDenseGetValueAndDerivativeThreaderType::New(); @@ -106,11 +103,8 @@ template void -ImpactImageToImageMetricv4::PrintSelf(std::ostream & os, Indent indent) const +ImpactImageToImageMetricv4:: + PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << "\nNormalizeLosses: " << (this->m_NormalizeLosses ? "on" : "off"); @@ -126,14 +120,15 @@ ImpactImageToImageMetricv4(this->GetSubsetFeatures()) - << "\nPCA: " << GetStringFromVector(this->GetPCA()) - << "\nLayersWeight: " << GetStringFromVector(this->GetLayersWeight()) - << "\nDistance: " << GetStringFromVector(this->GetDistance()) << "\nMode: " << this->GetMode() - << "\nDevice: " << this->GetDevice() << "\nSeed: " << this->GetSeed(); + << "\nPCA: " << GetStringFromVector(this->GetPCA()) + << "\nLayersWeight: " << GetStringFromVector(this->GetLayersWeight()) + << "\nDistance: " << GetStringFromVector(this->GetDistance()) << "\nMode: " << this->GetMode() + << "\nDevice: " << this->GetDevice() << "\nSeed: " << this->GetSeed(); if (this->GetMode() == "Static") { - os << "\nFeaturesMapUpdateInterval: " << this->GetFeaturesMapUpdateInterval() << "\nFeatureMapsPath: " << this->GetFeatureMapsPath(); + os << "\nFeaturesMapUpdateInterval: " << this->GetFeaturesMapUpdateInterval() + << "\nFeatureMapsPath: " << this->GetFeatureMapsPath(); } else { @@ -159,11 +154,12 @@ ImpactImageToImageMetricv4; std::vector featuresMaps; - + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetSplineOrder(3); - for(unsigned int i = 0; i < modelsConfiguration.size(); ++i){ + for (unsigned int i = 0; i < modelsConfiguration.size(); ++i) + { typename ImageToFeaturesMapType::Pointer imageToFeaturesMap = ImageToFeaturesMapType::New(); imageToFeaturesMap->SetModelConfiguration(modelsConfiguration[i]); imageToFeaturesMap->SetInterpolator(interpolator); @@ -196,7 +192,7 @@ ImpactImageToImageMetricv4GetFeatureMapsPath().empty()) { typename WriterType::Pointer writer = WriterType::New(); - std::string filename; + std::string filename; for (int it = 0; it < modelsConfiguration[i].GetVoxelSize().size(); ++it) { @@ -209,7 +205,7 @@ ImpactImageToImageMetricv4SetFileName(this->GetFeatureMapsPath() + (fct ? "/Moving_" : "/Fixed_")+ filename + ".mha"); + writer->SetFileName(this->GetFeatureMapsPath() + (fct ? "/Moving_" : "/Fixed_") + filename + ".mha"); typename TensorToImageFilterType::Pointer tensorToImageFilter = TensorToImageFilterType::New(); tensorToImageFilter->SetTensor(GetTensorInput(*imageToFeaturesMap, 0)); tensorToImageFilter->SetReferenceImage(image); @@ -259,11 +255,8 @@ template void -ImpactImageToImageMetricv4::UpdateMovingFeaturesMaps() +ImpactImageToImageMetricv4:: + UpdateMovingFeaturesMaps() { // Only Static holds a precomputed map; the online modes extract per point and have nothing // to refresh. @@ -271,14 +264,11 @@ ImpactImageToImageMetricv4m_Internals->m_MovingFeaturesMaps = - GetFeaturesMaps( - this->m_MovingImage.GetPointer(), - m_MovingModelsConfiguration, - std::function( - [this](const typename TMovingImage::PointType & point) { - return this->GetTransform()->TransformPoint(point); - })); + this->m_Internals->m_MovingFeaturesMaps = GetFeaturesMaps( + this->m_MovingImage.GetPointer(), + m_MovingModelsConfiguration, + std::function( + [this](const typename TMovingImage::PointType & point) { return this->GetTransform()->TransformPoint(point); })); // Keep the transform this map was built through. The map holds the features of M(T_k(p)), so // the threader reads it at the residual point x + (T(x) - T_k(x)): x itself at the moment of // the refresh, which makes that moment exact, and a point that moves with the parameters @@ -294,12 +284,9 @@ template void -ImpactImageToImageMetricv4::GetValueAndDerivative(typename Superclass::MeasureType & value, - typename Superclass::DerivativeType & derivative) const +ImpactImageToImageMetricv4:: + GetValueAndDerivative(typename Superclass::MeasureType & value, + typename Superclass::DerivativeType & derivative) const { // Under a Python interpreter libtorch's autograd engine refuses to run on a thread holding the // GIL, and the SWIG-wrapped optimizer holds it across this call; the Jacobian mode backpropagates @@ -323,11 +310,9 @@ template void -ImpactImageToImageMetricv4::Initialize(){ +ImpactImageToImageMetricv4:: + Initialize() +{ Superclass::Initialize(); this->m_LossNormalization.Reset(); // each level starts at 1 again this->m_features_indexes.clear(); @@ -367,8 +352,8 @@ ImpactImageToImageMetricv4m_PCA.size() != m_FixedModelsConfiguration.size()) { - itkExceptionMacro("PCA has " << this->m_PCA.size() << " entries but there are " - << m_FixedModelsConfiguration.size() << " models; it carries one entry per model."); + itkExceptionMacro("PCA has " << this->m_PCA.size() << " entries but there are " << m_FixedModelsConfiguration.size() + << " models; it carries one entry per model."); } if (this->GetMode() == "Static") @@ -379,7 +364,8 @@ ImpactImageToImageMetricv4m_Internals->m_Principal_components.clear(); - fixedFeaturesMaps = GetFeaturesMaps(this->m_FixedImage.GetPointer(), m_FixedModelsConfiguration); + fixedFeaturesMaps = GetFeaturesMaps(this->m_FixedImage.GetPointer(), + m_FixedModelsConfiguration); // The fixed map has to exist first: it is what fits the PCA basis the moving map reuses. this->m_CurrentIteration = 0; this->UpdateMovingFeaturesMaps(); @@ -396,10 +382,9 @@ ImpactImageToImageMetricv4GetNumberOfComponentsPerPixel() != movingFeaturesMaps[i].m_FeaturesMap->GetNumberOfComponentsPerPixel()) { - itkExceptionMacro( - "Mismatch in number of components per feature map at layer " - << i << ": fixed = " << fixedFeaturesMaps[i].m_FeaturesMap->GetNumberOfComponentsPerPixel() - << ", moving = " << movingFeaturesMaps[i].m_FeaturesMap->GetNumberOfComponentsPerPixel()); + itkExceptionMacro("Mismatch in number of components per feature map at layer " + << i << ": fixed = " << fixedFeaturesMaps[i].m_FeaturesMap->GetNumberOfComponentsPerPixel() + << ", moving = " << movingFeaturesMaps[i].m_FeaturesMap->GetNumberOfComponentsPerPixel()); } } @@ -455,8 +440,7 @@ ImpactImageToImageMetricv4 GetSubsetOfFeatures(const std::vector & features_index, - std::mt19937 & randomGenerator, - int n) const; + std::vector + GetSubsetOfFeatures(const std::vector & features_index, std::mt19937 & randomGenerator, int n) const; /** * \brief Per-work-unit accumulator of loss values and gradients, one loss object per @@ -111,13 +110,16 @@ class ITK_TEMPLATE_EXPORT ImpactImageToImageMetricv4GetValueAndDerivativeThreade struct LossPerThreadStruct { std::vector> m_losses; - std::vector m_layersWeight; - SizeValueType m_numberOfPixelsCounted; - int m_nb_parameters; - std::mt19937 m_randomGenerator; + std::vector m_layersWeight; + SizeValueType m_numberOfPixelsCounted; + int m_nb_parameters; + std::mt19937 m_randomGenerator; void - init(std::vector distance_name, std::vector layersWeight, unsigned int seed, unsigned int nb_parameters) + init(std::vector distance_name, + std::vector layersWeight, + unsigned int seed, + unsigned int nb_parameters) { if (seed > 0) { @@ -226,9 +228,8 @@ class ITK_TEMPLATE_EXPORT ImpactImageToImageMetricv4GetValueAndDerivativeThreade /** The plane a patch is cut on, as a matrix whose COLUMNS are the patch axes expressed in * the image's own frame. */ - using PatchPlaneType = Matrix; + using PatchPlaneType = + Matrix; /** Walk this work unit's share of the domain. Online mode walks it here so it can run the * models on a batch of points at a time; Static mode is left to the framework's own diff --git a/include/itkImpactImageToImageMetricv4GetValueAndDerivativeThreader.hxx b/include/itkImpactImageToImageMetricv4GetValueAndDerivativeThreader.hxx index e41d6a9..85de8dd 100644 --- a/include/itkImpactImageToImageMetricv4GetValueAndDerivativeThreader.hxx +++ b/include/itkImpactImageToImageMetricv4GetValueAndDerivativeThreader.hxx @@ -35,19 +35,16 @@ namespace itk { template -ImpactImageToImageMetricv4GetValueAndDerivativeThreader< - TDomainPartitioner, - TImageToImageMetric, - TImpactMetric>::ImpactImageToImageMetricv4GetValueAndDerivativeThreader() +ImpactImageToImageMetricv4GetValueAndDerivativeThreader:: + ImpactImageToImageMetricv4GetValueAndDerivativeThreader() : m_LossThreadStruct(nullptr) , m_ImpactAssociate(nullptr) {} template void -ImpactImageToImageMetricv4GetValueAndDerivativeThreader::BeforeThreadedExecution() +ImpactImageToImageMetricv4GetValueAndDerivativeThreader:: + BeforeThreadedExecution() { Superclass::BeforeThreadedExecution(); @@ -93,9 +90,8 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader void -ImpactImageToImageMetricv4GetValueAndDerivativeThreader::AfterThreadedExecution() +ImpactImageToImageMetricv4GetValueAndDerivativeThreader:: + AfterThreadedExecution() { // Let the base reduce the valid-point count (and per-point derivative contributions, // if any). It also sets m_Value from the per-point measures, which we leave at zero @@ -133,12 +129,8 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader std::vector -ImpactImageToImageMetricv4GetValueAndDerivativeThreader< - TDomainPartitioner, - TImageToImageMetric, - TImpactMetric>::GetSubsetOfFeatures(const std::vector & features_index, - std::mt19937 & randomGenerator, - int n) const +ImpactImageToImageMetricv4GetValueAndDerivativeThreader:: + GetSubsetOfFeatures(const std::vector & features_index, std::mt19937 & randomGenerator, int n) const { if (features_index.size() == static_cast(n)) return features_index; @@ -404,7 +396,8 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader(this->m_ImpactAssociate->m_MovingInterpolator->Evaluate(movingPoint)); + const double movingValue = + static_cast(this->m_ImpactAssociate->m_MovingInterpolator->Evaluate(movingPoint)); movingAccessor[row][flat] = static_cast(movingValue); if (!computeDerivative) { @@ -440,9 +433,9 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader(std::distance(first, last)); - const bool computeDerivative = this->GetComputeDerivative(); - const torch::Device device(this->m_ImpactAssociate->GetDevice()); + const auto batchSize = static_cast(std::distance(first, last)); + const bool computeDerivative = this->GetComputeDerivative(); + const torch::Device device(this->m_ImpactAssociate->GetDevice()); std::unique_lock deviceLock(Impact::DeviceForwardMutex(), std::defer_lock); // d(moving coordinate)/d(parameter) at every point of the batch, and the global parameter @@ -527,7 +520,7 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader(fixedConfig.GetNumberOfChannels()); + const auto channels = static_cast(fixedConfig.GetNumberOfChannels()); const std::vector & patchSize = fixedConfig.GetPatchSize(); const int64_t voxelCount = std::accumulate(patchSize.begin(), patchSize.end(), int64_t{ 1 }, std::multiplies<>()); @@ -675,20 +668,18 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader bool -ImpactImageToImageMetricv4GetValueAndDerivativeThreader< - TDomainPartitioner, - TImageToImageMetric, - TImpactMetric>::ProcessPoint(const VirtualIndexType & virtualIndex, - const VirtualPointType & virtualPoint, - const FixedImagePointType & mappedFixedPoint, - const FixedImagePixelType &, - const FixedImageGradientType &, - const MovingImagePointType & mappedMovingPoint, - const MovingImagePixelType &, - const MovingImageGradientType &, - MeasureType & metricValueReturn, - DerivativeType &, - const ThreadIdType threadId) const +ImpactImageToImageMetricv4GetValueAndDerivativeThreader:: + ProcessPoint(const VirtualIndexType & virtualIndex, + const VirtualPointType & virtualPoint, + const FixedImagePointType & mappedFixedPoint, + const FixedImagePixelType &, + const FixedImageGradientType &, + const MovingImagePointType & mappedMovingPoint, + const MovingImagePixelType &, + const MovingImageGradientType &, + MeasureType & metricValueReturn, + DerivativeType &, + const ThreadIdType threadId) const { // Static mode only: interpolate the (optionally subsampled) per-layer feature vector at the // mapped fixed/moving point and accumulate the per-layer loss. When a derivative is @@ -699,8 +690,8 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader< const bool computeDerivative = this->GetComputeDerivative(); - torch::Tensor transformJacobian; // [1, MovingDim, P] - torch::Tensor nonZeroJacobianIndices; // [1, P] + torch::Tensor transformJacobian; // [1, MovingDim, P] + torch::Tensor nonZeroJacobianIndices; // [1, P] if (computeDerivative) { using JacobianType = typename TImageToImageMetric::JacobianType; @@ -712,11 +703,10 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader< this->m_ImpactAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries( virtualPoint, jacobian, jacobianPositional); - constexpr unsigned int movingDimension = ImageToImageMetricv4Type::MovingImageDimension; + constexpr unsigned int movingDimension = ImageToImageMetricv4Type::MovingImageDimension; const NumberOfParametersType numberOfLocalParameters = this->GetCachedNumberOfLocalParameters(); - transformJacobian = - torch::empty({ 1, static_cast(movingDimension), static_cast(numberOfLocalParameters) }, - torch::kFloat32); + transformJacobian = torch::empty( + { 1, static_cast(movingDimension), static_cast(numberOfLocalParameters) }, torch::kFloat32); auto jacobianAccessor = transformJacobian.accessor(); for (unsigned int dimension = 0; dimension < movingDimension; ++dimension) { @@ -760,19 +750,19 @@ ImpactImageToImageMetricv4GetValueAndDerivativeThreader< loss.m_randomGenerator, this->m_ImpactAssociate->GetSubsetFeatures()[i]); - torch::Tensor fixedFeatures = - this->m_ImpactAssociate->m_Internals->m_FixedFeaturesMaps[i].m_FeaturesMapInterpolator->Evaluate(mappedFixedPoint, subsetOfFeatures) - .unsqueeze(0); - torch::Tensor movingFeatures = - this->m_ImpactAssociate->m_Internals->m_MovingFeaturesMaps[i].m_FeaturesMapInterpolator->Evaluate(movingReadPoint, subsetOfFeatures) - .unsqueeze(0); + torch::Tensor fixedFeatures = this->m_ImpactAssociate->m_Internals->m_FixedFeaturesMaps[i] + .m_FeaturesMapInterpolator->Evaluate(mappedFixedPoint, subsetOfFeatures) + .unsqueeze(0); + torch::Tensor movingFeatures = this->m_ImpactAssociate->m_Internals->m_MovingFeaturesMaps[i] + .m_FeaturesMapInterpolator->Evaluate(movingReadPoint, subsetOfFeatures) + .unsqueeze(0); if (computeDerivative) { // d(feature)/d(moving coordinate) : [1, C, MovingDim] torch::Tensor movingFeatureDerivative = - this->m_ImpactAssociate->m_Internals->m_MovingFeaturesMaps[i].m_FeaturesMapInterpolator - ->EvaluateDerivative(movingReadPoint, subsetOfFeatures) + this->m_ImpactAssociate->m_Internals->m_MovingFeaturesMaps[i] + .m_FeaturesMapInterpolator->EvaluateDerivative(movingReadPoint, subsetOfFeatures) .unsqueeze(0); // chain with d(moving coordinate)/d(parameter) : [1, C, P] torch::Tensor featureParameterJacobian = torch::bmm(movingFeatureDerivative, transformJacobian); diff --git a/include/itkImpactModelConfiguration.h b/include/itkImpactModelConfiguration.h index 6ba3486..0ec45f0 100644 --- a/include/itkImpactModelConfiguration.h +++ b/include/itkImpactModelConfiguration.h @@ -89,13 +89,13 @@ struct Impact_EXPORT ImpactModelConfiguration /** Load a TorchScript model and store its configuration. The model is loaded on the * CPU; move it with ModelTo() (itkImpactModelConfigurationDetail.h). */ ImpactModelConfiguration(std::string modelPath, - unsigned int dimension, - unsigned int numberOfChannels, - std::vector patchSize, - std::vector voxelSize, - std::vector overlap, - std::vector layersMask, - bool useMixedPrecision); + unsigned int dimension, + unsigned int numberOfChannels, + std::vector patchSize, + std::vector voxelSize, + std::vector overlap, + std::vector layersMask, + bool useMixedPrecision); /** An empty configuration with no associated model. */ ImpactModelConfiguration() = default; @@ -258,7 +258,7 @@ struct Impact_EXPORT ImpactModelConfiguration std::vector m_overlaps; std::string m_patchCombine{ "cosinus" }; std::string m_featureNormalization{ "none" }; - std::vector m_layersMask; + std::vector m_layersMask; // Precomputed physical patch offsets (torch-free); see GetPatchIndex(). std::vector> m_patchIndex; // shared_ptr keeps the special members usable with an incomplete impl type, so this diff --git a/include/itkImpactModelConfigurationDetail.h b/include/itkImpactModelConfigurationDetail.h index ff33418..d9eda42 100644 --- a/include/itkImpactModelConfigurationDetail.h +++ b/include/itkImpactModelConfigurationDetail.h @@ -97,7 +97,7 @@ inline std::vector Forward(const ImpactModelConfiguration & configuration, torch::Tensor inputPatch, const ImpactImageMetadata & metadata) { detail::ImpactModelConfigurationImpl * impl = configuration.GetImpl(); - std::vector args; + std::vector args; args.reserve(impl->nArgs); args.emplace_back(inputPatch); if (impl->nArgs >= 2) diff --git a/include/itkImpactOnlineInference.h b/include/itkImpactOnlineInference.h index 4825b14..b8918c5 100644 --- a/include/itkImpactOnlineInference.h +++ b/include/itkImpactOnlineInference.h @@ -294,30 +294,31 @@ GenerateOutputs(const std::vector & device, nbSample, [&](int64_t begin, int64_t end) { - torch::Tensor input = - patchValueTensor.narrow(0, begin, end - begin).to(device).repeat({ torch::IntArrayRef(resizeVector) }).clone(); - std::vector outputsList = itk::Forward(config, input); - std::vector layers; - layers.reserve(kept); - for (size_t it = 0; it < outputsList.size(); ++it) - { - if (mask[it]) + torch::Tensor input = patchValueTensor.narrow(0, begin, end - begin) + .to(device) + .repeat({ torch::IntArrayRef(resizeVector) }) + .clone(); + std::vector outputsList = itk::Forward(config, input); + std::vector layers; + layers.reserve(kept); + for (size_t it = 0; it < outputsList.size(); ++it) { - const size_t k = layers.size(); - // The point's feature vector normalized over all its channels, before the subset keeps a few. - layers.push_back(Impact::NormalizeFeatureChannels(outputsList[it] - .toTensor() - .index(itk::GetCentersIndexLayers(config)[a + k]) - .to(torch::kFloat32), - config.GetFeatureNormalization(), - 1) - .index_select(1, subsetsOfFeatures[a + k])); + if (mask[it]) + { + const size_t k = layers.size(); + // The point's feature vector normalized over all its channels, before the subset keeps a few. + layers.push_back( + Impact::NormalizeFeatureChannels( + outputsList[it].toTensor().index(itk::GetCentersIndexLayers(config)[a + k]).to(torch::kFloat32), + config.GetFeatureNormalization(), + 1) + .index_select(1, subsetsOfFeatures[a + k])); + } + } + for (size_t k = 0; k < layers.size(); ++k) + { + batches[k].push_back(layers[k]); } - } - for (size_t k = 0; k < layers.size(); ++k) - { - batches[k].push_back(layers[k]); - } }, [] {}); for (size_t k = 0; k < kept; ++k) @@ -366,9 +367,9 @@ GenerateOutputsAndJacobian( for (int64_t s = 0; s < nbSample; ++s) { - patchValueTensor[s] = imagesPatchValuesAndJacobiansEvaluator( - fixedPoints[s], imagesPatchesJacobians, patchIndex[i][s], tensorShape, s) - .to(itk::GetModelDtype(config)); + patchValueTensor[s] = + imagesPatchValuesAndJacobiansEvaluator(fixedPoints[s], imagesPatchesJacobians, patchIndex[i][s], tensorShape, s) + .to(itk::GetModelDtype(config)); } std::vector resizeVector(patchValueTensor.dim(), 1); @@ -413,15 +414,14 @@ GenerateOutputsAndJacobian( { continue; } - const size_t k = jacobians.size(); + const size_t k = jacobians.size(); // Normalized inside the graph, so the Jacobian below carries the normalization too. - torch::Tensor layer = Impact::NormalizeFeatureChannels(outputsList[it] - .toTensor() - .index(itk::GetCentersIndexLayers(config)[a + k]) - .to(torch::kFloat32), - config.GetFeatureNormalization(), - 1) - .index_select(1, subsetsOfFeatures[a + k]); + torch::Tensor layer = + Impact::NormalizeFeatureChannels( + outputsList[it].toTensor().index(itk::GetCentersIndexLayers(config)[a + k]).to(torch::kFloat32), + config.GetFeatureNormalization(), + 1) + .index_select(1, subsetsOfFeatures[a + k]); torch::Tensor fixedLayer = fixedOutputsTensor[a + k].narrow(0, begin, n); torch::Tensor gradientModulator = losses[a + k]->updateValueAndGetGradientModulator(fixedLayer, layer); jacobians.push_back( diff --git a/include/itkImpactPatchTiling.h b/include/itkImpactPatchTiling.h index c38377d..eb0f8a4 100644 --- a/include/itkImpactPatchTiling.h +++ b/include/itkImpactPatchTiling.h @@ -360,10 +360,9 @@ MakePatchGrid(const std::vector & shape, const int64_t stride = patchSize[d] - overlaps[d]; if (stride <= 0) { - itkGenericExceptionMacro("ImageToFeaturesMap: overlap (" << overlaps[d] - << ") must be smaller than the patch size (" - << patchSize[d] << ") along dimension " << d - << "; otherwise consecutive patches never advance."); + itkGenericExceptionMacro("ImageToFeaturesMap: overlap (" + << overlaps[d] << ") must be smaller than the patch size (" << patchSize[d] + << ") along dimension " << d << "; otherwise consecutive patches never advance."); } int64_t count = 1; if (extent > patchSize[d]) @@ -726,13 +725,14 @@ class Accumulator * the blend. A caller that needs the graph must run the model whole. */ inline std::vector -RunTiledModelOnce(const ImpactModelConfiguration & config, - const torch::Tensor & input, - const torch::Device & device, - const torch::Device & accumulateOn, - PatchCombineMode combine, - const std::function &)> & makeDestination, - const std::vector & patchOverride) // ITK order; a whole-extent axis cut down by RunTiledModel +RunTiledModelOnce( + const ImpactModelConfiguration & config, + const torch::Tensor & input, + const torch::Device & device, + const torch::Device & accumulateOn, + PatchCombineMode combine, + const std::function &)> & makeDestination, + const std::vector & patchOverride) // ITK order; a whole-extent axis cut down by RunTiledModel { torch::NoGradGuard noGrad; @@ -741,7 +741,8 @@ RunTiledModelOnce(const ImpactModelConfiguration & if (dimension == 0 || dimension > imageDimension) { itkGenericExceptionMacro("IMPACT: the model " << config.GetModelPath() << " is configured for " << dimension - << " dimension(s), which a " << imageDimension << "D image cannot be tiled for."); + << " dimension(s), which a " << imageDimension + << "D image cannot be tiled for."); } const unsigned int sweptAxes = imageDimension - dimension; @@ -826,16 +827,17 @@ RunTiledModelOnce(const ImpactModelConfiguration & } catch (const std::exception & e) { - itkGenericExceptionMacro( - "IMPACT: the model " << config.GetModelPath() - << " rejected its input. Check the number of channels (" << config.GetNumberOfChannels() - << "), the patch size " << GetStringFromVector(config.GetPatchSize()) << " and the dimension (" - << config.GetDimension() << ").\nDetails: " << e.what()); + itkGenericExceptionMacro("IMPACT: the model " + << config.GetModelPath() << " rejected its input. Check the number of channels (" + << config.GetNumberOfChannels() << "), the patch size " + << GetStringFromVector(config.GetPatchSize()) << " and the dimension (" + << config.GetDimension() << ").\nDetails: " << e.what()); } if (config.GetLayersMask().size() != outputsList.size()) { itkGenericExceptionMacro("IMPACT: " << config.GetModelPath() << " declares " << config.GetLayersMask().size() - << " layer mask entries but returned " << outputsList.size() << " layers."); + << " layer mask entries but returned " << outputsList.size() + << " layers."); } } else @@ -850,8 +852,7 @@ RunTiledModelOnce(const ImpactModelConfiguration & { continue; } - const torch::Tensor layerPatch = - outputsList[it].toTensor().squeeze(0).to(accumulateOn).to(torch::kFloat32); + const torch::Tensor layerPatch = outputsList[it].toTensor().squeeze(0).to(accumulateOn).to(torch::kFloat32); if (layerIndex == accumulators.size()) { @@ -874,11 +875,10 @@ RunTiledModelOnce(const ImpactModelConfiguration & layerOverlaps[d] = std::llround(overlaps[d] * scales[d]); if (layerOverlaps[d] >= layerPatchSize[d]) { - itkGenericExceptionMacro("IMPACT: the overlap (" << overlaps[d] - << ") scaled to this layer's resolution (" << layerOverlaps[d] - << ") is not smaller than the layer's patch extent (" - << layerPatchSize[d] << ") along dimension " << d - << ". Reduce the overlap or keep a finer layer."); + itkGenericExceptionMacro("IMPACT: the overlap (" + << overlaps[d] << ") scaled to this layer's resolution (" << layerOverlaps[d] + << ") is not smaller than the layer's patch extent (" << layerPatchSize[d] + << ") along dimension " << d << ". Reduce the overlap or keep a finer layer."); } } layerGrid = ScalePatchGrid(grid, scales); @@ -902,8 +902,7 @@ RunTiledModelOnce(const ImpactModelConfiguration & assembled.push_back(destination); } - accumulators.emplace_back( - layerShape, layerPatchSize, layerOverlaps, layerGrid, MakePatchCombine(combine)); + accumulators.emplace_back(layerShape, layerPatchSize, layerOverlaps, layerGrid, MakePatchCombine(combine)); // This slice's own sub-view of the layer's buffer. torch::Tensor destination = assembled[layerIndex]; @@ -924,7 +923,8 @@ RunTiledModelOnce(const ImpactModelConfiguration & // A layer that did not see every patch would be handed back as a valid map of size zero. if (!accumulators[i].IsFull()) { - itkGenericExceptionMacro("IMPACT: layer " << i << " of " << config.GetModelPath() + itkGenericExceptionMacro("IMPACT: layer " + << i << " of " << config.GetModelPath() << " was never assembled: the accumulator expected more patches than the tiling " "produced. This is an internal inconsistency, not a configuration error."); } diff --git a/include/itkImpactTorchRegistrationHelpers.h b/include/itkImpactTorchRegistrationHelpers.h index 961104b..5efa9f5 100644 --- a/include/itkImpactTorchRegistrationHelpers.h +++ b/include/itkImpactTorchRegistrationHelpers.h @@ -137,7 +137,8 @@ class PythonGilReleaseGuard } } PythonGilReleaseGuard(const PythonGilReleaseGuard &) = delete; - PythonGilReleaseGuard & operator=(const PythonGilReleaseGuard &) = delete; + PythonGilReleaseGuard & + operator=(const PythonGilReleaseGuard &) = delete; private: void * m_State{ nullptr }; @@ -181,7 +182,9 @@ ImageToBatchTensor(const TImage * image) * stay aligned, which is how the registration filters place a feature layer of any size over the image. */ template torch::Tensor -ResampleToVoxelSize(const torch::Tensor & image, const std::vector & spacing, const std::vector & voxelSize) +ResampleToVoxelSize(const torch::Tensor & image, + const std::vector & spacing, + const std::vector & voxelSize) { namespace F = torch::nn::functional; if (spacing.size() != Dim || voxelSize.size() != Dim) @@ -242,8 +245,7 @@ ExtractFeatureLayers(const std::vector & configs, if (!subset.empty()) { std::vector idx(subset.begin(), subset.end()); - subsetIndex = - torch::from_blob(idx.data(), { static_cast(idx.size()) }, torch::kLong).clone().to(device); + subsetIndex = torch::from_blob(idx.data(), { static_cast(idx.size()) }, torch::kLong).clone().to(device); } // Keep only the requested channels of a layer, validating the indices first. @@ -292,7 +294,8 @@ ExtractFeatureLayers(const std::vector & configs, sweepAxis)) { layers.push_back( - keepSubset(Impact::NormalizeFeatureChannels(map.unsqueeze(0), config.GetFeatureNormalization(), 1)).contiguous()); + keepSubset(Impact::NormalizeFeatureChannels(map.unsqueeze(0), config.GetFeatureNormalization(), 1)) + .contiguous()); } continue; } @@ -347,15 +350,15 @@ ExtractFeatureLayers(const std::vector & configs, { outputs = Forward(config, input); } - const std::vector & mask = config.GetLayersMask(); + const std::vector & mask = config.GetLayersMask(); for (size_t i = 0; i < outputs.size(); ++i) { if (i < mask.size() && !mask[i]) { continue; } - torch::Tensor layer = keepSubset( - Impact::NormalizeFeatureChannels(outputs[i].toTensor().to(torch::kFloat32), config.GetFeatureNormalization(), 1)); + torch::Tensor layer = keepSubset(Impact::NormalizeFeatureChannels( + outputs[i].toTensor().to(torch::kFloat32), config.GetFeatureNormalization(), 1)); // Keep the layer at its NATIVE resolution (see the function doc); the consumer resamples it. // Detach only in the no-grad path -- withGrad must preserve the graph back to `imageTensor`. layers.push_back(withGrad ? layer.contiguous() : layer.detach().contiguous()); @@ -595,10 +598,10 @@ SqueezedChannelLastPermutation() * where voxelOffset_xyz reverses the z,y,x components. This is the fixed->moving convention. */ template void -WriteVoxelFieldToDisplacement(const torch::Tensor & voxelField, - const typename Image, Dim>::SpacingType & spacing, - const typename Image, Dim>::DirectionType & direction, - Image, Dim> * output) +WriteVoxelFieldToDisplacement(const torch::Tensor & voxelField, + const typename Image, Dim>::SpacingType & spacing, + const typename Image, Dim>::DirectionType & direction, + Image, Dim> * output) { using FieldType = Image, Dim>; using VectorType = Vector; @@ -639,10 +642,10 @@ WriteVoxelFieldToDisplacement(const torch::Tensor & * orthonormal). Used to ingest an initial/warm-start field. */ template torch::Tensor -DisplacementToVoxelField(const Image, Dim> * field, - const typename Image, Dim>::SpacingType & spacing, - const typename Image, Dim>::DirectionType & direction, - const torch::Device & device) +DisplacementToVoxelField(const Image, Dim> * field, + const typename Image, Dim>::SpacingType & spacing, + const typename Image, Dim>::DirectionType & direction, + const torch::Device & device) { using FieldType = Image, Dim>; const auto size = field->GetLargestPossibleRegion().GetSize(); diff --git a/include/itkInterpolateVectorImageFunction.h b/include/itkInterpolateVectorImageFunction.h index 0672766..18ff2ed 100644 --- a/include/itkInterpolateVectorImageFunction.h +++ b/include/itkInterpolateVectorImageFunction.h @@ -38,7 +38,7 @@ namespace itk template -class ITK_TEMPLATE_EXPORT InterpolateVectorImageFunction +class ITK_TEMPLATE_EXPORT InterpolateVectorImageFunction : public InterpolateImageFunction { public: @@ -46,7 +46,8 @@ class ITK_TEMPLATE_EXPORT InterpolateVectorImageFunction /** Standard class type aliases. */ using Self = InterpolateVectorImageFunction; - using Superclass = InterpolateImageFunction; + using Superclass = + InterpolateImageFunction; using Pointer = SmartPointer; using ConstPointer = SmartPointer; @@ -65,9 +66,9 @@ class ITK_TEMPLATE_EXPORT InterpolateVectorImageFunction /** Dimension underlying input image. */ static constexpr unsigned int ImageDimension = Superclass::ImageDimension; - + using CovariantVectorType = itk::CovariantVector; - + InterpolateVectorImageFunction() = default; @@ -108,15 +109,17 @@ class ITK_TEMPLATE_EXPORT InterpolateVectorImageFunction EvaluateDerivative(ImagePointType point, std::vector subsetOfFeatures) const; using OutputType = typename Superclass::OutputType; -using ContinuousIndexType = typename Superclass::ContinuousIndexType; -using SizeType = typename Superclass::SizeType; + using ContinuousIndexType = typename Superclass::ContinuousIndexType; + using SizeType = typename Superclass::SizeType; - OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &) const override + OutputType + EvaluateAtContinuousIndex(const ContinuousIndexType &) const override { itkExceptionMacro("EvaluateAtContinuousIndex is not implemented."); } - SizeType GetRadius() const override + SizeType + GetRadius() const override { return SizeType{}; } diff --git a/include/itkInterpolateVectorImageFunction.hxx b/include/itkInterpolateVectorImageFunction.hxx index 94d036b..8dba9df 100644 --- a/include/itkInterpolateVectorImageFunction.hxx +++ b/include/itkInterpolateVectorImageFunction.hxx @@ -48,8 +48,8 @@ InterpolateVectorImageFunction::SetInputImage(ImagePointe */ template torch::Tensor -InterpolateVectorImageFunction::Evaluate(ImagePointType point, - std::vector subsetOfFeatures) const +InterpolateVectorImageFunction::Evaluate(ImagePointType point, + std::vector subsetOfFeatures) const { std::vector result; for (const unsigned int feature : subsetOfFeatures) @@ -65,8 +65,8 @@ InterpolateVectorImageFunction::Evaluate(ImagePointType p template torch::Tensor InterpolateVectorImageFunction::EvaluateDerivative( - ImagePointType point, - std::vector subsetOfFeatures) const + ImagePointType point, + std::vector subsetOfFeatures) const { diff --git a/include/itkTensorToImageFilter.h b/include/itkTensorToImageFilter.h index d294818..a05b805 100644 --- a/include/itkTensorToImageFilter.h +++ b/include/itkTensorToImageFilter.h @@ -74,7 +74,8 @@ class ITK_TEMPLATE_EXPORT TensorToImageFilter : public ImageSourceGetOrigin(); m_Spacing = referenceImage->GetSpacing(); @@ -132,7 +134,7 @@ class ITK_TEMPLATE_EXPORT TensorToImageFilter : public ImageSource; - + /** Set the direction of the image * \sa GetDirection() */ virtual void @@ -162,11 +164,10 @@ class ITK_TEMPLATE_EXPORT TensorToImageFilter : public ImageSource(1.0) }; - OriginType m_Origin{MakeFilled(0)}; + OriginType m_Origin{ MakeFilled(0) }; DirectionType m_Direction{ DirectionType::GetIdentity() }; SizeType m_Size{}; torch::Tensor m_tensor; - }; } // end namespace itk diff --git a/include/itkTensorToImageFilter.hxx b/include/itkTensorToImageFilter.hxx index 2e2788a..3307dae 100644 --- a/include/itkTensorToImageFilter.hxx +++ b/include/itkTensorToImageFilter.hxx @@ -77,9 +77,9 @@ TensorToImageFilter::GenerateOutputInformation() OutputImagePointer outputPtr = this->GetOutput(); SpacingType spacing; - RegionType region; - SizeType size; - + RegionType region; + SizeType size; + for (int s = 0; s < VImageDimension; ++s) { size[s] = m_tensor.size(VImageDimension - s); @@ -87,16 +87,19 @@ TensorToImageFilter::GenerateOutputInformation() region.SetSize(size); outputPtr->SetRegions(region); outputPtr->SetVectorLength(m_tensor.size(0)); - if (m_Size[0] > 0){ + if (m_Size[0] > 0) + { for (int i = 0; i < VImageDimension; ++i) { spacing[i] = m_Size[i] * m_Spacing[i] / size[i]; } outputPtr->SetSpacing(spacing); - } else { + } + else + { outputPtr->SetSpacing(m_Spacing); } - + outputPtr->SetOrigin(m_Origin); outputPtr->SetDirection(m_Direction); outputPtr->Allocate(); @@ -107,7 +110,7 @@ void TensorToImageFilter::GenerateData() { OutputImagePointer outputPtr = this->GetOutput(); - SizeType size = outputPtr->GetLargestPossibleRegion().GetSize(); + SizeType size = outputPtr->GetLargestPossibleRegion().GetSize(); std::vector dims; for (int64_t i = 1; i < m_tensor.dim(); ++i) @@ -118,8 +121,8 @@ TensorToImageFilter::GenerateData() torch::Tensor layers = m_tensor.permute(dims).contiguous().to(torch::kFloat32); - - unsigned int numberOfChannels = outputPtr->GetVectorLength(); + + unsigned int numberOfChannels = outputPtr->GetVectorLength(); const float * layersData = layers.data_ptr(); const unsigned int rowStride = size[0] * numberOfChannels; const unsigned int sliceStride = size[0] * size[1] * numberOfChannels; @@ -154,7 +157,7 @@ TensorToImageFilter::GenerateData() { for (int z = 0; z < size[0]; ++z) { - const float * pixelPtr = layersData + x * sliceStride + y * rowStride + z * numberOfChannels; + const float * pixelPtr = layersData + x * sliceStride + y * rowStride + z * numberOfChannels; VariableLengthVector variableLengthVector(numberOfChannels); for (unsigned int i = 0; i < numberOfChannels; ++i) { diff --git a/src/itkImpactModelConfiguration.cxx b/src/itkImpactModelConfiguration.cxx index 1572a00..fb9e9b2 100644 --- a/src/itkImpactModelConfiguration.cxx +++ b/src/itkImpactModelConfiguration.cxx @@ -26,13 +26,13 @@ namespace itk { ImpactModelConfiguration::ImpactModelConfiguration(std::string modelPath, - unsigned int dimension, - unsigned int numberOfChannels, - std::vector patchSize, - std::vector voxelSize, - std::vector overlap, - std::vector layersMask, - bool useMixedPrecision) + unsigned int dimension, + unsigned int numberOfChannels, + std::vector patchSize, + std::vector voxelSize, + std::vector overlap, + std::vector layersMask, + bool useMixedPrecision) : m_modelPath(modelPath) , m_dimension(dimension) , m_numberOfChannels(numberOfChannels) @@ -52,8 +52,8 @@ ImpactModelConfiguration::ImpactModelConfiguration(std::string mod // `prim::If`): such a graph can abort the profiling executor's TensorExpr fuser on some // LibTorch builds. Disabling the optimization globally to tolerate that slowed all // inference ~2x, so it is left on. - m_impl->model = std::make_shared( - torch::jit::load(m_modelPath, torch::Device(torch::kCPU))); + m_impl->model = + std::make_shared(torch::jit::load(m_modelPath, torch::Device(torch::kCPU))); m_impl->model->eval(); m_impl->model->to(m_impl->dtype); diff --git a/test/itkImageToFeaturesMapTest.cxx b/test/itkImageToFeaturesMapTest.cxx index 6fe4631..4c50b90 100644 --- a/test/itkImageToFeaturesMapTest.cxx +++ b/test/itkImageToFeaturesMapTest.cxx @@ -43,8 +43,8 @@ itkImageToFeaturesMapTest(int argc, char * argv[]) using InterpolatorType = itk::BSplineInterpolateImageFunction; using ImageToFeaturesMapType = itk::ImageToFeaturesMap; - auto image = ImageType::New(); - ImageType::SizeType size; + auto image = ImageType::New(); + ImageType::SizeType size; size.Fill(8); ImageType::RegionType region; region.SetSize(size); diff --git a/test/itkImageToTensorFilterTest.cxx b/test/itkImageToTensorFilterTest.cxx index fa8a32d..78713f8 100644 --- a/test/itkImageToTensorFilterTest.cxx +++ b/test/itkImageToTensorFilterTest.cxx @@ -44,8 +44,8 @@ itkImageToTensorFilterTest(int, char *[]) using InterpolatorType = itk::BSplineInterpolateImageFunction; using ImageToTensorFilterType = itk::ImageToTensorFilter; - auto image = ImageType::New(); - ImageType::SizeType size; + auto image = ImageType::New(); + ImageType::SizeType size; size.Fill(8); ImageType::RegionType region; region.SetSize(size); @@ -94,13 +94,12 @@ itkImageToTensorFilterTest(int, char *[]) } } } - const torch::Tensor expected = - torch::from_blob(expectedBuffer.data(), { 8, 8, 8 }, torch::kFloat32).clone(); + const torch::Tensor expected = torch::from_blob(expectedBuffer.data(), { 8, 8, 8 }, torch::kFloat32).clone(); if (!torch::allclose(t, expected, 1e-4, 1e-4)) { std::cerr << "The tensor does not reproduce the input ramp: max |difference| " - << (t - expected).abs().max().item() << ", non-zero voxels " - << t.count_nonzero().item() << " of " << t.numel() << std::endl; + << (t - expected).abs().max().item() << ", non-zero voxels " << t.count_nonzero().item() + << " of " << t.numel() << std::endl; return EXIT_FAILURE; } @@ -138,8 +137,8 @@ itkImageToTensorFilterTest(int, char *[]) const torch::Tensor rotatedTensor = rotatedFilter->GetTensor(); if (rotatedTensor.sizes() != t.sizes()) { - std::cerr << "Rotating the image changed the tensor shape: " << t.sizes() << " vs " - << rotatedTensor.sizes() << std::endl; + std::cerr << "Rotating the image changed the tensor shape: " << t.sizes() << " vs " << rotatedTensor.sizes() + << std::endl; return EXIT_FAILURE; } if (!torch::allclose(rotatedTensor, t, 1e-4, 1e-4)) @@ -148,8 +147,7 @@ itkImageToTensorFilterTest(int, char *[]) "rotated image gave a different tensor.\n" << " max |difference| : " << (rotatedTensor - t).abs().max().item() << '\n' << " non-zero voxels : " << rotatedTensor.count_nonzero().item() << " of " - << rotatedTensor.numel() << " (identity image: " << t.count_nonzero().item() << ')' - << std::endl; + << rotatedTensor.numel() << " (identity image: " << t.count_nonzero().item() << ')' << std::endl; return EXIT_FAILURE; } std::cout << "direction is honoured: the rotated image gives the same tensor" << std::endl; diff --git a/test/itkImpactBackendGTest.cxx b/test/itkImpactBackendGTest.cxx index 6ac1bda..1dbf42e 100644 --- a/test/itkImpactBackendGTest.cxx +++ b/test/itkImpactBackendGTest.cxx @@ -95,8 +95,8 @@ MakeRampImage(unsigned int n, int mode = 0) ImageType::Pointer MakeBlobImage(unsigned int n, double cx, double cy, double cz, double sigma) { - auto image = ImageType::New(); - ImageType::SizeType size; + auto image = ImageType::New(); + ImageType::SizeType size; size.Fill(n); ImageType::RegionType region; region.SetSize(size); @@ -109,8 +109,8 @@ MakeBlobImage(unsigned int n, double cx, double cy, double cz, double sigma) for (it.GoToBegin(); !it.IsAtEnd(); ++it) { const auto i = it.GetIndex(); - const double d2 = (i[0] - centerX) * (i[0] - centerX) + (i[1] - centerY) * (i[1] - centerY) + - (i[2] - centerZ) * (i[2] - centerZ); + const double d2 = + (i[0] - centerX) * (i[0] - centerX) + (i[1] - centerY) * (i[1] - centerY) + (i[2] - centerZ) * (i[2] - centerZ); image->SetPixel(i, static_cast(std::exp(-d2 / (2.0 * sigma * sigma)))); } return image; @@ -128,8 +128,8 @@ TEST(ImpactBackend, ImageToTensorFilterResampling) { using FilterType = itk::ImageToTensorFilter; const ImageType::SizeType size = { { 8, 6, 4 } }; // distinct dims to verify axis order - auto image = MakeRampImage(size); - auto interpolator = InterpolatorType::New(); + auto image = MakeRampImage(size); + auto interpolator = InterpolatorType::New(); interpolator->SetSplineOrder(3); auto filter = FilterType::New(); @@ -168,7 +168,7 @@ TEST(ImpactBackend, ImpactLossDiceNoMutationAndNaNSafe) { // Does not mutate its inputs. { - auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); + auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); dice->SetNumberOfParameters(1); torch::Tensor f = torch::rand({ 4, 3 }); torch::Tensor m = torch::rand({ 4, 3 }); @@ -180,7 +180,7 @@ TEST(ImpactBackend, ImpactLossDiceNoMutationAndNaNSafe) } // Identical inputs => perfect overlap (Dice == 1) => value == -1. { - auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); + auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); dice->SetNumberOfParameters(1); torch::Tensor f = torch::ones({ 4, 3 }); torch::Tensor m = f.clone(); @@ -189,7 +189,7 @@ TEST(ImpactBackend, ImpactLossDiceNoMutationAndNaNSafe) } // Empty/empty (all zeros) => finite value and gradient (no 0/0 NaN). { - auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); + auto dice = itk::Impact::LossFactory::Instance().Create("Dice"); dice->SetNumberOfParameters(1); torch::Tensor f = torch::zeros({ 4, 3 }); torch::Tensor m = torch::zeros({ 4, 3 }); @@ -238,7 +238,7 @@ TEST(ImpactBackend, ImageToFeaturesMapOverlapZeroKeepsSize) // overlap = 0 (6th ctor arg), patchSize {0,0,0} = full image -> single patch. itk::ImpactModelConfiguration config( ToyModelPath(), 3, 1, { 0, 0, 0 }, { 1.f, 1.f, 1.f }, { 0, 0, 0 }, { true, true }, false); - auto filter = FeatMapType::New(); + auto filter = FeatMapType::New(); filter->SetModelConfiguration(config); filter->SetInterpolator(interpolator); filter->AddInput(image); @@ -332,7 +332,7 @@ TEST(ImpactMetric, StaticFeaturesMapRefreshFollowsTheTransform) // The passthrough layer is the one to watch: its "features" ARE the intensities, so a map // that followed the transform is directly readable as the resampled image. auto build = [&](int interval, TransformType * transform) { - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -369,7 +369,7 @@ TEST(ImpactMetric, StaticFeaturesMapRefreshFollowsTheTransform) transformOff->SetIdentity(); auto off = build(0, transformOff); off->GetValueAndDerivative(value, derivative); - const double offBefore = static_cast(off->GetValue()); + const double offBefore = static_cast(off->GetValue()); MetricType::ParametersType shifted = transformOff->GetParameters(); shifted[0] += 3.0; transformOff->SetParameters(shifted); @@ -386,7 +386,7 @@ TEST(ImpactMetric, StaticFeaturesMapRefreshFollowsTheTransform) transformOn->SetIdentity(); auto on = build(1, transformOn); on->GetValueAndDerivative(value, derivative); - const double onBefore = static_cast(on->GetValue()); + const double onBefore = static_cast(on->GetValue()); MetricType::ParametersType shiftedOn = transformOn->GetParameters(); shiftedOn[0] += 3.0; transformOn->SetParameters(shiftedOn); @@ -486,7 +486,7 @@ TEST(ImpactMetric, StaticValueZeroForIdenticalPositiveForDifferent) auto fixed = MakeRampImage(8); auto buildAndEvaluate = [&](ImageType::Pointer moving) -> double { - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -497,9 +497,9 @@ TEST(ImpactMetric, StaticValueZeroForIdenticalPositiveForDifferent) std::vector{ true, false }, false); metric->SetModelsConfiguration(configs); - metric->SetDistance({ "Cosine" }); // non-normalized: value is the actual similarity + metric->SetDistance({ "Cosine" }); // non-normalized: value is the actual similarity metric->SetLayersWeight({ 1.f }); - metric->SetSubsetFeatures({ 4 }); // use all 4 conv feature channels deterministically + metric->SetSubsetFeatures({ 4 }); // use all 4 conv feature channels deterministically metric->SetPCA({ 0 }); metric->SetMode("Static"); metric->SetSeed(1); @@ -541,7 +541,7 @@ TEST(ImpactMetric, StaticDerivativeMatchesFiniteDifferences) auto checkLoss = [&](const char * lossName) { SCOPED_TRACE(std::string("loss ") + lossName); SCOPED_TRACE(lossName); - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -608,7 +608,7 @@ TEST(ImpactMetric, StaticDerivativeMatchesFiniteDifferences) }; checkLoss("Cosine"); // d(1 - cosine)/dm with the full dot-product cross term - checkLoss("NCC"); // cross-point statistics path (closed-form derivative) + checkLoss("NCC"); // cross-point statistics path (closed-form derivative) checkLoss("L1"); checkLoss("L2"); checkLoss("L1Cosine"); @@ -759,7 +759,7 @@ TEST(ImpactMetric, StaticDerivativeMatchesFiniteDifferencesDisplacementField) auto fixed = MakeRampImage(8, 0); auto moving = MakeRampImage(8, 1); - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -810,10 +810,10 @@ TEST(ImpactMetric, StaticDerivativeMatchesFiniteDifferencesDisplacementField) // The field has thousands of parameters; spot-check a few interior voxels. // Parameters are laid out [vox0_x, vox0_y, vox0_z, vox1_x, ...]. - const double h = 1e-3; - const MetricType::ParametersType params0 = displacement->GetParameters(); - const std::vector indices = { 3 * 100, 3 * 100 + 1, 3 * 100 + 2, 3 * 200 + 1, 3 * 250 + 2 }; - unsigned int nonZeroChecks = 0; + const double h = 1e-3; + const MetricType::ParametersType params0 = displacement->GetParameters(); + const std::vector indices = { 3 * 100, 3 * 100 + 1, 3 * 100 + 2, 3 * 200 + 1, 3 * 250 + 2 }; + unsigned int nonZeroChecks = 0; for (unsigned int j : indices) { if (j >= displacement->GetNumberOfParameters()) @@ -858,7 +858,7 @@ TEST(ImpactMetric, JacobianDerivativeMatchesFiniteDifferences) auto checkLoss = [&](const char * lossName) { SCOPED_TRACE(std::string("loss ") + lossName); SCOPED_TRACE(lossName); - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; // A strictly positive 3^3 patch at the image spacing: the toy conv (kernel 3, // padding 1) computes the exact neighborhood feature at the patch center. @@ -931,7 +931,7 @@ TEST(ImpactMetric, JacobianDerivativeMatchesFiniteDifferences) }; checkLoss("Cosine"); // normalized cross term - checkLoss("NCC"); // cross-point statistics path + checkLoss("NCC"); // cross-point statistics path checkLoss("L1"); checkLoss("L2"); checkLoss("L1Cosine"); @@ -956,7 +956,7 @@ TEST(ImpactMetric, JacobianResultIsIndependentOfBatchSize) auto identityFixed = itk::IdentityTransform::New(); auto evaluate = [&](const char * lossName, unsigned int batchSize, MetricType::DerivativeType & derivative) { - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -1133,7 +1133,7 @@ TEST(ImpactMetric, PerLayerVectorsMustMatchTheNumberOfComparedLayers) // One model keeping both layers, so two layers are compared: two distances, two weights, two // subsets -- and one PCA, which belongs to the model. auto build = [&](size_t distances, size_t weights, size_t subsets, size_t pca) { - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; configs.emplace_back(ToyModelPath(), 3, @@ -1192,7 +1192,7 @@ TEST(ImpactMetric, TwoDimensionalModelJacobianDerivativeMatchesFiniteDifferences transform->SetIdentity(); auto makeMetric = [&](unsigned int workUnits, unsigned int seed) { - auto metric = MetricType::New(); + auto metric = MetricType::New(); std::vector configs; // A 2D patch on a 3D image. The voxel size carries an entry per image axis because // ImpactModelConfiguration pairs the two when it precomputes its patch offsets; Jacobian mode @@ -1424,16 +1424,16 @@ TEST(ImpactMetric, TranslationRegistrationConvergesForEveryDistance) { for (const std::string distance : { "NCC", "L1", "L2", "Cosine", "L1Cosine" }) { - using VirtualImageType = itk::Image; - using MetricType = itk::ImpactImageToImageMetricv4; - using TransformType = itk::TranslationTransform; + using VirtualImageType = itk::Image; + using MetricType = itk::ImpactImageToImageMetricv4; + using TransformType = itk::TranslationTransform; - // Two sharp blobs straddling the image center, offset by +4 voxels in x (fixed at - // 6, moving at 10 in a size-16 volume). Both blobs sit well inside the image so the - // feature cross-correlation is not biased by boundary truncation. The optimal moving - // translation that brings the moving blob onto the fixed one is exactly (+4, 0, 0). - auto fixed = MakeBlobImage(16, -2.0, 0.0, 0.0, 2.0); - auto moving = MakeBlobImage(16, 2.0, 0.0, 0.0, 2.0); + // Two sharp blobs straddling the image center, offset by +4 voxels in x (fixed at + // 6, moving at 10 in a size-16 volume). Both blobs sit well inside the image so the + // feature cross-correlation is not biased by boundary truncation. The optimal moving + // translation that brings the moving blob onto the fixed one is exactly (+4, 0, 0). + auto fixed = MakeBlobImage(16, -2.0, 0.0, 0.0, 2.0); + auto moving = MakeBlobImage(16, 2.0, 0.0, 0.0, 2.0); auto metric = MetricType::New(); std::vector configs; @@ -1455,45 +1455,45 @@ TEST(ImpactMetric, TranslationRegistrationConvergesForEveryDistance) metric->SetFeaturesMapUpdateInterval(-1); metric->SetDevice("cpu"); - auto identityFixed = itk::IdentityTransform::New(); - auto transform = TransformType::New(); - transform->SetIdentity(); + auto identityFixed = itk::IdentityTransform::New(); + auto transform = TransformType::New(); + transform->SetIdentity(); - metric->SetFixedImage(fixed); - metric->SetMovingImage(moving); - metric->SetFixedTransform(identityFixed); - metric->SetMovingTransform(transform); - metric->SetUseFixedImageGradientFilter(false); - metric->SetUseMovingImageGradientFilter(false); - metric->SetMaximumNumberOfWorkUnits(1); - metric->Initialize(); + metric->SetFixedImage(fixed); + metric->SetMovingImage(moving); + metric->SetFixedTransform(identityFixed); + metric->SetMovingTransform(transform); + metric->SetUseFixedImageGradientFilter(false); + metric->SetUseMovingImageGradientFilter(false); + metric->SetMaximumNumberOfWorkUnits(1); + metric->Initialize(); - const double initialValue = static_cast(metric->GetValue()); - - // Drive a real ITK optimizer with the metric. - using OptimizerType = itk::RegularStepGradientDescentOptimizerv4; - using ScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift; - auto scales = ScalesEstimatorType::New(); - scales->SetMetric(metric); - auto optimizer = OptimizerType::New(); - optimizer->SetMetric(metric); - optimizer->SetScalesEstimator(scales); - optimizer->SetNumberOfIterations(300); - optimizer->SetLearningRate(2.0); - optimizer->SetMinimumStepLength(1e-4); - optimizer->SetRelaxationFactor(0.8); - optimizer->StartOptimization(); - - // The metric, through its analytic derivative, must drive the optimizer downhill. - const double finalValue = static_cast(metric->GetValue()); - EXPECT_LT(finalValue, initialValue) << distance << ": optimization should reduce the metric"; - - // The optimizer recovers the true geometric alignment: bringing the moving blob - // (x=10) onto the fixed one (x=6) is a +4 voxel translation in x, with no y/z motion. - const TransformType::ParametersType recovered = transform->GetParameters(); - EXPECT_NEAR(recovered[0], 4.0, 0.5) << distance << ": recovered x-translation should reach +4"; - EXPECT_LT(std::abs(recovered[1]), 0.2) << distance << ": motion should stay on the x axis"; - EXPECT_LT(std::abs(recovered[2]), 0.2) << distance << ": motion should stay on the x axis"; + const double initialValue = static_cast(metric->GetValue()); + + // Drive a real ITK optimizer with the metric. + using OptimizerType = itk::RegularStepGradientDescentOptimizerv4; + using ScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift; + auto scales = ScalesEstimatorType::New(); + scales->SetMetric(metric); + auto optimizer = OptimizerType::New(); + optimizer->SetMetric(metric); + optimizer->SetScalesEstimator(scales); + optimizer->SetNumberOfIterations(300); + optimizer->SetLearningRate(2.0); + optimizer->SetMinimumStepLength(1e-4); + optimizer->SetRelaxationFactor(0.8); + optimizer->StartOptimization(); + + // The metric, through its analytic derivative, must drive the optimizer downhill. + const double finalValue = static_cast(metric->GetValue()); + EXPECT_LT(finalValue, initialValue) << distance << ": optimization should reduce the metric"; + + // The optimizer recovers the true geometric alignment: bringing the moving blob + // (x=10) onto the fixed one (x=6) is a +4 voxel translation in x, with no y/z motion. + const TransformType::ParametersType recovered = transform->GetParameters(); + EXPECT_NEAR(recovered[0], 4.0, 0.5) << distance << ": recovered x-translation should reach +4"; + EXPECT_LT(std::abs(recovered[1]), 0.2) << distance << ": motion should stay on the x axis"; + EXPECT_LT(std::abs(recovered[2]), 0.2) << distance << ": motion should stay on the x axis"; } } @@ -1510,7 +1510,7 @@ TEST(ImpactBackend, MultiLayerFeatureMapsAreIndependent) auto interp = InterpolatorType::New(); interp->SetSplineOrder(3); itk::ImpactModelConfiguration cfg(ToyModelPath(), 3, 1, { 0, 0, 0 }, { 1.f, 1.f, 1.f }, { 2, 2, 2 }, mask, false); - auto f = FeatMapType::New(); + auto f = FeatMapType::New(); f->SetModelConfiguration(cfg); f->SetInterpolator(interp); f->AddInput(img); @@ -1586,7 +1586,7 @@ TEST(ImpactBackend, TiledFeatureMapReconstructsInputForEveryWindow) f->SetDevice("cpu"); ASSERT_NO_THROW(f->Update()) << window << " / overlap " << overlap; - auto map = f->GetOutput(1); // passthrough + auto map = f->GetOutput(1); // passthrough const std::string where = window + " / overlap " + std::to_string(overlap); ASSERT_EQ(map->GetLargestPossibleRegion().GetSize(), size) << where << ": map must lie on the input grid"; for (unsigned int d = 0; d < 3; ++d) @@ -1632,10 +1632,10 @@ TEST(ImpactBackend, PerAxisOverlapBlendsAndRejectsAnOverlapAsLargeAsThePatch) return f; }; - auto f = run({ 6, 0, 3 }); - auto map = f->GetOutput(1); + auto f = run({ 6, 0, 3 }); + auto map = f->GetOutput(1); itk::ImageRegionConstIteratorWithIndex it(img, img->GetLargestPossibleRegion()); - double maxDiff = 0.0; + double maxDiff = 0.0; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { maxDiff = std::max(maxDiff, std::abs(static_cast(map->GetPixel(it.GetIndex())[0]) - it.Get())); @@ -1659,7 +1659,7 @@ TEST(ImpactBackend, AnisotropicPatchIsTiledInItkAxisOrder) auto run = [&](const std::vector & patch) { itk::ImpactModelConfiguration cfg( ToyModelPath(), 3, 1, patch, { 1.f, 1.f, 1.f }, { 0, 0, 0 }, { true, true }, false); - auto interp = InterpolatorType::New(); + auto interp = InterpolatorType::New(); interp->SetSplineOrder(3); auto f = FeatMapType::New(); f->SetModelConfiguration(cfg); @@ -1710,14 +1710,19 @@ struct OnlineFixture std::vector configs; std::vector> points; std::vector>>> patchIndex{ 1 }; - std::vector subset{ torch::arange(4, torch::kInt64), - torch::arange(1, torch::kInt64) }; - torch::Device device{ torch::kCPU }; + std::vector subset{ torch::arange(4, torch::kInt64), torch::arange(1, torch::kInt64) }; + torch::Device device{ torch::kCPU }; explicit OnlineFixture(size_t samples) { - configs.emplace_back(ToyModelPath(), 3, 1, std::vector{ 3, 3, 3 }, std::vector{ 1.f, 1.f, 1.f }, - std::vector{ 0, 0, 0 }, std::vector{ true, true }, false); + configs.emplace_back(ToyModelPath(), + 3, + 1, + std::vector{ 3, 3, 3 }, + std::vector{ 1.f, 1.f, 1.f }, + std::vector{ 0, 0, 0 }, + std::vector{ true, true }, + false); itk::Impact::GetModelOutputsExample(configs, "fixed", device); for (size_t s = 0; s < samples; ++s) { @@ -1920,7 +1925,7 @@ TEST(ImpactBackend, OnlineInferenceShapesAnAnisotropicPatchForTheModel) std::vector{ true, false }, // keep the conv: it is not symmetric false); - const torch::Device device(torch::kCPU); + const torch::Device device(torch::kCPU); std::vector configs{ config }; // Probes the model and records where the center of each output layer is; GenerateOutputs // needs it, and it is the same call the metric makes before its first evaluation. @@ -1968,8 +1973,8 @@ TEST(ImpactBackend, OnlineInferenceShapesAnAnisotropicPatchForTheModel) ASSERT_EQ(patchIndex[0][0].size(), static_cast(px * py * pz)); const std::vector subset{ torch::arange(4, torch::kInt64) }; - std::vector outputs = itk::Impact::GenerateOutputs>( - configs, fixedPoints, patchIndex, subset, device, evaluator); + std::vector outputs = + itk::Impact::GenerateOutputs>(configs, fixedPoints, patchIndex, subset, device, evaluator); ASSERT_EQ(outputs.size(), 1u); ASSERT_EQ(outputs[0].sizes(), (std::vector{ 1, 4 })); @@ -1983,7 +1988,7 @@ TEST(ImpactBackend, OnlineInferenceShapesAnAnisotropicPatchForTheModel) TEST(ImpactBackend, TwoDimensionalModelIsSweptOverAVolume) { using FeatMapType = itk::ImageToFeaturesMap; - const std::string model2d = std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModel2D.pt"; + const std::string model2d = std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModel2D.pt"; const ImageType::SizeType size = { { 20, 17, 6 } }; // x, y tiled; z swept auto img = MakeRampImage(size); ImageType::SpacingType spacing; @@ -2303,7 +2308,7 @@ TEST(ImpactBackend, DownsamplingEncoderFeatureMapOverlaysInput) ASSERT_NO_THROW(f->Update()); // Layer 0 is an AvgPool (position-preserving), so its peak marks the true location. - auto map0 = f->GetOutput(0); + auto map0 = f->GetOutput(0); const auto size = map0->GetLargestPossibleRegion().GetSize(); const auto spacing = map0->GetSpacing(); const auto origin = map0->GetOrigin(); @@ -2314,10 +2319,19 @@ TEST(ImpactBackend, DownsamplingEncoderFeatureMapOverlaysInput) EXPECT_NEAR(origin[d], 0.0, 1e-6) << "origin should match the input (no half-voxel shift)"; } - itk::ImageRegionConstIteratorWithIndex> it( - map0, map0->GetLargestPossibleRegion()); - double best = -1; std::remove_reference_t::IndexType bi{}; - for (it.GoToBegin(); !it.IsAtEnd(); ++it) { const double v = std::abs(it.Get()[0]); if (v > best) { best = v; bi = it.GetIndex(); } } + itk::ImageRegionConstIteratorWithIndex> it(map0, + map0->GetLargestPossibleRegion()); + double best = -1; + std::remove_reference_t::IndexType bi{}; + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const double v = std::abs(it.Get()[0]); + if (v > best) + { + best = v; + bi = it.GetIndex(); + } + } std::remove_reference_t::PointType bp; map0->TransformIndexToPhysicalPoint(bi, bp); // The downsampled marker must still be at the input's physical position (within one @@ -2397,7 +2411,7 @@ TEST(ImpactLoss, ForwardValueGradientFiniteDifference) for (const std::string name : { "L1", "L2", "Cosine", "L1Cosine", "Dice", "NCC" }) { auto loss = itk::Impact::LossFactory::Instance().Create(name); - torch::Tensor fixed = torch::rand({ N, C }, opts) + 0.2; // strictly positive (Dice-safe) + torch::Tensor fixed = torch::rand({ N, C }, opts) + 0.2; // strictly positive (Dice-safe) torch::Tensor moving0 = torch::rand({ N, C }, opts) + 0.2; torch::Tensor moving = moving0.clone().set_requires_grad(true); @@ -2455,8 +2469,8 @@ MakeTorchAdamPattern(unsigned int n, const ImageType::SpacingType & spacing, const ImageType::DirectionType & direction) { - auto image = ImageType::New(); - ImageType::SizeType size; + auto image = ImageType::New(); + ImageType::SizeType size; size.Fill(n); ImageType::RegionType region; region.SetSize(size); @@ -2478,8 +2492,8 @@ InteriorMeanError(const TorchAdamFieldType * field, const itk::Vector { itk::Vector acc; acc.Fill(0.0); - const auto size = field->GetLargestPossibleRegion().GetSize(); - long count = 0; + const auto size = field->GetLargestPossibleRegion().GetSize(); + long count = 0; itk::ImageRegionConstIteratorWithIndex it(field, field->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { @@ -2539,7 +2553,7 @@ TEST(ImpactTorchAdam, ZeroIterationIdentity) filter->SetNumberOfIterations(0); filter->Update(); - double maxField = 0.0; + double maxField = 0.0; itk::ImageRegionConstIteratorWithIndex it( filter->GetDisplacementField(), filter->GetDisplacementField()->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) @@ -2547,12 +2561,11 @@ TEST(ImpactTorchAdam, ZeroIterationIdentity) maxField = std::max(maxField, static_cast(std::abs(it.Get()[d]))); EXPECT_LT(maxField, 1e-6) << "field must be exactly 0 with no optimization (variant " << variant << ")"; - double maxWarpDiff = 0.0; - const auto * warped = filter->GetWarpedMovingImage(); + double maxWarpDiff = 0.0; + const auto * warped = filter->GetWarpedMovingImage(); itk::ImageRegionConstIteratorWithIndex wi(warped, warped->GetLargestPossibleRegion()); for (wi.GoToBegin(); !wi.IsAtEnd(); ++wi) - maxWarpDiff = - std::max(maxWarpDiff, std::abs(static_cast(wi.Get() - fixed->GetPixel(wi.GetIndex())))); + maxWarpDiff = std::max(maxWarpDiff, std::abs(static_cast(wi.Get() - fixed->GetPixel(wi.GetIndex())))); EXPECT_LT(maxWarpDiff, 1e-4) << "identity grid_sample should reproduce the input (variant " << variant << ")"; } } @@ -2605,7 +2618,7 @@ TEST(ImpactTorchAdam, TranslationRecoveryIntensity) v += direction[r][c] * spacing[c] * ((c == 0) ? tx : (c == 1) ? ty : tz); expected[r] = v; } - const auto err = InteriorMeanError(filter->GetDisplacementField(), expected, 6); + const auto err = InteriorMeanError(filter->GetDisplacementField(), expected, 6); const auto & hist = filter->GetMetricValuesPerIteration(); EXPECT_LT(hist.back(), 0.1 * hist.front()) << "loss should drop (variant " << variant << ")"; @@ -2620,8 +2633,8 @@ TEST(ImpactTorchAdam, TranslationRecoveryIntensity) // (toy TorchScript model) via itk::Forward + ImpactLoss::forwardValue. TEST(ImpactTorchAdam, TranslationRecoveryFeatures) { - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -2660,7 +2673,7 @@ TEST(ImpactTorchAdam, TranslationRecoveryFeatures) // is 0 on the zero starting field), and it drops as the field moves. TEST(ImpactTorchAdam, NormalizedLossStartsAtTheSumOfTheWeights) { - ImageType::SpacingType spacing; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -2686,15 +2699,12 @@ TEST(ImpactTorchAdam, NormalizedLossStartsAtTheSumOfTheWeights) // normalization of the features; each still recovers a known translation. TEST(ImpactTorchAdam, SharedFeatureSettingsStillRecoverATranslation) { - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); - auto run = [&](const std::string & distance, - unsigned int subset, - float voxel, - const std::string & normalization) { + auto run = [&](const std::string & distance, unsigned int subset, float voxel, const std::string & normalization) { auto filter = TorchAdamFilterType::New(); filter->SetFixedImage(MakeTorchAdamPattern(20, 0, 0, 0, spacing, identity)); filter->SetMovingImage(MakeTorchAdamPattern(20, tx, ty, tz, spacing, identity)); @@ -2717,9 +2727,9 @@ TEST(ImpactTorchAdam, SharedFeatureSettingsStillRecoverATranslation) }; for (const auto & [label, error] : std::vector>>{ { "subset 2 of 4", run("L2", 2, 0.f, "none") }, - { "LNCC", run("LNCC", 0, 0.f, "none") }, - { "voxel size 2", run("L2", 0, 2.f, "none") }, - { "l2 normalization", run("L2", 0, 0.f, "l2") } }) + { "LNCC", run("LNCC", 0, 0.f, "none") }, + { "voxel size 2", run("L2", 0, 2.f, "none") }, + { "l2 normalization", run("L2", 0, 0.f, "l2") } }) { for (unsigned int d = 0; d < 3; ++d) EXPECT_LT(error[d], 0.5) << label << ": axis " << d << " off"; @@ -3086,15 +3096,22 @@ TEST(ImpactTorchAdam, SampledJacobianRecoversATranslation) // fixed features: a voxel size that resamples the image would hand back chunks of another grid, so it is refused. TEST(ImpactTorchAdam, ChunkedJacobianRefusesAResamplingModel) { - ImageType::SpacingType spacing; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); auto filter = TorchAdamFilterType::New(); filter->SetFixedImage(MakeTorchAdamPattern(16, 0, 0, 0, spacing, identity)); filter->SetMovingImage(MakeTorchAdamPattern(16, 1, 0, 0, spacing, identity)); - filter->AddModelConfiguration(itk::ImpactModelConfiguration(std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModel2D.pt", - 2, 1, { 0, 0 }, { 2.f, 2.f, 2.f }, { 0, 0 }, { true, false }, false)); + filter->AddModelConfiguration( + itk::ImpactModelConfiguration(std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModel2D.pt", + 2, + 1, + { 0, 0 }, + { 2.f, 2.f, 2.f }, + { 0, 0 }, + { true, false }, + false)); filter->SetDistance({ "L2" }); filter->SetMode("Jacobian"); filter->SetFeatureChunkSize(4); @@ -3195,8 +3212,8 @@ TEST(ImpactTorchAdam, TwoDimensionalJacobianSweepsARandomAxisEachIteration) TEST(ImpactTorchAdam, FeaturePCA) { - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3233,8 +3250,8 @@ TEST(ImpactTorchAdam, FeaturePCA) // the refresh path must run without error and a translation must still be recovered. TEST(ImpactTorchAdam, FeatureMapUpdateInterval) { - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3271,9 +3288,9 @@ TEST(ImpactTorchAdam, FeatureMapUpdateInterval) // the layers are upsampled back to the input resolution and a translation is still recovered. TEST(ImpactTorchAdam, DownsampledFeatureLayers) { - const std::string downModel = std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModelDown.pt"; - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const std::string downModel = std::string(IMPACT_TEST_DATA_DIR) + "/ImpactToyModelDown.pt"; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3310,8 +3327,8 @@ TEST(ImpactTorchAdam, DownsampledFeatureLayers) // units path and the control-grid smoothing. TEST(ImpactTorchAdam, LowResControlGrid) { - const double tx = 1.5, ty = -2.0, tz = 1.0; - ImageType::SpacingType spacing; + const double tx = 1.5, ty = -2.0, tz = 1.0; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3346,7 +3363,7 @@ TEST(ImpactTorchAdam, LowResControlGrid) // axis rather than averaging over zero elements (NaN). TEST(ImpactTorchAdam, DegenerateControlGridFinite) { - ImageType::SpacingType spacing; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3362,7 +3379,7 @@ TEST(ImpactTorchAdam, DegenerateControlGridFinite) filter->SetRegularizationWeight(1.0); filter->Update(); - bool finite = true; + bool finite = true; itk::ImageRegionConstIteratorWithIndex it( filter->GetDisplacementField(), filter->GetDisplacementField()->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) @@ -3382,7 +3399,7 @@ TEST(ImpactConvexAdam, CoarseInitThenAdamRefines) { using CoarseType = itk::ImpactCoarseRegistration; - ImageType::SpacingType spacing; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3446,7 +3463,7 @@ TEST(ImpactConvexAdam, CoarseInitThenAdamRefines) TEST(ImpactConvexAdam, CoarseInverseConsistency) { using CoarseType = itk::ImpactCoarseRegistration; - ImageType::SpacingType spacing; + ImageType::SpacingType spacing; spacing.Fill(1.0); ImageType::DirectionType identity; identity.SetIdentity(); @@ -3526,7 +3543,7 @@ TEST(ImpactConvexAdam, NormalizedCoarseCostIgnoresTheFeatureRange) ImageType::DirectionType identity; identity.SetIdentity(); auto scaled = [&](double factor, double tx, double ty, double tz) { - auto image = MakeTorchAdamPattern(24, tx, ty, tz, spacing, identity); + auto image = MakeTorchAdamPattern(24, tx, ty, tz, spacing, identity); itk::ImageRegionIteratorWithIndex it(image, image->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) it.Set(static_cast(it.Get() * factor)); @@ -3545,8 +3562,8 @@ TEST(ImpactConvexAdam, NormalizedCoarseCostIgnoresTheFeatureRange) return FieldType::Pointer(coarse->GetDisplacementField()); }; auto meanDifference = [](const FieldType * a, const FieldType * b) { - double sum = 0; - long n = 0; + double sum = 0; + long n = 0; itk::ImageRegionConstIteratorWithIndex it(a, a->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it, ++n) for (unsigned int d = 0; d < 3; ++d) @@ -4385,7 +4402,7 @@ TEST(ImpactConvexAdam, CoarseTakesPCAAndASeededChannelSubset) return FieldType::Pointer(coarse->GetDisplacementField()); }; auto maxDifference = [](const FieldType * a, const FieldType * b) { - double worst = 0; + double worst = 0; itk::ImageRegionConstIteratorWithIndex it(a, a->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) for (unsigned int d = 0; d < 3; ++d) @@ -4411,9 +4428,9 @@ namespace double InteriorNCC(const ImageType * a, const ImageType * b, int margin) { - const auto size = a->GetLargestPossibleRegion().GetSize(); - double sa = 0, sb = 0, saa = 0, sbb = 0, sab = 0; - long n = 0; + const auto size = a->GetLargestPossibleRegion().GetSize(); + double sa = 0, sb = 0, saa = 0, sbb = 0, sab = 0; + long n = 0; itk::ImageRegionConstIteratorWithIndex it(a, a->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { @@ -4458,9 +4475,9 @@ ResampleOnto(const ImageType * moving, const ImageType * reference) double FoldedFraction(const TorchAdamFieldType * field, int margin) { - const auto size = field->GetLargestPossibleRegion().GetSize(); - const auto spacing = field->GetSpacing(); - long folded = 0, total = 0; + const auto size = field->GetLargestPossibleRegion().GetSize(); + const auto spacing = field->GetSpacing(); + long folded = 0, total = 0; itk::ImageRegionConstIteratorWithIndex it(field, field->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { @@ -4578,8 +4595,7 @@ TEST(ImpactConvexAdam, RealLungCTCoarseInitThenAdamRefines) const double folded = FoldedFraction(fine->GetDisplacementField(), 2); std::cout << "[real-CT] device=" << device << " shrink=" << shrink << " loss " << hist.front() << " -> " - << hist.back() << " | NCC moving=" << s0 << " warped=" << s1 << " | foldedFraction=" << folded - << std::endl; + << hist.back() << " | NCC moving=" << s0 << " warped=" << s1 << " | foldedFraction=" << folded << std::endl; EXPECT_LT(hist.back(), 0.85 * hist.front()) << "feature loss should drop"; EXPECT_GT(s1, s0 + 0.02) << "warped moving should correlate better with fixed than the unregistered moving"; @@ -4596,8 +4612,8 @@ TEST(ImpactConvexAdam, Coarse2D) const double tx = 4.0, ty = -3.0; // index shift (x, y) auto make = [&](double sx, double sy) { - auto img = Image2D::New(); - Image2D::SizeType size; + auto img = Image2D::New(); + Image2D::SizeType size; size.Fill(n); Image2D::RegionType region; region.SetSize(size); @@ -4623,10 +4639,10 @@ TEST(ImpactConvexAdam, Coarse2D) coarse->SetDisplacementHalfWidth(4); coarse->Update(); - auto field = coarse->GetDisplacementField(); - const auto size = field->GetLargestPossibleRegion().GetSize(); - double ex = 0, ey = 0; - long cnt = 0; + auto field = coarse->GetDisplacementField(); + const auto size = field->GetLargestPossibleRegion().GetSize(); + double ex = 0, ey = 0; + long cnt = 0; itk::ImageRegionConstIteratorWithIndex it(field, field->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { diff --git a/test/itkInterpolateVectorImageFunctionTest.cxx b/test/itkInterpolateVectorImageFunctionTest.cxx index ba714da..fef54c9 100644 --- a/test/itkInterpolateVectorImageFunctionTest.cxx +++ b/test/itkInterpolateVectorImageFunctionTest.cxx @@ -38,8 +38,8 @@ itkInterpolateVectorImageFunctionTest(int, char *[]) using InterpolatorType = itk::BSplineInterpolateImageFunction, float, float>; using VectorInterpolatorType = itk::InterpolateVectorImageFunction; - auto image = FeaturesImageType::New(); - FeaturesImageType::SizeType size; + auto image = FeaturesImageType::New(); + FeaturesImageType::SizeType size; size.Fill(8); FeaturesImageType::RegionType region; region.SetSize(size); diff --git a/test/itkTensorToImageFilterTest.cxx b/test/itkTensorToImageFilterTest.cxx index 4c8c99c..afb4060 100644 --- a/test/itkTensorToImageFilterTest.cxx +++ b/test/itkTensorToImageFilterTest.cxx @@ -34,8 +34,8 @@ itkTensorToImageFilterTest(int, char *[]) using TensorToImageFilterType = itk::TensorToImageFilter; // Reference image only supplies geometry (origin/spacing/direction). - auto reference = ImageType::New(); - ImageType::SizeType size; + auto reference = ImageType::New(); + ImageType::SizeType size; size.Fill(8); ImageType::RegionType region; region.SetSize(size);