Skip to content

Commit

Permalink
ENH: Update style and use the test of ResampleInPlaceImageFilter
Browse files Browse the repository at this point in the history
Also use more complicated transform (further away from identity).
  • Loading branch information
dzenanz committed Apr 13, 2021
1 parent 53f8473 commit e08395a
Show file tree
Hide file tree
Showing 5 changed files with 103 additions and 137 deletions.
101 changes: 51 additions & 50 deletions Modules/Core/Transform/include/itkResampleInPlaceImageFilter.h
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
/*=========================================================================
*
* Copyright SINAPSE: Scalable Informatics for Neuroscience, Processing and Software Engineering
* The University of Iowa
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
Expand All @@ -16,16 +15,9 @@
* limitations under the License.
*
*=========================================================================*/
/*
* itkResampleInPlaceImageFilter.h
*
*
* Created by Wei Lu on 10/14/10.
*
*/

#ifndef __itkResampleInPlaceImageFilter_h
#define __itkResampleInPlaceImageFilter_h
#ifndef itkResampleInPlaceImageFilter_h
#define itkResampleInPlaceImageFilter_h

#include "itkImageToImageFilter.h"
#include "itkVersorRigid3DTransform.h"
Expand All @@ -36,22 +28,19 @@ namespace itk
* \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
* modified version of the input image if the input transform is not identity.
* The abuse use of the resample filter can be a source of problems.
* It can exhaust memory if the image is very large.
* It WILL reduce image quality when there are lots of transforms to be
* superimposed for the input image. We usually don't even care
* about those intermediate transformed images!
*
* If all the transforms are rigid, there is a far superior way to achieve a similar result.
* Updating image metadata in-place removes the accumulated resampling errors
* as well as eliminating the expense of accessing the physical memory of the image.
* We need to compose all the transforms beforehand to make use of this filter.
*
* \param RigidTransform -- Currently must be a \class 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
Expand All @@ -60,36 +49,49 @@ namespace itk
* 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.
* Given a set of PhysicalFixedImagePoints (i.e. from the fixedImage space)
* those points are converted to PhysicalMovingImagePoints = TfmF2M( PhysicalFixedImagePoints )
* and then MovingContinuousIndexPoints = movingImage->TransformPhysicalPointToContinuousIndex(
* PhysicalMovingImagePoints ) to get image values.
*
* We desire to change the moving image DirectionCosign [DC] and Origin O such that
* We desire to change the moving image DirectionCosine [DC] and Origin [O] such that
* we can compute the:
* MovingContinuousIndexPoints = newMovingImage->TransformPhysicalPointToContinuousIndex( PhysFixedImagePoints )
* MovingContinuousIndexPoints = newMovingImage->TransformPhysicalPointToContinuousIndex( PhysicalFixedImagePoints )
*
* Image Notations:
* \f$\mathbf{D}\f$: Direction cosine matrix
* \f$\mathbf{o}\f$: Origin vector
* \f$\mathbf{S}\f$: Spacing
* \f$\mathbf{ci}\f$: Continouos index
* \f$\mathbf{D}^{'}\f$: New direction cosine matrix
* \f$\mathbf{o}^{'}\f$: New origin vector
*
*Image Notations:
* DC-DirectionCosign
* O-Origin
* Rigid Transform Notations:
* \f$\mathbf{R}\f$: Rotation matrix
* \f$\mathbf{c}\f$: Center of rotation vector
* \f$\mathbf{t}\f$: Translation vector
*
*Rigid Transform Notations:
* R-Rotation
* C-Center Of Rotation
* T-Translation
* \f$\mathbf{f}_{p}\f$: fixed image points in physical space
* \f$\mathbf{m}_{p}\f$: moving image points in physical space
*
* TransformPhysicalPointToContinuousIndex:
* CI = [SP^-1][DC^-1]( PhysMovingImagePoints - O )
* PhysMovingImagePoints = [R](PhysFixedImagePoints - C) + C + T
* \f[
* \mathbf{ci} = \mathbf{S}^{-1}\mathbf{D}^{-1}( \mathbf{m}_{p} - \mathbf{o} ) \\
* \mathbf{m}_{p} = \mathbf{R}(\mathbf{f}_{p} - \mathbf{c}) + \mathbf{c} + \mathbf{T}
* \f]
*
* 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
* \f[
* \mathbf{m}_{c} = \underbrace{\mathbf{R}^{-1}\mathbf{D}}_\text{new cosine}\mathbf{S} * \mathbf{i} +
* \underbrace{\mathbf{R}^{-1} * \mathbf{o} - \mathbf{R}^{-1} * \mathbf{c} - \mathbf{R}^{-1}*T}_\text{new origin} +
* \mathbf{c} \\
* \\
* \mathbf{D}^{'} = \mathbf{R}^{-1}\mathbf{D} \\
* \mathbf{o}^{'} = \mathbf{R}^{-1} * ( \mathbf{o} - \mathbf{c} - \mathbf{t} ) + \mathbf{c}
* \f]
*
* \ingroup ITKTransform
*/
template <typename TInputImage, typename TOutputImage>
class ResampleInPlaceImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
Expand Down Expand Up @@ -146,9 +148,8 @@ class ResampleInPlaceImageFilter : public ImageToImageFilter<TInputImage, TOutpu
GetInputImage() const;

protected:
ResampleInPlaceImageFilter();
ResampleInPlaceImageFilter() = default;
~ResampleInPlaceImageFilter() override = default;
;

void
GenerateData() override;
Expand Down
32 changes: 4 additions & 28 deletions Modules/Core/Transform/include/itkResampleInPlaceImageFilter.hxx
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
/*=========================================================================
*
* Copyright SINAPSE: Scalable Informatics for Neuroscience, Processing and Software Engineering
* The University of Iowa
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
Expand All @@ -16,36 +15,15 @@
* limitations under the License.
*
*=========================================================================*/
/*
* itkResampleInPlaceImageFilter.hxx
*
*
* Created by Wei Lu on 10/14/10.
*
*/

#ifndef __itkResampleInPlaceImageFilter_hxx
#define __itkResampleInPlaceImageFilter_hxx
#ifndef itkResampleInPlaceImageFilter_hxx
#define itkResampleInPlaceImageFilter_hxx

#include "itkResampleInPlaceImageFilter.h"
#include "itkCastImageFilter.h"

namespace itk
{
/**
* Constructor
*/
template <typename TInputImage, typename TOutputImage>
ResampleInPlaceImageFilter<TInputImage, TOutputImage>::ResampleInPlaceImageFilter()
: m_OutputImage(nullptr)
, m_RigidTransform(nullptr)
{
this->SetNumberOfRequiredInputs(1);
}

/**
* Set/Get input image, required
*/
template <typename TInputImage, typename TOutputImage>
void
ResampleInPlaceImageFilter<TInputImage, TOutputImage>::SetInputImage(const InputImageType * image)
Expand All @@ -60,9 +38,7 @@ ResampleInPlaceImageFilter<TInputImage, TOutputImage>::GetInputImage() const
return this->GetInput(0);
}

/**
* GenerateData Performs the in-place resampling
*/

template <typename TInputImage, typename TOutputImage>
void
ResampleInPlaceImageFilter<TInputImage, TOutputImage>::GenerateData()
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
da19a074ef1405bd57ab5a4b1e03452a9604a2cc13c676a098df4203770bd17b528aaace5b09be1d6cd2267f743e23a573e076c8003a3beb16995954bc905efb
6 changes: 6 additions & 0 deletions Modules/Core/Transform/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ itkCompositeTransformTest.cxx
itkTransformCloneTest.cxx
itkMultiTransformTest.cxx
itkTestTransformGetInverse.cxx
itkResampleInPlaceImageFilterTest.cxx
)

CreateTestDriver(ITKTransform "${ITKTransform-Test_LIBRARIES}" "${ITKTransformTests}")
Expand Down Expand Up @@ -178,6 +179,11 @@ itk_add_test(NAME itkTestTransformGetInverse
COMMAND ITKTransformTestDriver itkTestTransformGetInverse)
set_property(TEST itkTestTransformGetInverse APPEND PROPERTY LABELS RUNS_LONG)
set_tests_properties( itkTestTransformGetInverse PROPERTIES COST 50 )
itk_add_test(NAME itkResampleInPlaceImageFilterTest
COMMAND ITKTransformTestDriver itkResampleInPlaceImageFilterTest
${ITK_EXAMPLE_DATA_ROOT}/BrainProtonDensity3Slices.mha
DATA{Baseline/BrainProtonDensity3SlicesHardened.mha}
${ITK_TEST_OUTPUT_DIR}/BrainProtonDensity3SlicesHardened.mha)

set(ITKTransformGTests
itkBSplineTransformGTest.cxx
Expand Down
Loading

0 comments on commit e08395a

Please sign in to comment.