Matlab sparse matrix

applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx

@@ -25,70 +25,7 @@
 #include "rtkJosephBackProjectionImageFilter.h"
 #include "rtkConfiguration.h"
 #include "rtkMatlabSparseMatrix.h"
-
-namespace rtk
-{
-namespace Functor
-{
-template <class TInput, class TCoordinateType, class TOutput = TCoordinateType>
-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<double> &
-  GetVnlSparseMatrix()
-  {
-    return m_SystemMatrix;
-  }
-  void
-  SetProjectionsBuffer(TInput * pb)
-  {
-    m_ProjectionsBuffer = pb;
-  }
-  void
-  SetVolumeBuffer(TOutput * vb)
-  {
-    m_VolumeBuffer = vb;
-  }
-
-private:
-  vnl_sparse_matrix<double> 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
@@ -183,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());
-  ofs << matlabSparseMatrix;
+  auto matlabSparseMatrix = rtk::MatlabSparseMatrix<OutputImageType>::New();
+  matlabSparseMatrix->SetMatrix(backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix());
+  matlabSparseMatrix->SetOutput(backProjection->GetOutput());
+  matlabSparseMatrix->Save(ofs);
   ofs.close();
   return EXIT_SUCCESS;
 }

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()

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})

examples/MatlabSparseMatrix/MatlabSparseMatrixExample.cxx

@@ -0,0 +1,81 @@
+#include "rtkMatlabSparseMatrix.h"
+#include "rtkThreeDCircularProjectionGeometry.h"
+#include "rtkConstantImageSource.h"
+#include "rtkJosephBackProjectionImageFilter.h"
+#include "rtkStoreSparseMatrixSplatWeightMultiplication.h"
+#include <vnl/vnl_sparse_matrix.h>
+#include <fstream>
+#include <itkImage.h>
+#include <iostream>
+
+int
+main()
+{
+  using OutputPixelType = float;
+  constexpr unsigned int Dimension = 3;
+  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
+
+  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<OutputImageType>::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<OutputImageType>::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<OutputPixelType, OutputPixelType>,
+    rtk::Functor::StoreSparseMatrixSplatWeightMultiplication<OutputPixelType, double, OutputPixelType>>;
+
+  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<double> & systemMatrix = backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix();
+  auto matlabMatrix = rtk::MatlabSparseMatrix<OutputImageType>::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;
+}

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")