From 06a8cfe044f91b5c2b0e0c660d13844c2540e961 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 16 Sep 2020 10:27:47 +0200 Subject: [PATCH 1/2] WIP: refactor RTK tests to follow ITK standards --- test/rtkfdktest.cxx | 225 -------------------------------------------- 1 file changed, 225 deletions(-) delete mode 100644 test/rtkfdktest.cxx diff --git a/test/rtkfdktest.cxx b/test/rtkfdktest.cxx deleted file mode 100644 index 8d987959f..000000000 --- a/test/rtkfdktest.cxx +++ /dev/null @@ -1,225 +0,0 @@ -#include -#include -#include - -#include "rtkTest.h" -#include "rtkSheppLoganPhantomFilter.h" -#include "rtkDrawSheppLoganFilter.h" -#include "rtkConstantImageSource.h" -#include "rtkFieldOfViewImageFilter.h" - -#ifdef USE_CUDA -# include "rtkCudaFDKConeBeamReconstructionFilter.h" -#else -# include "rtkFDKConeBeamReconstructionFilter.h" -#endif - -/** - * \file rtkfdktest.cxx - * - * \brief Functional test for classes performing FDK reconstructions - * - * This test generates the projections of a simulated Shepp-Logan phantom. - * A CT image is reconstructed from each set of generated projection images - * using the FDK algorithm and the reconstructed CT image is compared to the - * expected results which is analytically computed. - * - * \author Simon Rit and Marc Vila - */ - -int -main(int, char **) -{ - constexpr unsigned int Dimension = 3; - using OutputPixelType = float; -#ifdef USE_CUDA - using OutputImageType = itk::CudaImage; -#else - using OutputImageType = itk::Image; -#endif - -#if FAST_TESTS_NO_CHECKS - constexpr unsigned int NumberOfProjectionImages = 3; -#else - constexpr unsigned int NumberOfProjectionImages = 180; -#endif - - // Constant image sources - using ConstantImageSourceType = rtk::ConstantImageSource; - ConstantImageSourceType::PointType origin; - ConstantImageSourceType::SizeType size; - ConstantImageSourceType::SpacingType spacing; - - ConstantImageSourceType::Pointer tomographySource = ConstantImageSourceType::New(); - origin[0] = -127.; - origin[1] = -127.; - origin[2] = -127.; -#if FAST_TESTS_NO_CHECKS - size[0] = 32; - size[1] = 32; - size[2] = 32; - spacing[0] = 8.; - spacing[1] = 8.; - spacing[2] = 8.; -#else - size[0] = 128; - size[1] = 128; - size[2] = 128; - spacing[0] = 2.; - spacing[1] = 2.; - spacing[2] = 2.; -#endif - tomographySource->SetOrigin(origin); - tomographySource->SetSpacing(spacing); - tomographySource->SetSize(size); - tomographySource->SetConstant(0.); - - ConstantImageSourceType::Pointer projectionsSource = ConstantImageSourceType::New(); - origin[0] = -254.; - origin[1] = -254.; - origin[2] = -254.; -#if FAST_TESTS_NO_CHECKS - size[0] = 32; - size[1] = 32; - size[2] = NumberOfProjectionImages; - spacing[0] = 32.; - spacing[1] = 32.; - spacing[2] = 32.; -#else - size[0] = 128; - size[1] = 128; - size[2] = NumberOfProjectionImages; - spacing[0] = 4.; - spacing[1] = 4.; - spacing[2] = 4.; -#endif - projectionsSource->SetOrigin(origin); - projectionsSource->SetSpacing(spacing); - projectionsSource->SetSize(size); - projectionsSource->SetConstant(0.); - - std::cout << "\n\n****** Case 1: No streaming ******" << std::endl; - - // Geometry object - using GeometryType = rtk::ThreeDCircularProjectionGeometry; - GeometryType::Pointer geometry = GeometryType::New(); - for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++) - geometry->AddProjection(600., 1200., noProj * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); - - // Shepp Logan projections filter - using SLPType = rtk::SheppLoganPhantomFilter; - SLPType::Pointer slp = SLPType::New(); - slp->SetInput(projectionsSource->GetOutput()); - slp->SetGeometry(geometry); - slp->SetPhantomScale(116); - TRY_AND_EXIT_ON_ITK_EXCEPTION(slp->Update()); - - // Create a reference object (in this case a 3D phantom reference). - using DSLType = rtk::DrawSheppLoganFilter; - DSLType::Pointer dsl = DSLType::New(); - dsl->SetInput(tomographySource->GetOutput()); - dsl->SetPhantomScale(116); - TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) - - // FDK reconstruction filtering -#ifdef USE_CUDA - using FDKType = rtk::CudaFDKConeBeamReconstructionFilter; -#else - using FDKType = rtk::FDKConeBeamReconstructionFilter; -#endif - FDKType::Pointer feldkamp = FDKType::New(); - feldkamp->SetInput(0, tomographySource->GetOutput()); - feldkamp->SetInput(1, slp->GetOutput()); - feldkamp->SetGeometry(geometry); - TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update()); - - - // FOV - using FOVFilterType = rtk::FieldOfViewImageFilter; - FOVFilterType::Pointer fov = FOVFilterType::New(); - fov->SetInput(0, feldkamp->GetOutput()); - fov->SetProjectionsStack(slp->GetOutput()); - fov->SetGeometry(geometry); - TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); - - CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); - std::cout << "Test PASSED! " << std::endl; - - std::cout << "\n\n****** Case 2: perpendicular direction ******" << std::endl; - - ConstantImageSourceType::OutputImageType::DirectionType direction; - direction[0][0] = 0; - direction[0][1] = 1; - direction[0][2] = 0; - direction[1][0] = -1; - direction[1][1] = 0; - direction[1][2] = 0; - direction[2][0] = 0; - direction[2][1] = 0; - direction[2][2] = 1; - tomographySource->SetDirection(direction); - origin[0] = -127.; - origin[1] = 127.; - origin[2] = -127.; - tomographySource->SetOrigin(origin); - fov->GetOutput()->ResetPipeline(); - TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); - TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) - - CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); - std::cout << "Test PASSED! " << std::endl; - - std::cout << "\n\n****** Case 3: 45 degree tilt direction ******" << std::endl; - - direction[0][0] = 0.70710678118; - direction[0][1] = -0.70710678118; - direction[0][2] = 0.70710678118; - direction[1][0] = 0.70710678118; - direction[1][1] = 0; - direction[1][2] = 0; - direction[2][0] = 0; - direction[2][1] = 0; - direction[2][2] = 1; - tomographySource->SetDirection(direction); - origin[0] = -127.; - origin[1] = -127.; - origin[2] = -127.; - tomographySource->SetOrigin(origin); - TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) - TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); - - CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); - std::cout << "Test PASSED! " << std::endl; - - std::cout << "\n\n****** Case 4: streaming ******" << std::endl; - - // Make sure that the data will be recomputed by releasing them - fov->GetOutput()->ReleaseData(); - - using StreamingType = itk::StreamingImageFilter; - StreamingType::Pointer streamer = StreamingType::New(); - streamer->SetInput(0, fov->GetOutput()); - streamer->SetNumberOfStreamDivisions(8); - itk::ImageRegionSplitterDirection::Pointer splitter = itk::ImageRegionSplitterDirection::New(); - splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis - streamer->SetRegionSplitter(splitter); - TRY_AND_EXIT_ON_ITK_EXCEPTION(streamer->Update()); - - CheckImageQuality(streamer->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); - std::cout << "Test PASSED! " << std::endl; - - std::cout << "\n\n****** Case 5: small ROI ******" << std::endl; - origin[0] = -5.; - origin[1] = -13.; - origin[2] = -20.; - size[0] = 64; - size[1] = 64; - size[2] = 64; - tomographySource->SetOrigin(origin); - tomographySource->SetSize(size); - TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->UpdateLargestPossibleRegion()); - TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->UpdateLargestPossibleRegion()) - CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); - std::cout << "Test PASSED! " << std::endl; - return EXIT_SUCCESS; -} From 4d7aefed140c5644377bc61b2fc00e66f2748700 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Thu, 4 Jun 2020 00:10:32 +0200 Subject: [PATCH 2/2] ENH: add backprojection for super short scan --- ...kSuperShortScanBackProjectionImageFilter.h | 92 ++++++ ...uperShortScanBackProjectionImageFilter.hxx | 310 ++++++++++++++++++ test/rtkfdktest.cxx | 225 +++++++++++++ ...perShortScanBackProjectionImageFilter.wrap | 3 + 4 files changed, 630 insertions(+) create mode 100644 include/rtkSuperShortScanBackProjectionImageFilter.h create mode 100644 include/rtkSuperShortScanBackProjectionImageFilter.hxx create mode 100644 test/rtkfdktest.cxx create mode 100644 wrapping/rtkSuperShortScanBackProjectionImageFilter.wrap diff --git a/include/rtkSuperShortScanBackProjectionImageFilter.h b/include/rtkSuperShortScanBackProjectionImageFilter.h new file mode 100644 index 000000000..52e76f5d7 --- /dev/null +++ b/include/rtkSuperShortScanBackProjectionImageFilter.h @@ -0,0 +1,92 @@ +/*========================================================================= + * + * 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 + * + * http://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 rtkSuperShortScanBackProjectionImageFilter_h +#define rtkSuperShortScanBackProjectionImageFilter_h + +#include "rtkBackProjectionImageFilter.h" + +namespace rtk +{ + +/** \class SuperShortScanBackProjectionImageFilter + * \brief CPU version of the backprojection of the SuperShortScan algorithm. + * + * CPU implementation of the backprojection step of the + * [Feldkamp, Davis, Kress, 1984] algorithm for filtered backprojection + * reconstruction of cone-beam CT images with a circular source trajectory. + * + * \author Simon Rit + * + * \ingroup RTK Projector + */ +template +class ITK_EXPORT SuperShortScanBackProjectionImageFilter : public BackProjectionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(SuperShortScanBackProjectionImageFilter); + + /** Standard class type alias. */ + using Self = SuperShortScanBackProjectionImageFilter; + using Superclass = BackProjectionImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + using ProjectionMatrixType = typename Superclass::ProjectionMatrixType; + using OutputImageRegionType = typename TOutputImage::RegionType; + using ProjectionImageType = typename Superclass::ProjectionImageType; + using ProjectionImagePointer = typename ProjectionImageType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(SuperShortScanBackProjectionImageFilter, ImageToImageFilter); + +protected: + SuperShortScanBackProjectionImageFilter() = default; + ~SuperShortScanBackProjectionImageFilter() override = default; + + void + GenerateOutputInformation() override; + + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + + /** Optimized version when the rotation is parallel to X, i.e. matrix[1][0] + and matrix[2][0] are zeros. */ + void + OptimizedBackprojectionX(const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection) override; + + /** Optimized version when the rotation is parallel to Y, i.e. matrix[1][1] + and matrix[2][1] are zeros. */ + void + OptimizedBackprojectionY(const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection) override; +}; + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkSuperShortScanBackProjectionImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkSuperShortScanBackProjectionImageFilter.hxx b/include/rtkSuperShortScanBackProjectionImageFilter.hxx new file mode 100644 index 000000000..f9d9e8b3e --- /dev/null +++ b/include/rtkSuperShortScanBackProjectionImageFilter.hxx @@ -0,0 +1,310 @@ +/*========================================================================= + * + * 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 + * + * http://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 rtkSuperShortScanBackProjectionImageFilter_hxx +#define rtkSuperShortScanBackProjectionImageFilter_hxx + +#include "math.h" + +#include "rtkSuperShortScanBackProjectionImageFilter.h" + +#include +#include + +#define BILINEAR_BACKPROJECTION + +namespace rtk +{ + +template +void +SuperShortScanBackProjectionImageFilter::GenerateOutputInformation() +{ + // Check if detector is cylindrical + if (this->m_Geometry->GetRadiusCylindricalDetector() != 0) + { + itkGenericExceptionMacro( + << "Voxel-based SuperShortScan back projector can currently not handle cylindrical detectors"); + } + + // Run superclass' GenerateOutputInformation + Superclass::GenerateOutputInformation(); +} + +/** + * GenerateData performs the accumulation + */ +template +void +SuperShortScanBackProjectionImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + const unsigned int nProj = this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension - 1); + const unsigned int iFirstProj = this->GetInput(1)->GetLargestPossibleRegion().GetIndex(Dimension - 1); + + // Create interpolator, could be any interpolation + using InterpolatorType = itk::LinearInterpolateImageFunction; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + + // Iterators on volume input and output + using InputRegionIterator = itk::ImageRegionConstIterator; + InputRegionIterator itIn(this->GetInput(), outputRegionForThread); + using OutputRegionIterator = itk::ImageRegionIteratorWithIndex; + OutputRegionIterator itOut(this->GetOutput(), outputRegionForThread); + + // Initialize output region with input region in case the filter is not in + // place + if (this->GetInput() != this->GetOutput()) + { + itIn.GoToBegin(); + while (!itIn.IsAtEnd()) + { + itOut.Set(itIn.Get()); + ++itIn; + ++itOut; + } + } + + // Rotation center (assumed to be at 0 yet) + typename TInputImage::PointType rotCenterPoint; + rotCenterPoint.Fill(0.0); + itk::ContinuousIndex rotCenterIndex; + this->GetInput(0)->TransformPhysicalPointToContinuousIndex(rotCenterPoint, rotCenterIndex); + + // Continuous index at which we interpolate + itk::ContinuousIndex pointProj; + + // Go over each projection + for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) + { + // Extract the current slice + ProjectionImagePointer projection; + projection = this->template GetProjection(iProj); + interpolator->SetInputImage(projection); + + // Index to index matrix normalized to have a correct backprojection weight + // (1 / SID at the isocenter) + ProjectionMatrixType matrix = this->GetIndexToIndexProjectionMatrix(iProj); + double perspFactor = matrix[Dimension - 1][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + perspFactor += matrix[Dimension - 1][j] * rotCenterIndex[j]; + matrix *= this->GetGeometry()->GetSourceToIsocenterDistances()[iProj] / perspFactor; + + // Optimized version + if (fabs(matrix[1][0]) < 1e-10 && fabs(matrix[2][0]) < 1e-10) + { + OptimizedBackprojectionX(outputRegionForThread, matrix, projection); + continue; + } + if (fabs(matrix[1][1]) < 1e-10 && fabs(matrix[2][1]) < 1e-10) + { + OptimizedBackprojectionY(outputRegionForThread, matrix, projection); + continue; + } + + // Go over each voxel + itOut.GoToBegin(); + while (!itOut.IsAtEnd()) + { + // Compute projection index + for (unsigned int i = 0; i < Dimension - 1; i++) + { + pointProj[i] = matrix[i][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + pointProj[i] += matrix[i][j] * itOut.GetIndex()[j]; + } + + // Apply perspective + double perspFactor_local = matrix[Dimension - 1][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + perspFactor_local += matrix[Dimension - 1][j] * itOut.GetIndex()[j]; + perspFactor_local = 1 / perspFactor_local; + for (unsigned int i = 0; i < Dimension - 1; i++) + pointProj[i] = pointProj[i] * perspFactor_local; + + // Interpolate if in projection + if (interpolator->IsInsideBuffer(pointProj)) + { + itOut.Set(itOut.Get() + perspFactor_local * interpolator->EvaluateAtContinuousIndex(pointProj)); + } + + ++itOut; + } + } +} + +template +void +SuperShortScanBackProjectionImageFilter::OptimizedBackprojectionX( + const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection) +{ + typename ProjectionImageType::SizeType pSize = projection->GetBufferedRegion().GetSize(); + typename ProjectionImageType::IndexType pIndex = projection->GetBufferedRegion().GetIndex(); + typename TOutputImage::SizeType vBufferSize = this->GetOutput()->GetBufferedRegion().GetSize(); + typename TOutputImage::IndexType vBufferIndex = this->GetOutput()->GetBufferedRegion().GetIndex(); + typename TInputImage::PixelType * pProj = nullptr; + typename TOutputImage::PixelType * pVol = nullptr, *pVolZeroPointer = nullptr; + + // Pointers in memory to index (0,0,0) which do not necessarily exist + pVolZeroPointer = this->GetOutput()->GetBufferPointer(); + pVolZeroPointer -= vBufferIndex[0] + vBufferSize[0] * (vBufferIndex[1] + vBufferSize[1] * vBufferIndex[2]); + + // Continuous index at which we interpolate + double u = NAN, v = NAN, w = NAN; + int ui = 0, vi = 0; + double du = NAN; + + for (int k = region.GetIndex(2); k < region.GetIndex(2) + (int)region.GetSize(2); k++) + { + for (int j = region.GetIndex(1); j < region.GetIndex(1) + (int)region.GetSize(1); j++) + { + int i = region.GetIndex(0); + u = matrix[0][0] * i + matrix[0][1] * j + matrix[0][2] * k + matrix[0][3]; + v = matrix[1][1] * j + matrix[1][2] * k + matrix[1][3]; + w = matrix[2][1] * j + matrix[2][2] * k + matrix[2][3]; + + // Apply perspective + w = 1 / w; + u = u * w - pIndex[0]; + v = v * w - pIndex[1]; + du = w * matrix[0][0]; + +#ifdef BILINEAR_BACKPROJECTION + double u1 = NAN, u2 = NAN, v1 = NAN, v2 = NAN; + vi = itk::Math::floor(v); + if (vi >= 0 && vi < (int)pSize[1] - 1) + { + v1 = v - vi; + v2 = 1.0 - v1; +#else + vi = itk::Math::Round(v) - pIndex[1]; + if (vi >= 0 && vi < (int)pSize[1]) + { +#endif + + pProj = projection->GetBufferPointer() + vi * pSize[0]; + pVol = pVolZeroPointer + i + vBufferSize[0] * (j + k * vBufferSize[1]); + + // Innermost loop + for (; i < (region.GetIndex(0) + (int)region.GetSize(0)); i++, u += du, pVol++) + { +#ifdef BILINEAR_BACKPROJECTION + ui = itk::Math::floor(u); + if (ui >= 0 && ui < (int)pSize[0] - 1) + { + u1 = u - ui; + u2 = 1.0 - u1; + *pVol += w * (v2 * (u2 * *(pProj + ui) + u1 * *(pProj + ui + 1)) + + v1 * (u2 * *(pProj + ui + pSize[0]) + u1 * *(pProj + ui + pSize[0] + 1))); + } +#else + ui = itk::Math::Round(u); + if (ui >= 0 && ui < (int)pSize[0]) + { + *pVol += w * *(pProj + ui); + } +#endif + } // i + } + } // j + } // k +} + +template +void +SuperShortScanBackProjectionImageFilter::OptimizedBackprojectionY( + const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection) +{ + typename ProjectionImageType::SizeType pSize = projection->GetBufferedRegion().GetSize(); + typename ProjectionImageType::IndexType pIndex = projection->GetBufferedRegion().GetIndex(); + typename TOutputImage::SizeType vBufferSize = this->GetOutput()->GetBufferedRegion().GetSize(); + typename TOutputImage::IndexType vBufferIndex = this->GetOutput()->GetBufferedRegion().GetIndex(); + typename TInputImage::PixelType * pProj = nullptr; + typename TOutputImage::PixelType * pVol = nullptr, *pVolZeroPointer = nullptr; + + // Pointers in memory to index (0,0,0) which do not necessarily exist + pVolZeroPointer = this->GetOutput()->GetBufferPointer(); + pVolZeroPointer -= vBufferIndex[0] + vBufferSize[0] * (vBufferIndex[1] + vBufferSize[1] * vBufferIndex[2]); + + // Continuous index at which we interpolate + double u = NAN, v = NAN, w = NAN; + int ui = 0, vi = 0; + double du = NAN; + + for (int k = region.GetIndex(2); k < region.GetIndex(2) + (int)region.GetSize(2); k++) + { + for (int i = region.GetIndex(0); i < region.GetIndex(0) + (int)region.GetSize(0); i++) + { + int j = region.GetIndex(1); + u = matrix[0][0] * i + matrix[0][1] * j + matrix[0][2] * k + matrix[0][3]; + v = matrix[1][0] * i + matrix[1][2] * k + matrix[1][3]; + w = matrix[2][0] * i + matrix[2][2] * k + matrix[2][3]; + + // Apply perspective + w = 1 / w; + u = u * w - pIndex[0]; + v = v * w - pIndex[1]; + du = w * matrix[0][1]; + +#ifdef BILINEAR_BACKPROJECTION + vi = itk::Math::floor(v); + if (vi >= 0 && vi < (int)pSize[1] - 1) + { +#else + vi = itk::Math::Round(v); + if (vi >= 0 && vi < (int)pSize[1]) + { +#endif + pVol = pVolZeroPointer + i + vBufferSize[0] * (j + k * vBufferSize[1]); + for (; j < (region.GetIndex(1) + (int)region.GetSize(1)); j++, pVol += vBufferSize[0], u += du) + { +#ifdef BILINEAR_BACKPROJECTION + ui = itk::Math::floor(u); + if (ui >= 0 && ui < (int)pSize[0] - 1) + { + double u1 = NAN, u2 = NAN, v1 = NAN, v2 = NAN; + pProj = projection->GetBufferPointer() + vi * pSize[0] + ui; + v1 = v - vi; + v2 = 1.0 - v1; + u1 = u - ui; + u2 = 1.0 - u1; + *pVol += w * (v2 * (u2 * *(pProj) + u1 * *(pProj + 1)) + + v1 * (u2 * *(pProj + pSize[0]) + u1 * *(pProj + pSize[0] + 1))); + } +#else + ui = itk::Math::Round(u); + if (ui >= 0 && ui < (int)pSize[0]) + { + pProj = projection->GetBufferPointer() + vi * pSize[0]; + *pVol += w * *(pProj + ui); + } +#endif + } // j + } + } // i + } // k +} + +} // end namespace rtk + +#endif diff --git a/test/rtkfdktest.cxx b/test/rtkfdktest.cxx new file mode 100644 index 000000000..8d987959f --- /dev/null +++ b/test/rtkfdktest.cxx @@ -0,0 +1,225 @@ +#include +#include +#include + +#include "rtkTest.h" +#include "rtkSheppLoganPhantomFilter.h" +#include "rtkDrawSheppLoganFilter.h" +#include "rtkConstantImageSource.h" +#include "rtkFieldOfViewImageFilter.h" + +#ifdef USE_CUDA +# include "rtkCudaFDKConeBeamReconstructionFilter.h" +#else +# include "rtkFDKConeBeamReconstructionFilter.h" +#endif + +/** + * \file rtkfdktest.cxx + * + * \brief Functional test for classes performing FDK reconstructions + * + * This test generates the projections of a simulated Shepp-Logan phantom. + * A CT image is reconstructed from each set of generated projection images + * using the FDK algorithm and the reconstructed CT image is compared to the + * expected results which is analytically computed. + * + * \author Simon Rit and Marc Vila + */ + +int +main(int, char **) +{ + constexpr unsigned int Dimension = 3; + using OutputPixelType = float; +#ifdef USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = itk::Image; +#endif + +#if FAST_TESTS_NO_CHECKS + constexpr unsigned int NumberOfProjectionImages = 3; +#else + constexpr unsigned int NumberOfProjectionImages = 180; +#endif + + // Constant image sources + using ConstantImageSourceType = rtk::ConstantImageSource; + ConstantImageSourceType::PointType origin; + ConstantImageSourceType::SizeType size; + ConstantImageSourceType::SpacingType spacing; + + ConstantImageSourceType::Pointer tomographySource = ConstantImageSourceType::New(); + origin[0] = -127.; + origin[1] = -127.; + origin[2] = -127.; +#if FAST_TESTS_NO_CHECKS + size[0] = 32; + size[1] = 32; + size[2] = 32; + spacing[0] = 8.; + spacing[1] = 8.; + spacing[2] = 8.; +#else + size[0] = 128; + size[1] = 128; + size[2] = 128; + spacing[0] = 2.; + spacing[1] = 2.; + spacing[2] = 2.; +#endif + tomographySource->SetOrigin(origin); + tomographySource->SetSpacing(spacing); + tomographySource->SetSize(size); + tomographySource->SetConstant(0.); + + ConstantImageSourceType::Pointer projectionsSource = ConstantImageSourceType::New(); + origin[0] = -254.; + origin[1] = -254.; + origin[2] = -254.; +#if FAST_TESTS_NO_CHECKS + size[0] = 32; + size[1] = 32; + size[2] = NumberOfProjectionImages; + spacing[0] = 32.; + spacing[1] = 32.; + spacing[2] = 32.; +#else + size[0] = 128; + size[1] = 128; + size[2] = NumberOfProjectionImages; + spacing[0] = 4.; + spacing[1] = 4.; + spacing[2] = 4.; +#endif + projectionsSource->SetOrigin(origin); + projectionsSource->SetSpacing(spacing); + projectionsSource->SetSize(size); + projectionsSource->SetConstant(0.); + + std::cout << "\n\n****** Case 1: No streaming ******" << std::endl; + + // Geometry object + using GeometryType = rtk::ThreeDCircularProjectionGeometry; + GeometryType::Pointer geometry = GeometryType::New(); + for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++) + geometry->AddProjection(600., 1200., noProj * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + // Shepp Logan projections filter + using SLPType = rtk::SheppLoganPhantomFilter; + SLPType::Pointer slp = SLPType::New(); + slp->SetInput(projectionsSource->GetOutput()); + slp->SetGeometry(geometry); + slp->SetPhantomScale(116); + TRY_AND_EXIT_ON_ITK_EXCEPTION(slp->Update()); + + // Create a reference object (in this case a 3D phantom reference). + using DSLType = rtk::DrawSheppLoganFilter; + DSLType::Pointer dsl = DSLType::New(); + dsl->SetInput(tomographySource->GetOutput()); + dsl->SetPhantomScale(116); + TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) + + // FDK reconstruction filtering +#ifdef USE_CUDA + using FDKType = rtk::CudaFDKConeBeamReconstructionFilter; +#else + using FDKType = rtk::FDKConeBeamReconstructionFilter; +#endif + FDKType::Pointer feldkamp = FDKType::New(); + feldkamp->SetInput(0, tomographySource->GetOutput()); + feldkamp->SetInput(1, slp->GetOutput()); + feldkamp->SetGeometry(geometry); + TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update()); + + + // FOV + using FOVFilterType = rtk::FieldOfViewImageFilter; + FOVFilterType::Pointer fov = FOVFilterType::New(); + fov->SetInput(0, feldkamp->GetOutput()); + fov->SetProjectionsStack(slp->GetOutput()); + fov->SetGeometry(geometry); + TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); + + CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); + std::cout << "Test PASSED! " << std::endl; + + std::cout << "\n\n****** Case 2: perpendicular direction ******" << std::endl; + + ConstantImageSourceType::OutputImageType::DirectionType direction; + direction[0][0] = 0; + direction[0][1] = 1; + direction[0][2] = 0; + direction[1][0] = -1; + direction[1][1] = 0; + direction[1][2] = 0; + direction[2][0] = 0; + direction[2][1] = 0; + direction[2][2] = 1; + tomographySource->SetDirection(direction); + origin[0] = -127.; + origin[1] = 127.; + origin[2] = -127.; + tomographySource->SetOrigin(origin); + fov->GetOutput()->ResetPipeline(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) + + CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); + std::cout << "Test PASSED! " << std::endl; + + std::cout << "\n\n****** Case 3: 45 degree tilt direction ******" << std::endl; + + direction[0][0] = 0.70710678118; + direction[0][1] = -0.70710678118; + direction[0][2] = 0.70710678118; + direction[1][0] = 0.70710678118; + direction[1][1] = 0; + direction[1][2] = 0; + direction[2][0] = 0; + direction[2][1] = 0; + direction[2][2] = 1; + tomographySource->SetDirection(direction); + origin[0] = -127.; + origin[1] = -127.; + origin[2] = -127.; + tomographySource->SetOrigin(origin); + TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->Update()) + TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->Update()); + + CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); + std::cout << "Test PASSED! " << std::endl; + + std::cout << "\n\n****** Case 4: streaming ******" << std::endl; + + // Make sure that the data will be recomputed by releasing them + fov->GetOutput()->ReleaseData(); + + using StreamingType = itk::StreamingImageFilter; + StreamingType::Pointer streamer = StreamingType::New(); + streamer->SetInput(0, fov->GetOutput()); + streamer->SetNumberOfStreamDivisions(8); + itk::ImageRegionSplitterDirection::Pointer splitter = itk::ImageRegionSplitterDirection::New(); + splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis + streamer->SetRegionSplitter(splitter); + TRY_AND_EXIT_ON_ITK_EXCEPTION(streamer->Update()); + + CheckImageQuality(streamer->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); + std::cout << "Test PASSED! " << std::endl; + + std::cout << "\n\n****** Case 5: small ROI ******" << std::endl; + origin[0] = -5.; + origin[1] = -13.; + origin[2] = -20.; + size[0] = 64; + size[1] = 64; + size[2] = 64; + tomographySource->SetOrigin(origin); + tomographySource->SetSize(size); + TRY_AND_EXIT_ON_ITK_EXCEPTION(fov->UpdateLargestPossibleRegion()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->UpdateLargestPossibleRegion()) + CheckImageQuality(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0); + std::cout << "Test PASSED! " << std::endl; + return EXIT_SUCCESS; +} diff --git a/wrapping/rtkSuperShortScanBackProjectionImageFilter.wrap b/wrapping/rtkSuperShortScanBackProjectionImageFilter.wrap new file mode 100644 index 000000000..9ddc83c73 --- /dev/null +++ b/wrapping/rtkSuperShortScanBackProjectionImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("rtk::SuperShortScanBackProjectionImageFilter" POINTER) + itk_wrap_image_filter("${WRAP_ITK_REAL}" 2 3) +itk_end_wrap_class()