1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 // Software Guide : BeginLatex
20 //
21 // It is unfortunate that it is still very common to find medical image
22 // datasets that have been acquired with large inter-slice spacings that
23 // result in voxels with anisotropic shapes. In many cases these voxels have
24 // ratios of $[1:5]$ or even $[1:10]$ between the resolution in the plane $(x,y)$
25 // and the resolution along the $z$ axis. These datasets are close to
26 // \textbf{useless} for the purpose of computer-assisted image analysis. The
27 // abundance of datasets acquired with anisotropic voxel sizes bespeaks a
28 // dearth of understanding of the third dimension and its importance for
29 // medical image analysis in clinical settings and radiology reading rooms.
30 // Datasets acquired with large anisotropies bring with them the regressive
31 // message: \emph{``I do not think 3D is informative''}.
32 // They stubbornly insist: \emph{``all that you need to know, can be known
33 // by looking at individual slices, one by one''}. However, the fallacy of this
34 // statement is made evident by simply viewing the slices when
35 // reconstructed in any of the orthogonal planes. The rectangular pixel shape
36 // is ugly and distorted, and cripples any signal processing algorithm not
37 // designed specifically for this type of image.
38 //
39 // Image analysts have a long educational battle to fight in the radiological
40 // setting in order to bring the message that 3D datasets acquired with
41 // anisotropies larger than $[1:2]$ are simply dismissive of the most fundamental
42 // concept of digital signal processing: The Shannon Sampling
43 // Theorem~\cite{Shannon1948,Shannon1949}.
44 //
45 // Facing the inertia of many clinical imaging departments and their
46 // blithe insistence that these images are ``good enough''
47 // for image processing, some image analysts have stoically tried
48 // to deal with these poor datasets. These image analysts usually
49 // proceed to subsample the high in-plane resolution and to super-sample the
50 // inter-slice resolution with the purpose of faking the type of dataset that
51 // they should have received in the first place: an \textbf{isotropic} dataset.
52 // This example is an illustration of how such an operation can be performed using
53 // the filters available in the Insight Toolkit.
54 //
55 // Note that this example is not presented here as a \emph{solution} to the
56 // problem of anisotropic datasets. On the contrary, this is simply a
57 // \emph{dangerous palliative} which will only perpetuate the errant
58 // convictions of image acquisition departments. The real solution to the
59 // problem of the
60 // anisotropic dataset is to educate radiologists regarding the
61 // principles of image processing. If you really care about the technical
62 // decency of the medical image processing field, and you really care about
63 // providing your best effort to the patients who will receive health care
64 // directly or indirectly affected by your processed images, then it is your
65 // duty to reject anisotropic datasets and to patiently explain to your
66 // radiologist why anisotropic data are problematic for processing, and require
67 // crude workarounds which handicap your ability to draw accurate
68 // conclusions from the data and preclude his or her ability to provide
69 // quality care. Any barbarity such as a $[1:5]$ anisotropy ratio should be
70 // considered as a mere collection of slices, and not an authentic 3D dataset.
71 //
72 // Please, before employing the techniques covered in this section, do kindly
73 // invite your fellow radiologist to see the dataset in an orthogonal
74 // slice. Magnify that image in a viewer without any linear interpolation
75 // until you see the daunting reality of the rectangular pixels. Let her/him
76 // know how absurd it is to process digital data which have been sampled at
77 // ratios of $[1:5]$ or $[1:10]$. Then, inform them that your only option is
78 // to throw away all that high in-plane
79 // resolution and to \emph{make up} data between the slices in order to
80 // compensate for the low resolution. Only then will you be justified in
81 // using the following code.
82 //
83 // \index{Anisotropic data sets}
84 // \index{Subsampling}
85 // \index{Supersampling}
86 // \index{Resampling}
87 //
88 // Software Guide : EndLatex
89
90
91 // Software Guide : BeginLatex
92 //
93 // Let's now move into the code. It is appropriate for you to experience
94 // guilt\footnote{A feeling of regret or remorse for having committed some
95 // improper act; a recognition of one's own responsibility for doing something
96 // wrong.}, because your use the code below is the
97 // evidence that we have lost one more battle on the quest for real 3D dataset
98 // processing.
99 //
100 // This example performs subsampling on the in-plane resolution and performs
101 // super-sampling along the inter-slices resolution. The subsampling process
102 // requires that we preprocess the data with a smoothing filter in order to
103 // avoid the occurrence of aliasing effects due to overlap of the spectrum in
104 // the frequency domain~\cite{Shannon1948,Shannon1949}. The smoothing is
105 // performed here using the \code{RecursiveGaussian} filter, because it
106 // provides a convenient run-time performance.
107 //
108 // The first thing that you will need to do in order to resample this ugly
109 // anisotropic dataset is to include the header files for the
110 // \doxygen{ResampleImageFilter}, and the Gaussian smoothing filter.
111 //
112 // Software Guide : EndLatex
113
114 #include "itkImage.h"
115 #include "itkImageFileReader.h"
116 #include "itkImageFileWriter.h"
117
118 // Software Guide : BeginCodeSnippet
119 #include "itkResampleImageFilter.h"
120 #include "itkRecursiveGaussianImageFilter.h"
121 // Software Guide : EndCodeSnippet
122
123
124 // Software Guide : BeginLatex
125 //
126 // The resampling filter will need a Transform in order to map point
127 // coordinates and will need an interpolator in order to compute intensity
128 // values for the new resampled image. In this particular case we use the
129 // \doxygen{IdentityTransform} because the image is going to be resampled by
130 // preserving the physical extent of the sampled region. The Linear
131 // interpolator is used as a common trade-off\footnote{Although arguably we should use
132 // one type of interpolator for the in-plane subsampling process and another
133 // one for the inter-slice supersampling. But again, one should wonder why we
134 // apply any technical sophistication here, when we are covering up for
135 // an improper acquisition of medical data, trying
136 // to make it look as if it was correctly acquired.}.
137 //
138 // Software Guide : EndLatex
139
140
141 // Software Guide : BeginCodeSnippet
142 #include "itkIdentityTransform.h"
143 // Software Guide : EndCodeSnippet
144
145
146 // Software Guide : BeginLatex
147 //
148 // Note that, as part of the preprocessing of the image, in this example we are
149 // also rescaling the range of intensities. This operation has already been
150 // described as Intensity Windowing. In a real clinical application, this step
151 // requires careful consideration of the range of intensities that contain
152 // information about the anatomical structures that are of interest for the
153 // current clinical application. It practice you may want to remove this step
154 // of intensity rescaling.
155 //
156 // Software Guide : EndLatex
157
158 // Software Guide : BeginCodeSnippet
159 #include "itkIntensityWindowingImageFilter.h"
160 // Software Guide : EndCodeSnippet
161
162
main(int argc,char * argv[])163 int main( int argc, char * argv[] )
164 {
165 if( argc < 5 )
166 {
167 std::cerr << "Usage: " << std::endl;
168 std::cerr << argv[0] << " inputImageFile outputImageFile lower upper " << std::endl;
169 return EXIT_FAILURE;
170 }
171
172 // Software Guide : BeginLatex
173 //
174 // We make explicit now our choices for the pixel type and dimension of the
175 // input image to be processed, as well as the pixel type that we intend to use
176 // for the internal computation during the smoothing and resampling.
177 //
178 // Software Guide : EndLatex
179
180
181 // Software Guide : BeginCodeSnippet
182 constexpr unsigned int Dimension = 3;
183
184 using InputPixelType = unsigned short;
185 using InternalPixelType = float;
186
187 using InputImageType = itk::Image< InputPixelType, Dimension >;
188 using InternalImageType = itk::Image< InternalPixelType, Dimension >;
189 // Software Guide : EndCodeSnippet
190
191
192 using ReaderType = itk::ImageFileReader< InputImageType >;
193
194 ReaderType::Pointer reader = ReaderType::New();
195
196 reader->SetFileName( argv[1] );
197
198 try
199 {
200 reader->Update();
201 }
202 catch( itk::ExceptionObject & excep )
203 {
204 std::cerr << "Exception caught!" << std::endl;
205 std::cerr << excep << std::endl;
206 }
207
208 using IntensityFilterType = itk::IntensityWindowingImageFilter<
209 InputImageType,
210 InternalImageType >;
211
212 IntensityFilterType::Pointer intensityWindowing = IntensityFilterType::New();
213
214 intensityWindowing->SetWindowMinimum( std::stoi( argv[3] ) );
215 intensityWindowing->SetWindowMaximum( std::stoi( argv[4] ) );
216
217 intensityWindowing->SetOutputMinimum( 0.0 );
218 intensityWindowing->SetOutputMaximum( 255.0 ); // floats but in the range of chars.
219
220 intensityWindowing->SetInput( reader->GetOutput() );
221
222
223 // Software Guide : BeginLatex
224 //
225 // We instantiate the smoothing filter that will be used on the preprocessing
226 // for subsampling the in-plane resolution of the dataset.
227 //
228 // Software Guide : EndLatex
229
230 // Software Guide : BeginCodeSnippet
231 using GaussianFilterType = itk::RecursiveGaussianImageFilter<
232 InternalImageType,
233 InternalImageType >;
234 // Software Guide : EndCodeSnippet
235
236
237 // Software Guide : BeginLatex
238 //
239 // We create two instances of the smoothing filter: one will smooth along the
240 // $X$ direction while the other will smooth along the $Y$ direction. They are
241 // connected in a cascade in the pipeline, while taking their input from the
242 // intensity windowing filter. Note that you may want to skip the intensity
243 // windowing scale and simply take the input directly from the reader.
244 //
245 // Software Guide : EndLatex
246
247 // Software Guide : BeginCodeSnippet
248 GaussianFilterType::Pointer smootherX = GaussianFilterType::New();
249 GaussianFilterType::Pointer smootherY = GaussianFilterType::New();
250
251 smootherX->SetInput( intensityWindowing->GetOutput() );
252 smootherY->SetInput( smootherX->GetOutput() );
253 // Software Guide : EndCodeSnippet
254
255
256 // Software Guide : BeginLatex
257 //
258 // We must now provide the settings for the resampling itself. This is done by
259 // searching for a value of isotropic resolution that will provide a trade-off
260 // between the evil of subsampling and the evil of supersampling. We advance
261 // here the conjecture that the geometrical mean between the in-plane and the
262 // inter-slice resolutions should be a convenient isotropic resolution to use.
263 // This conjecture is supported on nothing other than intuition and common
264 // sense. You can rightfully argue that this choice deserves a more technical
265 // consideration, but then, if you are so concerned about the technical integrity
266 // of the image sampling process, you should not be using this code, and should
267 // discuss these issues with the radiologist who
268 // acquired this ugly anisotropic dataset.
269 //
270 // We take the image from the input and then request its array of pixel spacing
271 // values.
272 //
273 // Software Guide : EndLatex
274
275 // Software Guide : BeginCodeSnippet
276 InputImageType::ConstPointer inputImage = reader->GetOutput();
277
278 const InputImageType::SpacingType& inputSpacing = inputImage->GetSpacing();
279 // Software Guide : EndCodeSnippet
280
281
282 // Software Guide : BeginLatex
283 //
284 // and apply our ad-hoc conjecture that the correct anisotropic resolution
285 // to use is the geometrical mean of the in-plane and inter-slice resolutions.
286 // Then set this spacing as the Sigma value to be used for the Gaussian
287 // smoothing at the preprocessing stage.
288 //
289 // Software Guide : EndLatex
290
291 // Software Guide : BeginCodeSnippet
292 const double isoSpacing = std::sqrt( inputSpacing[2] * inputSpacing[0] );
293
294 smootherX->SetSigma( isoSpacing );
295 smootherY->SetSigma( isoSpacing );
296 // Software Guide : EndCodeSnippet
297
298
299 // Software Guide : BeginLatex
300 //
301 // We instruct the smoothing filters to act along the $X$ and $Y$ direction
302 // respectively.
303 //
304 // Software Guide : EndLatex
305
306
307 // Software Guide : BeginCodeSnippet
308 smootherX->SetDirection( 0 );
309 smootherY->SetDirection( 1 );
310 // Software Guide : EndCodeSnippet
311
312
313 // Software Guide : BeginLatex
314 //
315 // Now that we have taken care of the smoothing in-plane, we proceed to
316 // instantiate the resampling filter that will reconstruct an isotropic image.
317 // We start by declaring the pixel type to be used as the output of this filter,
318 // then instantiate the image type and the type for the resampling filter.
319 // Finally we construct an instantiation of the filter.
320 //
321 // Software Guide : EndLatex
322
323
324 // Software Guide : BeginCodeSnippet
325 using OutputPixelType = unsigned char;
326
327 using OutputImageType = itk::Image< OutputPixelType, Dimension >;
328
329 using ResampleFilterType = itk::ResampleImageFilter<
330 InternalImageType, OutputImageType >;
331
332 ResampleFilterType::Pointer resampler = ResampleFilterType::New();
333 // Software Guide : EndCodeSnippet
334
335
336 // Software Guide : BeginLatex
337 //
338 // The resampling filter requires that we provide a Transform, which in this
339 // particular case can simply be an identity transform.
340 //
341 // Software Guide : EndLatex
342
343 // Software Guide : BeginCodeSnippet
344 using TransformType = itk::IdentityTransform< double, Dimension >;
345
346 TransformType::Pointer transform = TransformType::New();
347 transform->SetIdentity();
348
349 resampler->SetTransform( transform );
350 // Software Guide : EndCodeSnippet
351
352
353 // Software Guide : BeginLatex
354 //
355 // The filter also requires an interpolator to be passed to it. In this case we
356 // chose to use a linear interpolator.
357 //
358 // Software Guide : EndLatex
359
360 // Software Guide : BeginCodeSnippet
361 using InterpolatorType = itk::LinearInterpolateImageFunction<
362 InternalImageType, double >;
363
364 InterpolatorType::Pointer interpolator = InterpolatorType::New();
365
366 resampler->SetInterpolator( interpolator );
367 // Software Guide : EndCodeSnippet
368
369
370 resampler->SetDefaultPixelValue( 255 ); // highlight regions without source
371
372
373 // Software Guide : BeginLatex
374 //
375 // The pixel spacing of the resampled dataset is loaded in a \code{SpacingType}
376 // and passed to the resampling filter.
377 //
378 // Software Guide : EndLatex
379
380 // Software Guide : BeginCodeSnippet
381 OutputImageType::SpacingType spacing;
382
383 spacing[0] = isoSpacing;
384 spacing[1] = isoSpacing;
385 spacing[2] = isoSpacing;
386
387 resampler->SetOutputSpacing( spacing );
388 // Software Guide : EndCodeSnippet
389
390
391 // Software Guide : BeginLatex
392 //
393 // The origin and orientation of the output image is maintained, since we
394 // decided to resample the image in the same physical extent of the input
395 // anisotropic image.
396 //
397 // Software Guide : EndLatex
398
399 // Software Guide : BeginCodeSnippet
400 resampler->SetOutputOrigin( inputImage->GetOrigin() );
401 resampler->SetOutputDirection( inputImage->GetDirection() );
402 // Software Guide : EndCodeSnippet
403
404
405 // Software Guide : BeginLatex
406 //
407 // The number of pixels to use along each dimension in the grid of the
408 // resampled image is computed using the ratio between the pixel spacings of the
409 // input image and those of the output image. Note that the computation of the
410 // number of pixels along the $Z$ direction is slightly different with the
411 // purpose of making sure that we don't attempt to compute pixels that are
412 // outside of the original anisotropic dataset.
413 //
414 // Software Guide : EndLatex
415
416 // Software Guide : BeginCodeSnippet
417 InputImageType::SizeType inputSize =
418 inputImage->GetLargestPossibleRegion().GetSize();
419
420 using SizeValueType = InputImageType::SizeType::SizeValueType;
421
422 const double dx = inputSize[0] * inputSpacing[0] / isoSpacing;
423 const double dy = inputSize[1] * inputSpacing[1] / isoSpacing;
424
425 const double dz = (inputSize[2] - 1 ) * inputSpacing[2] / isoSpacing;
426 // Software Guide : EndCodeSnippet
427
428
429 // Software Guide : BeginLatex
430 //
431 // Finally the values are stored in a \code{SizeType} and passed to the
432 // resampling filter. Note that this process requires a casting since the
433 // computations are performed in \code{double}, while the elements of the
434 // \code{SizeType} are integers.
435 //
436 // Software Guide : EndLatex
437
438 // Software Guide : BeginCodeSnippet
439 InputImageType::SizeType size;
440
441 size[0] = static_cast<SizeValueType>( dx );
442 size[1] = static_cast<SizeValueType>( dy );
443 size[2] = static_cast<SizeValueType>( dz );
444
445 resampler->SetSize( size );
446 // Software Guide : EndCodeSnippet
447
448
449 // Software Guide : BeginLatex
450 //
451 // Our last action is to take the input for the resampling image filter from
452 // the output of the cascade of smoothing filters, and then to trigger the
453 // execution of the pipeline by invoking the \code{Update()} method on the
454 // resampling filter.
455 //
456 // Software Guide : EndLatex
457
458
459 // Software Guide : BeginCodeSnippet
460 resampler->SetInput( smootherY->GetOutput() );
461
462 resampler->Update();
463 // Software Guide : EndCodeSnippet
464
465
466 // Software Guide : BeginLatex
467 //
468 // At this point we should take a moment in silence to reflect on the
469 // circumstances that have led us to accept this cover-up for the improper
470 // acquisition of medical data.
471 //
472 // Software Guide : EndLatex
473
474
475 using WriterType = itk::ImageFileWriter< OutputImageType >;
476
477 WriterType::Pointer writer = WriterType::New();
478
479 writer->SetFileName( argv[2] );
480 writer->SetInput( resampler->GetOutput() );
481
482 try
483 {
484 writer->Update();
485 }
486 catch( itk::ExceptionObject & excep )
487 {
488 std::cerr << "Exception caught !" << std::endl;
489 std::cerr << excep << std::endl;
490 }
491
492
493 return EXIT_SUCCESS;
494 }
495