Spectralvectorimage

applications/CMakeLists.txt

@@ -81,7 +81,6 @@ add_subdirectory(rtkspectralsimplexdecomposition)
 add_subdirectory(rtkspectralrooster)
 add_subdirectory(rtkspectralforwardmodel)
 add_subdirectory(rtkmaskcollimation)
-add_subdirectory(rtkspectraldenoiseprojections)
 add_subdirectory(rtkprojectionmatrix)
 add_subdirectory(rtkvectorconjugategradient)
 

applications/rtkinputprojections_group.py

@@ -113,7 +113,6 @@ def add_rtkinputprojections_group(parser):
         "--component",
         help="Vector component to extract, for multi-material projections",
         type=int,
-        default=0,
     )
     rtkinputprojections_group.add_argument(
         "--radius",

applications/rtkprojections/rtkprojections.cxx

@@ -21,32 +21,74 @@
 #include "rtkMacro.h"
 
 #include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+#include <itkVectorImage.h>
+
+#include <iostream>
 
 int
 main(int argc, char * argv[])
 {
   GGO(rtkprojections, args_info);
 
   constexpr unsigned int Dimension = 3;
-  using OutputImageType = itk::Image<float, Dimension>;
+  bool                   inputIsVectorImage = false;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION({
+    const auto fileNames = rtk::GetProjectionsFileNamesFromGgo(args_info);
+    if (!args_info.component_given && !fileNames.empty())
+    {
+      auto imageIO =
+        itk::ImageIOFactory::CreateImageIO(fileNames.front().c_str(), itk::ImageIOFactory::IOFileModeEnum::ReadMode);
+      if (imageIO.IsNotNull())
+      {
+        imageIO->SetFileName(fileNames.front());
+        imageIO->ReadImageInformation();
+        inputIsVectorImage = imageIO->GetPixelType() == itk::IOPixelEnum::VECTOR;
+      }
+    }
+  })
+  if (inputIsVectorImage)
+  {
+    using OutputImageType = itk::VectorImage<float, Dimension>;
+    if (args_info.poisson_given || args_info.gaussian_given)
+    {
+      std::cerr << "Noise addition is not supported for vector projections" << std::endl;
+      return EXIT_FAILURE;
+    }
+
+    using ReaderType = rtk::ProjectionsReader<OutputImageType>;
+    auto reader = ReaderType::New();
+    rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkprojections>(reader, args_info);
+
+    auto writer = itk::ImageFileWriter<OutputImageType>::New();
+    writer->SetFileName(args_info.output_arg);
+    writer->SetInput(reader->GetOutput());
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->UpdateOutputInformation())
+    writer->SetNumberOfStreamDivisions(1 + reader->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() /
+                                             (1024 * 1024 * 4));
+
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
+  }
+  else
+  {
+    using OutputImageType = itk::Image<float, Dimension>;
 
-  // Projections reader
-  using ReaderType = rtk::ProjectionsReader<OutputImageType>;
-  auto reader = ReaderType::New();
-  rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkprojections>(reader, args_info);
+    using ReaderType = rtk::ProjectionsReader<OutputImageType>;
+    auto reader = ReaderType::New();
+    rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkprojections>(reader, args_info);
 
-  OutputImageType::Pointer output =
-    rtk::AddNoiseFromGgo<OutputImageType, args_info_rtkprojections>(reader->GetOutput(), args_info);
+    OutputImageType::Pointer output =
+      rtk::AddNoiseFromGgo<OutputImageType, args_info_rtkprojections>(reader->GetOutput(), args_info);
 
-  // Write
-  auto writer = itk::ImageFileWriter<OutputImageType>::New();
-  writer->SetFileName(args_info.output_arg);
-  writer->SetInput(output);
-  TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->UpdateOutputInformation())
-  writer->SetNumberOfStreamDivisions(1 + reader->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() /
-                                           (1024 * 1024 * 4));
+    auto writer = itk::ImageFileWriter<OutputImageType>::New();
+    writer->SetFileName(args_info.output_arg);
+    writer->SetInput(output);
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->UpdateOutputInformation())
+    writer->SetNumberOfStreamDivisions(1 + reader->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() /
+                                             (1024 * 1024 * 4));
 
-  TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
+  }
 
   return EXIT_SUCCESS;
 }

applications/rtkprojections/rtkprojections.py

@@ -22,24 +22,55 @@ def process(args_info: argparse.Namespace):
 
     OutputPixelType = itk.F
     Dimension = 3
-    OutputImageType = itk.Image[OutputPixelType, Dimension]
-
-    # Projections reader
-    reader = rtk.ProjectionsReader[OutputImageType].New()
-    rtk.SetProjectionsReaderFromArgParse(reader, args_info)
-    if args_info.verbose:
-        print(f"Reading projections...")
-
-    # Add noise
-    output = rtk.AddNoiseFromArgParse(reader.GetOutput(), args_info)
-
-    # Write
-    writer = itk.ImageFileWriter[OutputImageType].New()
-    writer.SetFileName(args_info.output)
-    writer.SetInput(output)
-    if args_info.verbose:
-        print(f"Writing output to: {args_info.output}")
-    writer.Update()
+    fileNames = rtk.GetProjectionsFileNamesFromArgParse(args_info)
+
+    inputIsVectorImage = False
+    if args_info.component is None and fileNames:
+        imageio = itk.ImageIOFactory.CreateImageIO(
+            fileNames[0], itk.CommonEnums.IOFileMode_ReadMode
+        )
+        if imageio is not None:
+            imageio.SetFileName(fileNames[0])
+            imageio.ReadImageInformation()
+            inputIsVectorImage = (
+                imageio.GetPixelType() == itk.CommonEnums.IOPixel_VECTOR
+            )
+
+    if inputIsVectorImage:
+        if args_info.poisson is not None or args_info.gaussian is not None:
+            raise RuntimeError("Noise addition is not supported for vector projections")
+
+        OutputImageType = itk.VectorImage[OutputPixelType, Dimension]
+        reader = rtk.ProjectionsReader[OutputImageType].New()
+        rtk.SetProjectionsReaderFromArgParse(reader, args_info)
+        if args_info.verbose:
+            print(f"Reading projections...")
+
+        writer = itk.ImageFileWriter[OutputImageType].New()
+        writer.SetFileName(args_info.output)
+        writer.SetInput(reader.GetOutput())
+        if args_info.verbose:
+            print(f"Writing output to: {args_info.output}")
+        writer.Update()
+    else:
+        OutputImageType = itk.Image[OutputPixelType, Dimension]
+
+        # Projections reader
+        reader = rtk.ProjectionsReader[OutputImageType].New()
+        rtk.SetProjectionsReaderFromArgParse(reader, args_info)
+        if args_info.verbose:
+            print(f"Reading projections...")
+
+        # Add noise
+        output = rtk.AddNoiseFromArgParse(reader.GetOutput(), args_info)
+
+        # Write
+        writer = itk.ImageFileWriter[OutputImageType].New()
+        writer.SetFileName(args_info.output)
+        writer.SetInput(output)
+        if args_info.verbose:
+            print(f"Writing output to: {args_info.output}")
+        writer.Update()
 
 
 def main(argv=None):

applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx

@@ -24,6 +24,24 @@
 
 #include <itkImageFileReader.h>
 #include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+
+namespace
+{
+itk::ImageIOBase::Pointer
+GetFileHeader(const std::string & filename)
+{
+  itk::ImageIOBase::Pointer reader =
+    itk::ImageIOFactory::CreateImageIO(filename.c_str(), itk::ImageIOFactory::IOFileModeEnum::ReadMode);
+  if (!reader)
+  {
+    itkGenericExceptionMacro(<< "Could not read " << filename);
+  }
+  reader->SetFileName(filename);
+  reader->ReadImageInformation();
+  return reader;
+}
+} // namespace
 
 int
 main(int argc, char * argv[])
@@ -35,6 +53,7 @@ main(int argc, char * argv[])
   using DecomposedProjectionType = itk::VectorImage<PixelValueType, Dimension>;
   using MeasuredProjectionsType = itk::VectorImage<PixelValueType, Dimension>;
   using IncidentSpectrumImageType = itk::Image<PixelValueType, Dimension>;
+  using VectorSpectrumImageType = itk::VectorImage<PixelValueType, Dimension - 1>;
   using DetectorResponseImageType = itk::Image<PixelValueType, Dimension - 1>;
   using MaterialAttenuationsImageType = itk::Image<PixelValueType, Dimension - 1>;
 
@@ -43,7 +62,21 @@ main(int argc, char * argv[])
   TRY_AND_EXIT_ON_ITK_EXCEPTION(decomposedProjection = itk::ReadImage<DecomposedProjectionType>(args_info.input_arg))
 
   IncidentSpectrumImageType::Pointer incidentSpectrum;
-  TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage<IncidentSpectrumImageType>(args_info.incident_arg))
+  VectorSpectrumImageType::Pointer   vectorIncidentSpectrum;
+  unsigned int                       MaximumEnergy = 0;
+  itk::ImageIOBase::Pointer          incidentHeader;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentHeader = GetFileHeader(args_info.incident_arg))
+  if (incidentHeader->GetNumberOfDimensions() == Dimension - 1 && incidentHeader->GetNumberOfComponents() > 1)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(vectorIncidentSpectrum =
+                                    itk::ReadImage<VectorSpectrumImageType>(args_info.incident_arg))
+    MaximumEnergy = vectorIncidentSpectrum->GetVectorLength();
+  }
+  else
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage<IncidentSpectrumImageType>(args_info.incident_arg))
+    MaximumEnergy = incidentSpectrum->GetLargestPossibleRegion().GetSize()[0];
+  }
 
   DetectorResponseImageType::Pointer detectorResponse;
   TRY_AND_EXIT_ON_ITK_EXCEPTION(detectorResponse = itk::ReadImage<DetectorResponseImageType>(args_info.detector_arg))
@@ -55,7 +88,6 @@ main(int argc, char * argv[])
   // Get parameters from the images
   const unsigned int NumberOfMaterials = materialAttenuations->GetLargestPossibleRegion().GetSize()[0];
   const unsigned int NumberOfSpectralBins = args_info.thresholds_given;
-  const unsigned int MaximumEnergy = incidentSpectrum->GetLargestPossibleRegion().GetSize()[0];
 
   // Generate a set of zero-filled photon count projections
   auto measuredProjections = MeasuredProjectionsType::New();
@@ -97,7 +129,10 @@ main(int argc, char * argv[])
                                                       IncidentSpectrumImageType>::New();
   forward->SetInputDecomposedProjections(decomposedProjection);
   forward->SetInputMeasuredProjections(measuredProjections);
-  forward->SetInputIncidentSpectrum(incidentSpectrum);
+  if (vectorIncidentSpectrum.IsNotNull())
+    forward->SetInputIncidentSpectrum(vectorIncidentSpectrum);
+  else
+    forward->SetInputIncidentSpectrum(incidentSpectrum);
   forward->SetDetectorResponse(detectorResponse);
   forward->SetMaterialAttenuations(materialAttenuations);
   forward->SetThresholds(thresholds);

applications/rtkspectralonestep/rtkspectralonestep.cxx

@@ -31,6 +31,7 @@
 #include <vector> // std::vector
 
 #include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
 
 itk::ImageIOBase::Pointer
 GetFileHeader(const std::string & filename)
@@ -69,13 +70,24 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
   using DetectorResponseType = itk::Image<dataType, 2>;
   using MaterialAttenuationsType = itk::Image<dataType, 2>;
 #endif
+  using VectorSpectrumType = itk::VectorImage<dataType, Dimension - 1>;
 
   // Instantiate and update the readers
   typename MeasuredProjectionsType::Pointer mea;
   TRY_AND_EXIT_ON_ITK_EXCEPTION(mea = itk::ReadImage<MeasuredProjectionsType>(args_info.spectral_arg))
 
   IncidentSpectrumType::Pointer incidentSpectrum;
-  TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage<IncidentSpectrumType>(args_info.incident_arg))
+  VectorSpectrumType::Pointer   vectorIncidentSpectrum;
+  itk::ImageIOBase::Pointer     incidentHeader;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentHeader = GetFileHeader(args_info.incident_arg))
+  if (incidentHeader->GetNumberOfDimensions() == Dimension - 1 && incidentHeader->GetNumberOfComponents() > 1)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(vectorIncidentSpectrum = itk::ReadImage<VectorSpectrumType>(args_info.incident_arg))
+  }
+  else
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage<IncidentSpectrumType>(args_info.incident_arg))
+  }
 
   DetectorResponseType::Pointer detectorResponse;
   TRY_AND_EXIT_ON_ITK_EXCEPTION(detectorResponse = itk::ReadImage<DetectorResponseType>(args_info.detector_arg))
@@ -188,7 +200,10 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
   SetForwardProjectionFromGgo(args_info, mechlemOneStep.GetPointer());
   SetBackProjectionFromGgo(args_info, mechlemOneStep.GetPointer());
   mechlemOneStep->SetInputMaterialVolumes(input);
-  mechlemOneStep->SetInputIncidentSpectrum(incidentSpectrum);
+  if (vectorIncidentSpectrum.IsNotNull())
+    mechlemOneStep->SetInputIncidentSpectrum(vectorIncidentSpectrum);
+  else
+    mechlemOneStep->SetInputIncidentSpectrum(incidentSpectrum);
   mechlemOneStep->SetBinnedDetectorResponse(drm);
   mechlemOneStep->SetMaterialAttenuations(materialAttenuationsMatrix);
   mechlemOneStep->SetNumberOfIterations(args_info.niterations_arg);