diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 72730bf95..05bf8983a 100644 --- a/applications/CMakeLists.txt +++ b/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) diff --git a/applications/rtkinputprojections_group.py b/applications/rtkinputprojections_group.py index ced898bdb..23da63dba 100644 --- a/applications/rtkinputprojections_group.py +++ b/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", diff --git a/applications/rtkprojections/rtkprojections.cxx b/applications/rtkprojections/rtkprojections.cxx index ce56d4e74..f7fb9dc55 100644 --- a/applications/rtkprojections/rtkprojections.cxx +++ b/applications/rtkprojections/rtkprojections.cxx @@ -21,6 +21,10 @@ #include "rtkMacro.h" #include +#include +#include + +#include int main(int argc, char * argv[]) @@ -28,25 +32,63 @@ main(int argc, char * argv[]) GGO(rtkprojections, args_info); constexpr unsigned int Dimension = 3; - using OutputImageType = itk::Image; + 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; + 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; + auto reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(reader, args_info); + + auto writer = itk::ImageFileWriter::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; - // Projections reader - using ReaderType = rtk::ProjectionsReader; - auto reader = ReaderType::New(); - rtk::SetProjectionsReaderFromGgo(reader, args_info); + using ReaderType = rtk::ProjectionsReader; + auto reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(reader, args_info); - OutputImageType::Pointer output = - rtk::AddNoiseFromGgo(reader->GetOutput(), args_info); + OutputImageType::Pointer output = + rtk::AddNoiseFromGgo(reader->GetOutput(), args_info); - // Write - auto writer = itk::ImageFileWriter::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::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; } diff --git a/applications/rtkprojections/rtkprojections.py b/applications/rtkprojections/rtkprojections.py index 208837434..f074d8af2 100644 --- a/applications/rtkprojections/rtkprojections.py +++ b/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): diff --git a/applications/rtkspectraldenoiseprojections/CMakeLists.txt b/applications/rtkspectraldenoiseprojections/CMakeLists.txt deleted file mode 100644 index c64473237..000000000 --- a/applications/rtkspectraldenoiseprojections/CMakeLists.txt +++ /dev/null @@ -1,31 +0,0 @@ -wrap_ggo(rtkspectraldenoiseprojections_GGO_C rtkspectraldenoiseprojections.ggo) -add_executable( - rtkspectraldenoiseprojections - rtkspectraldenoiseprojections.cxx - ${rtkspectraldenoiseprojections_GGO_C} -) -target_link_libraries(rtkspectraldenoiseprojections ${RTK_APPLICATION_TARGETS}) - -set_target_properties( - rtkspectraldenoiseprojections - PROPERTIES - RUNTIME_OUTPUT_DIRECTORY - "${CMAKE_BINARY_DIR}/bin" -) - -# Installation code -if(NOT RTK_INSTALL_NO_EXECUTABLES) - install( - TARGETS - rtkspectraldenoiseprojections - RUNTIME - DESTINATION ${RTK_INSTALL_RUNTIME_DIR} - COMPONENT Runtime - LIBRARY - DESTINATION ${RTK_INSTALL_LIB_DIR} - COMPONENT RuntimeLibraries - ARCHIVE - DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} - COMPONENT Development - ) -endif() diff --git a/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.cxx b/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.cxx deleted file mode 100644 index a1f990cd1..000000000 --- a/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.cxx +++ /dev/null @@ -1,57 +0,0 @@ -/*========================================================================= - * - * 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. - * - *=========================================================================*/ - -#include "rtkspectraldenoiseprojections_ggo.h" -#include "rtkConditionalMedianImageFilter.h" -#include "rtkGgoFunctions.h" -#include "rtkMacro.h" - -#include -#include - -int -main(int argc, char * argv[]) -{ - GGO(rtkspectraldenoiseprojections, args_info); - - constexpr unsigned int Dimension = 3; - using OutputImageType = itk::VectorImage; - - // Reader - OutputImageType::Pointer input; - TRY_AND_EXIT_ON_ITK_EXCEPTION(input = itk::ReadImage(args_info.input_arg)) - - // Remove aberrant pixels - using MedianType = rtk::ConditionalMedianImageFilter; - auto median = MedianType::New(); - median->SetThresholdMultiplier(args_info.multiplier_arg); - MedianType::MedianRadiusType radius; - if (args_info.radius_given) - { - radius.Fill(args_info.radius_arg[0]); - for (unsigned int i = 0; i < args_info.radius_given; i++) - radius[i] = args_info.radius_arg[i]; - } - median->SetRadius(radius); - median->SetInput(input); - - // Write - TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(median->GetOutput(), args_info.output_arg)) - - return EXIT_SUCCESS; -} diff --git a/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.ggo b/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.ggo deleted file mode 100644 index 9c7f1331d..000000000 --- a/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.ggo +++ /dev/null @@ -1,6 +0,0 @@ -purpose "Replaces aberrant pixels by the median in a small neighborhood around them. Pixels are aberrant if the difference between their value and the median is larger that threshold multiplier * the standard deviation in the neighborhood" - -option "input" i "Input file name" string yes -option "output" o "Output file name" string yes -option "multiplier" m "Threshold multiplier (actual threshold is obtained by multiplying by standard dev. of neighborhood)" double no default="1" -option "radius" r "Radius of neighborhood in each direction (actual radius is 2r+1)" int multiple no default="1" diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx index d5c68bffa..565ac2aa2 100644 --- a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx +++ b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx @@ -24,6 +24,24 @@ #include #include +#include + +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; using MeasuredProjectionsType = itk::VectorImage; using IncidentSpectrumImageType = itk::Image; + using VectorSpectrumImageType = itk::VectorImage; using DetectorResponseImageType = itk::Image; using MaterialAttenuationsImageType = itk::Image; @@ -43,7 +62,21 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(decomposedProjection = itk::ReadImage(args_info.input_arg)) IncidentSpectrumImageType::Pointer incidentSpectrum; - TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(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(args_info.incident_arg)) + MaximumEnergy = vectorIncidentSpectrum->GetVectorLength(); + } + else + { + TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(args_info.incident_arg)) + MaximumEnergy = incidentSpectrum->GetLargestPossibleRegion().GetSize()[0]; + } DetectorResponseImageType::Pointer detectorResponse; TRY_AND_EXIT_ON_ITK_EXCEPTION(detectorResponse = itk::ReadImage(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); diff --git a/applications/rtkspectralonestep/rtkspectralonestep.cxx b/applications/rtkspectralonestep/rtkspectralonestep.cxx index f15d3a780..8d191eb7d 100644 --- a/applications/rtkspectralonestep/rtkspectralonestep.cxx +++ b/applications/rtkspectralonestep/rtkspectralonestep.cxx @@ -31,6 +31,7 @@ #include // std::vector #include +#include itk::ImageIOBase::Pointer GetFileHeader(const std::string & filename) @@ -69,13 +70,24 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info) using DetectorResponseType = itk::Image; using MaterialAttenuationsType = itk::Image; #endif + using VectorSpectrumType = itk::VectorImage; // Instantiate and update the readers typename MeasuredProjectionsType::Pointer mea; TRY_AND_EXIT_ON_ITK_EXCEPTION(mea = itk::ReadImage(args_info.spectral_arg)) IncidentSpectrumType::Pointer incidentSpectrum; - TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(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(args_info.incident_arg)) + } + else + { + TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(args_info.incident_arg)) + } DetectorResponseType::Pointer detectorResponse; TRY_AND_EXIT_ON_ITK_EXCEPTION(detectorResponse = itk::ReadImage(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); diff --git a/applications/rtkspectralsimplexdecomposition/rtkspectralsimplexdecomposition.cxx b/applications/rtkspectralsimplexdecomposition/rtkspectralsimplexdecomposition.cxx index a9331eb9d..2a0607b2f 100644 --- a/applications/rtkspectralsimplexdecomposition/rtkspectralsimplexdecomposition.cxx +++ b/applications/rtkspectralsimplexdecomposition/rtkspectralsimplexdecomposition.cxx @@ -24,6 +24,24 @@ #include #include +#include + +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[]) @@ -38,6 +56,7 @@ main(int argc, char * argv[]) using SpectralProjectionsType = itk::VectorImage; using IncidentSpectrumImageType = itk::Image; + using VectorSpectrumImageType = itk::VectorImage; using DetectorResponseImageType = itk::Image; @@ -51,7 +70,21 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(spectralProjection = itk::ReadImage(args_info.spectral_arg)) IncidentSpectrumImageType::Pointer incidentSpectrum; - TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(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(args_info.incident_arg)) + MaximumEnergy = vectorIncidentSpectrum->GetVectorLength(); + } + else + { + TRY_AND_EXIT_ON_ITK_EXCEPTION(incidentSpectrum = itk::ReadImage(args_info.incident_arg)) + MaximumEnergy = incidentSpectrum->GetLargestPossibleRegion().GetSize()[0]; + } DetectorResponseImageType::Pointer detectorResponse; TRY_AND_EXIT_ON_ITK_EXCEPTION(detectorResponse = itk::ReadImage(args_info.detector_arg)) @@ -63,7 +96,6 @@ main(int argc, char * argv[]) // Get parameters from the images const unsigned int NumberOfMaterials = materialAttenuations->GetLargestPossibleRegion().GetSize()[0]; const unsigned int NumberOfSpectralBins = spectralProjection->GetVectorLength(); - const unsigned int MaximumEnergy = incidentSpectrum->GetLargestPossibleRegion().GetSize()[0]; // Read the thresholds on command line and check their number itk::VariableLengthVector thresholds; @@ -106,7 +138,10 @@ main(int argc, char * argv[]) simplex->SetInputDecomposedProjections(decomposedProjection); simplex->SetGuessInitialization(args_info.guess_flag); simplex->SetInputMeasuredProjections(spectralProjection); - simplex->SetInputIncidentSpectrum(incidentSpectrum); + if (vectorIncidentSpectrum.IsNotNull()) + simplex->SetInputIncidentSpectrum(vectorIncidentSpectrum); + else + simplex->SetInputIncidentSpectrum(incidentSpectrum); simplex->SetDetectorResponse(detectorResponse); simplex->SetMaterialAttenuations(materialAttenuations); simplex->SetThresholds(thresholds); diff --git a/include/rtkMechlemOneStepSpectralReconstructionFilter.h b/include/rtkMechlemOneStepSpectralReconstructionFilter.h index e5cbb2ee8..1b57bd937 100644 --- a/include/rtkMechlemOneStepSpectralReconstructionFilter.h +++ b/include/rtkMechlemOneStepSpectralReconstructionFilter.h @@ -27,16 +27,19 @@ #include "rtkReorderProjectionsImageFilter.h" #include "rtkSeparableQuadraticSurrogateRegularizationImageFilter.h" #include "rtkThreeDCircularProjectionGeometry.h" +#include "rtkVectorImageToImageFilter.h" #include "rtkWeidingerForwardModelImageFilter.h" #include #include #include +#include #include #ifdef RTK_USE_CUDA # include "rtkCudaWeidingerForwardModelImageFilter.h" +# include #endif namespace rtk @@ -185,6 +188,25 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter using SingleComponentImageType = typename TOutputImage::template RebindImageType; + // Vector incident spectra are flattened with CPU iterators first. When the reconstruction uses CUDA, the transfer to + // TIncidentSpectrum must happen after this CPU-only conversion step. + using CPUIncidentSpectrumType = itk::Image; + using VectorSpectrumImageType = itk::VectorImage; + using FlattenVectorSpectrumFilterType = + rtk::VectorImageToImageFilter; + using PermuteSpectrumFilterType = itk::PermuteAxesImageFilter; +#ifdef RTK_USE_CUDA + // The VectorImage overload first converts the spectrum to a scalar CPU image. When the reconstruction uses CUDA, that + // CPU image must then be adapted to itk::CudaImage; a plain CastImageFilter would leave the CUDA data manager + // uninitialized and CudaWeidingerForwardModelImageFilter would read an invalid GPU buffer. + using CastIncidentSpectrumFilterType = + std::conditional_t, + itk::CastImageFilter, + itk::CudaImageFromImageFilter>; +#else + using CastIncidentSpectrumFilterType = itk::CastImageFilter; +#endif + #if !defined(ITK_WRAPPING_PARSER) # ifdef RTK_USE_CUDA using WeidingerForwardModelType = typename std::conditional_t< @@ -264,6 +286,8 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter SetInputMeasuredProjections(const VectorImageType * measuredProjections); void SetInputIncidentSpectrum(const TIncidentSpectrum * incidentSpectrum); + void + SetInputIncidentSpectrum(const VectorSpectrumImageType * incidentSpectrum); #ifndef ITK_FUTURE_LEGACY_REMOVE void SetInputPhotonCounts(const TMeasuredProjections * measuredProjections); @@ -332,6 +356,9 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter typename MultiplyGradientFilterType::Pointer m_MultiplyGradientToBeBackprojectedFilter; typename ReorderMeasuredProjectionsFilterType::Pointer m_ReorderMeasuredProjectionsFilter; typename ReorderProjectionsWeightsFilterType::Pointer m_ReorderProjectionsWeightsFilter; + typename FlattenVectorSpectrumFilterType::Pointer m_FlattenVectorSpectrumFilter; + typename PermuteSpectrumFilterType::Pointer m_PermuteSpectrumFilter; + typename CastIncidentSpectrumFilterType::Pointer m_CastIncidentSpectrumFilter; #endif /** The inputs of this filter have the same type but not the same meaning diff --git a/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx b/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx index 810349c16..aa58781a1 100644 --- a/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx +++ b/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx @@ -61,6 +61,9 @@ MechlemOneStepSpectralReconstructionFilterSetConstant(itk::NumericTraits::ZeroValue()); @@ -119,6 +122,22 @@ MechlemOneStepSpectralReconstructionFilterSetNthInput(2, const_cast(incidentSpectrum)); } +template +void +MechlemOneStepSpectralReconstructionFilter:: + SetInputIncidentSpectrum(const VectorSpectrumImageType * incidentSpectrum) +{ + m_FlattenVectorSpectrumFilter->SetInput(incidentSpectrum); + m_PermuteSpectrumFilter->SetInput(m_FlattenVectorSpectrumFilter->GetOutput()); + typename PermuteSpectrumFilterType::PermuteOrderArrayType order; + order[0] = 2; + order[1] = 0; + order[2] = 1; + m_PermuteSpectrumFilter->SetOrder(order); + m_CastIncidentSpectrumFilter->SetInput(m_PermuteSpectrumFilter->GetOutput()); + this->SetInputIncidentSpectrum(m_CastIncidentSpectrumFilter->GetOutput()); +} + #ifndef ITK_FUTURE_LEGACY_REMOVE template void diff --git a/include/rtkProjectionsReader.h b/include/rtkProjectionsReader.h index 0a4a12f22..b7f485886 100644 --- a/include/rtkProjectionsReader.h +++ b/include/rtkProjectionsReader.h @@ -313,8 +313,8 @@ class ITK_TEMPLATE_EXPORT ProjectionsReader : public itk::ImageSource::Pointer m_RawCastFilter; /** Pointers for post-processing filters that are created only when required. */ - typename WaterPrecorrectionType::Pointer m_WaterPrecorrectionFilter; - typename StreamingType::Pointer m_StreamingFilter; + itk::ProcessObject::Pointer m_WaterPrecorrectionFilter; + typename StreamingType::Pointer m_StreamingFilter; /** Image IO object which is stored to create the pipe only when required */ itk::ImageIOBase::Pointer m_ImageIO{ nullptr }; diff --git a/include/rtkProjectionsReader.hxx b/include/rtkProjectionsReader.hxx index 84909c27b..e164cfb93 100644 --- a/include/rtkProjectionsReader.hxx +++ b/include/rtkProjectionsReader.hxx @@ -28,7 +28,10 @@ #include #include #include +#include #include +#include +#include // RTK #include "rtkIOFactories.h" @@ -36,6 +39,8 @@ #include "rtkLUTbasedVariableI0RawToAttenuationImageFilter.h" #include "rtkConditionalMedianImageFilter.h" +#include + // Varian Obi includes #include "rtkHndImageIOFactory.h" #include "rtkXimImageIOFactory.h" @@ -108,7 +113,11 @@ template ProjectionsReader::ProjectionsReader() { // Filters common to all input types and that do not depend on the input image type. - m_WaterPrecorrectionFilter = WaterPrecorrectionType::New(); + if constexpr (std::is_same_v>) + m_WaterPrecorrectionFilter = nullptr; + else + m_WaterPrecorrectionFilter = WaterPrecorrectionType::New(); m_StreamingFilter = StreamingType::New(); // Default values of parameters @@ -182,9 +191,46 @@ ProjectionsReader::GenerateOutputInformation() m_RawCastFilter = nullptr; // Start creation - if ((!strcmp(imageIO->GetNameOfClass(), "EdfImageIO") && - imageIO->GetComponentType() == itk::ImageIOBase::IOComponentEnum::USHORT) || - !strcmp(imageIO->GetNameOfClass(), "XRadImageIO")) + if constexpr (std::is_same_v>) + { + if (imageIO->GetPixelType() != itk::IOPixelEnum::VECTOR) + { + itkGenericExceptionMacro(<< "Vector image output is only supported for vector pixel input files"); + } + + if (m_FileNames.size() == 1) + { + using ReaderType = itk::ImageFileReader; + auto reader = ReaderType::New(); + m_RawDataReader = reader; + } + else + { + using ReaderType = itk::ImageSeriesReader; + auto reader = ReaderType::New(); + m_RawDataReader = reader; + } + + using ChangeInfoType = itk::ChangeInformationImageFilter; + auto cif = ChangeInfoType::New(); + m_ChangeInformationFilter = cif; + + using CropType = itk::CropImageFilter; + auto crop = CropType::New(); + m_CropFilter = crop; + + using ConditionalMedianType = rtk::ConditionalMedianImageFilter; + auto cond = ConditionalMedianType::New(); + m_ConditionalMedianFilter = cond; + + using BinType = itk::BinShrinkImageFilter; + auto bin = BinType::New(); + m_BinningFilter = bin; + } + else if ((!strcmp(imageIO->GetNameOfClass(), "EdfImageIO") && + imageIO->GetComponentType() == itk::ImageIOBase::IOComponentEnum::USHORT) || + !strcmp(imageIO->GetNameOfClass(), "XRadImageIO")) { using InputPixelType = unsigned short; using InputImageType = itk::Image; @@ -362,7 +408,7 @@ ProjectionsReader::GenerateOutputInformation() else { // If ImageIO has vector pixels, extract one component from it - if (!strcmp(imageIO->GetPixelTypeAsString(imageIO->GetPixelType()).c_str(), "vector")) + if (imageIO->GetPixelType() == itk::IOPixelEnum::VECTOR) { SET_INPUT_IMAGE_VECTOR_TYPE(float, 1) SET_INPUT_IMAGE_VECTOR_TYPE(float, 2) @@ -412,7 +458,10 @@ ProjectionsReader::GenerateOutputInformation() } // Parameter propagation - if (imageIO->GetComponentType() == itk::ImageIOBase::IOComponentEnum::USHORT) + if constexpr (std::is_same_v>) + PropagateParametersToMiniPipeline(); + else if (imageIO->GetComponentType() == itk::ImageIOBase::IOComponentEnum::USHORT) PropagateParametersToMiniPipeline>(); else if (!strcmp(imageIO->GetNameOfClass(), "HndImageIO") || !strcmp(imageIO->GetNameOfClass(), "XimImageIO")) PropagateParametersToMiniPipeline>(); @@ -426,6 +475,11 @@ ProjectionsReader::GenerateOutputInformation() output->SetSpacing(m_StreamingFilter->GetOutput()->GetSpacing()); output->SetDirection(m_StreamingFilter->GetOutput()->GetDirection()); output->SetLargestPossibleRegion(m_StreamingFilter->GetOutput()->GetLargestPossibleRegion()); + if constexpr (std::is_same_v>) + { + output->SetVectorLength(m_StreamingFilter->GetOutput()->GetVectorLength()); + } } //-------------------------------------------------------------------- @@ -468,12 +522,36 @@ ProjectionsReader::PropagateParametersToMiniPipeline() else // Regular case { // Raw - using RawType = typename itk::ImageSeriesReader; - auto * raw = dynamic_cast(m_RawDataReader.GetPointer()); - assert(raw != nullptr); - raw->SetFileNames(this->GetFileNames()); - raw->SetImageIO(m_ImageIO); - nextInput = raw->GetOutput(); + if constexpr (std::is_same_v>) + { + using RawFileType = typename itk::ImageFileReader; + auto * rawFile = dynamic_cast(m_RawDataReader.GetPointer()); + if (rawFile != nullptr) + { + rawFile->SetFileName(this->GetFileNames().front()); + rawFile->SetImageIO(m_ImageIO); + nextInput = rawFile->GetOutput(); + } + else + { + using RawSeriesType = typename itk::ImageSeriesReader; + auto * rawSeries = dynamic_cast(m_RawDataReader.GetPointer()); + assert(rawSeries != nullptr); + rawSeries->SetFileNames(this->GetFileNames()); + rawSeries->SetImageIO(m_ImageIO); + nextInput = rawSeries->GetOutput(); + } + } + else + { + using RawType = typename itk::ImageSeriesReader; + auto * raw = dynamic_cast(m_RawDataReader.GetPointer()); + assert(raw != nullptr); + raw->SetFileNames(this->GetFileNames()); + raw->SetImageIO(m_ImageIO); + nextInput = raw->GetOutput(); + } // Image information OutputImageSpacingType defaultSpacing; @@ -584,90 +662,109 @@ ProjectionsReader::PropagateParametersToMiniPipeline() } } - // Boellaard scatter correction - if (m_NonNegativityConstraintThreshold != itk::NumericTraits::NonpositiveMin() || - m_ScatterToPrimaryRatio != 0.) + if constexpr (std::is_same_v>) { - if (m_ScatterFilter.GetPointer() == nullptr) - { - itkGenericExceptionMacro(<< "Can not use Boellaard scatter correction with this input (not implemented)"); - } - else + if (m_NonNegativityConstraintThreshold != itk::NumericTraits::NonpositiveMin() || + m_ScatterToPrimaryRatio != 0. || m_I0 != itk::NumericTraits::NonpositiveMin() || + !m_WaterPrecorrectionCoefficients.empty()) { - using ScatterFilterType = rtk::BoellaardScatterCorrectionImageFilter; - auto * scatter = dynamic_cast(m_ScatterFilter.GetPointer()); - assert(scatter != nullptr); - scatter->SetAirThreshold(m_AirThreshold); - scatter->SetScatterToPrimaryRatio(m_ScatterToPrimaryRatio); - if (m_NonNegativityConstraintThreshold != itk::NumericTraits::NonpositiveMin()) - scatter->SetNonNegativityConstraintThreshold(m_NonNegativityConstraintThreshold); - scatter->SetInput(nextInput); - nextInput = scatter->GetOutput(); + itkGenericExceptionMacro( + << "Raw projection conversion, scatter correction, and water precorrection are not supported for vector " + "projections"); } - } - // LUTbasedVariableI0RawToAttenuationImageFilter - if (m_I0 != itk::NumericTraits::NonpositiveMin()) + output = dynamic_cast(nextInput); + assert(output != nullptr); + } + else { - if (m_RawToAttenuationFilter.GetPointer() == nullptr) + // Boellaard scatter correction + if (m_NonNegativityConstraintThreshold != itk::NumericTraits::NonpositiveMin() || + m_ScatterToPrimaryRatio != 0.) { - itkGenericExceptionMacro( - << "Can not use I0 in LUTbasedVariableI0RawToAttenuationImageFilter with this input (not implemented)"); + if (m_ScatterFilter.GetPointer() == nullptr) + { + itkGenericExceptionMacro(<< "Can not use Boellaard scatter correction with this input (not implemented)"); + } + else + { + using ScatterFilterType = rtk::BoellaardScatterCorrectionImageFilter; + auto * scatter = dynamic_cast(m_ScatterFilter.GetPointer()); + assert(scatter != nullptr); + scatter->SetAirThreshold(m_AirThreshold); + scatter->SetScatterToPrimaryRatio(m_ScatterToPrimaryRatio); + if (m_NonNegativityConstraintThreshold != itk::NumericTraits::NonpositiveMin()) + scatter->SetNonNegativityConstraintThreshold(m_NonNegativityConstraintThreshold); + scatter->SetInput(nextInput); + nextInput = scatter->GetOutput(); + } } - else + + // LUTbasedVariableI0RawToAttenuationImageFilter + if (m_I0 != itk::NumericTraits::NonpositiveMin()) { - nextInputBase = dynamic_cast *>(nextInput); - assert(nextInputBase != nullptr); - PropagateI0(&nextInputBase); - nextInput = dynamic_cast(nextInputBase); - assert(nextInput != nullptr); + if (m_RawToAttenuationFilter.GetPointer() == nullptr) + { + itkGenericExceptionMacro( + << "Can not use I0 in LUTbasedVariableI0RawToAttenuationImageFilter with this input (not implemented)"); + } + else + { + nextInputBase = dynamic_cast *>(nextInput); + assert(nextInputBase != nullptr); + PropagateI0(&nextInputBase); + nextInput = dynamic_cast(nextInputBase); + assert(nextInput != nullptr); + } } - } - // Raw to attenuation or cast filter, change of type - if (m_RawToAttenuationFilter.GetPointer() != nullptr) - { - // Check if Ora pointer - using OraRawType = rtk::OraLookupTableImageFilter; - auto * oraraw = dynamic_cast(m_RawToAttenuationFilter.GetPointer()); - if (oraraw != nullptr) + // Raw to attenuation or cast filter, change of type + if (m_RawToAttenuationFilter.GetPointer() != nullptr) { - oraraw->SetComputeLineIntegral(m_ComputeLineIntegral); - oraraw->SetFileNames(m_FileNames); - } + // Check if Ora pointer + using OraRawType = rtk::OraLookupTableImageFilter; + auto * oraraw = dynamic_cast(m_RawToAttenuationFilter.GetPointer()); + if (oraraw != nullptr) + { + oraraw->SetComputeLineIntegral(m_ComputeLineIntegral); + oraraw->SetFileNames(m_FileNames); + } - // Cast or convert to line integral depending on m_ComputeLineIntegral - using IToIFilterType = itk::ImageToImageFilter; - IToIFilterType * itoi = nullptr; - if (m_ComputeLineIntegral || oraraw != nullptr) - itoi = dynamic_cast(m_RawToAttenuationFilter.GetPointer()); + // Cast or convert to line integral depending on m_ComputeLineIntegral + using IToIFilterType = itk::ImageToImageFilter; + IToIFilterType * itoi = nullptr; + if (m_ComputeLineIntegral || oraraw != nullptr) + itoi = dynamic_cast(m_RawToAttenuationFilter.GetPointer()); + else + itoi = dynamic_cast(m_RawCastFilter.GetPointer()); + assert(itoi != nullptr); + itoi->SetInput(nextInput); + output = itoi->GetOutput(); + + // Release output data of m_RawDataReader if conversion occurs + itoi->ReleaseDataFlagOn(); + } else - itoi = dynamic_cast(m_RawCastFilter.GetPointer()); - assert(itoi != nullptr); - itoi->SetInput(nextInput); - output = itoi->GetOutput(); - - // Release output data of m_RawDataReader if conversion occurs - itoi->ReleaseDataFlagOn(); - } - else - { - output = dynamic_cast(nextInput); - assert(output != nullptr); - } + { + output = dynamic_cast(nextInput); + assert(output != nullptr); + } - // ESRF raw to attenuation converter also needs the filenames - auto * edf = dynamic_cast *>( - m_RawToAttenuationFilter.GetPointer()); - if (edf) - edf->SetFileNames(this->GetFileNames()); + // ESRF raw to attenuation converter also needs the filenames + auto * edf = dynamic_cast *>( + m_RawToAttenuationFilter.GetPointer()); + if (edf) + edf->SetFileNames(this->GetFileNames()); - // Water coefficients - if (!m_WaterPrecorrectionCoefficients.empty()) - { - m_WaterPrecorrectionFilter->SetCoefficients(m_WaterPrecorrectionCoefficients); - m_WaterPrecorrectionFilter->SetInput(output); - output = m_WaterPrecorrectionFilter->GetOutput(); + // Water coefficients + if (!m_WaterPrecorrectionCoefficients.empty()) + { + auto * wpc = static_cast(m_WaterPrecorrectionFilter.GetPointer()); + wpc->SetCoefficients(m_WaterPrecorrectionCoefficients); + wpc->SetInput(output); + output = wpc->GetOutput(); + } } } // Streaming image filter diff --git a/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.h b/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.h index a49c4a2f4..69ede1dcf 100644 --- a/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.h +++ b/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.h @@ -69,12 +69,10 @@ class ITK_TEMPLATE_EXPORT SimplexSpectralProjectionsDecompositionImageFilter using DecomposedProjectionsDataType = typename DecomposedProjectionsType::PixelType::ValueType; using MeasuredProjectionsDataType = typename MeasuredProjectionsType::PixelType::ValueType; -#ifndef ITK_FUTURE_LEGACY_REMOVE /** Additional types to overload SetInputIncidentSpectrum */ using VectorSpectrumImageType = itk::VectorImage; using FlattenVectorFilterType = rtk::VectorImageToImageFilter; using PermuteFilterType = itk::PermuteAxesImageFilter; -#endif /** Standard New method. */ itkNewMacro(Self); @@ -122,10 +120,8 @@ class ITK_TEMPLATE_EXPORT SimplexSpectralProjectionsDecompositionImageFilter /** Set/Get the incident spectrum input images */ void SetInputIncidentSpectrum(const IncidentSpectrumImageType * IncidentSpectrum); -#ifndef ITK_FUTURE_LEGACY_REMOVE void SetInputIncidentSpectrum(const VectorSpectrumImageType * IncidentSpectrum); -#endif typename IncidentSpectrumImageType::ConstPointer GetInputIncidentSpectrum(); @@ -209,11 +205,9 @@ class ITK_TEMPLATE_EXPORT SimplexSpectralProjectionsDecompositionImageFilter typename itk::ImageSource::Pointer m_CastDecomposedProjectionsPointer; typename itk::ImageSource::Pointer m_CastMeasuredProjectionsPointer; -#ifndef ITK_FUTURE_LEGACY_REMOVE /** Filters required for the overload of SetInputIncidentSpectrum */ typename FlattenVectorFilterType::Pointer m_FlattenFilter; typename PermuteFilterType::Pointer m_PermuteFilter; -#endif }; // end of class } // end namespace rtk diff --git a/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.hxx b/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.hxx index 45a33b2d1..397ececfa 100644 --- a/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.hxx +++ b/include/rtkSimplexSpectralProjectionsDecompositionImageFilter.hxx @@ -66,11 +66,9 @@ SimplexSpectralProjectionsDecompositionImageFilter< m_LogTransformEachBin = false; m_GuessInitialization = false; -#ifndef ITK_FUTURE_LEGACY_REMOVE // Instantiate the filters required in the overload of SetInputIncidentSpectrum m_FlattenFilter = FlattenVectorFilterType::New(); m_PermuteFilter = PermuteFilterType::New(); -#endif } template SetInput("IncidentSpectrum", const_cast(IncidentSpectrum)); } -#ifndef ITK_FUTURE_LEGACY_REMOVE template m_PermuteFilter->SetOrder(order); this->SetInputIncidentSpectrum(m_PermuteFilter->GetOutput()); } -#endif template ; using FlattenVectorFilterType = rtk::VectorImageToImageFilter; using PermuteFilterType = itk::PermuteAxesImageFilter; -#endif /** Standard New method. */ itkNewMacro(Self); @@ -85,10 +83,8 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter /** Set/Get the incident spectrum input images */ void SetInputIncidentSpectrum(const IncidentSpectrumImageType * IncidentSpectrum); -#ifndef ITK_FUTURE_LEGACY_REMOVE void SetInputIncidentSpectrum(const VectorSpectrumImageType * IncidentSpectrum); -#endif typename IncidentSpectrumImageType::ConstPointer GetInputIncidentSpectrum(); @@ -200,11 +196,9 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter typename itk::ImageSource::Pointer m_CastDecomposedProjectionsPointer; typename itk::ImageSource::Pointer m_CastMeasuredProjectionsPointer; -#ifndef ITK_FUTURE_LEGACY_REMOVE /** Filters required for the overload of SetInputIncidentSpectrum */ typename FlattenVectorFilterType::Pointer m_FlattenFilter; typename PermuteFilterType::Pointer m_PermuteFilter; -#endif }; // end of class // Function to bin a detector response matrix according to given energy thresholds diff --git a/include/rtkSpectralForwardModelImageFilter.hxx b/include/rtkSpectralForwardModelImageFilter.hxx index 19bf40a80..cda93d206 100644 --- a/include/rtkSpectralForwardModelImageFilter.hxx +++ b/include/rtkSpectralForwardModelImageFilter.hxx @@ -52,11 +52,9 @@ SpectralForwardModelImageFilterSetNthOutput(2, this->MakeOutput(2)); -#ifndef ITK_FUTURE_LEGACY_REMOVE // Instantiate the filters required in the overload of SetInputIncidentSpectrum m_FlattenFilter = FlattenVectorFilterType::New(); m_PermuteFilter = PermuteFilterType::New(); -#endif } template SetInput("IncidentSpectrum", const_cast(IncidentSpectrum)); } -#ifndef ITK_FUTURE_LEGACY_REMOVE template m_PermuteFilter->SetOrder(order); this->SetInputIncidentSpectrum(m_PermuteFilter->GetOutput()); } -#endif template +#include + +#include + +/** + * \file rtkprojectionsreadervectorimagetest.cxx + * + * \brief Test rtk::ProjectionsReader with an itk::VectorImage output. + */ + +int +rtkprojectionsreadervectorimagetest(int argc, char * argv[]) +{ + if (argc < 2) + { + std::cerr << "Usage: " << std::endl; + std::cerr << argv[0] << " vectorProjectionImage " << std::endl; + return EXIT_FAILURE; + } + + using VectorImageType = itk::VectorImage; + + auto referenceReader = itk::ImageFileReader::New(); + referenceReader->SetFileName(argv[1]); + TRY_AND_EXIT_ON_ITK_EXCEPTION(referenceReader->Update()); + + auto projectionsReader = rtk::ProjectionsReader::New(); + projectionsReader->SetFileNames(std::vector{ argv[1] }); + TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update()); + + CheckVariableLengthVectorImageQuality( + projectionsReader->GetOutput(), referenceReader->GetOutput(), 1e-5, 100, 2000.0); + + using MedianType = rtk::ConditionalMedianImageFilter; + auto directMedian = MedianType::New(); + MedianType::MedianRadiusType radius; + radius[0] = 1; + radius[1] = 1; + radius[2] = 0; + directMedian->SetRadius(radius); + directMedian->SetThresholdMultiplier(1.); + directMedian->SetInput(referenceReader->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(directMedian->Update()); + + auto denoisingProjectionsReader = rtk::ProjectionsReader::New(); + denoisingProjectionsReader->SetFileNames(std::vector{ argv[1] }); + denoisingProjectionsReader->SetMedianRadius(radius); + denoisingProjectionsReader->SetConditionalMedianThresholdMultiplier(1.); + TRY_AND_EXIT_ON_ITK_EXCEPTION(denoisingProjectionsReader->Update()); + + CheckVariableLengthVectorImageQuality( + denoisingProjectionsReader->GetOutput(), directMedian->GetOutput(), 0.04, 60, 2000.0); + + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/rtkspectralonesteptest.cxx b/test/rtkspectralonesteptest.cxx index 4fe513af9..81e35b720 100644 --- a/test/rtkspectralonesteptest.cxx +++ b/test/rtkspectralonesteptest.cxx @@ -1,17 +1,22 @@ #include "rtkTest.h" #include "rtkConstantImageSource.h" #include "rtkDrawEllipsoidImageFilter.h" +#include "rtkImageToVectorImageFilter.h" #include "rtkMechlemOneStepSpectralReconstructionFilter.h" #include "rtkRayEllipsoidIntersectionImageFilter.h" #include "rtkSpectralForwardModelImageFilter.h" +#include "rtkVectorImageToImageFilter.h" #ifdef USE_CUDA # include "itkCudaImage.h" #endif #include +#include +#include #include #include +#include #include /** @@ -61,6 +66,14 @@ rtkspectralonesteptest(int argc, char * argv[]) using MaterialAttenuationsImageType = itk::Image; using SingleComponentImageType = itk::Image; #endif + using ScalarSpectrumImageType = itk::Image; + using VectorSpectrumImageType = itk::VectorImage; + using PermuteIncidentSpectrumFilterType = itk::PermuteAxesImageFilter; + using ChangeIncidentSpectrumInformationFilterType = itk::ChangeInformationImageFilter; + using IncidentSpectrumToVectorFilterType = + rtk::ImageToVectorImageFilter; + using VectorSpectrumToScalarFilterType = + rtk::VectorImageToImageFilter; // Cast filters to convert between vector image types @@ -71,12 +84,16 @@ rtkspectralonesteptest(int argc, char * argv[]) incidentSpectrumReader->SetFileName(argv[1]); incidentSpectrumReader->Update(); + auto scalarIncidentSpectrumReader = itk::ImageFileReader::New(); + scalarIncidentSpectrumReader->SetFileName(argv[1]); + scalarIncidentSpectrumReader->Update(); + auto detectorResponseReader = itk::ImageFileReader::New(); - detectorResponseReader->SetFileName(argv[3]); + detectorResponseReader->SetFileName(argv[2]); detectorResponseReader->Update(); auto materialAttenuationsReader = itk::ImageFileReader::New(); - materialAttenuationsReader->SetFileName(argv[4]); + materialAttenuationsReader->SetFileName(argv[3]); materialAttenuationsReader->Update(); #if FAST_TESTS_NO_CHECKS @@ -306,5 +323,60 @@ rtkspectralonesteptest(int argc, char * argv[]) CheckVectorImageQuality(mechlemOneStep->GetOutput(), composeVols->GetOutput(), 0.08, 23, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; + auto scalarIncidentSpectrumResult = itk::ImageDuplicator::New(); + scalarIncidentSpectrumResult->SetInputImage(mechlemOneStep->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(scalarIncidentSpectrumResult->Update()) + + auto permuteIncidentSpectrum = PermuteIncidentSpectrumFilterType::New(); + permuteIncidentSpectrum->SetInput(scalarIncidentSpectrumReader->GetOutput()); + typename PermuteIncidentSpectrumFilterType::PermuteOrderArrayType order; + order[0] = 1; + order[1] = 2; + order[2] = 0; + permuteIncidentSpectrum->SetOrder(order); + + auto changeIncidentSpectrumInformation = ChangeIncidentSpectrumInformationFilterType::New(); + changeIncidentSpectrumInformation->SetInput(permuteIncidentSpectrum->GetOutput()); + typename ScalarSpectrumImageType::DirectionType identityDirection; + identityDirection.SetIdentity(); + changeIncidentSpectrumInformation->SetOutputDirection(identityDirection); + changeIncidentSpectrumInformation->ChangeDirectionOn(); + + auto vectorIncidentSpectrum = IncidentSpectrumToVectorFilterType::New(); + vectorIncidentSpectrum->SetInput(changeIncidentSpectrumInformation->GetOutput()); + + // Validate the test input itself: the vector spectrum must round-trip to the original scalar spectrum before it is + // used to compare the overloaded SetInputIncidentSpectrum path. + auto scalarIncidentSpectrumFromVector = VectorSpectrumToScalarFilterType::New(); + scalarIncidentSpectrumFromVector->SetInput(vectorIncidentSpectrum->GetOutput()); + + auto inversePermuteIncidentSpectrum = PermuteIncidentSpectrumFilterType::New(); + inversePermuteIncidentSpectrum->SetInput(scalarIncidentSpectrumFromVector->GetOutput()); + order[0] = 2; + order[1] = 0; + order[2] = 1; + inversePermuteIncidentSpectrum->SetOrder(order); + TRY_AND_EXIT_ON_ITK_EXCEPTION(inversePermuteIncidentSpectrum->Update()) + CheckImageQuality( + inversePermuteIncidentSpectrum->GetOutput(), scalarIncidentSpectrumReader->GetOutput(), 1e-7, 120, 1000.0); + +#ifdef RTK_USE_CUDA + std::cout << "\n\n****** Case 6: CUDA voxel-based Backprojector, 4 subsets, with regularization, " + "itkVectorImage incident spectrum converted from scalar input ******" + << std::endl; +#else + std::cout << "\n\n****** Case 5: Voxel-based Backprojector, 4 subsets, with regularization, " + "itkVectorImage incident spectrum converted from scalar input ******" + << std::endl; +#endif + + mechlemOneStep->SetInputIncidentSpectrum(vectorIncidentSpectrum->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(mechlemOneStep->Update()); + + CheckVectorImageQuality( + mechlemOneStep->GetOutput(), scalarIncidentSpectrumResult->GetOutput(), 1e-7, 10, 1.0); + CheckVectorImageQuality(mechlemOneStep->GetOutput(), composeVols->GetOutput(), 0.08, 23, 2.0); + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; } diff --git a/wrapping/rtkConditionalMedianImageFilter.wrap b/wrapping/rtkConditionalMedianImageFilter.wrap index 2bb30818f..e7ddc11a1 100644 --- a/wrapping/rtkConditionalMedianImageFilter.wrap +++ b/wrapping/rtkConditionalMedianImageFilter.wrap @@ -1,3 +1,4 @@ itk_wrap_class("rtk::ConditionalMedianImageFilter" POINTER) itk_wrap_image_filter("${WRAP_ITK_REAL}" 1) + itk_wrap_template("VI${ITKM_F}3" "itk::VectorImage<${ITKT_F}, 3>") itk_end_wrap_class() diff --git a/wrapping/rtkProjectionsReader.wrap b/wrapping/rtkProjectionsReader.wrap index ad7649256..1a259d3e0 100644 --- a/wrapping/rtkProjectionsReader.wrap +++ b/wrapping/rtkProjectionsReader.wrap @@ -6,5 +6,6 @@ itk_wrap_class("rtk::ProjectionsReader" POINTER) #WARNING: "UC;${WRAP_ITK_REAL}" is too restrictive. # Types should be similar to vnl wrapped types : # "D;F;SI;SL;LD;SC;UC;US;UI;UL" - #Removing VectorImage<> (VI) and WRAP_ITK_SCALAR + #Removing WRAP_ITK_SCALAR + itk_wrap_template("VI${ITKM_F}3" "itk::VectorImage<${ITKT_F}, 3>") itk_end_wrap_class()