diff --git a/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.h b/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.h new file mode 100644 index 00000000000..2218eb26adf --- /dev/null +++ b/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.h @@ -0,0 +1,169 @@ +/*========================================================================= + * + * Copyright SINAPSE: Scalable Informatics for Neuroscience, Processing and Software Engineering + * The University of Iowa + * + * 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. + * + *=========================================================================*/ +/* + * itkResampleInPlaceImageFilter.h + * + * + * Created by Wei Lu on 10/14/10. + * + */ + +#ifndef __itkResampleInPlaceImageFilter_h +#define __itkResampleInPlaceImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkVersorRigid3DTransform.h" + +namespace itk +{ +/** \class ResampleInPlaceImageFilter + * \brief Resample an image in place. + * + * The current ITK resample image filter will generate a physical memory- + * modified version of the input image if the input transform is not identity. The + * abuse use of the filter can be cumbersome in a situation that the image is very + * large, and there are lots of transform to be superimposed for the input image, and + * we even don't care about those intermediate transformed images. + * + * If all the transforms are rigid, a far superior way to achieve a similar result + * while minimizing the accumulative resampling errors as well as eliminating the expense + * on accessing the physical memory of the image is to compose all the + * transforms before hand if it is possible, and then we only need to resample the + * input image with the final composed transform once. + * + * Here we present a more compact alternative, all information is stored in the header + * of the image and there is no need to maintain the final transform any longer. ITK + * image class has innate support for doing this. + * + * \param RigidTransform -- Currently must be a VersorRigid3D + * \param InputImage -- The image to be duplicated and modified to incorporate the + * rigid transform + * \return -- An image with the same voxels values as the input, but with differnt + * physical space representation affected by the rigid transform. + * + * The purpose of this code is to generate the new origin and direction + * that will remove the need for using the transform. + * + * Given a set of PhysFixedImagePoints (i.e. from the fixedImage space) + * those points are converted to PhysMovingImagePoints = TfmF2M( PhysFixedImagePoints ) + * and then MovingContinuousIndexPoints = movingImage->TransformPhysicalPointToContinuousIndex( PhysMovingImagePoints ) + * to get image values. + * + * We desire to change the moving image DirectionCosign [DC] and Origin O such that + * we can compute the: + * MovingContinuousIndexPoints = newMovingImage->TransformPhysicalPointToContinuousIndex( PhysFixedImagePoints ) + * + *Image Notations: + * DC-DirectionCosign + * O-Origin + * + *Rigid Transform Notations: + * R-Rotation + * C-Center Of Rotation + * T-Translation + * + * TransformPhysicalPointToContinuousIndex: + * CI = [SP^-1][DC^-1]( PhysMovingImagePoints - O ) + * PhysMovingImagePoints = [R](PhysFixedImagePoints - C) + C + T + * + * After substitutions: + * MovingContinuousIndexPoints = [R^-1][DC][SP] * CI + [R^-1] * O - [R^-1] * C - [R^-1]*T + C + * ---------- -------------------------------------- + * NewDC NewOrigin + * NewDC = [R^-1][DC] + * NewOrigin = [R^-1] * ( O - C - T ) + C + * + * \ingroup GeometricTransforms + */ +template +class ResampleInPlaceImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(ResampleInPlaceImageFilter); + + /** Standard class type alias */ + using Self = ResampleInPlaceImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory */ + itkNewMacro(Self); + + /** Run-time type information (and related methods) */ + itkTypeMacro(ResampleInPlaceImageFilter, ImageToImageFilter); + + /** input/output image type alias */ + using InputImageType = TInputImage; + using InputImagePointer = typename InputImageType::Pointer; + using InputImageRegionType = typename InputImageType::RegionType; + using InputImagePixelType = typename InputImageType::PixelType; + using InputImagePointType = typename InputImageType::PointType; + + using OutputImageType = TOutputImage; + using OutputImagePointer = typename OutputImageType::Pointer; + using OutputImageRegionType = typename OutputImageType::RegionType; + using OutputImagePixelType = typename OutputImageType::PixelType; + + /** ImageDimension constants */ + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + +#ifdef ITK_USE_CONCEPT_CHECKING + itkConceptMacro(SameDimensionCheck, (Concept::SameDimension)); + itkConceptMacro(InputConvertibleToOutputCheck, (Concept::Convertible)); +#endif + + /** Transform type alias */ + using RigidTransformType = VersorRigid3DTransform; + using RigidTransformConstPointer = typename RigidTransformType::ConstPointer; + + /** Set/Get rigid transform. The default is an identity transform */ + itkSetConstObjectMacro(RigidTransform, RigidTransformType); + itkGetConstObjectMacro(RigidTransform, RigidTransformType); + + /** Set/Get required input image. (A wrapper to this->Set/GetInput()) */ + void + SetInputImage(const InputImageType * image); + + const InputImageType * + GetInputImage() const; + +protected: + ResampleInPlaceImageFilter(); + ~ResampleInPlaceImageFilter() override = default; + ; + + void + GenerateData() override; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + +private: + OutputImagePointer m_OutputImage; + RigidTransformConstPointer m_RigidTransform; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkResampleInPlaceImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.hxx b/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.hxx new file mode 100644 index 00000000000..e6d5d455e91 --- /dev/null +++ b/Modules/Core/Transform/include/itkResampleInPlaceImageFilter.hxx @@ -0,0 +1,117 @@ +/*========================================================================= + * + * Copyright SINAPSE: Scalable Informatics for Neuroscience, Processing and Software Engineering + * The University of Iowa + * + * 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. + * + *=========================================================================*/ +/* + * itkResampleInPlaceImageFilter.hxx + * + * + * Created by Wei Lu on 10/14/10. + * + */ + +#ifndef __itkResampleInPlaceImageFilter_hxx +#define __itkResampleInPlaceImageFilter_hxx + +#include "itkResampleInPlaceImageFilter.h" +#include "itkCastImageFilter.h" + +namespace itk +{ +/** + * Constructor + */ +template +ResampleInPlaceImageFilter::ResampleInPlaceImageFilter() + : m_OutputImage(nullptr) + , m_RigidTransform(nullptr) +{ + this->SetNumberOfRequiredInputs(1); +} + +/** + * Set/Get input image, required + */ +template +void +ResampleInPlaceImageFilter::SetInputImage(const InputImageType * image) +{ + this->SetInput(0, image); +} + +template +const typename ResampleInPlaceImageFilter::InputImageType * +ResampleInPlaceImageFilter::GetInputImage() const +{ + return this->GetInput(0); +} + +/** + * GenerateData Performs the in-place resampling + */ +template +void +ResampleInPlaceImageFilter::GenerateData() +{ + if (!this->GetInput()) + { + itkExceptionMacro(<< "Input image has not been connected"); + return; + } + + // SEE HEADER FILE FOR MATH DESCRIPTION + + { + /** make a cast copied version of the input image **/ + using DuplicatorType = CastImageFilter; + typename DuplicatorType::Pointer CastFilter = DuplicatorType::New(); + CastFilter->SetInput(this->GetInput()); + CastFilter->Update(); + m_OutputImage = CastFilter->GetOutput(); + } + + RigidTransformConstPointer FMTxfm = this->m_RigidTransform.GetPointer(); + const typename RigidTransformType::MatrixType inverseRotation(FMTxfm->GetMatrix().GetInverse()); + + // Modify the origin and direction info of the image to reflect the transform. + itk::Vector newOriginVector = + inverseRotation * (this->GetInput()->GetOrigin().GetVectorFromOrigin() - FMTxfm->GetCenter().GetVectorFromOrigin() - + FMTxfm->GetTranslation()) // NewOrigin = [R^-1] * ( O - C - T ) + C + + FMTxfm->GetCenter().GetVectorFromOrigin(); + itk::Point newOriginPoint; + for (int i = 0; i < 3; ++i) + { + newOriginPoint[i] = newOriginVector[i]; + } + m_OutputImage->SetOrigin(newOriginPoint); + m_OutputImage->SetDirection(inverseRotation * this->GetInput()->GetDirection()); // NewDC = [R^-1][DC] + + this->GraftOutput(m_OutputImage); +} + +template +void +ResampleInPlaceImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "Input transform: " << m_RigidTransform << std::endl; + os << indent << "Output image: " << m_OutputImage << std::endl; +} +} // end namespace itk + +#endif diff --git a/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx b/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx new file mode 100644 index 00000000000..77df227081a --- /dev/null +++ b/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx @@ -0,0 +1,159 @@ +/* ========================================================================= + * + * Copyright Insight Software 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. + * + *=========================================================================*/ +/* ========================================================================= + * Created by Wei Lu on 10/14/10. + * Copyright 2010 The University of Iowa + *=========================================================================*/ + +#include "itkResampleInPlaceImageFilter.h" + +#include "itkImage.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkVersorRigid3DTransform.h" +#include "itkImageRegionConstIterator.h" + +#include + +// return true if it fails the validation +inline bool +Validate(double input, double desired, double tolerance) +{ + return std::abs(input - desired) > tolerance * std::abs(desired); +} + +int +main(int argc, char * argv[]) +{ + // Simple parameter check + if (argc < 3) + { + std::cerr << "Wrong arguments!" << std::endl; + std::cerr << "Usage: ./" << argv[0] << " inputImage baselineImage outputImage" << std::endl; + return EXIT_FAILURE; + } + + bool result = false; // test result default = no failure + double tol = 1.e-3; // tolerance + + // Image, filter, transform typedef + constexpr unsigned int LocalImageDimension = 3; + using PixelType = short; + + using ImageType = itk::Image; + using ImagePointer = ImageType::Pointer; + using ImagePointType = ImageType::PointType; + using ImageDirectionType = ImageType::DirectionType; + using ImageSpacingType = ImageType::SpacingType; + + using ImageConstIterator = itk::ImageRegionConstIterator; + using ReaderType = itk::ImageFileReader; + using TransformType = itk::VersorRigid3DTransform; + + using FilterType = itk::ResampleInPlaceImageFilter; + + // Read in input test image + ImagePointer inputImage; + { + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & err) + { + std::cerr << " Error while reading image file(s) with ITK:\n " << err << std::endl; + } + inputImage = reader->GetOutput(); + } + + // Set up transforms + itk::Vector rotationAxis; + rotationAxis.Fill(0.); + rotationAxis[0] = 1.; + double rotationAngle = .5; // in rad + itk::Vector translation; + translation.Fill(0.); + translation[1] = 300.; // in mm along P-axis + TransformType::Pointer transform = TransformType::New(); + transform->SetIdentity(); + transform->SetRotation(rotationAxis, rotationAngle); + transform->Translate(translation, true); + + // Set up the resample filter + FilterType::Pointer filter = FilterType::New(); + filter->SetInputImage(inputImage); + filter->SetRigidTransform(transform); + filter->Update(); + ImagePointer outputImage = filter->GetOutput(); + + // Get image info + ImagePointType origin = outputImage->GetOrigin(); + ImageDirectionType direction = outputImage->GetDirection(); + ImageSpacingType spacing = outputImage->GetSpacing(); + + // Read in baseline image + ImagePointer baselineImage = nullptr; + { + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[2]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & err) + { + std::cerr << " Error while reading image file(s) with ITK:\n " << err << std::endl; + } + baselineImage = reader->GetOutput(); + } + ImagePointType origin_d = baselineImage->GetOrigin(); + ImageDirectionType direction_d = baselineImage->GetDirection(); + ImageSpacingType spacing_d = baselineImage->GetSpacing(); + // Image info validation + for (unsigned int i = 0; i < LocalImageDimension; ++i) + { + result = (result || Validate(origin[i], origin_d[i], tol)); + result = (result || Validate(spacing[i], spacing_d[i], tol)); + for (unsigned int j = 0; j < LocalImageDimension; ++j) + { + result = (result || Validate(direction(i, j), direction_d(i, j), tol)); + } + } + + // Voxel contents validation + ImageConstIterator it1(outputImage, outputImage->GetRequestedRegion()); + ImageConstIterator it2(baselineImage, baselineImage->GetRequestedRegion()); + it1.GoToBegin(); + it2.GoToBegin(); + while (!it1.IsAtEnd()) + { + result = (result || Validate(it1.Get(), it2.Get(), tol)); + ++it1; + ++it2; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(argv[3]); + writer->SetInput(outputImage); + writer->Update(); + + return result; +}