forked from InsightSoftwareConsortium/ITK
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ImageSliceIteratorWithIndex.cxx
322 lines (285 loc) · 11.2 KB
/
ImageSliceIteratorWithIndex.cxx
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
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: ImageSliceIteratorWithIndex.cxx
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
// Software Guide : BeginLatex
//
// \index{Iterators!and image slices}
//
// The \doxygen{ImageSliceIteratorWithIndex} class is an extension of
// \doxygen{ImageLinearIteratorWithIndex} from iteration along lines to
// iteration along both lines \emph{and planes} in an image.
// A \emph{slice} is a 2D
// plane spanned by two vectors pointing along orthogonal coordinate axes. The
// slice orientation of the slice iterator is defined by specifying its two
// spanning axes.
//
// \begin{itemize}
// \index{itk::Image\-Slice\-Iterator\-With\-Index!SetFirstDirection()}
// \item \textbf{\code{SetFirstDirection()}}
// Specifies the first coordinate axis
// direction of the slice plane.
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!SetSecondDirection()}
// \item \textbf{\code{SetSecondDirection()}}
// Specifies the second coordinate axis
// direction of the slice plane.
// \end{itemize}
//
// Several new methods control movement from slice to slice.
//
// \begin{itemize}
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!NextSlice()}
// \item \textbf{\code{NextSlice()}} Moves the iterator to the beginning pixel
// location of the next slice in the image. The origin of the next slice is
// calculated by incrementing the current origin index along the fastest
// increasing dimension of the image subspace which excludes the first and
// second dimensions of the iterator.
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!PreviousSlice()}
// \item \textbf{\code{PreviousSlice()}} Moves the iterator to the \emph{last
// valid pixel location} in the previous slice. The origin of the previous
// slice is calculated by decrementing the current origin index along the
// fastest increasing dimension of the image subspace which excludes the first
// and second dimensions of the iterator.
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!IsAtReverseEndOfSlice()}
// \item \textbf{\code{IsAtReverseEndOfSlice()}} Returns true if the iterator
// points to \emph{one position before} the beginning pixel of the current
// slice.
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!IsAtEndOfSlice()}
// \item \textbf{\code{IsAtEndOfSlice()}} Returns true if the iterator points
// to \emph{one position past} the last valid pixel of the current slice.
//
// \end{itemize}
//
// The slice iterator moves line by line using \code{NextLine()} and
// \code{PreviousLine()}. The line direction is parallel to the \emph{second}
// coordinate axis direction of the slice plane (see also
// Section~\ref{sec:itkImageLinearIteratorWithIndex}).
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!example of using|(}
// The next code example calculates the maximum intensity projection along one
// of the coordinate axes of an image volume. The algorithm is straightforward
// using ImageSliceIteratorWithIndex because we can coordinate
// movement through a slice of the 3D input image with movement through the 2D
// planar output.
//
// Here is how the algorithm works. For each 2D slice of the input, iterate
// through all the pixels line by line. Copy a pixel value to the corresponding
// position in the 2D output image if it is larger than the value already
// contained there. When all slices have been processed, the output image is
// the desired maximum intensity projection.
//
// We include a header for the const version of the slice iterator. For writing
// values to the 2D projection image, we use the linear iterator from the
// previous section. The linear iterator is chosen because it can be set to
// follow the same path in its underlying 2D image that the slice iterator
// follows over each slice of the 3D image.
//
// Software Guide : EndLatex
#include "itkImage.h"
#include "vnl/vnl_math.h"
// Software Guide : BeginCodeSnippet
#include "itkImageSliceConstIteratorWithIndex.h"
#include "itkImageLinearIteratorWithIndex.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
int main( int argc, char *argv[] )
{
// Verify the number of parameters on the command line.
if ( argc < 4 )
{
std::cerr << "Missing parameters. " << std::endl;
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0]
<< " inputImageFile outputImageFile projectionDirection"
<< std::endl;
return -1;
}
// Software Guide : BeginLatex
//
// The pixel type is defined as \code{unsigned short}. For this application,
// we need two image types, a 3D image for the input, and a 2D image for the
// intensity projection.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef unsigned short PixelType;
typedef itk::Image< PixelType, 2 > ImageType2D;
typedef itk::Image< PixelType, 3 > ImageType3D;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// A slice iterator type is defined to walk the input image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageLinearIteratorWithIndex< ImageType2D > LinearIteratorType;
typedef itk::ImageSliceConstIteratorWithIndex< ImageType3D > SliceIteratorType;
// Software Guide : EndCodeSnippet
typedef itk::ImageFileReader< ImageType3D > ReaderType;
typedef itk::ImageFileWriter< ImageType2D > WriterType;
ImageType3D::ConstPointer inputImage;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( argv[1] );
try
{
reader->Update();
inputImage = reader->GetOutput();
}
catch ( itk::ExceptionObject &err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
// Software Guide : BeginLatex
//
// The projection direction is read from the command line. The projection image
// will be the size of the 2D plane orthogonal to the projection direction.
// Its spanning vectors are the two remaining coordinate axes in the volume.
// These axes are recorded in the \code{direction} array.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
unsigned int projectionDirection =
static_cast<unsigned int>( ::atoi( argv[3] ) );
unsigned int i, j;
unsigned int direction[2];
for (i = 0, j = 0; i < 3; ++i )
{
if (i != projectionDirection)
{
direction[j] = i;
j++;
}
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The \code{direction} array is now used to define the projection image size
// based on the input image size. The output image is created so that its
// common dimension(s) with the input image are the same size. For example,
// if we project along the $x$ axis of the input, the size and origin of the
// $y$ axes of the input and output will match. This makes the code slightly
// more complicated, but prevents a counter-intuitive rotation of the output.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ImageType2D::RegionType region;
ImageType2D::RegionType::SizeType size;
ImageType2D::RegionType::IndexType index;
ImageType3D::RegionType requestedRegion = inputImage->GetRequestedRegion();
index[ direction[0] ] = requestedRegion.GetIndex()[ direction[0] ];
index[ 1- direction[0] ] = requestedRegion.GetIndex()[ direction[1] ];
size[ direction[0] ] = requestedRegion.GetSize()[ direction[0] ];
size[ 1- direction[0] ] = requestedRegion.GetSize()[ direction[1] ];
region.SetSize( size );
region.SetIndex( index );
ImageType2D::Pointer outputImage = ImageType2D::New();
outputImage->SetRegions( region );
outputImage->Allocate();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Next we create the necessary iterators. The const slice iterator walks
// the 3D input image, and the non-const linear iterator walks the 2D output
// image. The iterators are initialized to walk the same linear path through
// a slice. Remember that the \emph{second} direction of the slice iterator
// defines the direction that linear iteration walks within a slice.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
SliceIteratorType inputIt( inputImage, inputImage->GetRequestedRegion() );
LinearIteratorType outputIt( outputImage, outputImage->GetRequestedRegion() );
inputIt.SetFirstDirection( direction[1] );
inputIt.SetSecondDirection( direction[0] );
outputIt.SetDirection( 1 - direction[0] );
// Software Guide : EndCodeSnippet
// Software Guide: BeginLatex
//
// Now we are ready to compute the projection. The first step is to initialize
// all of the projection values to their nonpositive minimum value. The
// projection values are then updated row by row from the first slice of the
// input. At the end of the first slice, the input iterator steps to the first
// row in the next slice, while the output iterator, whose underlying image
// consists of only one slice, rewinds to its first row. The process repeats
// until the last slice of the input is processed.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
outputIt.GoToBegin();
while ( ! outputIt.IsAtEnd() )
{
while ( ! outputIt.IsAtEndOfLine() )
{
outputIt.Set( itk::NumericTraits<unsigned short>::NonpositiveMin() );
++outputIt;
}
outputIt.NextLine();
}
inputIt.GoToBegin();
outputIt.GoToBegin();
while( !inputIt.IsAtEnd() )
{
while ( !inputIt.IsAtEndOfSlice() )
{
while ( !inputIt.IsAtEndOfLine() )
{
outputIt.Set( vnl_math_max( outputIt.Get(), inputIt.Get() ));
++inputIt;
++outputIt;
}
outputIt.NextLine();
inputIt.NextLine();
}
outputIt.GoToBegin();
inputIt.NextSlice();
}
// Software Guide : EndCodeSnippet
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( argv[2] );
writer->SetInput(outputImage);
try
{
writer->Update();
}
catch ( itk::ExceptionObject &err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
// Software Guide : BeginLatex
//
// Running this example code on the 3D image
// \code{Examples/Data/BrainProtonDensity3Slices.mha} using the $z$-axis as
// the axis of projection gives the image shown in
// Figure~\ref{fig:ImageSliceIteratorWithIndexOutput}.
//
// \begin{figure}
// \centering
// \includegraphics[width=0.4\textwidth]{ImageSliceIteratorWithIndexOutput.eps}
// \itkcaption[Maximum intensity projection using ImageSliceIteratorWithIndex]{The
// maximum intensity projection through three slices of a volume.}
// \protect\label{fig:ImageSliceIteratorWithIndexOutput}
// \end{figure}
//
//
// \index{itk::Image\-Slice\-Iterator\-With\-Index!example of using|)}
//
// Software Guide : EndLatex
return 0;
}