Skip to content

Commit

Permalink
STYLE: Rename variables and types in spectral applications
Browse files Browse the repository at this point in the history
This commit makes variable and type names in spectral applications consistent
with the recent changes in spectral filters
  • Loading branch information
cyrilmory authored and SimonRit committed Feb 11, 2025
1 parent faf4926 commit cd0c841
Show file tree
Hide file tree
Showing 9 changed files with 212 additions and 168 deletions.
21 changes: 11 additions & 10 deletions applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ main(int argc, char * argv[])
using PixelValueType = float;
constexpr unsigned int Dimension = 3;
using DecomposedProjectionType = itk::VectorImage<PixelValueType, Dimension>;
using SpectralProjectionsType = itk::VectorImage<PixelValueType, Dimension>;
using MeasuredProjectionsType = itk::VectorImage<PixelValueType, Dimension>;
using IncidentSpectrumImageType = itk::VectorImage<PixelValueType, Dimension - 1>;
using DetectorResponseImageType = itk::Image<PixelValueType, Dimension - 1>;
using MaterialAttenuationsImageType = itk::Image<PixelValueType, Dimension - 1>;
Expand All @@ -58,10 +58,10 @@ main(int argc, char * argv[])
const unsigned int MaximumEnergy = incidentSpectrum->GetVectorLength();

// Generate a set of zero-filled photon count projections
SpectralProjectionsType::Pointer photonCounts = SpectralProjectionsType::New();
photonCounts->CopyInformation(decomposedProjection);
photonCounts->SetVectorLength(NumberOfSpectralBins);
photonCounts->Allocate();
MeasuredProjectionsType::Pointer measuredProjections = MeasuredProjectionsType::New();
measuredProjections->CopyInformation(decomposedProjection);
measuredProjections->SetVectorLength(NumberOfSpectralBins);
measuredProjections->Allocate();

// Read the thresholds on command line
itk::VariableLengthVector<unsigned int> thresholds;
Expand All @@ -81,11 +81,12 @@ main(int argc, char * argv[])
<< decomposedProjection->GetPixel(indexDecomp).Size() << ", should be "
<< NumberOfMaterials);

SpectralProjectionsType::IndexType indexSpect;
MeasuredProjectionsType::IndexType indexSpect;
indexSpect.Fill(0);
if (photonCounts->GetPixel(indexSpect).Size() != NumberOfSpectralBins)
if (measuredProjections->GetPixel(indexSpect).Size() != NumberOfSpectralBins)
itkGenericExceptionMacro(<< "Spectral projections (i.e. photon count data) image has vector size "
<< photonCounts->GetPixel(indexSpect).Size() << ", should be " << NumberOfSpectralBins);
<< measuredProjections->GetPixel(indexSpect).Size() << ", should be "
<< NumberOfSpectralBins);

IncidentSpectrumImageType::IndexType indexIncident;
indexIncident.Fill(0);
Expand All @@ -100,10 +101,10 @@ main(int argc, char * argv[])

// Create and set the filter
using ForwardModelFilterType =
rtk::SpectralForwardModelImageFilter<DecomposedProjectionType, SpectralProjectionsType, IncidentSpectrumImageType>;
rtk::SpectralForwardModelImageFilter<DecomposedProjectionType, MeasuredProjectionsType, IncidentSpectrumImageType>;
ForwardModelFilterType::Pointer forward = ForwardModelFilterType::New();
forward->SetInputDecomposedProjections(decomposedProjection);
forward->SetInputMeasuredProjections(photonCounts);
forward->SetInputMeasuredProjections(measuredProjections);
forward->SetInputIncidentSpectrum(incidentSpectrum);
forward->SetDetectorResponse(detectorResponse);
forward->SetMaterialAttenuations(materialAttenuations);
Expand Down
16 changes: 8 additions & 8 deletions applications/rtkspectralonestep/rtkspectralonestep.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -58,21 +58,21 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
// Define types for the input images
#ifdef RTK_USE_CUDA
using MaterialVolumesType = typename itk::CudaImage<itk::Vector<dataType, VNumberOfMaterials>, Dimension>;
using PhotonCountsType = typename itk::CudaImage<itk::Vector<dataType, VNumberOfBins>, Dimension>;
using MeasuredProjectionsType = typename itk::CudaImage<itk::Vector<dataType, VNumberOfBins>, Dimension>;
using IncidentSpectrumType = itk::CudaImage<dataType, Dimension>;
using DetectorResponseType = itk::CudaImage<dataType, 2>;
using MaterialAttenuationsType = itk::CudaImage<dataType, 2>;
#else
using MaterialVolumesType = typename itk::Image<itk::Vector<dataType, VNumberOfMaterials>, Dimension>;
using PhotonCountsType = typename itk::Image<itk::Vector<dataType, VNumberOfBins>, Dimension>;
using MeasuredProjectionsType = typename itk::Image<itk::Vector<dataType, VNumberOfBins>, Dimension>;
using IncidentSpectrumType = itk::Image<dataType, Dimension>;
using DetectorResponseType = itk::Image<dataType, 2>;
using MaterialAttenuationsType = itk::Image<dataType, 2>;
#endif

// Instantiate and update the readers
typename PhotonCountsType::Pointer photonCounts;
TRY_AND_EXIT_ON_ITK_EXCEPTION(photonCounts = itk::ReadImage<PhotonCountsType>(args_info.spectral_arg))
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))
Expand Down Expand Up @@ -183,7 +183,7 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)

// Set the forward and back projection filters to be used
using MechlemFilterType = typename rtk::
MechlemOneStepSpectralReconstructionFilter<MaterialVolumesType, PhotonCountsType, IncidentSpectrumType>;
MechlemOneStepSpectralReconstructionFilter<MaterialVolumesType, MeasuredProjectionsType, IncidentSpectrumType>;
typename MechlemFilterType::Pointer mechlemOneStep = MechlemFilterType::New();
SetForwardProjectionFromGgo(args_info, mechlemOneStep.GetPointer());
SetBackProjectionFromGgo(args_info, mechlemOneStep.GetPointer());
Expand All @@ -201,7 +201,7 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
mechlemOneStep->SetSupportMask(supportmask);
if (args_info.regul_spatial_weights_given)
mechlemOneStep->SetSpatialRegularizationWeights(spatialRegulWeighs);
mechlemOneStep->SetInputMeasuredProjections(photonCounts);
mechlemOneStep->SetInputMeasuredProjections(mea);
mechlemOneStep->SetGeometry(geometry);
if (args_info.projection_weights_given)
mechlemOneStep->SetProjectionWeights(projectionWeights);
Expand All @@ -221,8 +221,8 @@ main(int argc, char * argv[])
GGO(rtkspectralonestep, args_info);
try
{
itk::ImageIOBase::Pointer headerInputPhotonCounts = GetFileHeader(args_info.spectral_arg);
unsigned int nBins = headerInputPhotonCounts->GetNumberOfComponents();
itk::ImageIOBase::Pointer headerInputMeasuredProjections = GetFileHeader(args_info.spectral_arg);
unsigned int nBins = headerInputMeasuredProjections->GetNumberOfComponents();
itk::ImageIOBase::Pointer headerAttenuations = GetFileHeader(args_info.attenuations_arg);
unsigned int nMaterials = headerAttenuations->GetDimensions(0);

Expand Down
38 changes: 22 additions & 16 deletions include/rtkCudaWeidingerForwardModelImageFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjection

template <class TDecomposedProjections, class TMeasuredProjections, class TIncidentSpectrum, class TProjections>
void
CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjections, TIncidentSpectrum, TProjections>::GPUGenerateData()
CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjections, TIncidentSpectrum, TProjections>::
GPUGenerateData()
{
this->AllocateOutputs();

Expand Down Expand Up @@ -71,21 +72,26 @@ CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjection
# endif

// Run the forward projection with a slab of SLAB_SIZE or less projections
CUDA_WeidingerForwardModel(
projectionSize,
this->m_MaterialAttenuations.data_block(),
this->m_BinnedDetectorResponse.data_block(),
pMatProj,
pPhoCount,
pSpectrum,
pProjOnes,
pOut1,
pOut2,
CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjections, TIncidentSpectrum, TProjections>::nBins,
nEnergies,
CudaWeidingerForwardModelImageFilter<TDecomposedProjections, TMeasuredProjections, TIncidentSpectrum, TProjections>::nMaterials,
this->m_NumberOfProjectionsInIncidentSpectrum,
this->GetInputMeasuredProjections()->GetBufferedRegion().GetIndex(Dimension - 1));
CUDA_WeidingerForwardModel(projectionSize,
this->m_MaterialAttenuations.data_block(),
this->m_BinnedDetectorResponse.data_block(),
pMatProj,
pPhoCount,
pSpectrum,
pProjOnes,
pOut1,
pOut2,
CudaWeidingerForwardModelImageFilter<TDecomposedProjections,
TMeasuredProjections,
TIncidentSpectrum,
TProjections>::nBins,
nEnergies,
CudaWeidingerForwardModelImageFilter<TDecomposedProjections,
TMeasuredProjections,
TIncidentSpectrum,
TProjections>::nMaterials,
this->m_NumberOfProjectionsInIncidentSpectrum,
this->GetInputMeasuredProjections()->GetBufferedRegion().GetIndex(Dimension - 1));
}

} // end namespace rtk
Expand Down
54 changes: 28 additions & 26 deletions include/rtkMechlemOneStepSpectralReconstructionFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -192,9 +192,10 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter

#if !defined(ITK_WRAPPING_PARSER)
# ifdef RTK_USE_CUDA
typedef typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
WeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>,
CudaWeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>>::type
typedef typename std::conditional<
std::is_same<TOutputImage, CPUOutputImageType>::value,
WeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>,
CudaWeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>>::type
WeidingerForwardModelType;
typedef
typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
Expand All @@ -206,7 +207,8 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter
CudaBackProjectionImageFilter<HessiansImageType>>::type
CudaHessiansBackProjectionImageFilterType;
# else
using WeidingerForwardModelType = WeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>;
using WeidingerForwardModelType =
WeidingerForwardModelImageFilter<TOutputImage, TMeasuredProjections, TIncidentSpectrum>;
using CudaSingleComponentForwardProjectionImageFilterType =
JosephForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>;
using CudaHessiansBackProjectionImageFilterType = BackProjectionImageFilter<HessiansImageType, HessiansImageType>;
Expand Down Expand Up @@ -307,28 +309,28 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter

#if !defined(ITK_WRAPPING_PARSER)
/** Member pointers to the filters used internally (for convenience)*/
typename ExtractMeasuredProjectionsFilterType::Pointer m_ExtractMeasuredProjectionsFilter;
typename AddFilterType::Pointer m_AddGradients;
typename SingleComponentForwardProjectionFilterType::Pointer m_SingleComponentForwardProjectionFilter;
typename MaterialProjectionsSourceType::Pointer m_ProjectionsSource;
typename SingleComponentImageSourceType::Pointer m_SingleComponentProjectionsSource;
typename SingleComponentImageSourceType::Pointer m_SingleComponentVolumeSource;
typename GradientsSourceType::Pointer m_GradientsSource;
typename HessiansSourceType::Pointer m_HessiansSource;
typename WeidingerForwardModelType::Pointer m_WeidingerForward;
typename SQSRegularizationType::Pointer m_SQSRegul;
typename AddMatrixAndDiagonalFilterType::Pointer m_AddHessians;
typename NewtonFilterType::Pointer m_NewtonFilter;
typename NesterovFilterType::Pointer m_NesterovFilter;
typename ForwardProjectionFilterType::Pointer m_ForwardProjectionFilter;
typename GradientsBackProjectionFilterType::Pointer m_GradientsBackProjectionFilter;
typename HessiansBackProjectionFilterType::Pointer m_HessiansBackProjectionFilter;
typename MultiplyFilterType::Pointer m_MultiplySupportFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyRegulGradientsFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyRegulHessiansFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyGradientToBeBackprojectedFilter;
typename ReorderMeasuredProjectionsFilterType::Pointer m_ReorderMeasuredProjectionsFilter;
typename ReorderProjectionsWeightsFilterType::Pointer m_ReorderProjectionsWeightsFilter;
typename ExtractMeasuredProjectionsFilterType::Pointer m_ExtractMeasuredProjectionsFilter;
typename AddFilterType::Pointer m_AddGradients;
typename SingleComponentForwardProjectionFilterType::Pointer m_SingleComponentForwardProjectionFilter;
typename MaterialProjectionsSourceType::Pointer m_ProjectionsSource;
typename SingleComponentImageSourceType::Pointer m_SingleComponentProjectionsSource;
typename SingleComponentImageSourceType::Pointer m_SingleComponentVolumeSource;
typename GradientsSourceType::Pointer m_GradientsSource;
typename HessiansSourceType::Pointer m_HessiansSource;
typename WeidingerForwardModelType::Pointer m_WeidingerForward;
typename SQSRegularizationType::Pointer m_SQSRegul;
typename AddMatrixAndDiagonalFilterType::Pointer m_AddHessians;
typename NewtonFilterType::Pointer m_NewtonFilter;
typename NesterovFilterType::Pointer m_NesterovFilter;
typename ForwardProjectionFilterType::Pointer m_ForwardProjectionFilter;
typename GradientsBackProjectionFilterType::Pointer m_GradientsBackProjectionFilter;
typename HessiansBackProjectionFilterType::Pointer m_HessiansBackProjectionFilter;
typename MultiplyFilterType::Pointer m_MultiplySupportFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyRegulGradientsFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyRegulHessiansFilter;
typename MultiplyGradientFilterType::Pointer m_MultiplyGradientToBeBackprojectedFilter;
typename ReorderMeasuredProjectionsFilterType::Pointer m_ReorderMeasuredProjectionsFilter;
typename ReorderProjectionsWeightsFilterType::Pointer m_ReorderProjectionsWeightsFilter;
#endif

/** The inputs of this filter have the same type but not the same meaning
Expand Down
Loading

0 comments on commit cd0c841

Please sign in to comment.