diff --git a/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx b/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx index caddfda70..04cd0e0ff 100644 --- a/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx +++ b/include/rtkMechlemOneStepSpectralReconstructionFilter.hxx @@ -341,6 +341,10 @@ MechlemOneStepSpectralReconstructionFilterSetInputGeometry(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 { @@ -449,16 +453,6 @@ MechlemOneStepSpectralReconstructionFilterUpdate(); - 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++) diff --git a/include/rtkReorderProjectionsImageFilter.h b/include/rtkReorderProjectionsImageFilter.h index 07405fc61..9e651a5ce 100644 --- a/include/rtkReorderProjectionsImageFilter.h +++ b/include/rtkReorderProjectionsImageFilter.h @@ -96,6 +96,9 @@ class ITK_TEMPLATE_EXPORT ReorderProjectionsImageFilter : public itk::ImageToIma void VerifyPreconditions() const override; + void + GenerateOutputInformation() override; + void GenerateData() override; @@ -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 m_NewIndices; }; // end of class } // end namespace rtk diff --git a/include/rtkReorderProjectionsImageFilter.hxx b/include/rtkReorderProjectionsImageFilter.hxx index 8513c07e6..5e015ad87 100644 --- a/include/rtkReorderProjectionsImageFilter.hxx +++ b/include/rtkReorderProjectionsImageFilter.hxx @@ -65,17 +65,19 @@ ReorderProjectionsImageFilter::VerifyPreconditions() template void -ReorderProjectionsImageFilter::GenerateData() +ReorderProjectionsImageFilter::GenerateOutputInformation() { + Superclass::GenerateOutputInformation(); + unsigned int NumberOfProjections = this->GetInput()->GetLargestPossibleRegion().GetSize()[TInputImage::ImageDimension - 1]; - std::vector 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): @@ -92,25 +94,58 @@ ReorderProjectionsImageFilter::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 +void +ReorderProjectionsImageFilter::GenerateData() +{ + // Allocate the pixels of the output this->GetOutput()->SetBufferedRegion(this->GetOutput()->GetRequestedRegion()); this->GetOutput()->Allocate(); - this->GetOutput()->FillBuffer(itk::NumericTraits::ZeroValue()); // Declare regions used in the loop typename TInputImage::RegionType inputRegion = this->GetOutput()->GetRequestedRegion(); @@ -118,18 +153,13 @@ ReorderProjectionsImageFilter::GenerateData() 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 inputProjsIt(this->GetInput(), inputRegion); @@ -142,26 +172,6 @@ ReorderProjectionsImageFilter::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]]); } }