ANIMA  4.0
animaAverageImages.cxx
Go to the documentation of this file.
1 #include <animaReadWriteFunctions.h>
2 #include <itkImageRegionIterator.h>
3 #include <itkMultiplyImageFilter.h>
4 #include <itkDivideImageFilter.h>
5 #include <itkMaskImageFilter.h>
6 #include <itkCastImageFilter.h>
7 
8 #include <tclap/CmdLine.h>
9 
10 int main(int argc, char **argv)
11 {
12  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
13 
14  TCLAP::ValueArg<std::string> inArg("i","inputfiles","Input image list in text file",true,"","input image list",cmd);
15  TCLAP::ValueArg<std::string> maskArg("m","maskfiles","Input masks list in text file (mask images should contain only zeros or ones)",false,"","input masks list",cmd);
16  TCLAP::ValueArg<std::string> weightsArg("w","weights","Weights list in text file",false,"","input weights list",cmd);
17  TCLAP::ValueArg<std::string> outArg("o","outputfile","Output image",true,"","output image",cmd);
18 
19  try
20  {
21  cmd.parse(argc,argv);
22  }
23  catch (TCLAP::ArgException& e)
24  {
25  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
26  return EXIT_FAILURE;
27  }
28 
29  std::ifstream imageIn(inArg.getValue());
30 
31  char refN[2048];
32  imageIn.getline(refN,2048);
33  while ((strcmp(refN,"") == 0)&&(!imageIn.eof()))
34  imageIn.getline(refN,2048);
35 
36  itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(refN, itk::ImageIOFactory::ReadMode);
37 
38  if (!imageIO)
39  {
40  std::cerr << "Unable to read input image " << inArg.getValue() << std::endl;
41  return EXIT_FAILURE;
42  }
43 
44  imageIO->SetFileName(refN);
45  imageIO->ReadImageInformation();
46 
47  // Return to file start
48  imageIn.clear();
49  imageIn.seekg(0, std::ios::beg);
50 
51  bool vectorImage = (imageIO->GetNumberOfComponents() > 1);
52 
53  using ImageType = itk::Image <long double, 3>;
54  using SaveImageType = itk::Image <double, 3>;
55  using VectorImageType = itk::VectorImage <long double, 3>;
56  using SaveVectorImageType = itk::VectorImage <double, 3>;
57  typedef itk::MultiplyImageFilter <ImageType, ImageType, ImageType> MultiplyFilterType;
58  typedef itk::DivideImageFilter <ImageType, ImageType, ImageType> DivideFilterType;
59  typedef itk::MultiplyImageFilter <VectorImageType, ImageType, VectorImageType> MultiplyVectorFilterType;
60  typedef itk::DivideImageFilter <VectorImageType, ImageType, VectorImageType> DivideVectorFilterType;
61 
62  using CastFilterType = itk::CastImageFilter <ImageType, SaveImageType>;
63  using CastVectorFilterType = itk::CastImageFilter <VectorImageType, SaveVectorImageType>;
64 
65  std::ifstream masksIn;
66  if (maskArg.getValue() != "")
67  masksIn.open(maskArg.getValue());
68 
69  std::ifstream weightsIn;
70  if (weightsArg.getValue() != "")
71  weightsIn.open(weightsArg.getValue());
72 
73  double sumWeights = 0;
74  char maskN[2048];
75 
76  if (!vectorImage)
77  {
78  ImageType::Pointer tmpOutput;
79  ImageType::Pointer tmpSumMasks;
80 
81  while (tmpOutput.IsNull())
82  {
83  imageIn.getline(refN,2048);
84 
85  if (masksIn.is_open())
86  masksIn.getline(maskN,2048);
87  if (strcmp(refN,"") == 0)
88  continue;
89 
90  tmpOutput = anima::readImage <ImageType> (refN);
91  double imageWeight = 1.0;
92  if (weightsIn.is_open())
93  weightsIn >> imageWeight;
94 
95  std::cout << "Adding image " << refN << " with weight " << imageWeight << "..." << std::endl;
96 
97  sumWeights += imageWeight;
98 
99  if (imageWeight != 1.0)
100  {
101  MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
102  multiplyFilter->SetInput1(tmpOutput);
103  multiplyFilter->SetConstant(imageWeight);
104  multiplyFilter->Update();
105 
106  tmpOutput = multiplyFilter->GetOutput();
107  tmpOutput->DisconnectPipeline();
108  }
109 
110  if (masksIn.is_open())
111  {
112  tmpSumMasks = anima::readImage <ImageType> (maskN);
113 
114  MultiplyFilterType::Pointer maskFilter = MultiplyFilterType::New();
115  maskFilter->SetInput1(tmpOutput);
116  maskFilter->SetInput2(tmpSumMasks);
117  maskFilter->Update();
118 
119  tmpOutput = maskFilter->GetOutput();
120  tmpOutput->DisconnectPipeline();
121 
122  if (imageWeight != 1.0)
123  {
124  MultiplyFilterType::Pointer multiplyMaskFilter = MultiplyFilterType::New();
125  multiplyMaskFilter->SetInput1(tmpSumMasks);
126  multiplyMaskFilter->SetConstant(imageWeight);
127  multiplyMaskFilter->Update();
128 
129  tmpSumMasks = multiplyMaskFilter->GetOutput();
130  tmpSumMasks->DisconnectPipeline();
131  }
132  }
133  }
134 
135  while(!imageIn.eof())
136  {
137  imageIn.getline(refN,2048);
138 
139  if (masksIn.is_open())
140  masksIn.getline(maskN,2048);
141  if (strcmp(refN,"") == 0)
142  continue;
143 
144  double imageWeight = 1.0;
145  if (weightsIn.is_open())
146  weightsIn >> imageWeight;
147 
148  sumWeights += imageWeight;
149 
150  std::cout << "Adding image " << refN << " with weight " << imageWeight << "..." << std::endl;
151 
152  ImageType::Pointer tmpImage = anima::readImage <ImageType> (refN);
153  itk::ImageRegionIterator <ImageType> tmpIt(tmpImage,tmpImage->GetLargestPossibleRegion());
154  itk::ImageRegionIterator <ImageType> resIt(tmpOutput,tmpImage->GetLargestPossibleRegion());
155 
156  if (masksIn.is_open())
157  {
158  ImageType::Pointer tmpMask = anima::readImage <ImageType> (maskN);
159 
160  itk::ImageRegionIterator <ImageType> tmpMaskIt(tmpMask,tmpImage->GetLargestPossibleRegion());
161  itk::ImageRegionIterator <ImageType> sumMasksIt(tmpSumMasks,tmpImage->GetLargestPossibleRegion());
162 
163  while (!sumMasksIt.IsAtEnd())
164  {
165  if (tmpMaskIt.Get() == 0)
166  {
167  ++tmpMaskIt;
168  ++sumMasksIt;
169  ++tmpIt;
170  ++resIt;
171  continue;
172  }
173 
174  sumMasksIt.Set(sumMasksIt.Get() + imageWeight * tmpMaskIt.Get());
175  resIt.Set(resIt.Get() + imageWeight * tmpIt.Get());
176 
177  ++tmpMaskIt;
178  ++sumMasksIt;
179  ++tmpIt;
180  ++resIt;
181  }
182  }
183  else
184  {
185  while (!resIt.IsAtEnd())
186  {
187  resIt.Set(resIt.Get() + imageWeight * tmpIt.Get());
188  ++tmpIt;
189  ++resIt;
190  }
191  }
192  }
193 
194  if (masksIn.is_open())
195  masksIn.close();
196 
197  if (weightsIn.is_open())
198  weightsIn.close();
199 
200  if (tmpSumMasks.IsNotNull())
201  {
202  DivideFilterType::Pointer divideFilter = DivideFilterType::New();
203  divideFilter->SetInput1(tmpOutput);
204  divideFilter->SetInput2(tmpSumMasks);
205  divideFilter->Update();
206 
207  typedef itk::MaskImageFilter <ImageType, ImageType, ImageType> MaskFilterType;
208  MaskFilterType::Pointer maskFilter = MaskFilterType::New();
209  maskFilter->SetInput(divideFilter->GetOutput());
210  maskFilter->SetMaskImage(tmpSumMasks);
211  maskFilter->Update();
212 
213  tmpOutput = maskFilter->GetOutput();
214  tmpOutput->DisconnectPipeline();
215  }
216  else
217  {
218  DivideFilterType::Pointer divideFilter = DivideFilterType::New();
219  divideFilter->SetInput1(tmpOutput);
220  divideFilter->SetConstant(sumWeights);
221  divideFilter->Update();
222 
223  tmpOutput = divideFilter->GetOutput();
224  tmpOutput->DisconnectPipeline();
225  }
226 
227  CastFilterType::Pointer castFilter = CastFilterType::New();
228  castFilter->SetInput(tmpOutput);
229  castFilter->Update();
230 
231  anima::writeImage <SaveImageType> (outArg.getValue(),castFilter->GetOutput());
232  }
233  else // vecArg is activated
234  {
235  VectorImageType::Pointer outputData;
236  ImageType::Pointer tmpSumMasks;
237 
238  while (outputData.IsNull())
239  {
240  imageIn.getline(refN,2048);
241 
242  if (masksIn.is_open())
243  masksIn.getline(maskN,2048);
244  if (strcmp(refN,"") == 0)
245  continue;
246 
247  outputData = anima::readImage <VectorImageType> (refN);
248  double imageWeight = 1.0;
249  if (weightsIn.is_open())
250  weightsIn >> imageWeight;
251 
252  std::cout << "Adding image " << refN << " with weight " << imageWeight << "..." << std::endl;
253 
254  sumWeights += imageWeight;
255 
256  if (imageWeight != 1.0)
257  {
258  MultiplyVectorFilterType::Pointer multiplyFilter = MultiplyVectorFilterType::New();
259  multiplyFilter->SetInput1(outputData);
260  multiplyFilter->SetConstant(imageWeight);
261  multiplyFilter->Update();
262 
263  outputData = multiplyFilter->GetOutput();
264  outputData->DisconnectPipeline();
265  }
266 
267  if (masksIn.is_open())
268  {
269  tmpSumMasks = anima::readImage <ImageType> (maskN);
270 
271  MultiplyVectorFilterType::Pointer maskFilter = MultiplyVectorFilterType::New();
272  maskFilter->SetInput1(outputData);
273  maskFilter->SetInput2(tmpSumMasks);
274  maskFilter->Update();
275 
276  outputData = maskFilter->GetOutput();
277  outputData->DisconnectPipeline();
278 
279  if (imageWeight != 1.0)
280  {
281  MultiplyFilterType::Pointer multiplyMaskFilter = MultiplyFilterType::New();
282  multiplyMaskFilter->SetInput1(tmpSumMasks);
283  multiplyMaskFilter->SetConstant(imageWeight);
284  multiplyMaskFilter->Update();
285 
286  tmpSumMasks = multiplyMaskFilter->GetOutput();
287  tmpSumMasks->DisconnectPipeline();
288  }
289  }
290  }
291 
292  while(!imageIn.eof())
293  {
294  imageIn.getline(refN,2048);
295 
296  if (masksIn.is_open())
297  masksIn.getline(maskN,2048);
298  if (strcmp(refN,"") == 0)
299  continue;
300 
301  double imageWeight = 1.0;
302  if (weightsIn.is_open())
303  weightsIn >> imageWeight;
304 
305  sumWeights += imageWeight;
306 
307  std::cout << "Adding image " << refN << " with weight " << imageWeight << "..." << std::endl;
308 
309  VectorImageType::Pointer tmpImage = anima::readImage <VectorImageType> (refN);
310 
311  itk::ImageRegionIterator <VectorImageType> tmpIt(tmpImage,tmpImage->GetLargestPossibleRegion());
312  itk::ImageRegionIterator <VectorImageType> resIt(outputData,tmpImage->GetLargestPossibleRegion());
313 
314  if (masksIn.is_open())
315  {
316  ImageType::Pointer tmpMask = anima::readImage <ImageType> (maskN);
317 
318  itk::ImageRegionIterator <ImageType> tmpMaskIt(tmpMask,tmpImage->GetLargestPossibleRegion());
319  itk::ImageRegionIterator <ImageType> sumMasksIt(tmpSumMasks,tmpImage->GetLargestPossibleRegion());
320 
321  while (!sumMasksIt.IsAtEnd())
322  {
323  if (tmpMaskIt.Get() == 0)
324  {
325  ++tmpMaskIt;
326  ++sumMasksIt;
327  ++tmpIt;
328  ++resIt;
329  continue;
330  }
331 
332  sumMasksIt.Set(sumMasksIt.Get() + imageWeight * tmpMaskIt.Get());
333  resIt.Set(resIt.Get() + imageWeight * tmpIt.Get());
334 
335  ++tmpMaskIt;
336  ++sumMasksIt;
337  ++tmpIt;
338  ++resIt;
339  }
340  }
341  else
342  {
343  while (!resIt.IsAtEnd())
344  {
345  resIt.Set(resIt.Get() + imageWeight * tmpIt.Get());
346  ++tmpIt;
347  ++resIt;
348  }
349  }
350  }
351 
352  if (masksIn.is_open())
353  masksIn.close();
354 
355  if (weightsIn.is_open())
356  weightsIn.close();
357 
358  if (tmpSumMasks.IsNotNull())
359  {
360  DivideVectorFilterType::Pointer divideFilter = DivideVectorFilterType::New();
361  divideFilter->SetInput1(outputData);
362  divideFilter->SetInput2(tmpSumMasks);
363  divideFilter->Update();
364 
365  typedef itk::MaskImageFilter <VectorImageType, ImageType, VectorImageType> MaskFilterType;
366  MaskFilterType::Pointer maskFilter = MaskFilterType::New();
367  maskFilter->SetInput(divideFilter->GetOutput());
368  maskFilter->SetMaskImage(tmpSumMasks);
369  maskFilter->Update();
370 
371  outputData = maskFilter->GetOutput();
372  outputData->DisconnectPipeline();
373  }
374  else
375  {
376  DivideVectorFilterType::Pointer divideFilter = DivideVectorFilterType::New();
377  divideFilter->SetInput1(outputData);
378  divideFilter->SetConstant(sumWeights);
379  divideFilter->Update();
380 
381  outputData = divideFilter->GetOutput();
382  outputData->DisconnectPipeline();
383  }
384 
385  CastVectorFilterType::Pointer castFilter = CastVectorFilterType::New();
386  castFilter->SetInput(outputData);
387  castFilter->Update();
388 
389  anima::writeImage <SaveVectorImageType> (outArg.getValue(),castFilter->GetOutput());
390  }
391 
392  imageIn.close();
393 
394  return EXIT_SUCCESS;
395 }
main
int main(int argc, char **argv)
Definition: animaAverageImages.cxx:10
Generated by   doxygen 1.8.13