Skip to content

Commit

Permalink
ENH: Update geometry and signal in GenerateOutputInformation
Browse files Browse the repository at this point in the history 
A filter generating a geometry should do it in GenerateOutputInformation
  • Loading branch information
@SimonRit
Simon Rit authored and SimonRit committed Sep 30, 2024
1 parent 9cbea00 commit 2c8effc
Show file tree
Hide file tree
Showing 3 changed files with 56 additions and 47 deletions.
14 changes: 4 additions & 10 deletions include/rtkMechlemOneStepSpectralReconstructionFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,6 +341,10 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
m_ReorderPhotonCountsFilter->SetInputGeometry(this->m_Geometry);
m_ReorderPhotonCountsFilter->SetPermutation(ReorderProjectionsFilterPhotonCountsType::SHUFFLE);
m_ExtractPhotonCountsFilter->SetInput(m_ReorderPhotonCountsFilter->GetOutput());
m_ForwardProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_SingleComponentForwardProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_GradientsBackProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
m_HessiansBackProjectionFilter->SetGeometry(m_ReorderPhotonCountsFilter->GetOutputGeometry());
}
else
{
Expand Down Expand Up @@ -449,16 +453,6 @@ 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
5 changes: 5 additions & 0 deletions include/rtkReorderProjectionsImageFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,9 @@ class ITK_TEMPLATE_EXPORT ReorderProjectionsImageFilter : public itk::ImageToIma
void
VerifyPreconditions() const override;

void
GenerateOutputInformation() override;

void
GenerateData() override;

Expand All @@ -111,6 +114,8 @@ class ITK_TEMPLATE_EXPORT ReorderProjectionsImageFilter : public itk::ImageToIma
/** Permutation type */
PermutationType m_Permutation;

/** Indices in the original stack for the new stack */
std::vector<unsigned int> m_NewIndices;
}; // end of class

} // end namespace rtk
Expand Down
84 changes: 47 additions & 37 deletions include/rtkReorderProjectionsImageFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,17 +65,19 @@ ReorderProjectionsImageFilter<TInputImage, TOutputImage>::VerifyPreconditions()

template <class TInputImage, class TOutputImage>
void
ReorderProjectionsImageFilter<TInputImage, TOutputImage>::GenerateData()
ReorderProjectionsImageFilter<TInputImage, TOutputImage>::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();

unsigned int NumberOfProjections =
this->GetInput()->GetLargestPossibleRegion().GetSize()[TInputImage::ImageDimension - 1];
std::vector<unsigned int> permutation;
m_NewIndices.clear();
switch (m_Permutation)
{
case (NONE):
{
for (unsigned int i = 0; i < NumberOfProjections; i++)
permutation.push_back(i);
m_NewIndices.push_back(i);
break;
}
case (SORT):
Expand All @@ -92,44 +94,72 @@ ReorderProjectionsImageFilter<TInputImage, TOutputImage>::GenerateData()

// Extract the permutated indices
for (unsigned int i = 0; i < NumberOfProjections; i++)
permutation.push_back(pairsVector[i].second);
m_NewIndices.push_back(pairsVector[i].second);
break;
}
case (SHUFFLE):
{
for (unsigned int i = 0; i < NumberOfProjections; i++)
permutation.push_back(i);
m_NewIndices.push_back(i);
std::default_random_engine randomGenerator(0); // The seed is hard-coded to 0 to make the behavior reproducible
std::shuffle(permutation.begin(), permutation.end(), randomGenerator);
std::shuffle(m_NewIndices.begin(), m_NewIndices.end(), randomGenerator);
break;
}
default:
itkGenericExceptionMacro(<< "Unhandled projection reordering method");
}

// Allocate the pixels of the output, and at first fill them with zeros
// Initialize objects (otherwise, if the filter runs several times,
// the outputs become incorrect)
m_OutputGeometry->Clear();
m_OutputSignal.clear();

// Perform the copies
for (unsigned int proj = 0; proj < NumberOfProjections; proj++)
{
// Copy the geometry
m_OutputGeometry->SetRadiusCylindricalDetector(m_InputGeometry->GetRadiusCylindricalDetector());
m_OutputGeometry->AddProjectionInRadians(m_InputGeometry->GetSourceToIsocenterDistances()[m_NewIndices[proj]],
m_InputGeometry->GetSourceToDetectorDistances()[m_NewIndices[proj]],
m_InputGeometry->GetGantryAngles()[m_NewIndices[proj]],
m_InputGeometry->GetProjectionOffsetsX()[m_NewIndices[proj]],
m_InputGeometry->GetProjectionOffsetsY()[m_NewIndices[proj]],
m_InputGeometry->GetOutOfPlaneAngles()[m_NewIndices[proj]],
m_InputGeometry->GetInPlaneAngles()[m_NewIndices[proj]],
m_InputGeometry->GetSourceOffsetsX()[m_NewIndices[proj]],
m_InputGeometry->GetSourceOffsetsY()[m_NewIndices[proj]]);
m_OutputGeometry->SetCollimationOfLastProjection(m_InputGeometry->GetCollimationUInf()[m_NewIndices[proj]],
m_InputGeometry->GetCollimationUSup()[m_NewIndices[proj]],
m_InputGeometry->GetCollimationVInf()[m_NewIndices[proj]],
m_InputGeometry->GetCollimationVSup()[m_NewIndices[proj]]);

// Copy the signal, if any
if (m_Permutation == SORT)
m_OutputSignal.push_back(m_InputSignal[m_NewIndices[proj]]);
}
}

template <class TInputImage, class TOutputImage>
void
ReorderProjectionsImageFilter<TInputImage, TOutputImage>::GenerateData()
{
// Allocate the pixels of the output
this->GetOutput()->SetBufferedRegion(this->GetOutput()->GetRequestedRegion());
this->GetOutput()->Allocate();
this->GetOutput()->FillBuffer(itk::NumericTraits<typename TInputImage::PixelType>::ZeroValue());

// Declare regions used in the loop
typename TInputImage::RegionType inputRegion = this->GetOutput()->GetRequestedRegion();
typename TInputImage::RegionType outputRegion = this->GetOutput()->GetRequestedRegion();
inputRegion.SetSize(2, 1);
outputRegion.SetSize(2, 1);

// Initialize objects (otherwise, if the filter runs several times,
// the outputs become incorrect)
m_OutputGeometry->Clear();
m_OutputSignal.clear();

// Perform the copies
for (unsigned int proj = 0; proj < this->GetOutput()->GetRequestedRegion().GetSize()[2]; proj++)
// Copy the projection data
for (unsigned int i = 0; i < this->GetOutput()->GetRequestedRegion().GetSize()[2]; i++)
{
// Copy the projection data
const unsigned int proj = i + this->GetOutput()->GetRequestedRegion().GetIndex()[2];

// Regions
inputRegion.SetIndex(2, permutation[proj]);
inputRegion.SetIndex(2, m_NewIndices[proj]);
outputRegion.SetIndex(2, proj);

itk::ImageRegionConstIterator<TInputImage> inputProjsIt(this->GetInput(), inputRegion);
Expand All @@ -142,26 +172,6 @@ ReorderProjectionsImageFilter<TInputImage, TOutputImage>::GenerateData()
++outputProjsIt;
++inputProjsIt;
}

// Copy the geometry
m_OutputGeometry->SetRadiusCylindricalDetector(m_InputGeometry->GetRadiusCylindricalDetector());
m_OutputGeometry->AddProjectionInRadians(m_InputGeometry->GetSourceToIsocenterDistances()[permutation[proj]],
m_InputGeometry->GetSourceToDetectorDistances()[permutation[proj]],
m_InputGeometry->GetGantryAngles()[permutation[proj]],
m_InputGeometry->GetProjectionOffsetsX()[permutation[proj]],
m_InputGeometry->GetProjectionOffsetsY()[permutation[proj]],
m_InputGeometry->GetOutOfPlaneAngles()[permutation[proj]],
m_InputGeometry->GetInPlaneAngles()[permutation[proj]],
m_InputGeometry->GetSourceOffsetsX()[permutation[proj]],
m_InputGeometry->GetSourceOffsetsY()[permutation[proj]]);
m_OutputGeometry->SetCollimationOfLastProjection(m_InputGeometry->GetCollimationUInf()[permutation[proj]],
m_InputGeometry->GetCollimationUSup()[permutation[proj]],
m_InputGeometry->GetCollimationVInf()[permutation[proj]],
m_InputGeometry->GetCollimationVSup()[permutation[proj]]);

// Copy the signal, if any
if (m_Permutation == SORT)
m_OutputSignal.push_back(m_InputSignal[permutation[proj]]);
}
}

Expand Down

0 comments on commit 2c8effc

Please sign in to comment.