From ccee041e6885b9e1d478946d5e89067815eb54b7 Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Sat, 20 Jun 2026 19:45:13 -0500 Subject: [PATCH 1/2] PERF: Hoist constant tensor-basis SVD out of the per-voxel DTI loop DiffusionTensor3DReconstructionImageFilter constructed a vnl_svd of the constant m_TensorBasis once per voxel (~16.7M times for a 256^3 DWI) to solve the Stejskal-Tanner equations. The pseudo-inverse (the dual tensor basis) is the same for every voxel, so precompute it once in BeforeThreadedGenerateData and apply it per voxel as a matrix-vector product. vnl_svd::solve(X) is pinverse()*X, so the reconstructed tensors are unchanged. Removing the per-voxel SVD also frees the loop of the thread-unsafe netlib dsvdc call, so the stale single-thread warning is dropped. --- ...ffusionTensor3DReconstructionImageFilter.h | 4 ++++ ...usionTensor3DReconstructionImageFilter.hxx | 20 +++++++++---------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h index 5f1ebde8e1b..9d2283c994b 100644 --- a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h +++ b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h @@ -330,6 +330,10 @@ class ITK_TEMPLATE_EXPORT DiffusionTensor3DReconstructionImageFilter /* Tensor basis coeffs */ TensorBasisMatrixType m_TensorBasis{}; + /* Pseudo-inverse (dual tensor basis) of m_TensorBasis, precomputed once so the + * per-voxel loop is a matrix-vector product rather than a fresh SVD per voxel. */ + vnl_matrix m_TensorBasisInverse{}; + CoefficientMatrixType m_BMatrix{}; /** container to hold gradient directions */ diff --git a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx index 57d5cc73366..4432b883de1 100644 --- a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx +++ b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx @@ -154,10 +154,6 @@ DiffusionTensor3DReconstructionImageFilter pseudoInverseSolver{ m_TensorBasis.as_matrix() }; if (m_NumberOfGradientDirections > 6) { - D = pseudoInverseSolver.solve(m_BMatrix * B); + D = m_TensorBasisInverse * (m_BMatrix * B); } else { - D = pseudoInverseSolver.solve(B); + D = m_TensorBasisInverse * B; } tensor(0, 0) = D[0]; @@ -376,14 +371,13 @@ DiffusionTensor3DReconstructionImageFilter pseudoInverseSolver{ m_TensorBasis.as_matrix() }; if (m_NumberOfGradientDirections > 6) { - D = pseudoInverseSolver.solve(m_BMatrix * B); + D = m_TensorBasisInverse * (m_BMatrix * B); } else { - D = pseudoInverseSolver.solve(B); + D = m_TensorBasisInverse * B; } tensor(0, 0) = D[0]; @@ -451,6 +445,12 @@ DiffusionTensor3DReconstructionImageFilter{ m_TensorBasis.as_matrix() }.pinverse(); + m_BMatrix.inplace_transpose(); } From ee656f3d61c0f1bc71d14ffb686be1aeafe6e0d1 Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Sat, 20 Jun 2026 21:42:04 -0500 Subject: [PATCH 2/2] ENH: Add single-image-path test for DTI tensor reconstruction The existing DiffusionTensor3DReconstructionImageFilter test only exercises the multi-image gradient path (AddGradientImage). Add a test that reconstructs the same data through the single-VectorImage path (SetGradientImage) and verifies it produces identical tensors, covering the previously-untested code path. --- .../DiffusionTensorImage/test/CMakeLists.txt | 7 + ...nTensor3DReconstructionSingleImageTest.cxx | 139 ++++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100644 Modules/Filtering/DiffusionTensorImage/test/itkDiffusionTensor3DReconstructionSingleImageTest.cxx diff --git a/Modules/Filtering/DiffusionTensorImage/test/CMakeLists.txt b/Modules/Filtering/DiffusionTensorImage/test/CMakeLists.txt index c0091016bb2..66f889917db 100644 --- a/Modules/Filtering/DiffusionTensorImage/test/CMakeLists.txt +++ b/Modules/Filtering/DiffusionTensorImage/test/CMakeLists.txt @@ -2,6 +2,7 @@ itk_module_test() set(ITKDiffusionTensorImageTests itkDiffusionTensor3DTest.cxx itkDiffusionTensor3DReconstructionImageFilterTest.cxx + itkDiffusionTensor3DReconstructionSingleImageTest.cxx itkTensorRelativeAnisotropyImageFilterTest.cxx itkTensorFractionalAnisotropyImageFilterTest.cxx) @@ -13,6 +14,12 @@ itk_add_test( COMMAND ITKDiffusionTensorImageTestDriver itkDiffusionTensor3DTest) +itk_add_test( + NAME + itkDiffusionTensor3DReconstructionSingleImageTest + COMMAND + ITKDiffusionTensorImageTestDriver + itkDiffusionTensor3DReconstructionSingleImageTest) itk_add_test( NAME itkDiffusionTensor3DReconstructionImageFilterTest diff --git a/Modules/Filtering/DiffusionTensorImage/test/itkDiffusionTensor3DReconstructionSingleImageTest.cxx b/Modules/Filtering/DiffusionTensorImage/test/itkDiffusionTensor3DReconstructionSingleImageTest.cxx new file mode 100644 index 00000000000..d4207c3fbe5 --- /dev/null +++ b/Modules/Filtering/DiffusionTensorImage/test/itkDiffusionTensor3DReconstructionSingleImageTest.cxx @@ -0,0 +1,139 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +// Exercises the single-VectorImage gradient path (SetGradientImage), which the +// existing reconstruction test does not cover, and verifies it produces the same +// tensors as the multi-image path (AddGradientImage) on identical data. + +#include "itkDiffusionTensor3DReconstructionImageFilter.h" +#include "itkImageRegionIterator.h" +#include +#include "itkTestingMacros.h" + +int +itkDiffusionTensor3DReconstructionSingleImageTest(int, char *[]) +{ + using ReferencePixelType = short; + using GradientPixelType = short; + using TensorPrecisionType = double; + using FilterType = + itk::DiffusionTensor3DReconstructionImageFilter; + + constexpr unsigned int numberOfGradientImages = 6; + const double gradientDirections[numberOfGradientImages][3] = { + { -1.000000, 0.000000, 0.000000 }, { -0.166000, 0.986000, 0.000000 }, { 0.110000, 0.664000, 0.740000 }, + { -0.901000, -0.419000, -0.110000 }, { 0.169000, -0.601000, 0.781000 }, { 0.815000, -0.386000, 0.433000 } + }; + constexpr ReferencePixelType b0 = 100; + auto gradientValue = [](unsigned int i) { return static_cast((i + 1) * (i + 1) * (i + 1)); }; + + using ReferenceImageType = FilterType::ReferenceImageType; + using RegionType = ReferenceImageType::RegionType; + constexpr RegionType::SizeType imageSize{ 4, 4, 4 }; + const RegionType region{ imageSize }; + + // --- Path A: multi-image (AddGradientImage) --- + auto multiFilter = FilterType::New(); + { + auto reference = ReferenceImageType::New(); + reference->SetRegions(region); + reference->Allocate(); + reference->FillBuffer(b0); + multiFilter->SetReferenceImage(reference); + for (unsigned int i = 0; i < numberOfGradientImages; ++i) + { + using GradientImageType = FilterType::GradientImageType; + auto gradientImage = GradientImageType::New(); + gradientImage->SetRegions(region); + gradientImage->Allocate(); + gradientImage->FillBuffer(gradientValue(i)); + FilterType::GradientDirectionType d; + d[0] = gradientDirections[i][0]; + d[1] = gradientDirections[i][1]; + d[2] = gradientDirections[i][2]; + multiFilter->AddGradientImage(d, gradientImage); + } + } + multiFilter->SetBValue(1000.0); + multiFilter->Update(); + + // --- Path B: single VectorImage (SetGradientImage) with the SAME data --- + // Component 0 is the b0 baseline (direction (0,0,0)); components 1..n are the + // diffusion-weighted measurements. + auto singleFilter = FilterType::New(); + { + using VectorImageType = itk::VectorImage; + auto vectorImage = VectorImageType::New(); + vectorImage->SetRegions(region); + vectorImage->SetVectorLength(numberOfGradientImages + 1); + vectorImage->Allocate(); + itk::VariableLengthVector value(numberOfGradientImages + 1); + value[0] = b0; + for (unsigned int i = 0; i < numberOfGradientImages; ++i) + { + value[i + 1] = gradientValue(i); + } + vectorImage->FillBuffer(value); + + auto directions = FilterType::GradientDirectionContainerType::New(); + FilterType::GradientDirectionType zero{}; + zero.fill(0.0); + directions->InsertElement(0, zero); + for (unsigned int i = 0; i < numberOfGradientImages; ++i) + { + FilterType::GradientDirectionType d; + d[0] = gradientDirections[i][0]; + d[1] = gradientDirections[i][1]; + d[2] = gradientDirections[i][2]; + directions->InsertElement(i + 1, d); + } + singleFilter->SetGradientImage(directions, vectorImage); + } + singleFilter->SetBValue(1000.0); + singleFilter->Update(); + + // --- The two paths must reconstruct identical tensors --- + using TensorImageType = FilterType::TensorImageType; + itk::ImageRegionConstIterator mit(multiFilter->GetOutput(), + multiFilter->GetOutput()->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator sit(singleFilter->GetOutput(), + singleFilter->GetOutput()->GetLargestPossibleRegion()); + double maxDiff = 0.0; + constexpr double tolerance = 1e-8; + bool anyNonZero = false; + for (mit.GoToBegin(), sit.GoToBegin(); !mit.IsAtEnd(); ++mit, ++sit) + { + for (unsigned int k = 0; k < 6; ++k) + { + maxDiff = std::max(maxDiff, std::abs(mit.Get()[k] - sit.Get()[k])); + anyNonZero = anyNonZero || (std::abs(sit.Get()[k]) > 1e-12); + } + } + + std::cout << "SingleImage vs MultiImage max tensor difference: " << maxDiff << std::endl; + ITK_TEST_EXPECT_TRUE(anyNonZero); + if (maxDiff > tolerance) + { + std::cerr << "Single-image path disagrees with multi-image path (max diff " << maxDiff << " > " << tolerance << ')' + << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +}