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/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()