diff --git a/include/rtkExtractImageSubRegion.h b/include/rtkExtractImageSubRegion.h new file mode 100644 index 000000000..7bfae3c7d --- /dev/null +++ b/include/rtkExtractImageSubRegion.h @@ -0,0 +1,141 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef rtkExtractImageSubRegion_h +#define rtkExtractImageSubRegion_h + +#include +#include +#include + +namespace rtk +{ + +/** \class ExtractImageSubRegion + * \brief Create an image that is a view of a sub-region of another image, + * without copying pixel data when the region is contiguous. + * + * This is a lightweight alternative to itk::ExtractImageFilter for the case + * where input and output image types are the same dimension (no dimension + * collapse). It avoids the overhead of the filter pipeline machinery by + * directly creating an image that shares the same pixel buffer as the input, + * with adjusted metadata (origin, region). + * + * The zero-copy optimization only applies when the extraction region is + * contiguous in memory, i.e. all dimensions except the last span the full + * input extent. When the region is not contiguous, falls back to + * itk::ExtractImageFilter. + * + * When the input buffer is not yet allocated (e.g., during + * GenerateOutputInformation), the output image is created with the correct + * metadata only. When the buffer is available (e.g., during GenerateData), + * the pixel buffer is shared via SetImportPointer for zero-copy access. + * + * Warning: since the output shares the input's pixel buffer, downstream + * filters that operate in-place (InPlaceOn) will corrupt the source data. + * Callers must ensure that no in-place filter modifies this image's buffer. + * + * This is useful in mini-pipelines where a sub-stack of projections is + * repeatedly extracted from a projection stack (e.g., FDK, SART, OSEM). + * + * \author Axel Garcia + * + * \ingroup RTK + */ +/** Check if a sub-region is contiguous in the input buffer (zero-copy possible). + * True when all non-last dimensions match the input exactly. */ +template +bool +IsContiguousSubRegion(const TImage * input, const itk::ImageRegion & region) +{ + constexpr unsigned int Dimension = TImage::ImageDimension; + const auto & inputRegion = input->GetLargestPossibleRegion(); + for (unsigned int d = 0; d < Dimension - 1; ++d) + { + if (region.GetIndex()[d] != inputRegion.GetIndex()[d] || region.GetSize()[d] != inputRegion.GetSize()[d]) + return false; + } + return true; +} + +template +typename TImage::Pointer +ExtractImageSubRegion(const TImage * input, const itk::ImageRegion & extractionRegion) +{ + constexpr unsigned int Dimension = TImage::ImageDimension; + using PixelType = typename TImage::PixelType; + using RegionType = itk::ImageRegion; + using SizeType = itk::Size; + using IndexType = itk::Index; + using SpacingType = typename TImage::SpacingType; + using PointType = typename TImage::PointType; + using DirectionType = typename TImage::DirectionType; + + const RegionType & inputRegion = input->GetLargestPossibleRegion(); + const IndexType & inputIndex = inputRegion.GetIndex(); + + if (!IsContiguousSubRegion(input, extractionRegion)) + { + using ExtractFilterType = itk::ExtractImageFilter; + typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New(); + extractFilter->SetInput(input); + extractFilter->SetExtractionRegion(extractionRegion); + extractFilter->SetDirectionCollapseToSubmatrix(); + extractFilter->Update(); + return extractFilter->GetOutput(); + } + + const SpacingType & spacing = input->GetSpacing(); + const PointType & inputOrigin = input->GetOrigin(); + const DirectionType & direction = input->GetDirection(); + + // Create output with correct metadata (skips Allocate for CudaImage). + typename TImage::Pointer output = TImage::New(); + output->SetRegions(extractionRegion); + output->SetSpacing(spacing); + output->SetOrigin(inputOrigin); + output->SetDirection(direction); + + // If the input buffer is available, share it (zero-copy). + // Otherwise, return a metadata-only image. + if (input->GetBufferPointer()) + { + const IndexType & extractIndex = extractionRegion.GetIndex(); + + // Pixels per slice: product of all input sizes except the last. + typename SizeType::SizeValueType sliceSize = 1; + for (unsigned int d = 0; d < Dimension - 1; ++d) + sliceSize *= inputRegion.GetSize()[d]; + + const long sliceOffset = extractIndex[Dimension - 1] - inputIndex[Dimension - 1]; + const PixelType * bufferPtr = input->GetBufferPointer() + sliceOffset * sliceSize; + + const typename SizeType::SizeValueType numPixels = extractionRegion.GetNumberOfPixels(); + output->GetPixelContainer()->SetImportPointer(const_cast(bufferPtr), numPixels, false); + + // Re-assign pixel container to sync subclass containers (e.g. + // CudaDataManager reads the CPU pointer and marks GPU dirty). + output->SetPixelContainer(output->GetPixelContainer()); + } + + return output; +} + +} // namespace rtk + +#endif // rtkExtractImageSubRegion_h diff --git a/include/rtkFDKConeBeamReconstructionFilter.h b/include/rtkFDKConeBeamReconstructionFilter.h index 85bc1104f..9ca467b79 100644 --- a/include/rtkFDKConeBeamReconstructionFilter.h +++ b/include/rtkFDKConeBeamReconstructionFilter.h @@ -23,8 +23,7 @@ #include "rtkConfiguration.h" #include "rtkFDKBackProjectionImageFilter.h" #include "rtkFFTRampImageFilter.h" - -#include +#include "rtkExtractImageSubRegion.h" namespace rtk { @@ -38,8 +37,8 @@ namespace rtk * - rtk::FFTRampImageFilter for ramp filtering, * - rtk::FDKBackProjectionImageFilter for backprojection. * The input stack of projections is processed piece by piece (the size is - * controlled with ProjectionSubsetSize) via the use of itk::ExtractImageFilter - * to extract sub-stacks. + * controlled with ProjectionSubsetSize) by extracting sub-stacks directly + * from the input buffer pointer (zero-copy). * * \dot * digraph FDKConeBeamReconstructionFilter { @@ -76,7 +75,6 @@ class ITK_TEMPLATE_EXPORT FDKConeBeamReconstructionFilter : public itk::InPlaceI using OutputImageType = TOutputImage; /** Typedefs of each subfilter of this composite filter */ - using ExtractFilterType = itk::ExtractImageFilter; using WeightFilterType = rtk::FDKWeightProjectionFilter; using RampFilterType = rtk::FFTRampImageFilter; using BackProjectionFilterType = rtk::FDKBackProjectionImageFilter; @@ -142,10 +140,9 @@ class ITK_TEMPLATE_EXPORT FDKConeBeamReconstructionFilter : public itk::InPlaceI {} /** Pointers to each subfilter of this composite filter */ - typename ExtractFilterType::Pointer m_ExtractFilter; - typename WeightFilterType::Pointer m_WeightFilter; - typename RampFilterType::Pointer m_RampFilter; - BackProjectionFilterPointer m_BackProjectionFilter; + typename WeightFilterType::Pointer m_WeightFilter; + typename RampFilterType::Pointer m_RampFilter; + BackProjectionFilterPointer m_BackProjectionFilter; private: /** Number of projections processed at a time. */ diff --git a/include/rtkFDKConeBeamReconstructionFilter.hxx b/include/rtkFDKConeBeamReconstructionFilter.hxx index aac8d6d25..06be4e7a5 100644 --- a/include/rtkFDKConeBeamReconstructionFilter.hxx +++ b/include/rtkFDKConeBeamReconstructionFilter.hxx @@ -31,17 +31,14 @@ FDKConeBeamReconstructionFilter::FDKCo this->SetNumberOfRequiredInputs(2); // Create each filter of the composite filter - m_ExtractFilter = ExtractFilterType::New(); m_WeightFilter = WeightFilterType::New(); m_RampFilter = RampFilterType::New(); this->SetBackProjectionFilter(BackProjectionFilterType::New()); // Permanent internal connections - m_WeightFilter->SetInput(m_ExtractFilter->GetOutput()); m_RampFilter->SetInput(m_WeightFilter->GetOutput()); // Default parameters - m_ExtractFilter->SetDirectionCollapseToSubmatrix(); m_WeightFilter->InPlaceOn(); // Default to one projection per subset when FFTW is not available @@ -76,9 +73,11 @@ FDKConeBeamReconstructionFilter::Gener // SR: is this useful? m_BackProjectionFilter->SetInput(0, this->GetInput(0)); m_BackProjectionFilter->SetInPlace(this->GetInPlace()); - m_ExtractFilter->SetInput(this->GetInput(1)); m_BackProjectionFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); + + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + inputPtr1->SetRequestedRegion(this->GetInput(1)->GetLargestPossibleRegion()); } template @@ -87,21 +86,29 @@ FDKConeBeamReconstructionFilter::Gener { const unsigned int Dimension = this->InputImageDimension; + // Trigger upstream update, filters like DisplacedDetector can change regions. + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + inputPtr1->UpdateOutputInformation(); + m_WeightFilter->SetGeometry(m_Geometry); m_BackProjectionFilter->SetGeometry(m_Geometry); // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function - typename ExtractFilterType::InputImageRegionType projRegion; + typename InputImageType::RegionType projRegion; projRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int firstStackSize = std::min(m_ProjectionSubsetSize, (unsigned int)projRegion.GetSize(Dimension - 1)); projRegion.SetSize(Dimension - 1, firstStackSize); - m_ExtractFilter->SetExtractionRegion(projRegion); + + // Create a zero-copy view of the first sub-stack + typename InputImageType::Pointer subStack = rtk::ExtractImageSubRegion(this->GetInput(1), projRegion); + m_WeightFilter->SetInput(subStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion)) + m_WeightFilter->InPlaceOff(); // Run composite filter update m_BackProjectionFilter->SetInput(0, this->GetInput(0)); m_BackProjectionFilter->SetInPlace(this->GetInPlace()); - m_ExtractFilter->SetInput(this->GetInput(1)); m_BackProjectionFilter->UpdateOutputInformation(); // Update output information @@ -117,13 +124,14 @@ FDKConeBeamReconstructionFilter::Gener { const unsigned int Dimension = this->InputImageDimension; - // The backprojection works on a small stack of projections, not the full stack - typename ExtractFilterType::InputImageRegionType subsetRegion; + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + inputPtr1->Update(); + + typename InputImageType::RegionType subsetRegion; subsetRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int nProj = subsetRegion.GetSize(Dimension - 1); + unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1); - // The progress accumulator tracks the progress of the pipeline - // Each filter is equally weighted across all iterations of the stack auto progress = itk::ProgressAccumulator::New(); progress->SetMiniPipelineFilter(this); auto frac = (1.0f / 3) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize); @@ -133,22 +141,23 @@ FDKConeBeamReconstructionFilter::Gener for (unsigned int i = 0; i < nProj; i += m_ProjectionSubsetSize) { - // After the first bp update, we need to use its output as input. + subsetRegion.SetIndex(Dimension - 1, baseIndex + i); + subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i)); + typename InputImageType::Pointer subStack = rtk::ExtractImageSubRegion(this->GetInput(1), subsetRegion); + m_WeightFilter->SetInput(subStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion)) + m_WeightFilter->InPlaceOff(); + if (i) { typename TInputImage::Pointer pimg = m_BackProjectionFilter->GetOutput(); pimg->DisconnectPipeline(); m_BackProjectionFilter->SetInput(pimg); - // Change projection subset - subsetRegion.SetIndex(Dimension - 1, i); - subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i)); - m_ExtractFilter->SetExtractionRegion(subsetRegion); - - // This is required to reset the full pipeline m_BackProjectionFilter->GetOutput()->UpdateOutputInformation(); m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); } + m_BackProjectionFilter->Update(); } diff --git a/include/rtkFDKVarianceReconstructionFilter.h b/include/rtkFDKVarianceReconstructionFilter.h index 4e3b7a8b9..29549db63 100644 --- a/include/rtkFDKVarianceReconstructionFilter.h +++ b/include/rtkFDKVarianceReconstructionFilter.h @@ -23,8 +23,7 @@ #include "rtkConfiguration.h" #include "rtkFDKBackProjectionImageFilter.h" #include "rtkFFTVarianceRampImageFilter.h" - -#include +#include "rtkExtractImageSubRegion.h" namespace rtk { @@ -68,7 +67,6 @@ class ITK_TEMPLATE_EXPORT FDKVarianceReconstructionFilter : public itk::InPlaceI using OutputImageType = TOutputImage; /** Typedefs of each subfilter of this composite filter */ - using ExtractFilterType = itk::ExtractImageFilter; using WeightFilterType = rtk::FDKWeightProjectionFilter; using VarianceRampFilterType = rtk::FFTVarianceRampImageFilter; using BackProjectionFilterType = rtk::FDKBackProjectionImageFilter; @@ -127,7 +125,6 @@ class ITK_TEMPLATE_EXPORT FDKVarianceReconstructionFilter : public itk::InPlaceI {} /** Pointers to each subfilter of this composite filter */ - typename ExtractFilterType::Pointer m_ExtractFilter; typename WeightFilterType::Pointer m_WeightFilter1; typename WeightFilterType::Pointer m_WeightFilter2; typename VarianceRampFilterType::Pointer m_VarianceRampFilter; diff --git a/include/rtkFDKVarianceReconstructionFilter.hxx b/include/rtkFDKVarianceReconstructionFilter.hxx index 0dde3f13f..9f95a75ec 100644 --- a/include/rtkFDKVarianceReconstructionFilter.hxx +++ b/include/rtkFDKVarianceReconstructionFilter.hxx @@ -32,19 +32,16 @@ FDKVarianceReconstructionFilter::FDKVa this->SetNumberOfRequiredInputs(2); // Create each filter of the composite filter - m_ExtractFilter = ExtractFilterType::New(); m_WeightFilter1 = WeightFilterType::New(); m_WeightFilter2 = WeightFilterType::New(); m_VarianceRampFilter = VarianceRampFilterType::New(); this->SetBackProjectionFilter(BackProjectionFilterType::New()); // Permanent internal connections - m_WeightFilter1->SetInput(m_ExtractFilter->GetOutput()); m_WeightFilter2->SetInput(m_WeightFilter1->GetOutput()); m_VarianceRampFilter->SetInput(m_WeightFilter2->GetOutput()); // Default parameters - m_ExtractFilter->SetDirectionCollapseToSubmatrix(); m_WeightFilter1->InPlaceOn(); m_WeightFilter2->InPlaceOn(); @@ -80,9 +77,11 @@ FDKVarianceReconstructionFilter::Gener // SR: is this useful? m_BackProjectionFilter->SetInput(0, this->GetInput(0)); m_BackProjectionFilter->SetInPlace(this->GetInPlace()); - m_ExtractFilter->SetInput(this->GetInput(1)); m_BackProjectionFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); + + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + inputPtr1->SetRequestedRegion(this->GetInput(1)->GetLargestPossibleRegion()); } template @@ -91,22 +90,29 @@ FDKVarianceReconstructionFilter::Gener { const unsigned int Dimension = this->InputImageDimension; + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + inputPtr1->UpdateOutputInformation(); + m_WeightFilter1->SetGeometry(m_Geometry); m_WeightFilter2->SetGeometry(m_Geometry); m_BackProjectionFilter->SetGeometry(m_Geometry); // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function - typename ExtractFilterType::InputImageRegionType projRegion; + typename InputImageType::RegionType projRegion; projRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int firstStackSize = std::min(m_ProjectionSubsetSize, (unsigned int)projRegion.GetSize(Dimension - 1)); projRegion.SetSize(Dimension - 1, firstStackSize); - m_ExtractFilter->SetExtractionRegion(projRegion); + + // Create a zero-copy view of the first sub-stack + typename InputImageType::Pointer subStack = rtk::ExtractImageSubRegion(this->GetInput(1), projRegion); + m_WeightFilter1->SetInput(subStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion)) + m_WeightFilter1->InPlaceOff(); // Run composite filter update m_BackProjectionFilter->SetInput(0, this->GetInput(0)); m_BackProjectionFilter->SetInPlace(this->GetInPlace()); - m_ExtractFilter->SetInput(this->GetInput(1)); m_BackProjectionFilter->UpdateOutputInformation(); // Update output information @@ -123,9 +129,10 @@ FDKVarianceReconstructionFilter::Gener const unsigned int Dimension = this->InputImageDimension; // The backprojection works on a small stack of projections, not the full stack - typename ExtractFilterType::InputImageRegionType subsetRegion; + typename InputImageType::RegionType subsetRegion; subsetRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int nProj = subsetRegion.GetSize(Dimension - 1); + unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1); // The progress accumulator tracks the progress of the pipeline // Each filter is equally weighted across all iterations of the stack @@ -146,15 +153,19 @@ FDKVarianceReconstructionFilter::Gener pimg->DisconnectPipeline(); m_BackProjectionFilter->SetInput(pimg); - // Change projection subset - subsetRegion.SetIndex(Dimension - 1, i); - subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i)); - m_ExtractFilter->SetExtractionRegion(subsetRegion); - // This is required to reset the full pipeline m_BackProjectionFilter->GetOutput()->UpdateOutputInformation(); m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); } + + // Always create the substack for the current subset + subsetRegion.SetIndex(Dimension - 1, baseIndex + i); + subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i)); + typename InputImageType::Pointer subStack = rtk::ExtractImageSubRegion(this->GetInput(1), subsetRegion); + m_WeightFilter1->SetInput(subStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion)) + m_WeightFilter1->InPlaceOff(); + m_BackProjectionFilter->Update(); } diff --git a/include/rtkOSEMConeBeamReconstructionFilter.h b/include/rtkOSEMConeBeamReconstructionFilter.h index 55c346d5f..024eff888 100644 --- a/include/rtkOSEMConeBeamReconstructionFilter.h +++ b/include/rtkOSEMConeBeamReconstructionFilter.h @@ -26,10 +26,10 @@ #include #include #include -#include #include #include "rtkConstantImageSource.h" +#include "rtkExtractImageSubRegion.h" #include "rtkIterativeConeBeamReconstructionFilter.h" namespace rtk @@ -132,7 +132,6 @@ class ITK_TEMPLATE_EXPORT OSEMConeBeamReconstructionFilter using ProjectionType = TProjectionImage; /** Typedefs of each subfilter of this composite filter */ - using ExtractFilterType = itk::ExtractImageFilter; using MultiplyFilterType = itk::MultiplyImageFilter; using ForwardProjectionFilterType = rtk::ForwardProjectionImageFilter; using BackProjectionFilterType = rtk::BackProjectionImageFilter; @@ -199,7 +198,6 @@ class ITK_TEMPLATE_EXPORT OSEMConeBeamReconstructionFilter {} /** Pointers to each subfilter of this composite filter */ - typename ExtractFilterType::Pointer m_ExtractFilter; typename ForwardProjectionFilterType::Pointer m_ForwardProjectionFilter; typename MultiplyFilterType::Pointer m_MultiplyFilter; typename BackProjectionFilterType::Pointer m_BackProjectionFilter; diff --git a/include/rtkOSEMConeBeamReconstructionFilter.hxx b/include/rtkOSEMConeBeamReconstructionFilter.hxx index 080d0e8ae..c2a5591a4 100644 --- a/include/rtkOSEMConeBeamReconstructionFilter.hxx +++ b/include/rtkOSEMConeBeamReconstructionFilter.hxx @@ -36,7 +36,6 @@ OSEMConeBeamReconstructionFilter::OSEMConeBeamRe this->SetNumberOfRequiredInputs(2); // Create each filter of the composite filter - m_ExtractFilter = ExtractFilterType::New(); m_MultiplyFilter = MultiplyFilterType::New(); m_ConstantImageSource = ConstantImageSourceType::New(); m_ZeroConstantProjectionStackSource = ConstantProjectionSourceType::New(); @@ -49,11 +48,9 @@ OSEMConeBeamReconstructionFilter::OSEMConeBeamRe m_DivideVolumeFilter = DivideVolumeFilterType::New(); // Permanent internal connections - m_DivideProjectionFilter->SetInput1(m_ExtractFilter->GetOutput()); m_DivideVolumeFilter->SetInput1(m_MultiplyFilter->GetOutput()); // Default parameters - m_ExtractFilter->SetDirectionCollapseToSubmatrix(); } template @@ -86,7 +83,7 @@ OSEMConeBeamReconstructionFilter::GenerateOutput // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function - typename ExtractFilterType::InputImageRegionType projRegion; + typename ProjectionType::RegionType projRegion; // Set forward projection filter m_ForwardProjectionFilter = this->InstantiateForwardProjectionFilter(this->m_CurrentForwardProjectionConfiguration); @@ -97,22 +94,21 @@ OSEMConeBeamReconstructionFilter::GenerateOutput this->InstantiateBackProjectionFilter(this->m_CurrentBackProjectionConfiguration); projRegion = this->GetInput(1)->GetLargestPossibleRegion(); - m_ExtractFilter->SetExtractionRegion(projRegion); - - m_ExtractFilter->SetInput(this->GetInput(1)); - m_ExtractFilter->UpdateOutputInformation(); // Links with the forward and back projection filters should be set here // and not in the constructor, as these filters are set at runtime m_ConstantImageSource->SetInformationFromImage(const_cast(this->GetInput(0))); m_ConstantImageSource->SetConstant(0); - m_OneConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + // Create a zero-copy sub-region view for metadata setup + typename ProjectionType::Pointer projSubStack = rtk::ExtractImageSubRegion(this->GetInput(1), projRegion); + m_DivideProjectionFilter->SetInput1(projSubStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion)) + m_DivideProjectionFilter->InPlaceOff(); + m_OneConstantProjectionStackSource->SetInformationFromImage(projSubStack.GetPointer()); m_OneConstantProjectionStackSource->SetConstant(1); - m_ZeroConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + m_ZeroConstantProjectionStackSource->SetInformationFromImage(projSubStack.GetPointer()); m_ZeroConstantProjectionStackSource->SetConstant(0); m_BackProjectionFilter->SetInput(0, m_ConstantImageSource->GetOutput()); @@ -163,9 +159,10 @@ OSEMConeBeamReconstructionFilter::GenerateData() const unsigned int Dimension = this->InputImageDimension; // The backprojection works on one projection at a time - typename ExtractFilterType::InputImageRegionType subsetRegion; + typename ProjectionType::RegionType subsetRegion; subsetRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int nProj = subsetRegion.GetSize(Dimension - 1); + unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1); subsetRegion.SetSize(Dimension - 1, 1); // Fill and shuffle randomly the projection order. @@ -195,12 +192,24 @@ OSEMConeBeamReconstructionFilter::GenerateData() for (unsigned int i = 0; i < nProj; i++) { // Change projection subset - subsetRegion.SetIndex(Dimension - 1, projOrder[i]); - m_ExtractFilter->SetExtractionRegion(subsetRegion); - m_ExtractFilter->UpdateOutputInformation(); + subsetRegion.SetIndex(Dimension - 1, baseIndex + projOrder[i]); + typename ProjectionType::Pointer projSubStack = rtk::ExtractImageSubRegion(this->GetInput(1), subsetRegion); + + // Constant projection for normalization backprojection + typename ConstantProjectionSourceType::Pointer oneNormProj = ConstantProjectionSourceType::New(); + oneNormProj->SetInformationFromImage(projSubStack.GetPointer()); + oneNormProj->SetConstant(1); + + // Constant projection for forward projection input + typename ConstantProjectionSourceType::Pointer zeroProj = ConstantProjectionSourceType::New(); + zeroProj->SetInformationFromImage(projSubStack.GetPointer()); + zeroProj->SetConstant(0); + zeroProj->Update(); - m_ZeroConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + m_DivideProjectionFilter->SetInput1(projSubStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion)) + m_DivideProjectionFilter->InPlaceOff(); + m_ForwardProjectionFilter->SetInput(0, zeroProj->GetOutput()); // This is required to reset the full pipeline m_BackProjectionFilter->GetOutput()->UpdateOutputInformation(); @@ -209,8 +218,8 @@ OSEMConeBeamReconstructionFilter::GenerateData() m_BackProjectionFilter->Update(); if (iter == 0 || !m_StoreNormalizationImages) { - m_OneConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + oneNormProj->SetInformationFromImage(projSubStack.GetPointer()); + m_BackProjectionNormalizationFilter->SetInput(1, oneNormProj->GetOutput()); m_BackProjectionNormalizationFilter->GetOutput()->UpdateOutputInformation(); m_BackProjectionNormalizationFilter->GetOutput()->PropagateRequestedRegion(); m_BackProjectionNormalizationFilter->Update(); diff --git a/include/rtkSARTConeBeamReconstructionFilter.h b/include/rtkSARTConeBeamReconstructionFilter.h index 5ca4b9822..0795fcf16 100644 --- a/include/rtkSARTConeBeamReconstructionFilter.h +++ b/include/rtkSARTConeBeamReconstructionFilter.h @@ -30,11 +30,12 @@ #include #include #include -#include #include #include #include +#include "rtkExtractImageSubRegion.h" + namespace rtk { @@ -48,8 +49,8 @@ namespace rtk * - SubtractImageFilter, * - BackProjectionImageFilter. * The input stack of projections is processed piece by piece (the size is - * controlled with ProjectionSubsetSize) via the use of itk::ExtractImageFilter - * to extract sub-stacks. + * controlled with ProjectionSubsetSize) by extracting sub-stacks directly + * from the input buffer pointer (zero-copy). * * Two weighting steps must be applied when processing a given projection: * - each pixel of the forward projection must be divided by the total length of the @@ -156,7 +157,6 @@ class ITK_TEMPLATE_EXPORT SARTConeBeamReconstructionFilter using ProjectionPixelType = typename ProjectionType::PixelType; /** Typedefs of each subfilter of this composite filter */ - using ExtractFilterType = itk::ExtractImageFilter; using MultiplyFilterType = itk::MultiplyImageFilter; using ForwardProjectionFilterType = rtk::ForwardProjectionImageFilter; using SubtractFilterType = itk::SubtractImageFilter; @@ -246,8 +246,6 @@ class ITK_TEMPLATE_EXPORT SARTConeBeamReconstructionFilter {} /** Pointers to each subfilter of this composite filter */ - typename ExtractFilterType::Pointer m_ExtractFilter; - typename ExtractFilterType::Pointer m_ExtractFilterRayBox; typename MultiplyFilterType::Pointer m_ZeroMultiplyFilter; typename ForwardProjectionFilterType::Pointer m_ForwardProjectionFilter; typename SubtractFilterType::Pointer m_SubtractFilter; diff --git a/include/rtkSARTConeBeamReconstructionFilter.hxx b/include/rtkSARTConeBeamReconstructionFilter.hxx index f28c30836..c25364901 100644 --- a/include/rtkSARTConeBeamReconstructionFilter.hxx +++ b/include/rtkSARTConeBeamReconstructionFilter.hxx @@ -37,7 +37,6 @@ SARTConeBeamReconstructionFilter::SARTConeBeamRe m_Lambda = 0.3; // Create each filter of the composite filter - m_ExtractFilter = ExtractFilterType::New(); m_ZeroMultiplyFilter = MultiplyFilterType::New(); m_SubtractFilter = SubtractFilterType::New(); m_AddFilter = AddFilterType::New(); @@ -49,7 +48,6 @@ SARTConeBeamReconstructionFilter::SARTConeBeamRe // Create the filters required for correct weighting of the difference // projection - m_ExtractFilterRayBox = ExtractFilterType::New(); m_RayBoxFilter = RayBoxIntersectionFilterType::New(); m_DivideProjectionFilter = DivideProjectionFilterType::New(); m_ConstantProjectionStackSource = ConstantProjectionSourceType::New(); @@ -64,23 +62,18 @@ SARTConeBeamReconstructionFilter::SARTConeBeamRe m_DivisionThreshold = m_DivideVolumeFilter->GetThreshold(); // Permanent internal connections + m_SubtractFilter->InPlaceOn(); m_ZeroMultiplyFilter->SetInput1(itk::NumericTraits::ZeroValue()); - m_ZeroMultiplyFilter->SetInput2(m_ExtractFilter->GetOutput()); - - m_SubtractFilter->SetInput(0, m_ExtractFilter->GetOutput()); m_MultiplyFilter->SetInput1(m_Lambda); m_MultiplyFilter->SetInput2(m_SubtractFilter->GetOutput()); - m_ExtractFilterRayBox->SetInput(m_ConstantProjectionStackSource->GetOutput()); - m_RayBoxFilter->SetInput(m_ExtractFilterRayBox->GetOutput()); + m_RayBoxFilter->SetInput(m_ConstantProjectionStackSource->GetOutput()); m_DivideProjectionFilter->SetInput1(m_MultiplyFilter->GetOutput()); m_DivideProjectionFilter->SetInput2(m_RayBoxFilter->GetOutput()); m_DisplacedDetectorFilter->SetInput(m_DivideProjectionFilter->GetOutput()); // Default parameters - m_ExtractFilter->SetDirectionCollapseToSubmatrix(); - m_ExtractFilterRayBox->SetDirectionCollapseToSubmatrix(); m_DisplacedDetectorFilter->SetPadOnTruncatedSide(false); } @@ -132,11 +125,19 @@ SARTConeBeamReconstructionFilter::GenerateOutput // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function - typename ExtractFilterType::InputImageRegionType projRegion; + typename ProjectionType::RegionType projRegion; projRegion = this->GetInput(1)->GetLargestPossibleRegion(); - m_ExtractFilter->SetExtractionRegion(projRegion); - m_ExtractFilterRayBox->SetExtractionRegion(projRegion); + + // Create zero-copy view of the full projection stack + typename ProjectionType::Pointer projSubStack = rtk::ExtractImageSubRegion(this->GetInput(1), projRegion); + m_ZeroMultiplyFilter->SetInput2(projSubStack); + m_SubtractFilter->SetInput(0, projSubStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion)) + { + m_ZeroMultiplyFilter->InPlaceOff(); + m_SubtractFilter->InPlaceOff(); + } // Set forward projection filter m_ForwardProjectionFilter = this->InstantiateForwardProjectionFilter(this->m_CurrentForwardProjectionConfiguration); @@ -152,8 +153,7 @@ SARTConeBeamReconstructionFilter::GenerateOutput m_ConstantImageSource->SetConstant(0); m_ConstantImageSource->UpdateOutputInformation(); - m_OneConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + m_OneConstantProjectionStackSource->SetInformationFromImage(projSubStack.GetPointer()); m_OneConstantProjectionStackSource->SetConstant(1); m_BackProjectionFilter->SetInput(0, m_ConstantImageSource->GetOutput()); @@ -186,7 +186,6 @@ SARTConeBeamReconstructionFilter::GenerateOutput m_ForwardProjectionFilter->SetInput(0, m_ZeroMultiplyFilter->GetOutput()); m_ForwardProjectionFilter->SetInput(1, this->GetInput(0)); - m_ExtractFilter->SetInput(this->GetInput(1)); m_SubtractFilter->SetInput(1, m_ForwardProjectionFilter->GetOutput()); m_ForwardProjectionFilter->SetGeometry(this->m_Geometry); @@ -250,9 +249,10 @@ SARTConeBeamReconstructionFilter::GenerateData() const unsigned int Dimension = this->InputImageDimension; // The backprojection works on one projection at a time - typename ExtractFilterType::InputImageRegionType subsetRegion; + typename ProjectionType::RegionType subsetRegion; subsetRegion = this->GetInput(1)->GetLargestPossibleRegion(); unsigned int nProj = subsetRegion.GetSize(Dimension - 1); + unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1); subsetRegion.SetSize(Dimension - 1, 1); // Fill and shuffle randomly the projection order. @@ -278,14 +278,22 @@ SARTConeBeamReconstructionFilter::GenerateData() unsigned int projectionsProcessedInSubset = 0; for (unsigned int i = 0; i < nProj; i++) { - // Change projection subset - subsetRegion.SetIndex(Dimension - 1, projOrder[i]); - m_ExtractFilter->SetExtractionRegion(subsetRegion); - m_ExtractFilterRayBox->SetExtractionRegion(subsetRegion); - m_ExtractFilter->UpdateOutputInformation(); - - m_OneConstantProjectionStackSource->SetInformationFromImage( - const_cast(m_ExtractFilter->GetOutput())); + // Change projection subset and create zero-copy view for forward projection path + subsetRegion.SetIndex(Dimension - 1, baseIndex + projOrder[i]); + typename ProjectionType::Pointer projSubStack = rtk::ExtractImageSubRegion(this->GetInput(1), subsetRegion); + m_ZeroMultiplyFilter->SetInput2(projSubStack); + m_SubtractFilter->SetInput(0, projSubStack); + if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion)) + { + m_ZeroMultiplyFilter->InPlaceOff(); + m_SubtractFilter->InPlaceOff(); + } + + // Update constant sources with current projection's metadata + m_OneConstantProjectionStackSource->SetInformationFromImage(projSubStack.GetPointer()); + m_OneConstantProjectionStackSource->SetConstant(1); + m_ConstantProjectionStackSource->SetInformationFromImage(projSubStack.GetPointer()); + m_ConstantProjectionStackSource->SetConstant(0); // Set gating weight for the current projection if (m_IsGated) diff --git a/src/rtkCudaFDKConeBeamReconstructionFilter.cxx b/src/rtkCudaFDKConeBeamReconstructionFilter.cxx index aef8e2bcc..02f47aa09 100644 --- a/src/rtkCudaFDKConeBeamReconstructionFilter.cxx +++ b/src/rtkCudaFDKConeBeamReconstructionFilter.cxx @@ -29,7 +29,6 @@ CudaFDKConeBeamReconstructionFilter ::CudaFDKConeBeamReconstructionFilter() m_BackProjectionFilter = BackProjectionFilterType::New(); // Permanent internal connections - m_WeightFilter->SetInput(m_ExtractFilter->GetOutput()); m_RampFilter->SetInput(m_WeightFilter->GetOutput()); m_BackProjectionFilter->SetInput(1, m_RampFilter->GetOutput()); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 18da9b27c..ef31cc040 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -113,6 +113,9 @@ endif() rtk_add_test(rtkFDKTest rtkfdktest.cxx) rtk_add_cuda_test(rtkFDKCudaTest rtkfdktest.cxx) +rtk_add_test(rtkExtractImageSubRegionTest rtkextractimagesubregiontest.cxx) +rtk_add_cuda_test(rtkExtractImageSubRegionCudaTest rtkextractimagesubregiontest.cxx) + rtk_add_cuda_test(rtkFDKProjWeightCompCudaTest rtkfdkprojweightcompcudatest.cxx) rtk_add_test(rtkFBPParallelTest rtkfbpparalleltest.cxx) diff --git a/test/rtkextractimagesubregiontest.cxx b/test/rtkextractimagesubregiontest.cxx new file mode 100644 index 000000000..a981ebccc --- /dev/null +++ b/test/rtkextractimagesubregiontest.cxx @@ -0,0 +1,190 @@ +#include + +#include "rtkConstantImageSource.h" +#include "rtkExtractImageSubRegion.h" +#include "rtkTest.h" + +#ifdef USE_CUDA +# include +#endif + +/** + * \file rtkextractimagesubregiontest.cxx + * + * \brief Functional test for rtk::ExtractImageSubRegion + * + * This test verifies that rtk::ExtractImageSubRegion produces a zero-copy + * view when the requested region is contiguous along the last dimension, + * and falls back to a copy when it is not. + * + * \author Axel Garcia + */ + +int +rtkextractimagesubregiontest(int, char *[]) +{ + constexpr unsigned int Dimension = 3; + using PixelType = float; +#ifdef USE_CUDA + using ImageType = itk::CudaImage; +#else + using ImageType = itk::Image; +#endif + using RegionType = itk::ImageRegion; + using IndexType = itk::Index; + using SizeType = itk::Size; + + // Create a 4x5x6 image with known pixel values + auto source = rtk::ConstantImageSource::New(); + source->SetOrigin(itk::MakePoint(0., 0., 0.)); + source->SetSpacing(itk::MakeVector(1., 1., 1.)); + source->SetSize(itk::MakeSize(4, 5, 6)); + source->SetConstant(3.14f); + TRY_AND_EXIT_ON_ITK_EXCEPTION(source->UpdateLargestPossibleRegion()); + + ImageType::Pointer input = source->GetOutput(); + + // Fill with a gradient so each pixel is unique: value = x + y*10 + z*100 + itk::ImageRegionIterator it(input, input->GetLargestPossibleRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + IndexType idx = it.GetIndex(); + it.Set(static_cast(idx[0] + idx[1] * 10 + idx[2] * 100)); + } + + // ===== Case 1: Contiguous extraction (dims 0,1 span full input) ===== + std::cout << "\n\n****** Case 1: contiguous extraction (zero-copy) ******" << std::endl; + + RegionType contiguousRegion; + contiguousRegion.SetIndex(itk::MakeIndex(0, 0, 2)); + contiguousRegion.SetSize(itk::MakeSize(4, 5, 3)); + + ImageType::Pointer subRegion = rtk::ExtractImageSubRegion(input.GetPointer(), contiguousRegion); + + // Verify metadata + if (subRegion->GetLargestPossibleRegion() != contiguousRegion) + { + std::cerr << "Region mismatch!" << std::endl; + return EXIT_FAILURE; + } + if (subRegion->GetSpacing() != input->GetSpacing()) + { + std::cerr << "Spacing mismatch!" << std::endl; + return EXIT_FAILURE; + } + + // Verify zero-copy: buffer pointer should point into the input's buffer + const PixelType * inputBuf = input->GetBufferPointer(); + const PixelType * outputBuf = subRegion->GetBufferPointer(); + ptrdiff_t offset = outputBuf - inputBuf; + // For a contiguous extraction starting at z=2 with slice size 4*5=20: + // offset should be 2*20 = 40 + ptrdiff_t expectedOffset = 2 * 4 * 5; + if (offset != expectedOffset) + { + std::cerr << "Zero-copy FAILED: expected offset " << expectedOffset << ", got " << offset << std::endl; + return EXIT_FAILURE; + } + std::cout << "Zero-copy: buffer offset = " << offset << " (correct)" << std::endl; + + // Verify pixel values are accessible and correct + itk::ImageRegionIterator outIt(subRegion, subRegion->GetLargestPossibleRegion()); + bool valuesCorrect = true; + for (outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt) + { + IndexType idx = outIt.GetIndex(); + PixelType expected = static_cast(idx[0] + idx[1] * 10 + idx[2] * 100); + if (std::abs(outIt.Get() - expected) > 1e-6f) + { + std::cerr << "Value mismatch at " << idx << ": got " << outIt.Get() << ", expected " << expected << std::endl; + valuesCorrect = false; + break; + } + } + if (!valuesCorrect) + return EXIT_FAILURE; + std::cout << "Pixel values: correct" << std::endl; + + // ===== Case 2: Non-contiguous extraction (middle slice) ===== + std::cout << "\n\n****** Case 2: non-contiguous extraction (copy) ******" << std::endl; + + RegionType nonContiguousRegion; + nonContiguousRegion.SetIndex(itk::MakeIndex(1, 1, 3)); + nonContiguousRegion.SetSize(itk::MakeSize(2, 2, 2)); + + ImageType::Pointer subRegion2 = rtk::ExtractImageSubRegion(input.GetPointer(), nonContiguousRegion); + + // Verify metadata + if (subRegion2->GetLargestPossibleRegion() != nonContiguousRegion) + { + std::cerr << "Region mismatch!" << std::endl; + return EXIT_FAILURE; + } + + // Verify zero-copy does NOT apply: buffer pointers should differ + const PixelType * outputBuf2 = subRegion2->GetBufferPointer(); + if (outputBuf2 >= inputBuf && outputBuf2 < inputBuf + input->GetLargestPossibleRegion().GetNumberOfPixels()) + { + std::cerr << "Non-contiguous extraction should NOT share the input buffer!" << std::endl; + return EXIT_FAILURE; + } + std::cout << "Buffer: independent copy (correct)" << std::endl; + + // Verify pixel values via reference image + auto reference = rtk::ConstantImageSource::New(); + reference->SetOrigin(itk::MakePoint(0., 0., 0.)); + reference->SetSpacing(itk::MakeVector(1., 1., 1.)); + reference->SetSize(itk::MakeSize(2, 2, 2)); + reference->SetConstant(0.f); + reference->UpdateLargestPossibleRegion(); + + // Fill reference with expected values + itk::ImageRegionIterator refIt(reference->GetOutput(), reference->GetOutput()->GetLargestPossibleRegion()); + for (refIt.GoToBegin(); !refIt.IsAtEnd(); ++refIt) + { + IndexType idx = refIt.GetIndex(); + // Map back to input coordinates + PixelType expected = static_cast((idx[0] + 1) + (idx[1] + 1) * 10 + (idx[2] + 3) * 100); + refIt.Set(expected); + } + + CheckImageQuality(subRegion2, reference->GetOutput(), 0.001, 120, 432.f); + std::cout << "Pixel values: correct" << std::endl; + + // ===== Case 3: No buffer yet (GenerateOutputInformation scenario) ===== + std::cout << "\n\n****** Case 3: metadata-only (no buffer) ******" << std::endl; + + auto noBufferSource = rtk::ConstantImageSource::New(); + noBufferSource->SetOrigin(itk::MakePoint(0., 0., 0.)); + noBufferSource->SetSpacing(itk::MakeVector(1., 1., 1.)); + noBufferSource->SetSize(itk::MakeSize(4, 5, 6)); + noBufferSource->SetConstant(0.f); + // Don't call Update — buffer is not allocated + noBufferSource->UpdateOutputInformation(); + + RegionType metaRegion; + metaRegion.SetIndex(itk::MakeIndex(0, 0, 1)); + metaRegion.SetSize(itk::MakeSize(4, 5, 2)); + + // Should not crash — returns metadata-only image + ImageType::Pointer metaSub = rtk::ExtractImageSubRegion(noBufferSource->GetOutput(), metaRegion); + if (metaSub->GetLargestPossibleRegion() != metaRegion) + { + std::cerr << "Metadata region mismatch!" << std::endl; + return EXIT_FAILURE; + } +#ifndef USE_CUDA + // For CPU images, buffer should be null when input has no data. + // For CudaImage, GetBufferPointer() has side effects on CudaDataManager + // that prevent a clean nullptr check, so we skip this for CUDA. + if (metaSub->GetBufferPointer() != nullptr) + { + std::cerr << "Expected null buffer for metadata-only image!" << std::endl; + return EXIT_FAILURE; + } +#endif + std::cout << "Metadata-only: correct" << std::endl; + + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; +}