Skip to content

Commit

Permalink
ENH: Reorder the photon counts inside one-step filter
Browse files Browse the repository at this point in the history 
When subsets are used, the photon counts are shuffle before doing the
reconstruction. This shuffle was done in the application and not in the
filter which was not in adequation with what is done for the others
reconstruction filters.
  • Loading branch information
@arobert01
arobert01 committed Sep 24, 2024
1 parent 5cd1f78 commit e65e2db
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 58 deletions.
38 changes: 4 additions & 34 deletions applications/rtkspectralonestep/rtkspectralonestep.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,40 +201,10 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
mechlemOneStep->SetSupportMask(supportmask);
if (args_info.regul_spatial_weights_given)
mechlemOneStep->SetSpatialRegularizationWeights(spatialRegulWeighs);

// If subsets are used, reorder projections and geometry according to
// a random permutation
if (args_info.subsets_arg != 1)
{
using ReorderProjectionsFilterType = rtk::ReorderProjectionsImageFilter<PhotonCountsType>;
typename ReorderProjectionsFilterType::Pointer reorder_PhotonsCounts = ReorderProjectionsFilterType::New();
reorder_PhotonsCounts->SetInput(photonCounts);
reorder_PhotonsCounts->SetInputGeometry(geometry);
reorder_PhotonsCounts->SetPermutation(rtk::ReorderProjectionsImageFilter<PhotonCountsType>::SHUFFLE);
reorder_PhotonsCounts->Update();
mechlemOneStep->SetInputPhotonCounts(reorder_PhotonsCounts->GetOutput());
mechlemOneStep->SetGeometry(reorder_PhotonsCounts->GetOutputGeometry());
if (args_info.projection_weights_given)
{
// N.B: The permutation of photon counts and weights correspond because we use the SHUFFLE case mode initialized
// to zero (cf case (SHUFFLE): in ReorderProjectionsImageFilter.hxx)
using ReorderWeightsFilterType = rtk::ReorderProjectionsImageFilter<IncidentSpectrumType>;
typename ReorderWeightsFilterType::Pointer reorder_ProjectionWeights = ReorderWeightsFilterType::New();
reorder_ProjectionWeights->SetInput(projectionWeights);
// Here we get geometry which has not been shuffled
reorder_ProjectionWeights->SetInputGeometry(geometry);
reorder_ProjectionWeights->SetPermutation(rtk::ReorderProjectionsImageFilter<IncidentSpectrumType>::SHUFFLE);
reorder_ProjectionWeights->Update();
mechlemOneStep->SetProjectionWeights(reorder_ProjectionWeights->GetOutput());
}
}
else
{
mechlemOneStep->SetInputPhotonCounts(photonCounts);
mechlemOneStep->SetGeometry(geometry);
if (args_info.projection_weights_given)
mechlemOneStep->SetProjectionWeights(projectionWeights);
}
mechlemOneStep->SetInputPhotonCounts(photonCounts);
mechlemOneStep->SetGeometry(geometry);
if (args_info.projection_weights_given)
mechlemOneStep->SetProjectionWeights(projectionWeights);

REPORT_ITERATIONS(mechlemOneStep, MechlemFilterType, MaterialVolumesType);

Expand Down
49 changes: 27 additions & 22 deletions include/rtkMechlemOneStepSpectralReconstructionFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "rtkNesterovUpdateImageFilter.h"
#include "rtkSeparableQuadraticSurrogateRegularizationImageFilter.h"
#include "rtkAddMatrixAndDiagonalImageFilter.h"
#include "rtkReorderProjectionsImageFilter.h"

#include <itkExtractImageFilter.h>
#include <itkAddImageFilter.h>
Expand Down Expand Up @@ -233,10 +234,12 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter
using NewtonFilterType = rtk::GetNewtonUpdateImageFilter<GradientsImageType, HessiansImageType>;
using MultiplyFilterType = itk::MultiplyImageFilter<TOutputImage, SingleComponentImageType>;
using MultiplyGradientFilterType = itk::MultiplyImageFilter<GradientsImageType, SingleComponentImageType>;
using ReorderProjectionsFilterPhotonCountsType = rtk::ReorderProjectionsImageFilter<TPhotonCounts>;
using ReorderProjectionsFilterProjectionsWeightsType = rtk::ReorderProjectionsImageFilter<SingleComponentImageType>;
#endif

/** Pass the geometry to all filters needing it */
itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);
itkSetObjectMacro(Geometry, ThreeDCircularProjectionGeometry);

itkSetMacro(NumberOfIterations, int);
itkGetMacro(NumberOfIterations, int);
Expand Down Expand Up @@ -296,26 +299,28 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter

#if !defined(ITK_WRAPPING_PARSER)
/** Member pointers to the filters used internally (for convenience)*/
typename ExtractPhotonCountsFilterType::Pointer m_ExtractPhotonCountsFilter;
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 ExtractPhotonCountsFilterType::Pointer m_ExtractPhotonCountsFilter;
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 ReorderProjectionsFilterPhotonCountsType::Pointer m_ReorderPhotonCountsFilter;
typename ReorderProjectionsFilterProjectionsWeightsType::Pointer m_ReorderProjectionsWeightsFilter;
#endif

/** The inputs of this filter have the same type but not the same meaning
Expand Down Expand Up @@ -355,7 +360,7 @@ class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter
InstantiateHessiansBackProjectionFilter(int bptype);
#endif

ThreeDCircularProjectionGeometry::ConstPointer m_Geometry;
ThreeDCircularProjectionGeometry::Pointer m_Geometry;

int m_NumberOfIterations;
int m_NumberOfProjectionsPerSubset;
Expand Down
37 changes: 35 additions & 2 deletions include/rtkMechlemOneStepSpectralReconstructionFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,8 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
m_NewtonFilter = NewtonFilterType::New();
m_NesterovFilter = NesterovFilterType::New();
m_MultiplySupportFilter = MultiplyFilterType::New();
m_ReorderPhotonCountsFilter = ReorderProjectionsFilterPhotonCountsType::New();
m_ReorderProjectionsWeightsFilter = ReorderProjectionsFilterProjectionsWeightsType::New();

// Set permanent parameters
m_ProjectionsSource->SetConstant(itk::NumericTraits<typename TOutputImage::PixelType>::ZeroValue());
Expand Down Expand Up @@ -333,7 +335,17 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
// Set runtime connections. Links with the forward and back projection filters should be set here,
// since those filters are not instantiated by the constructor, but by
// a call to SetForwardProjectionFilter() and SetBackProjectionFilter()
m_ExtractPhotonCountsFilter->SetInput(this->GetInputPhotonCounts());
if (m_NumberOfSubsets != 1)
{
m_ReorderPhotonCountsFilter->SetInput(this->GetInputPhotonCounts());
m_ReorderPhotonCountsFilter->SetInputGeometry(this->m_Geometry);
m_ReorderPhotonCountsFilter->SetPermutation(ReorderProjectionsFilterPhotonCountsType::SHUFFLE);
m_ExtractPhotonCountsFilter->SetInput(m_ReorderPhotonCountsFilter->GetOutput());
}
else
{
m_ExtractPhotonCountsFilter->SetInput(this->GetInputPhotonCounts());
}

m_ForwardProjectionFilter->SetInput(0, m_ProjectionsSource->GetOutput());
m_ForwardProjectionFilter->SetInput(1, this->GetInputMaterialVolumes());
Expand All @@ -350,7 +362,17 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
if (this->GetProjectionWeights().GetPointer() != nullptr)
{
m_MultiplyGradientToBeBackprojectedFilter->SetInput1(m_WeidingerForward->GetOutput1());
m_MultiplyGradientToBeBackprojectedFilter->SetInput2(this->GetProjectionWeights());
if (m_NumberOfSubsets != 1)
{
m_ReorderProjectionsWeightsFilter->SetInput(this->GetProjectionWeights());
m_ReorderProjectionsWeightsFilter->SetInputGeometry(this->m_Geometry);
m_ReorderProjectionsWeightsFilter->SetPermutation(ReorderProjectionsFilterProjectionsWeightsType::SHUFFLE);
m_MultiplyGradientToBeBackprojectedFilter->SetInput2(m_ReorderProjectionsWeightsFilter->GetOutput());
}
else
{
m_MultiplyGradientToBeBackprojectedFilter->SetInput2(this->GetProjectionWeights());
}
m_GradientsBackProjectionFilter->SetInput(1, m_MultiplyGradientToBeBackprojectedFilter->GetOutput());
}
else
Expand Down Expand Up @@ -427,6 +449,17 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
{
itk::IterationReporter iterationReporter(this, 0, 1);

// If subsets are used the photons counts and geometry are shuffled, then geometry need to be set again.
if (m_NumberOfSubsets != 1)
{
m_ReorderPhotonCountsFilter->Update();
this->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_ForwardProjectionFilter->SetGeometry(this->m_Geometry);
m_SingleComponentForwardProjectionFilter->SetGeometry(this->m_Geometry);
m_GradientsBackProjectionFilter->SetGeometry(this->m_Geometry.GetPointer());
m_HessiansBackProjectionFilter->SetGeometry(this->m_Geometry.GetPointer());
}

// Run the iteration loop
typename TOutputImage::Pointer Next_Zk;
for (int iter = 0; iter < m_NumberOfIterations; iter++)
Expand Down

0 comments on commit e65e2db

Please sign in to comment.