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 shuffled before doing the
reconstruction. The shuffle was done in the application and not in the
filter which is not how it is done for the others reconstruction filters.
  • Loading branch information
@arobert01 @SimonRit
arobert01 authored and SimonRit committed Sep 26, 2024
1 parent 85e92df commit 9cbea00
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 60 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
59 changes: 35 additions & 24 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 @@ -97,12 +98,17 @@ namespace rtk
* MultiplySupport [ label="itk::MultiplyImageFilter" URL="\ref itk::MultiplyImageFilter" style=dashed];
* Alphak [ label="", fixedsize="false", width=0, height=0, shape=none];
* NextAlphak [ label="", fixedsize="false", width=0, height=0, shape=none];
* ReorderProjections [label="rtk::ReorderProjectionsImageFilter (Photon Counts)"
* URL="\ref ReorderProjectionsImageFilter" style=dashed];
* ReorderProjectionsWeights [label="rtk::ReorderProjectionsImageFilter (Projection Weights)"
* URL="\ref ReorderProjectionsImageFilter" style=dashed];
*
* Input0 -> Alphak [arrowhead=none];
* Alphak -> ForwardProjection;
* Alphak -> SQSRegul;
* ProjectionsSource -> ForwardProjection;
* Input1 -> Extract;
* Input1 -> ReorderProjections;
* ReorderProjections -> Extract;
* Extract -> Weidinger;
* Input2 -> Weidinger;
* ForwardProjection -> Weidinger;
Expand All @@ -116,7 +122,8 @@ namespace rtk
* Input4 -> MultiplyRegulGradient;
* SQSRegul -> MultiplyRegulGradient;
* MultiplyRegulGradient -> AddGradients;
* Input5 -> MultiplyGradientToBeBackProjected;
* Input5 -> ReorderProjectionsWeights;
* ReorderProjectionsWeights-> MultiplyGradientToBeBackProjected;
* Weidinger -> MultiplyGradientToBeBackProjected
* MultiplyGradientToBeBackProjected -> AddGradients;
* AddGradients -> Newton;
Expand Down Expand Up @@ -233,10 +240,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 +305,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 +366,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
36 changes: 34 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,16 @@ 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();
m_ForwardProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_SingleComponentForwardProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_GradientsBackProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_HessiansBackProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
}

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

0 comments on commit 9cbea00

Please sign in to comment.