From 6611471c5907e41f372ca6d518983bffeb670f54 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 6 Jan 2026 14:12:15 +0100 Subject: [PATCH 1/5] ENH: Move rtkMatlabSparseMatrix to include directory --- .../rtkMatlabSparseMatrix.h | 16 +- .../rtkMatlabSparseMatrix.hxx | 148 ++---------------- 2 files changed, 28 insertions(+), 136 deletions(-) rename {applications/rtkprojectionmatrix => include}/rtkMatlabSparseMatrix.h (91%) rename {applications/rtkprojectionmatrix => include}/rtkMatlabSparseMatrix.hxx (60%) diff --git a/applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.h b/include/rtkMatlabSparseMatrix.h similarity index 91% rename from applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.h rename to include/rtkMatlabSparseMatrix.h index de3dab3ee..a27b08489 100644 --- a/applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.h +++ b/include/rtkMatlabSparseMatrix.h @@ -19,25 +19,32 @@ #ifndef rtkMatlabSparseMatrix_h #define rtkMatlabSparseMatrix_h +#include "rtkConfiguration.h" #include namespace rtk { + /** \class MatlabSparseMatrix * - * sparse matrix in Matlab format - * Initilaize it with vnl_sparse_matrix - * Save it into .mat format + * \brief Sparse matrix in Matlab format + * + * Converts a vnl_sparse_matrix to Matlab .mat format. + * Initialize it with vnl_sparse_matrix and save it into .mat format. * + * \author Thomas Baudier + * + * \ingroup RTK */ - class MatlabSparseMatrix { public: template MatlabSparseMatrix(const vnl_sparse_matrix & sparseMatrix, TOutputImage * output); + void Save(std::ostream & out); + void Print(); @@ -79,6 +86,7 @@ class MatlabSparseMatrix protected: MatlabSparseMatrixStruct m_MatlabSparseMatrix; }; + } // namespace rtk #ifndef ITK_MANUAL_INSTANTIATION diff --git a/applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.hxx b/include/rtkMatlabSparseMatrix.hxx similarity index 60% rename from applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.hxx rename to include/rtkMatlabSparseMatrix.hxx index f2e9febb3..83c8c9ae6 100644 --- a/applications/rtkprojectionmatrix/rtkMatlabSparseMatrix.hxx +++ b/include/rtkMatlabSparseMatrix.hxx @@ -16,12 +16,19 @@ * *=========================================================================*/ +#ifndef rtkMatlabSparseMatrix_hxx +#define rtkMatlabSparseMatrix_hxx + #include #include +#include + +namespace rtk +{ //==================================================================== template -rtk::MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseMatrix, TOutputImage * output) +MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseMatrix, TOutputImage * output) { unsigned int nbColumn = output->GetLargestPossibleRegion().GetSize()[0] * output->GetLargestPossibleRegion().GetSize()[2]; // it's not sparseMatrix.columns() @@ -106,10 +113,12 @@ rtk::MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sp } } m_MatlabSparseMatrix.s_columnIndex[m_MatlabSparseMatrix.s_dimensionNbColumn] = m_MatlabSparseMatrix.s_arrayNzmax; + + delete[] columnsVector; } -void -rtk::MatlabSparseMatrix::Save(std::ostream & out) +inline void +MatlabSparseMatrix::Save(std::ostream & out) { // Write data out.write(reinterpret_cast(&m_MatlabSparseMatrix.s_headerMatlab), 116); @@ -150,8 +159,8 @@ rtk::MatlabSparseMatrix::Save(std::ostream & out) out.write(reinterpret_cast(&m_MatlabSparseMatrix.s_value[i]), 8); } -void -rtk::MatlabSparseMatrix::Print() +inline void +MatlabSparseMatrix::Print() { std::cout << "m_MatlabSparseMatrix.s_headerMatlab : \"" << m_MatlabSparseMatrix.s_headerMatlab << "\"" << std::endl; std::cout << "m_MatlabSparseMatrix.s_headerOffset : \"" << m_MatlabSparseMatrix.s_headerOffset << "\"" << std::endl; @@ -209,131 +218,6 @@ rtk::MatlabSparseMatrix::Print() std::cout << "\"" << std::endl; } -/* -std::istream& operator>>(std::istream& in) { - char temp1(0); - unsigned short int temp2(0); - unsigned long int temp4(0); - double temp8(0); - - //header - for (unsigned int i=0; i<116; ++i) { - in.read(reinterpret_cast(&temp1),1); - m_MatlabSparseMatrix.s_headerMatlab[i] = temp1; - } - std::cout << "m_MatlabSparseMatrix.s_headerMatlab : \"" << m_MatlabSparseMatrix.s_headerMatlab << "\"" << std::endl; - for (unsigned int i=0; i<8; ++i) { - in.read(reinterpret_cast(&temp1),1); - m_MatlabSparseMatrix.s_headerOffset[i] = temp1; - } - std::cout << "m_MatlabSparseMatrix.s_headerOffset : \"" << m_MatlabSparseMatrix.s_headerOffset << "\"" << std::endl; - in.read(reinterpret_cast(&temp2), 2); - m_MatlabSparseMatrix.s_headerVersion = temp2; - std::cout << "m_MatlabSparseMatrix.s_headerVersion : \"" << std::hex << +m_MatlabSparseMatrix.s_headerVersion << "\"" -<< std::endl; for (unsigned int i=0; i<2; ++i) { in.read(reinterpret_cast(&temp1),1); - m_MatlabSparseMatrix.s_headerEndian[i] = temp1; - } - std::cout << "m_MatlabSparseMatrix.s_headerEndian : \"" << m_MatlabSparseMatrix.s_headerEndian << "\"" << std::endl; - //data - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_mainTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_mainTag : \"" << m_MatlabSparseMatrix.s_mainTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_dataLength = temp4; - std::cout << "m_MatlabSparseMatrix.s_dataLength : \"" << m_MatlabSparseMatrix.s_dataLength << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_arrayTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_arrayTag : \"" << m_MatlabSparseMatrix.s_arrayTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_arrayLength = temp4; - std::cout << "m_MatlabSparseMatrix.s_arrayLength : \"" << m_MatlabSparseMatrix.s_arrayLength << "\"" << std::endl; - in.read(reinterpret_cast(&temp1), 1); - m_MatlabSparseMatrix.s_arrayFlags = temp1; - std::cout << "m_MatlabSparseMatrix.s_arrayFlags : \"" << (unsigned int)m_MatlabSparseMatrix.s_arrayFlags << "\"" << -std::endl; in.read(reinterpret_cast(&temp1),1); m_MatlabSparseMatrix.s_arrayClass = temp1; std::cout << -"m_MatlabSparseMatrix.s_arrayClass : \"" << (unsigned int)m_MatlabSparseMatrix.s_arrayClass << "\"" << std::endl; - in.read(reinterpret_cast(&temp2), 2); - m_MatlabSparseMatrix.s_arrayUndefined = temp2; - std::cout << "m_MatlabSparseMatrix.s_arrayUndefined : \"" << m_MatlabSparseMatrix.s_arrayUndefined << "\"" << -std::endl; in.read(reinterpret_cast(&temp4), 4); m_MatlabSparseMatrix.s_arrayNzmax = temp4; std::cout << -"m_MatlabSparseMatrix.s_arrayNzmax : \"" << m_MatlabSparseMatrix.s_arrayNzmax << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_dimensionTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_dimensionTag : \"" << m_MatlabSparseMatrix.s_dimensionTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_dimensionLength = temp4; - std::cout << "m_MatlabSparseMatrix.s_dimensionLength : \"" << m_MatlabSparseMatrix.s_dimensionLength << "\"" << -std::endl; in.read(reinterpret_cast(&temp4), 4); m_MatlabSparseMatrix.s_dimensionNbRow = temp4; std::cout << -"m_MatlabSparseMatrix.s_dimensionNbRow : \"" << m_MatlabSparseMatrix.s_dimensionNbRow << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_dimensionNbColumn = temp4; - std::cout << "m_MatlabSparseMatrix.s_dimensionNbColumn : \"" << m_MatlabSparseMatrix.s_dimensionNbColumn << "\"" << -std::endl; in.read(reinterpret_cast(&temp2), 2); m_MatlabSparseMatrix.s_nameLength = temp2; std::cout << -"m_MatlabSparseMatrix.s_nameLength : \"" << m_MatlabSparseMatrix.s_nameLength << "\"" << std::endl; - in.read(reinterpret_cast(&temp2), 2); - m_MatlabSparseMatrix.s_nameTag = temp2; - std::cout << "m_MatlabSparseMatrix.s_nameTag : \"" << m_MatlabSparseMatrix.s_nameTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp1),1); - m_MatlabSparseMatrix.s_nameChar = temp1; - std::cout << "m_MatlabSparseMatrix.s_nameChar : \"" << m_MatlabSparseMatrix.s_nameChar << "\"" << std::endl; - for (unsigned int i=0; i<3; ++i) { - in.read(reinterpret_cast(&temp1),1); - m_MatlabSparseMatrix.s_namePadding[i] = temp1; - } - std::cout << "m_MatlabSparseMatrix.s_namePadding : \"" << m_MatlabSparseMatrix.s_namePadding << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_rowIndexTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_rowIndexTag : \"" << m_MatlabSparseMatrix.s_rowIndexTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_rowIndexLength = temp4; - std::cout << "m_MatlabSparseMatrix.s_rowIndexLength : \"" << m_MatlabSparseMatrix.s_rowIndexLength << "\"" << -std::endl; m_MatlabSparseMatrix.s_rowIndex = new unsigned long int[m_MatlabSparseMatrix.s_arrayNzmax]; std::cout << -"m_MatlabSparseMatrix.s_rowIndex : \""; for (unsigned int i=0; i(&temp4),4); - m_MatlabSparseMatrix.s_rowIndex[i] = temp4; - std::cout << m_MatlabSparseMatrix.s_rowIndex[i] << " "; - } - std::cout << "\"" << std::endl; - if (m_MatlabSparseMatrix.s_arrayNzmax*4%8 != 0) { - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_rowIndexPadding = temp4; - std::cout << "m_MatlabSparseMatrix.s_rowIndexPadding : \"" << m_MatlabSparseMatrix.s_rowIndexPadding << "\"" << -std::endl; - } - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_columnIndexTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_columnIndexTag : \"" << m_MatlabSparseMatrix.s_columnIndexTag << "\"" << -std::endl; in.read(reinterpret_cast(&temp4), 4); m_MatlabSparseMatrix.s_columnIndexLength = temp4; std::cout << -"m_MatlabSparseMatrix.s_columnIndexLength : \"" << m_MatlabSparseMatrix.s_columnIndexLength << "\"" << std::endl; - m_MatlabSparseMatrix.s_columnIndex = new unsigned long int[m_MatlabSparseMatrix.s_dimensionNbColumn+1]; - std::cout << "m_MatlabSparseMatrix.s_columnIndex : \""; - for (unsigned int i=0; i(&temp4),4); - m_MatlabSparseMatrix.s_columnIndex[i] = temp4; - std::cout << m_MatlabSparseMatrix.s_columnIndex[i] << " "; - } - std::cout << "\"" << std::endl; - if ((m_MatlabSparseMatrix.s_dimensionNbColumn+1)*4%8 != 0) { - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_columnIndexPadding = temp4; - std::cout << "m_MatlabSparseMatrix.s_columnIndexPadding : \"" << m_MatlabSparseMatrix.s_columnIndexPadding << "\"" -<< std::endl; - } - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_valueTag = temp4; - std::cout << "m_MatlabSparseMatrix.s_valueTag : \"" << m_MatlabSparseMatrix.s_valueTag << "\"" << std::endl; - in.read(reinterpret_cast(&temp4), 4); - m_MatlabSparseMatrix.s_valueLength = temp4; - std::cout << "m_MatlabSparseMatrix.s_valueLength : \"" << m_MatlabSparseMatrix.s_valueLength << "\"" << std::endl; - m_MatlabSparseMatrix.s_value = new double[m_MatlabSparseMatrix.s_arrayNzmax]; - std::cout << "m_MatlabSparseMatrix.s_value : \""; - for (unsigned int i=0; i(&temp8),8); - m_MatlabSparseMatrix.s_value[i] = temp8; - std::cout << m_MatlabSparseMatrix.s_value[i] << " "; - } - std::cout << "\"" << std::endl; +} // namespace rtk - return in; -} -*/ +#endif From 9774c402fb862b171422865f37267789d118f222 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 6 Jan 2026 14:13:44 +0100 Subject: [PATCH 2/5] ENH: Move rtkStoreSparseMatrixSplatWeightMultiplication to include --- .../rtkprojectionmatrix.cxx | 67 +--------- ...oreSparseMatrixSplatWeightMultiplication.h | 119 ++++++++++++++++++ 2 files changed, 121 insertions(+), 65 deletions(-) create mode 100644 include/rtkStoreSparseMatrixSplatWeightMultiplication.h diff --git a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx index 5482cca97..14d86c27e 100644 --- a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx +++ b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx @@ -25,70 +25,7 @@ #include "rtkJosephBackProjectionImageFilter.h" #include "rtkConfiguration.h" #include "rtkMatlabSparseMatrix.h" - -namespace rtk -{ -namespace Functor -{ -template -class StoreSparseMatrixSplatWeightMultiplication -{ -public: - StoreSparseMatrixSplatWeightMultiplication() = default; - ~StoreSparseMatrixSplatWeightMultiplication() = default; - - bool - operator!=(const StoreSparseMatrixSplatWeightMultiplication &) const - { - return false; - } - bool - operator==(const StoreSparseMatrixSplatWeightMultiplication & other) const - { - return !(*this != other); - } - - inline void - operator()(const TInput & rayValue, - TOutput & output, - const double stepLengthInVoxel, - const double voxelSize, - const TCoordinateType weight) - { - // One row of the matrix is one ray, it should be thread safe - m_SystemMatrix.put( - &rayValue - m_ProjectionsBuffer, &output - m_VolumeBuffer, weight * voxelSize * stepLengthInVoxel); - } - vnl_sparse_matrix & - GetVnlSparseMatrix() - { - return m_SystemMatrix; - } - void - SetProjectionsBuffer(TInput * pb) - { - m_ProjectionsBuffer = pb; - } - void - SetVolumeBuffer(TOutput * vb) - { - m_VolumeBuffer = vb; - } - -private: - vnl_sparse_matrix m_SystemMatrix; - TInput * m_ProjectionsBuffer; - TOutput * m_VolumeBuffer; -}; -} // namespace Functor -} // namespace rtk - -std::ostream & -operator<<(std::ostream & out, rtk::MatlabSparseMatrix & matlabSparseMatrix) -{ - matlabSparseMatrix.Save(out); - return out; -} +#include "rtkStoreSparseMatrixSplatWeightMultiplication.h" int @@ -185,7 +122,7 @@ main(int argc, char * argv[]) } rtk::MatlabSparseMatrix matlabSparseMatrix(backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix(), backProjection->GetOutput()); - ofs << matlabSparseMatrix; + matlabSparseMatrix.Save(ofs); ofs.close(); return EXIT_SUCCESS; } diff --git a/include/rtkStoreSparseMatrixSplatWeightMultiplication.h b/include/rtkStoreSparseMatrixSplatWeightMultiplication.h new file mode 100644 index 000000000..012ded761 --- /dev/null +++ b/include/rtkStoreSparseMatrixSplatWeightMultiplication.h @@ -0,0 +1,119 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkStoreSparseMatrixSplatWeightMultiplication_h +#define rtkStoreSparseMatrixSplatWeightMultiplication_h + +#include +#include "rtkConfiguration.h" + +namespace rtk +{ +namespace Functor +{ +/** + * \class StoreSparseMatrixSplatWeightMultiplication + * + * \brief Functor to capture and store the back-projection system matrix. + * + * This functor is used with JosephBackProjectionImageFilter to capture + * the sparse system matrix entries during back-projection computation. + * Each matrix entry A[i,j] represents the contribution of volume voxel j + * to projection pixel i during the back-projection operation. + * + * The matrix entry is computed as: weight * voxelSize * stepLengthInVoxel + * + * \ingroup RTK + */ +template +class StoreSparseMatrixSplatWeightMultiplication +{ +public: + StoreSparseMatrixSplatWeightMultiplication() = default; + ~StoreSparseMatrixSplatWeightMultiplication() = default; + + bool + operator!=(const StoreSparseMatrixSplatWeightMultiplication &) const + { + return false; + } + + bool + operator==(const StoreSparseMatrixSplatWeightMultiplication & other) const + { + return !(*this != other); + } + + /** + * \brief Store matrix entry for this voxel-ray intersection. + * + * \param rayValue Pointer to projection pixel value + * \param output Pointer to volume voxel value + * \param stepLengthInVoxel Distance traveled through voxel + * \param voxelSize Physical size of voxel + * \param weight Interpolation weight from back-projection filter + */ + inline void + operator()(const TInput & rayValue, + TOutput & output, + const double stepLengthInVoxel, + const double voxelSize, + const TCoordinateType weight) + { + // One row of the matrix is one ray, it should be thread safe + m_SystemMatrix.put( + &rayValue - m_ProjectionsBuffer, &output - m_VolumeBuffer, weight * voxelSize * stepLengthInVoxel); + } + + /** + * \brief Get reference to the sparse matrix. + */ + vnl_sparse_matrix & + GetVnlSparseMatrix() + { + return m_SystemMatrix; + } + + /** + * \brief Set pointer to projection data buffer for index computation. + */ + void + SetProjectionsBuffer(TInput * pb) + { + m_ProjectionsBuffer = pb; + } + + /** + * \brief Set pointer to volume data buffer for index computation. + */ + void + SetVolumeBuffer(TOutput * vb) + { + m_VolumeBuffer = vb; + } + +private: + vnl_sparse_matrix m_SystemMatrix; + TInput * m_ProjectionsBuffer; + TOutput * m_VolumeBuffer; +}; + +} // namespace Functor +} // namespace rtk + +#endif From ddf031b557a13e8fc55df6e7a7b78a7aa1579580 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 6 Jan 2026 14:24:28 +0100 Subject: [PATCH 3/5] ENH : Add rtkMatlabSparseMatrix example and test --- .../rtkprojectionmatrix.cxx | 7 +- examples/MatlabSparseMatrix/CMakeLists.txt | 12 +++ .../MatlabSparseMatrixExample.cxx | 81 +++++++++++++++++++ examples/MatlabSparseMatrix/README.md | 26 ++++++ examples/README.md | 1 + include/rtkMatlabSparseMatrix.h | 55 +++++++++++-- include/rtkMatlabSparseMatrix.hxx | 48 ++++++++--- test/CMakeLists.txt | 3 + 8 files changed, 209 insertions(+), 24 deletions(-) create mode 100644 examples/MatlabSparseMatrix/CMakeLists.txt create mode 100644 examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx create mode 100644 examples/MatlabSparseMatrix/README.md diff --git a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx index 14d86c27e..c3c5aab6d 100644 --- a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx +++ b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx @@ -120,9 +120,10 @@ main(int argc, char * argv[]) std::cerr << "Failed to open " << args_info.output_arg << std::endl; return EXIT_FAILURE; } - rtk::MatlabSparseMatrix matlabSparseMatrix(backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix(), - backProjection->GetOutput()); - matlabSparseMatrix.Save(ofs); + auto matlabSparseMatrix = rtk::MatlabSparseMatrix::New(); + matlabSparseMatrix->SetMatrix(backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix()); + matlabSparseMatrix->SetOutput(backProjection->GetOutput()); + matlabSparseMatrix->Save(ofs); ofs.close(); return EXIT_SUCCESS; } diff --git a/examples/MatlabSparseMatrix/CMakeLists.txt b/examples/MatlabSparseMatrix/CMakeLists.txt new file mode 100644 index 000000000..824d3ed30 --- /dev/null +++ b/examples/MatlabSparseMatrix/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(MatlabSparseMatrixExample) + +# Find ITK with RTK +find_package(ITK REQUIRED COMPONENTS RTK) +include(${ITK_USE_FILE}) + +# Executable(s) +add_executable(MatlabSparseMatrixExample MatlabSparseMatrixExample.cxx) +target_link_libraries(MatlabSparseMatrixExample ${ITK_LIBRARIES}) diff --git a/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx b/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx new file mode 100644 index 000000000..74c46a5d3 --- /dev/null +++ b/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx @@ -0,0 +1,81 @@ +#include "rtkMatlabSparseMatrix.h" +#include "rtkThreeDCircularProjectionGeometry.h" +#include "rtkConstantImageSource.h" +#include "rtkJosephBackProjectionImageFilter.h" +#include "rtkStoreSparseMatrixSplatWeightMultiplication.h" +#include +#include +#include +#include + +int +main() +{ + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + using OutputImageType = itk::Image; + + constexpr unsigned int numberOfProjections = 10; + constexpr unsigned int sid = 600; // source to isocenter distance + constexpr unsigned int sdd = 1200; // source to detector distance + + // Set up the geometry + auto geometry = rtk::ThreeDCircularProjectionGeometry::New(); + for (unsigned int i = 0; i < numberOfProjections; ++i) + { + double angle = (360.0 * i) / numberOfProjections; + geometry->AddProjectionInRadians(0, 0, angle * itk::Math::pi / 180.0, sid, sdd, 0, 0); + } + + // Create projection images + auto projectionsSource = rtk::ConstantImageSource::New(); + projectionsSource->SetSize(itk::MakeSize(512, 1, numberOfProjections)); + projectionsSource->SetSpacing(itk::MakeVector(1.0, 1.0, 1.0)); + projectionsSource->SetOrigin(itk::MakePoint(-256.0, 0.0, 0.0)); + projectionsSource->SetConstant(1.0); + projectionsSource->Update(); + + // Create volume + auto volumeSource = rtk::ConstantImageSource::New(); + volumeSource->SetSize(itk::MakeSize(32, 32, 32)); + volumeSource->SetSpacing(itk::MakeVector(2.0, 2.0, 2.0)); + volumeSource->SetOrigin(itk::MakePoint(-32.0, -32.0, -32.0)); + volumeSource->SetConstant(0.0); + volumeSource->Update(); + + // Back-project with matrix capture + using BackProjectionType = rtk::JosephBackProjectionImageFilter< + OutputImageType, + OutputImageType, + rtk::Functor::InterpolationWeightMultiplicationBackProjection, + rtk::Functor::StoreSparseMatrixSplatWeightMultiplication>; + + auto backProjection = BackProjectionType::New(); + backProjection->SetInput(volumeSource->GetOutput()); + backProjection->SetInput(1, projectionsSource->GetOutput()); + backProjection->SetGeometry(geometry); + + // Initialize and configure matrix capture + backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix().resize( + projectionsSource->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels(), + volumeSource->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels()); + backProjection->GetSplatWeightMultiplication().SetProjectionsBuffer( + projectionsSource->GetOutput()->GetBufferPointer()); + backProjection->GetSplatWeightMultiplication().SetVolumeBuffer(volumeSource->GetOutput()->GetBufferPointer()); + + backProjection->Update(); + + // Export to Matlab format + vnl_sparse_matrix & systemMatrix = backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix(); + auto matlabMatrix = rtk::MatlabSparseMatrix::New(); + matlabMatrix->SetMatrix(systemMatrix); + matlabMatrix->SetOutput(volumeSource->GetOutput()); + + std::ofstream outputFile("backprojection_matrix.mat", std::ios::binary); + matlabMatrix->Save(outputFile); + + // Print matrix information + matlabMatrix->Print(); + + return 0; +} diff --git a/examples/MatlabSparseMatrix/README.md b/examples/MatlabSparseMatrix/README.md new file mode 100644 index 000000000..f8163bcec --- /dev/null +++ b/examples/MatlabSparseMatrix/README.md @@ -0,0 +1,26 @@ +# Matlab Sparse Matrix + +This example demonstrates how to compute the tomographic back-projection system matrix and export it in Matlab .mat format using RTK's [JosephBackProjectionImageFilter](../../documentation/docs/Projectors.md) + +## Using in Matlab + +Load and analyze the generated sparse matrix: + +```matlab +% Load the sparse matrix (stored as variable 'A') +load('backprojection_matrix.mat') + +% Visualize the sparsity pattern +spy(A) + +% Display matrix statistics +fprintf('Matrix size: %d x %d\n', size(A,1), size(A,2)) +fprintf('Non-zero elements: %d\n', nnz(A)) +fprintf('Sparsity: %.2f%%\n', 100 * (1 - nnz(A)/(size(A,1)*size(A,2)))) +``` + +## Code + +```{literalinclude} MatlabSparseMatrixExample.cxx +:language: c++ +``` diff --git a/examples/README.md b/examples/README.md index 48aeaf70a..7e6d167df 100644 --- a/examples/README.md +++ b/examples/README.md @@ -23,4 +23,5 @@ providing efficient implementations for high-performance applications. ./InlineReconstruction/README.md ./AddNoise/README.md ./GeometricPhantom/README.md +./MatlabSparseMatrix/README.md ``` diff --git a/include/rtkMatlabSparseMatrix.h b/include/rtkMatlabSparseMatrix.h index a27b08489..df32cf2f9 100644 --- a/include/rtkMatlabSparseMatrix.h +++ b/include/rtkMatlabSparseMatrix.h @@ -19,9 +19,17 @@ #ifndef rtkMatlabSparseMatrix_h #define rtkMatlabSparseMatrix_h +#include "RTKExport.h" #include "rtkConfiguration.h" + +#include +#include + #include +#include +#include + namespace rtk { @@ -36,11 +44,33 @@ namespace rtk * * \ingroup RTK */ -class MatlabSparseMatrix +template +class ITK_TEMPLATE_EXPORT MatlabSparseMatrix : public itk::Object { public: - template - MatlabSparseMatrix(const vnl_sparse_matrix & sparseMatrix, TOutputImage * output); + ITK_DISALLOW_COPY_AND_MOVE(MatlabSparseMatrix); + + using Self = MatlabSparseMatrix; + using Superclass = itk::Object; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + using OutputImageType = TOutputImage; + using MatrixType = vnl_sparse_matrix; + + itkOverrideGetNameOfClassMacro(MatlabSparseMatrix); + itkFactorylessNewMacro(Self); + + // Custom setter needed because vnl_sparse_matrix doesn't support operator<< for itkSetMacro + void + SetMatrix(const MatrixType & matrix) + { + m_Matrix = matrix; + this->Modified(); + } + itkGetMacro(Matrix, MatrixType); + + itkSetConstObjectMacro(Output, OutputImageType); + itkGetConstObjectMacro(Output, OutputImageType); void Save(std::ostream & out); @@ -48,6 +78,10 @@ class MatlabSparseMatrix void Print(); +protected: + MatlabSparseMatrix() = default; + ~MatlabSparseMatrix() override; + struct MatlabSparseMatrixStruct { unsigned char s_headerMatlab[116]; @@ -72,19 +106,24 @@ class MatlabSparseMatrix unsigned char s_namePadding[3]; unsigned long int s_rowIndexTag; unsigned long int s_rowIndexLength; - unsigned long int * s_rowIndex; + unsigned long int * s_rowIndex = nullptr; unsigned long int s_rowIndexPadding; unsigned long int s_columnIndexTag; unsigned long int s_columnIndexLength; - unsigned long int * s_columnIndex; + unsigned long int * s_columnIndex = nullptr; unsigned long int s_columnIndexPadding; unsigned long int s_valueTag; unsigned long int s_valueLength; - double * s_value; + double * s_value = nullptr; }; -protected: - MatlabSparseMatrixStruct m_MatlabSparseMatrix; +private: + void + BuildMatlabMatrix(); + + MatlabSparseMatrixStruct m_MatlabSparseMatrix; + MatrixType m_Matrix{}; + typename OutputImageType::ConstPointer m_Output; }; } // namespace rtk diff --git a/include/rtkMatlabSparseMatrix.hxx b/include/rtkMatlabSparseMatrix.hxx index 83c8c9ae6..d806a456d 100644 --- a/include/rtkMatlabSparseMatrix.hxx +++ b/include/rtkMatlabSparseMatrix.hxx @@ -21,6 +21,7 @@ #include #include +#include #include namespace rtk @@ -28,10 +29,27 @@ namespace rtk //==================================================================== template -MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseMatrix, TOutputImage * output) +MatlabSparseMatrix::~MatlabSparseMatrix() { - unsigned int nbColumn = output->GetLargestPossibleRegion().GetSize()[0] * - output->GetLargestPossibleRegion().GetSize()[2]; // it's not sparseMatrix.columns() + if (m_MatlabSparseMatrix.s_rowIndex) + delete[] m_MatlabSparseMatrix.s_rowIndex; + if (m_MatlabSparseMatrix.s_columnIndex) + delete[] m_MatlabSparseMatrix.s_columnIndex; + if (m_MatlabSparseMatrix.s_value) + delete[] m_MatlabSparseMatrix.s_value; +} + +template +void +MatlabSparseMatrix::BuildMatlabMatrix() +{ + if (!m_Output) + { + itkExceptionMacro("Output image not set"); + } + + unsigned int nbColumn = m_Output->GetLargestPossibleRegion().GetSize()[0] * + m_Output->GetLargestPossibleRegion().GetSize()[2]; // it's not m_Matrix.columns() // Take the number of non-zero elements: // Compute the column index @@ -39,21 +57,21 @@ MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseM unsigned int nonZeroElement(0); using sparseMatrixColumn = std::vector>; auto * columnsVector = new sparseMatrixColumn[nbColumn]; - sparseMatrix.reset(); - while (sparseMatrix.next()) + m_Matrix.reset(); + while (m_Matrix.next()) { - typename TOutputImage::IndexType idx = output->ComputeIndex(sparseMatrix.getcolumn()); + typename TOutputImage::IndexType idx = m_Output->ComputeIndex(m_Matrix.getcolumn()); if (idx[1] != 1) continue; - unsigned int indexColumn = idx[0] + idx[2] * output->GetLargestPossibleRegion().GetSize()[2]; + unsigned int indexColumn = idx[0] + idx[2] * m_Output->GetLargestPossibleRegion().GetSize()[2]; auto it = columnsVector[indexColumn].begin(); while (it != columnsVector[indexColumn].end()) { - if ((unsigned int)sparseMatrix.getrow() < it->first) + if ((unsigned int)m_Matrix.getrow() < it->first) break; ++it; } - columnsVector[indexColumn].insert(it, std::make_pair(sparseMatrix.getrow(), sparseMatrix.value())); + columnsVector[indexColumn].insert(it, std::make_pair(m_Matrix.getrow(), m_Matrix.value())); ++nonZeroElement; } @@ -76,7 +94,7 @@ MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseM m_MatlabSparseMatrix.s_arrayNzmax = nonZeroElement; m_MatlabSparseMatrix.s_dimensionTag = 5; m_MatlabSparseMatrix.s_dimensionLength = 8; - m_MatlabSparseMatrix.s_dimensionNbRow = sparseMatrix.rows(); + m_MatlabSparseMatrix.s_dimensionNbRow = m_Matrix.rows(); m_MatlabSparseMatrix.s_dimensionNbColumn = nbColumn; m_MatlabSparseMatrix.s_nameLength = 1; m_MatlabSparseMatrix.s_nameTag = 1; @@ -117,9 +135,12 @@ MatlabSparseMatrix::MatlabSparseMatrix(const vnl_sparse_matrix & sparseM delete[] columnsVector; } -inline void -MatlabSparseMatrix::Save(std::ostream & out) +template +void +MatlabSparseMatrix::Save(std::ostream & out) { + BuildMatlabMatrix(); + // Write data out.write(reinterpret_cast(&m_MatlabSparseMatrix.s_headerMatlab), 116); out.write(reinterpret_cast(&m_MatlabSparseMatrix.s_headerOffset), 8); @@ -159,8 +180,9 @@ MatlabSparseMatrix::Save(std::ostream & out) out.write(reinterpret_cast(&m_MatlabSparseMatrix.s_value[i]), 8); } +template inline void -MatlabSparseMatrix::Print() +MatlabSparseMatrix::Print() { std::cout << "m_MatlabSparseMatrix.s_headerMatlab : \"" << m_MatlabSparseMatrix.s_headerMatlab << "\"" << std::endl; std::cout << "m_MatlabSparseMatrix.s_headerOffset : \"" << m_MatlabSparseMatrix.s_headerOffset << "\"" << std::endl; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 107d7016a..4b3fca1e3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -376,6 +376,9 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkMatlabSparseMatrixExampleTest + ../examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx +) if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py) From 1789f213d2f0006f663ae4d099ac61aa3aff8f0a Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Mon, 12 Jan 2026 10:59:41 +0100 Subject: [PATCH 4/5] ENH: Wrap rtkprojectionmatrix and add python application and example --- .../rtkprojectionmatrix.py | 115 ++++++++++++++++++ .../MatlabSparseMatrixExample.py | 77 ++++++++++++ test/rtktestexamples.py | 4 + wrapping/rtkMatlabSparseMatrix.wrap | 5 + 4 files changed, 201 insertions(+) create mode 100644 applications/rtkprojectionmatrix/rtkprojectionmatrix.py create mode 100644 examples/MatlabSparseMatrix/MatlabSparseMatrixExample.py create mode 100644 wrapping/rtkMatlabSparseMatrix.wrap diff --git a/applications/rtkprojectionmatrix/rtkprojectionmatrix.py b/applications/rtkprojectionmatrix/rtkprojectionmatrix.py new file mode 100644 index 000000000..44e4eda7b --- /dev/null +++ b/applications/rtkprojectionmatrix/rtkprojectionmatrix.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python +import argparse +import sys +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description="Saves the sparse system matrix of rtk::JosephBackProjectionImageFilter in a file. " + "Only works in 2D." + ) + + # General options + parser.add_argument( + "--geometry", + "-g", + help="XML geometry file name", + type=str, + required=True, + ) + parser.add_argument( + "--output", + "-o", + help="Output file name for the sparse matrix", + type=str, + required=True, + ) + + # RTK specific groups + rtk.add_rtkinputprojections_group(parser) + rtk.add_rtk3Doutputimage_group(parser) + + return parser + + +def process(args_info: argparse.Namespace): + OutputPixelType = itk.F + Dimension = 3 + OutputImageType = itk.Image[OutputPixelType, Dimension] + + reader = rtk.ProjectionsReader[OutputImageType].New() + rtk.SetProjectionsReaderFromArgParse(reader, args_info) + + if args_info.verbose: + print("Reading projections...") + reader.Update() + + # Create back projection image filter + if reader.GetOutput().GetLargestPossibleRegion().GetSize()[1] != 1: + print( + "This tool has been designed for 2D, i.e., with one row in the sinogram only." + ) + sys.exit(1) + + direction = reader.GetOutput().GetDirection() + if abs(direction[0, 0]) != 1.0 or abs(direction[1, 1]) != 1.0: + print("Projections with non-diagonal Direction is not handled.") + sys.exit(1) + + if args_info.verbose: + print(f"Reading geometry information from {args_info.geometry}...") + geometry = rtk.ReadGeometry(args_info.geometry) + + constant_image_source = rtk.ConstantImageSource[OutputImageType].New() + rtk.SetConstantImageSourceFromArgParse(constant_image_source, args_info) + constant_image_source.Update() + + if constant_image_source.GetOutput().GetLargestPossibleRegion().GetSize()[1] != 3: + print( + "This tool has been designed for 2D with Joseph project. " + "Joseph requires at least 2 slices in the y direction for bilinear interpolation. " + "To have one slice exactly in front of the row, use 3 slices in the volume " + "with the central slice in front of the single projection row." + ) + sys.exit(1) + + direction = constant_image_source.GetOutput().GetDirection() + if not direction.GetVnlMatrix().is_identity(): + print("Volume with non-identity Direction is not handled.") + sys.exit(1) + + if reader.GetOutput().GetLargestPossibleRegion().GetSize()[2] != len( + geometry.GetGantryAngles() + ): + print("Number of projections in the geometry and in the stack do not match.") + sys.exit(1) + + if args_info.verbose: + print("Backprojecting volume and recording matrix values...") + + backProjection = rtk.JosephBackProjectionImageFilter[ + OutputImageType, + OutputImageType, + ].New() + backProjection.SetInput(constant_image_source.GetOutput()) + backProjection.SetInput(1, reader.GetOutput()) + backProjection.SetGeometry(geometry) + backProjection.Update() + + if args_info.verbose: + print("Writing matrix to disk...") + + matlab_matrix = rtk.MatlabSparseMatrix[OutputImageType].New() + matlab_matrix.Save(args_info.output) + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.py b/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.py new file mode 100644 index 000000000..8f5aa83c0 --- /dev/null +++ b/examples/MatlabSparseMatrix/MatlabSparseMatrixExample.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python + +import itk +import numpy as np + +# Define types +OutputPixelType = itk.F +Dimension = 3 +OutputImageType = itk.Image[OutputPixelType, Dimension] + +# Parameters +numberOfProjections = 10 +sid = 600 # source to isocenter distance +sdd = 1200 # source to detector distance + +# Set up the geometry +GeometryType = itk.ThreeDCircularProjectionGeometry +geometry = GeometryType.New() +for i in range(numberOfProjections): + angle = (360.0 * i) / numberOfProjections + geometry.AddProjectionInRadians(0, 0, angle * np.pi / 180.0, sid, sdd, 0, 0) + +# Create projection images +projectionsSource = itk.ConstantImageSource[OutputImageType].New() +projectionsSource.SetSize([512, 1, numberOfProjections]) +projectionsSource.SetSpacing([1.0, 1.0, 1.0]) +projectionsSource.SetOrigin([-256.0, 0.0, 0.0]) +projectionsSource.SetConstant(1.0) +projectionsSource.Update() + +# Create volume +volumeSource = itk.ConstantImageSource[OutputImageType].New() +volumeSource.SetSize([32, 32, 32]) +volumeSource.SetSpacing([2.0, 2.0, 2.0]) +volumeSource.SetOrigin([-32.0, -32.0, -32.0]) +volumeSource.SetConstant(0.0) +volumeSource.Update() + +# Back-project with matrix capture +BackProjectionType = itk.JosephBackProjectionImageFilter[ + OutputImageType, + OutputImageType, + itk.Functor.InterpolationWeightMultiplicationBackProjection[OutputPixelType, OutputPixelType], + itk.Functor.StoreSparseMatrixSplatWeightMultiplication[OutputPixelType, itk.D, OutputPixelType] +] +backProjection = BackProjectionType.New() +backProjection.SetInput(volumeSource.GetOutput()) +backProjection.SetInput(1, projectionsSource.GetOutput()) +backProjection.SetGeometry(geometry) + +# Initialize and configure matrix capture +backProjection.GetSplatWeightMultiplication().GetVnlSparseMatrix().resize( + projectionsSource.GetOutput().GetLargestPossibleRegion().GetNumberOfPixels(), + volumeSource.GetOutput().GetLargestPossibleRegion().GetNumberOfPixels() +) +backProjection.GetSplatWeightMultiplication().SetProjectionsBuffer( + projectionsSource.GetOutput().GetBufferPointer() +) +backProjection.GetSplatWeightMultiplication().SetVolumeBuffer( + volumeSource.GetOutput().GetBufferPointer() +) + +backProjection.Update() + +# Export to Matlab format +systemMatrix = backProjection.GetSplatWeightMultiplication().GetVnlSparseMatrix() +matlabMatrix = itk.MatlabSparseMatrix[OutputImageType].New() +matlabMatrix.SetMatrix(systemMatrix) +matlabMatrix.SetOutput(volumeSource.GetOutput()) + +with open("backprojection_matrix.mat", "wb") as outputFile: + matlabMatrix.Save(outputFile) + +# Print matrix information +matlabMatrix.Print() + +print("Matrix successfully saved to backprojection_matrix.mat") diff --git a/test/rtktestexamples.py b/test/rtktestexamples.py index e694646b9..98d73d9c4 100644 --- a/test/rtktestexamples.py +++ b/test/rtktestexamples.py @@ -68,3 +68,7 @@ def test_ConjugateGradient(tmp_path, thorax_file): run_example( tmp_path, "ConjugateGradient/ConjugateGradient.py", thorax_file, out_img ) + + +def test_MatlabSparseMatrixExample(tmp_path): + run_example(tmp_path, "MatlabSparseMatrix/MatlabSparseMatrixExample.py") diff --git a/wrapping/rtkMatlabSparseMatrix.wrap b/wrapping/rtkMatlabSparseMatrix.wrap new file mode 100644 index 000000000..cf457f081 --- /dev/null +++ b/wrapping/rtkMatlabSparseMatrix.wrap @@ -0,0 +1,5 @@ +itk_wrap_class("rtk::MatlabSparseMatrix" POINTER) + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("I${ITKM_${t}}3" "itk::Image<${ITKT_${t}}, 3>") + endforeach() +itk_end_wrap_class() From 33844d36259d9381f2d87d969d06e14f68ca97ec Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 20 Jan 2026 10:16:22 +0100 Subject: [PATCH 5/5] ENH: Add functors in rtkJosephBackProjectionImageFilter wrapper --- .../rtkJosephBackProjectionImageFilter.wrap | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/wrapping/rtkJosephBackProjectionImageFilter.wrap b/wrapping/rtkJosephBackProjectionImageFilter.wrap index d1d437fe2..df8618304 100644 --- a/wrapping/rtkJosephBackProjectionImageFilter.wrap +++ b/wrapping/rtkJosephBackProjectionImageFilter.wrap @@ -1,14 +1,32 @@ set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) +itk_wrap_include("rtkStoreSparseMatrixSplatWeightMultiplication.h") + itk_wrap_named_class("rtk::Functor::SplatWeightMultiplication" "rtkFunctorSplatWeightMultiplication") foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("${ITKM_${t}}D${ITKM_${t}}" "${ITKT_${t}}, double, ${ITKT_${t}}") endforeach() itk_end_wrap_class() + +itk_wrap_named_class("rtk::Functor::StoreSparseMatrixSplatWeightMultiplication" "rtkFunctorStoreSparseMatrixSplatWeightMultiplication") + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("${ITKM_${t}}D${ITKM_${t}}" "${ITKT_${t}}, double, ${ITKT_${t}}") + endforeach() +itk_end_wrap_class() + +itk_wrap_named_class("rtk::Functor::InterpolationWeightMultiplicationBackProjection" "rtkFunctorInterpolationWeightMultiplicationBackProjection") + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") + endforeach() +itk_end_wrap_class() set(WRAPPER_AUTO_INCLUDE_HEADERS ON) itk_wrap_class("rtk::JosephBackProjectionImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) + # Standard version itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3SWM${ITKM_${t}}D${ITKM_${t}}" "itk::Image<${ITKT_${t}}, 3>, itk::Image< ${ITKT_${t}}, 3>") + # With matrix capture functors + itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3IWMBP${ITKM_${t}}${ITKM_${t}}SSMSW${ITKM_${t}}D${ITKM_${t}}" + "itk::Image<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 3>, rtk::Functor::InterpolationWeightMultiplicationBackProjection<${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::StoreSparseMatrixSplatWeightMultiplication<${ITKT_${t}}, double, ${ITKT_${t}}>") endforeach() itk_end_wrap_class()