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 // This example shows how to use the \doxygen{ImageLinearIteratorWithIndex} for
22 // computing the mean across time of a 4D image where the first three
23 // dimensions correspond to spatial coordinates and the fourth dimension
24 // corresponds to time. The result of the mean across time is to be stored in a
25 // 3D image.
26 //
27 // \index{Iterators!and 4D images}
28 // \index{ImageLinearIteratorWithIndex!4D images}
29 //
30 // Software Guide : EndLatex
31
32 #include "itkImage.h"
33 // Software Guide : BeginCodeSnippet
34 #include "itkImageLinearConstIteratorWithIndex.h"
35 // Software Guide : EndCodeSnippet
36 #include "itkImageFileReader.h"
37 #include "itkImageFileWriter.h"
38
main(int argc,char * argv[])39 int main( int argc, char *argv[] )
40 {
41 // Verify the number of parameters on the command line.
42 if ( argc < 3 )
43 {
44 std::cerr << "Missing parameters. " << std::endl;
45 std::cerr << "Usage: " << std::endl;
46 std::cerr << argv[0]
47 << " input4DImageFile output3DImageFile"
48 << std::endl;
49 return EXIT_FAILURE;
50 }
51
52 // Software Guide : BeginLatex
53 //
54 // First we declare the types of the images, the 3D and 4D readers.
55 //
56 // Software Guide : EndLatex
57
58
59 // Software Guide : BeginCodeSnippet
60 using PixelType = unsigned char;
61 using Image3DType = itk::Image< PixelType, 3 >;
62 using Image4DType = itk::Image< PixelType, 4 >;
63
64 using Reader4DType = itk::ImageFileReader< Image4DType >;
65 using Writer3DType = itk::ImageFileWriter< Image3DType >;
66 // Software Guide : EndCodeSnippet
67
68 Reader4DType::Pointer reader4D = Reader4DType::New();
69 reader4D->SetFileName( argv[1] );
70
71 try
72 {
73 reader4D->Update();
74 }
75 catch( itk::ExceptionObject & excp )
76 {
77 std::cerr << "Error reading the image" << std::endl;
78 std::cerr << excp << std::endl;
79 return EXIT_FAILURE;
80 }
81
82 Image4DType::ConstPointer image4D = reader4D->GetOutput();
83
84 // Software Guide : BeginLatex
85 //
86 // Next, define the necessary types for indices, points, spacings, and size.
87 //
88 // Software Guide : EndLatex
89
90 // Software Guide : BeginCodeSnippet
91 Image3DType::Pointer image3D = Image3DType::New();
92 using Index3DType = Image3DType::IndexType;
93 using Size3DType = Image3DType::SizeType;
94 using Region3DType = Image3DType::RegionType;
95 using Spacing3DType = Image3DType::SpacingType;
96 using Origin3DType = Image3DType::PointType;
97
98 using Index4DType = Image4DType::IndexType;
99 using Size4DType = Image4DType::SizeType;
100 using Spacing4DType = Image4DType::SpacingType;
101 using Origin4DType = Image4DType::PointType;
102 // Software Guide : EndCodeSnippet
103
104 Index3DType index3D;
105 Size3DType size3D;
106 Spacing3DType spacing3D;
107 Origin3DType origin3D;
108
109 Image4DType::RegionType region4D = image4D->GetBufferedRegion();
110
111 Index4DType index4D = region4D.GetIndex();
112 Size4DType size4D = region4D.GetSize();
113 Spacing4DType spacing4D = image4D->GetSpacing();
114 Origin4DType origin4D = image4D->GetOrigin();
115
116 // Software Guide : BeginLatex
117 //
118 // Here we make sure that the values for our resultant 3D mean image
119 // match up with the input 4D image.
120 //
121 // Software Guide : EndLatex
122
123 // Software Guide : BeginCodeSnippet
124 for( unsigned int i=0; i < 3; i++)
125 {
126 size3D[i] = size4D[i];
127 index3D[i] = index4D[i];
128 spacing3D[i] = spacing4D[i];
129 origin3D[i] = origin4D[i];
130 }
131
132 image3D->SetSpacing( spacing3D );
133 image3D->SetOrigin( origin3D );
134
135 Region3DType region3D;
136 region3D.SetIndex( index3D );
137 region3D.SetSize( size3D );
138
139 image3D->SetRegions( region3D );
140 image3D->Allocate();
141 // Software Guide : EndCodeSnippet
142
143 using SumType = itk::NumericTraits< PixelType >::AccumulateType;
144 using MeanType = itk::NumericTraits< SumType >::RealType;
145
146 const unsigned int timeLength = region4D.GetSize()[3];
147
148 using IteratorType = itk::ImageLinearConstIteratorWithIndex<Image4DType>;
149
150 // Software Guide : BeginLatex
151 //
152 // Next we iterate over time in the input image series, compute the average,
153 // and store that value in the corresponding pixel of the output 3D image.
154 //
155 // Software Guide : EndLatex
156
157 // Software Guide : BeginCodeSnippet
158 IteratorType it( image4D, region4D );
159 it.SetDirection( 3 ); // Walk along time dimension
160 it.GoToBegin();
161 while( !it.IsAtEnd() )
162 {
163 SumType sum = itk::NumericTraits< SumType >::ZeroValue();
164 it.GoToBeginOfLine();
165 index4D = it.GetIndex();
166 while( !it.IsAtEndOfLine() )
167 {
168 sum += it.Get();
169 ++it;
170 }
171 MeanType mean = static_cast< MeanType >( sum ) /
172 static_cast< MeanType >( timeLength );
173
174 index3D[0] = index4D[0];
175 index3D[1] = index4D[1];
176 index3D[2] = index4D[2];
177
178 image3D->SetPixel( index3D, static_cast< PixelType >( mean ) );
179 it.NextLine();
180 }
181 // Software Guide : EndCodeSnippet
182
183
184 // Software Guide : BeginLatex
185 //
186 // As you can see, we avoid to use a 3D iterator to walk
187 // over the mean image. The reason is that there is no
188 // guarantee that the 3D iterator will walk in the same
189 // order as the 4D. Iterators just adhere to their contract
190 // of visiting every pixel, but do not enforce any particular
191 // order for the visits. The linear iterator guarantees it will
192 // visit the pixels along a line of the image in the order
193 // in which they are placed in the line, but does not state
194 // in what order one line will be visited with respect to
195 // other lines. Here we simply take advantage of knowing
196 // the first three components of the 4D iterator index,
197 // and use them to place the resulting mean value in the
198 // output 3D image.
199 //
200 // Software Guide : EndLatex
201
202 Writer3DType::Pointer writer3D = Writer3DType::New();
203 writer3D->SetFileName( argv[2] );
204 writer3D->SetInput( image3D );
205
206 try
207 {
208 writer3D->Update();
209 }
210 catch( itk::ExceptionObject & excp )
211 {
212 std::cerr << "Error writing the image" << std::endl;
213 std::cerr << excp << std::endl;
214 return EXIT_FAILURE;
215 }
216
217 return EXIT_SUCCESS;
218 }
219