1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
|
/*=========================================================================
*
* 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.
*
*=========================================================================*/
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ResampleImageFilterOutput1.png}
// ARGUMENTS: 0
// Software Guide : EndCommandLineArgs
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ResampleImageFilterOutput2.png}
// ARGUMENTS: 1
// Software Guide : EndCommandLineArgs
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ResampleImageFilterOutput3.png}
// ARGUMENTS: 2
// Software Guide : EndCommandLineArgs
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ResampleImageFilterOutput4.png}
// ARGUMENTS: 3
// Software Guide : EndCommandLineArgs
//
// Software Guide : BeginLatex
//
// Resampling an image is a very important task in image analysis. It is
// especially important in the frame of image registration. The
// \doxygen{ResampleImageFilter} implements image resampling through the use
// of \doxygen{Transform}s. The inputs expected by this filter are an image,
// a transform and an interpolator. The space coordinates of the image are
// mapped through the transform in order to generate a new image. The extent
// and spacing of the resulting image are selected by the user. Resampling
// is performed in space coordinates, not pixel/grid coordinates. It is
// quite important to ensure that image spacing is properly set on the
// images involved. The interpolator is required since the mapping from one
// space to the other will often require evaluation of the intensity of the
// image at non-grid positions.
//
// \index{itk::ResampleImageFilter}
//
// Software Guide : EndLatex
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
// Software Guide : BeginLatex
//
// The header file corresponding to this filter should be included first.
//
// \index{itk::ResampleImageFilter!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkResampleImageFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The header files corresponding to the transform and interpolator must
// also be included.
//
// \index{itk::AffineTransform!header}
// \index{itk::Nearest\-Neighbor\-Interpolate\-Image\-Function!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkAffineTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
// Software Guide : EndCodeSnippet
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,2,3}]" << std::endl;
return EXIT_FAILURE;
}
int exampleAction = 0;
if( argc >= 4 )
{
exampleAction = atoi( argv[3] );
}
// Software Guide : BeginLatex
//
// The dimension and pixel types for input and output image must be
// defined and with them the image types can be instantiated.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 2;
typedef unsigned char InputPixelType;
typedef unsigned char OutputPixelType;
typedef itk::Image< InputPixelType, Dimension > InputImageType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
// Software Guide : EndCodeSnippet
typedef itk::ImageFileReader< InputImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName( argv[1] );
writer->SetFileName( argv[2] );
// Software Guide : BeginLatex
//
// Using the image and transform types it is now possible to instantiate
// the filter type and create the filter object.
//
// \index{itk::ResampleImageFilter!instantiation}
// \index{itk::ResampleImageFilter!New()}
// \index{itk::ResampleImageFilter!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ResampleImageFilter<InputImageType,OutputImageType> FilterType;
FilterType::Pointer filter = FilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The transform type is typically defined using the image dimension
// and the type used for representing space coordinates.
//
// \index{itk::AffineTransform!instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::AffineTransform< double, Dimension > TransformType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// An instance of the transform object is instantiated and passed to the
// resample filter. By default, the parameters of the transform are set to
// represent the identity transform.
//
// \index{itk::ResampleImageFilter!SetTransform()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::Pointer transform = TransformType::New();
filter->SetTransform( transform );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The interpolator type is defined using the full image type and the type
// used for representing space coordinates.
//
// \index{itk::Nearest\-Neighbor\-Interpolate\-Image\-Function!instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::NearestNeighborInterpolateImageFunction<
InputImageType, double > InterpolatorType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// An instance of the interpolator object is instantiated and passed to
// the resample filter.
//
// \index{itk::ResampleImageFilter!SetInterpolator()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator( interpolator );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Given that some pixels of the output image may end up being mapped
// outside the extent of the input image it is necessary to decide what
// values to assign to them. This is done by invoking the
// \code{SetDefaultPixelValue()} method.
//
// \index{itk::ResampleImageFilter!SetDefaultPixelValue()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
filter->SetDefaultPixelValue( 0 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The sampling grid of the output space is specified with the spacing along
// each dimension and the origin.
//
// \index{itk::ResampleImageFilter!SetOutputOrigin()}
// \index{itk::ResampleImageFilter!SetOutputSpacing()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// pixel spacing in millimeters along X and Y
const double spacing[ Dimension ] = { 1.0, 1.0 };
filter->SetOutputSpacing( spacing );
// Physical space coordinate of origin for X and Y
const double origin[ Dimension ] = { 0.0, 0.0 };
filter->SetOutputOrigin( origin );
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
InputImageType::DirectionType direction;
direction.SetIdentity();
filter->SetOutputDirection( direction );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The extent of the sampling grid on the output image is defined by a
// \code{SizeType} and is set using the \code{SetSize()} method.
//
// \index{itk::ResampleImageFilter!SetSize()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
InputImageType::SizeType size;
size[0] = 300; // number of pixels along X
size[1] = 300; // number of pixels along Y
filter->SetSize( size );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The input to the filter can be taken from any other filter, for example
// a reader. The output can be passed down the pipeline to other filters,
// for example a writer. An update call on any downstream filter will
// trigger the execution of the resampling filter.
//
// \index{itk::ResampleImageFilter!SetInput()}
// \index{itk::ResampleImageFilter!GetOutput()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
filter->SetInput( reader->GetOutput() );
writer->SetInput( filter->GetOutput() );
writer->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{BrainProtonDensitySlice}
// \includegraphics[width=0.44\textwidth]{ResampleImageFilterOutput1}
// \itkcaption[Effect of the Resample filter]{Effect of the resample filter.}
// \label{fig:ResampleImageFilterOutput1}
// \end{figure}
//
// \begin{figure}
// \center
// \includegraphics[width=\textwidth]{ResampleImageFilterOutput1Analysis}
// \itkcaption[Analysis of resampling in common coordinate system]{Analysis of
// the resample image done in a common coordinate system.}
// \label{fig:ResampleImageFilterOutput1Analysis}
// \end{figure}
//
// Figure \ref{fig:ResampleImageFilterOutput1} illustrates the effect of
// this filter on a slice of MRI brain image using an affine transform
// containing an identity transform. Note that any analysis of the
// behavior of this filter must be done on the space coordinate system in
// millimeters, not with respect to the sampling grid in pixels. The
// figure shows the resulting image in the lower left quarter of the
// extent. This may seem odd if analyzed in terms of the image grid but is
// quite clear when seen with respect to space coordinates. Figure
// \ref{fig:ResampleImageFilterOutput1} is particularly misleading
// because the images are rescaled to fit nicely on the text of this book.
// Figure \ref{fig:ResampleImageFilterOutput1Analysis} clarifies the
// situation. It shows the two same images placed on an equally-scaled
// coordinate system. It becomes clear here that an identity transform is
// being used to map the image data, and that simply, we have requested to
// resample additional empty space around the image. The input image is
// $181 \times 217$ pixels in size and we have requested an output of $300
// \times 300$ pixels. In this case, the input and output images both have
// spacing of $1mm \times 1mm$ and origin of $(0.0,0.0)$.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Let's now set values on the transform. Note that the supplied transform
// represents the mapping of points from the output space to the input
// space. The following code sets up a translation.
//
// \index{itk::AffineTransform!Translate()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::OutputVectorType translation;
translation[0] = -30; // X translation in millimeters
translation[1] = -50; // Y translation in millimeters
transform->Translate( translation );
// Software Guide : EndCodeSnippet
if( exampleAction == 1 )
{
writer->Update();
}
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{BrainProtonDensitySlice}
// \includegraphics[width=0.44\textwidth]{ResampleImageFilterOutput2}
// \itkcaption[ResampleImageFilter with a translation by
// $(-30,-50)$]{ResampleImageFilter with a translation by $(-30,-50)$.}
// \label{fig:ResampleImageFilterOutput2}
// \end{figure}
//
// \begin{figure}
// \center
// \includegraphics[width=\textwidth]{ResampleImageFilterOutput2Analysis}
// \itkcaption[ResampleImageFilter. Analysis of a translation by
// $(-30,-50)$]{ResampleImageFilter. Analysis of a translation by
// $(-30,-50)$.}
// \label{fig:ResampleImageFilterOutput2Analysis}
// \end{figure}
//
// The output image resulting from the translation can be seen in Figure
// \ref{fig:ResampleImageFilterOutput2}. Again, it is better to interpret
// the result in a common coordinate system as illustrated in Figure
// \ref{fig:ResampleImageFilterOutput2Analysis}.
//
// Probably the most important thing to keep in mind when resampling images
// is that the transform is used to map points from the \textbf{output}
// image space into the \textbf{input} image space. In this case, Figure
// \ref{fig:ResampleImageFilterOutput2Analysis} shows that the translation
// is applied to every point of the output image and the resulting position
// is used to read the intensity from the input image. In this way, the
// gray level of the point $P$ in the output image is taken from the point
// $T(P)$ in the input image. Where $T$ is the transformation. In the
// specific case of the Figure
// \ref{fig:ResampleImageFilterOutput2Analysis}, the value of point
// $(105,188)$ in the output image is taken from the point $(75,138)$ of
// the input image because the transformation applied was a translation of
// $(-30,-50)$.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// It is sometimes useful to intentionally set the default output value to
// a distinct gray value in order to highlight the mapping of the image
// borders. For example, the following code sets the default external
// value of $100$. The result is shown in the right side of Figure
// \ref{fig:ResampleImageFilterOutput3Analysis}.
//
// \index{itk::ResampleImageFilter!SetDefaultPixelValue()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
filter->SetDefaultPixelValue( 100 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=\textwidth]{ResampleImageFilterOutput3Analysis}
// \itkcaption[ResampleImageFilter highlighting image
// borders]{ResampleImageFilter highlighting image borders with
// SetDefaultPixelValue().}
// \label{fig:ResampleImageFilterOutput3Analysis}
// \end{figure}
//
// With this change we can better appreciate the effect of the previous
// translation transform on the image resampling. Figure
// \ref{fig:ResampleImageFilterOutput3Analysis} illustrates how the point
// $(30,50)$ of the output image gets its gray value from the point $(0,0)$
// of the input image.
//
// Software Guide : EndLatex
if( exampleAction == 2 )
{
writer->Update();
}
return EXIT_SUCCESS;
}
|