ITK  6.0.0
Insight Toolkit
Examples/Filtering/ResampleImageFilter3.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* 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.
*
*=========================================================================*/
// Software Guide : BeginLatex
//
// Previous examples have described the basic principles behind the
// \doxygen{ResampleImageFilter}. Now it's time to have some fun with it.
//
// Figure \ref{fig:ResampleImageFilterTransformComposition6} illustrates the
// general case of the resampling process. The origin and spacing of the
// output image has been selected to be different from those of the input
// image. The circles represent the \emph{center} of pixels. They are
// inscribed in a rectangle representing the \emph{coverage} of this pixel.
// The spacing specifies the distance between pixel centers along every
// dimension.
//
// The transform applied is a rotation of $30$ degrees. It is important to
// note here that the transform supplied to the
// \doxygen{ResampleImageFilter} is a \emph{clockwise} rotation. This
// transform rotates the \emph{coordinate system} of the output image 30
// degrees clockwise. When the two images are relocated in a common
// coordinate system---as in Figure
// \ref{fig:ResampleImageFilterTransformComposition6}---the result is that
// the frame of the output image appears rotated 30 degrees
// \emph{clockwise}. If the output image is seen with its coordinate system
// vertically aligned---as in Figure
// \ref{fig:ResampleImageFilterOutput9}---the image content appears rotated
// 30 degrees \emph{counter-clockwise}. Before continuing to read this
// section, you may want to meditate a bit on this fact while enjoying a cup
// of (Colombian) coffee.
//
// \begin{figure}
// \center
// \includegraphics[height=6cm]{ResampleImageFilterInput2x3}
// \includegraphics[height=4cm]{ResampleImageFilterOutput9}
// \itkcaption[Effect of a rotation on the resampling filter.]{Effect of a
// rotation on the resampling filter. Input image at left, output image at
// right.}
// \label{fig:ResampleImageFilterOutput9}
// \end{figure}
//
// \begin{figure}
// \center
// \includegraphics[width=\textwidth]{ResampleImageFilterTransformComposition6}
// \itkcaption[Input and output image placed in a common reference
// system]{Input and output image placed in a common reference system.}
// \label{fig:ResampleImageFilterTransformComposition6}
// \end{figure}
//
// The following code implements the conditions illustrated in Figure
// \ref{fig:ResampleImageFilterTransformComposition6} with two differences:
// the output spacing is $40$ times smaller and there are $40$ times more
// pixels in both dimensions. Without these changes, few
// details will be recognizable in the images. Note that the spacing and
// origin of the input image should be prepared in advance by using other
// means since this filter cannot alter the actual content of the
// input image in any way.
//
// Software Guide : EndLatex
#include "itkImage.h"
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " inputImageFile outputImageFile";
std::cerr << " [exampleAction={0,1}]" << std::endl;
return EXIT_FAILURE;
}
int exampleAction = 0;
if (argc >= 4)
{
exampleAction = std::stoi(argv[3]);
}
constexpr unsigned int Dimension = 2;
using InputPixelType = unsigned char;
using OutputPixelType = unsigned char;
using InputImageType = itk::Image<InputPixelType, Dimension>;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
auto reader = ReaderType::New();
auto writer = WriterType::New();
reader->SetFileName(argv[1]);
writer->SetFileName(argv[2]);
using FilterType =
auto filter = FilterType::New();
auto transform = TransformType::New();
using InterpolatorType =
auto interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
// Software Guide : BeginLatex
//
// In order to facilitate the interpretation of the transform we set the
// default pixel value to value be distinct from the image background.
//
// \index{itk::ResampleImageFilter!SetDefaultPixelValue()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
filter->SetDefaultPixelValue(100);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The spacing is selected here to be 40 times smaller than the one
// illustrated in Figure
// \ref{fig:ResampleImageFilterTransformComposition6}.
//
// \index{itk::ResampleImageFilter!SetOutputSpacing()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double spacing[Dimension];
spacing[0] = 40.0 / 40.0; // pixel spacing in millimeters along X
spacing[1] = 30.0 / 40.0; // pixel spacing in millimeters along Y
filter->SetOutputSpacing(spacing);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We will preserve the orientation of the input image by using the
// following call.
//
// \index{itk::ResampleImageFilter!SetOutputOrigin()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
filter->SetOutputDirection(reader->GetOutput()->GetDirection());
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Let us now set up the origin of the output image. Note that the values
// provided here will be those of the space coordinates for the output
// image pixel of index $(0,0)$.
//
// \index{itk::ResampleImageFilter!SetOutputOrigin()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double origin[Dimension];
origin[0] = 50.0; // X space coordinate of origin
origin[1] = 130.0; // Y space coordinate of origin
filter->SetOutputOrigin(origin);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The output image size is defined to be $40$ times the one illustrated
// on the Figure \ref{fig:ResampleImageFilterTransformComposition6}.
//
// \index{itk::ResampleImageFilter!SetSize()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
size[0] = 5 * 40; // number of pixels along X
size[1] = 4 * 40; // number of pixels along Y
filter->SetSize(size);
// Software Guide : EndCodeSnippet
filter->SetInput(reader->GetOutput());
writer->SetInput(filter->GetOutput());
// Software Guide : BeginLatex
//
// Rotations are performed around the origin of physical coordinates---not
// the image origin nor the image center. Hence, the process of
// positioning the output image frame as it is shown in Figure
// \ref{fig:ResampleImageFilterTransformComposition6} requires three
// steps. First, the image origin must be moved to the origin of the
// coordinate system. This is done by applying a translation equal to the
// negative values of the image origin.
//
// \index{itk::AffineTransform!Translate()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::OutputVectorType translation1;
translation1[0] = -origin[0];
translation1[1] = -origin[1];
transform->Translate(translation1);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In a second step, a rotation of $30$ degrees is performed. In the
// \doxygen{AffineTransform}, angles are specified in
// \emph{radians}. Also, a second boolean argument is used to specify if
// the current modification of the transform should be pre-composed or
// post-composed with the current transform content. In this case the
// argument is set to \code{false} to indicate that the rotation should be
// applied \emph{after} the current transform content.
//
// \index{itk::AffineTransform!Rotate2D()}
// \index{itk::AffineTransform!Composition}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const double degreesToRadians = std::atan(1.0) / 45.0;
transform->Rotate2D(-30.0 * degreesToRadians, false);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The third and final step implies translating the image origin back to
// its previous location. This is be done by applying a translation equal
// to the origin values.
//
// \index{itk::AffineTransform!Translate()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::OutputVectorType translation2;
translation2[0] = origin[0];
translation2[1] = origin[1];
transform->Translate(translation2, false);
filter->SetTransform(transform);
// Software Guide : EndCodeSnippet
if (exampleAction == 0)
{
try
{
writer->Update();
}
catch (const itk::ExceptionObject & excep)
{
std::cerr << "Exception caught !" << std::endl;
std::cerr << excep << std::endl;
}
}
// Software Guide : BeginLatex
//
// Figure \ref{fig:ResampleImageFilterOutput9} presents the actual input
// and output images of this example as shown by a correct viewer which
// takes spacing into account. Note the \emph{clockwise} versus
// \emph{counter-clockwise} effect discussed previously between the
// representation in Figure
// \ref{fig:ResampleImageFilterTransformComposition6} and Figure
// \ref{fig:ResampleImageFilterOutput9}.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// As a final exercise, let's track the mapping of an individual pixel.
// Keep in mind that the transformation is initiated by walking through
// the pixels of the \emph{output} image. This is the only way to ensure
// that the image will be generated without holes or redundant
// values. When you think about transformation it is always useful to
// analyze things from the output image towards the input image.
//
// Let's take the pixel with index $I=(1,2)$ from the output image. The
// physical coordinates of this point in the output image reference system
// are $P=( 1 \times 40.0 + 50.0, 2 \times 30.0 + 130.0 ) = (90.0,190.0)$
// millimeters.
//
// This point $P$ is now mapped through the \doxygen{AffineTransform} into
// the input image space. The operation subtracts the origin,
// applies a $30$ degrees rotation and adds the origin back. Let's follow
// those steps. Subtracting the origin from $P$ leads to
// $P1=(40.0,60.0)$, the rotation maps $P1$ to $P2=( 40.0 \times cos
// (30.0) + 60.0 \times sin (30.0), 40.0 \times sin(30.0) - 60.0 \times
// cos(30.0)) = (64.64,31.96)$. Finally this point is translated back by
// the amount of the image origin. This moves $P2$ to
// $P3=(114.64,161.96)$.
//
// The point $P3$ is now in the coordinate system of the input image. The
// pixel of the input image associated with this physical position is
// computed using the origin and spacing of the input image. $I=( ( 114.64
// - 60.0 )/ 20.0 , ( 161 - 70.0 ) / 30.0 )$ which results in
// $I=(2.7,3.0)$. Note that this is a non-grid position since the values
// are non-integers. This means that the gray value to be assigned to the
// output image pixel $I=(1,2)$ must be computed by interpolation of the
// input image values.
//
// In this particular code the interpolator used is simply a\newline
// \doxygen{NearestNeighborInterpolateImageFunction} which will assign the
// value of the closest pixel. This ends up being the pixel of index
// $I=(3,3)$ and can be seen from Figure
// \ref{fig:ResampleImageFilterTransformComposition6}.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
itkImageFileReader.h
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itkImage.h
itkAffineTransform.h
itk::AffineTransform
Definition: itkAffineTransform.h:101
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:90
itkNearestNeighborInterpolateImageFunction.h
itkImageFileWriter.h
itk::Size::SetSize
void SetSize(const SizeValueType val[VDimension])
Definition: itkSize.h:179
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
New
static Pointer New()
itk::NearestNeighborInterpolateImageFunction
Nearest neighbor interpolation of a scalar image.
Definition: itkNearestNeighborInterpolateImageFunction.h:39
itkResampleImageFilter.h
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44