DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
4DImageMath.cxx
Go to the documentation of this file.
1 
11 #include "4DImageMathCLP.h"
12 
13 #include "itkRegionOfInterestImageFilter.h"
14 #include "itkPermuteAxesImageFilter.h"
15 
18 #include "utl.h"
19 #include "itkFunctors.h"
22 
25 
27 
28 typedef double ScalarType;
29 typedef itk::Image<ScalarType, 4> NDImageType;
30 typedef itk::VectorImage<ScalarType, 3> VectorImageType;
31 
32 template <class IndexType>
33 bool
34 IsOutsideBox ( const IndexType& index, const std::vector<int>& box )
35 {
36  return index[0]<box[0] || index[0]>box[1] || index[1]<box[2] || index[1]>box[3] || index[2]<box[4] || index[2]>box[5];
37 }
38 
39 
40 enum
41 {
42  // no op
44  // binary op (scalar)
52  // unary op (scalar)
59  // unary op (vector)
66  // binary op (vector)
68  // multiple op
70  // complex op
73 };
74 
76 std::vector<int> _realBox;
77 std::vector<int> _size;
79 std::string _axis;
80 int _axisNumber=-1;
81 
82 void
83 SetOperationWithChecking( int &op, int value )
84 {
85  if ( op==OP_NULL )
86  op = value;
87  else if ( op==value )
88  ;
89  else
90  utlGlobalException(true, "Only one type of operation is allowed!");
91 }
92 
93 #define _SetOperationWithChecking(op, value, numberOfInputsCond) \
94 do \
95  { \
96  SetOperationWithChecking(op, value); \
97  utlSAGlobalException(!(numberOfInputsCond))(numberOfInputsCond).msg("wrong number of inputs!"); \
98  } while ( 0 );
99 
100 
101 void
102 GetImageFiles(const std::vector<std::string>& imageVec, std::string& inImage0, std::string& outImage)
103 {
104  utlGlobalException(imageVec.size()<2, "need to set at least 2 image files (at least one input image, and one output image)");
105  inImage0 = imageVec[0];
106  outImage = imageVec.back();
107 }
108 
109 int GetNumberFromAxis(std::string axis)
110 {
111  if (axis=="X" || axis=="x" || axis=="0") return 0;
112  else if (axis=="Y" || axis=="y" || axis=="1" ) return 1;
113  else if (axis=="Z" || axis=="z" || axis=="2" ) return 2;
114  else if (axis=="T" || axis=="t" || axis=="3" ) return 3;
115  else
116  utlGlobalException(true, "wrong logic");
117 }
118 
119 template <class ImageType, class Image2Type>
120 void
121 ConvertImage ( const itk::SmartPointer<ImageType>& image, itk::SmartPointer<Image2Type>& imageOut )
122 {
123  imageOut = image;
124 }
125 
126 template<>
127 void
128 ConvertImage<NDImageType, VectorImageType>( const itk::SmartPointer<NDImageType>& image, itk::SmartPointer<VectorImageType>& imageOut )
129 {
130  itk::MultiVolumeToVectorImage(image, imageOut);
131 }
132 
133 template<>
134 void
135 ConvertImage<VectorImageType, NDImageType>( const itk::SmartPointer<VectorImageType>& image, itk::SmartPointer<NDImageType>& imageOut )
136 {
137  itk::VectorToMultiVolumeImage(image, imageOut);
138 }
139 
140 
141 // template <class ImageType>
142 // void
143 // PermuteImage ( const itk::SmartPointer<ImageType>& image, const std::string& axis, itk::SmartPointer<NDImageType>& imageOut )
144 // {
145 // imageOut = NDImageType::New();
146 // typedef itk::PermuteAxesImageFilter <NDImageType> PermuteAxesImageFilterType;
147 // PermuteAxesImageFilterType::Pointer permuteAxesFilter = PermuteAxesImageFilterType::New();
148 // itk::FixedArray<unsigned int, 4> order;
149 // if (axis=="X")
150 // order[0]=3, order[1]=1, order[2]=2, order[3]=0;
151 // else if (axis=="Y")
152 // order[0]=0, order[1]=3, order[2]=2, order[3]=1;
153 // else if (axis=="Z")
154 // order[0]=0, order[1]=1, order[2]=3, order[3]=2;
155 // else
156 // utlGlobalException(true, "axis must be X or Y or Z");
157 // permuteAxesFilter->SetOrder(order);
158 // auto imageTmp = NDImageType::New();
159 // ConvertImage(image, imageTmp);
160 // permuteAxesFilter->SetInput(imageTmp);
161 // permuteAxesFilter->Update();
162 // imageOut = permuteAxesFilter->GetOutput();
163 // }
164 
165 #define __SetImage(ImageType, fileName, imageName, imageItName, axis) \
166  typename ImageType::Pointer imageName = ImageType::New(); \
167  itk::VectorImageRegionIteratorWithIndex<ImageType> imageItName; \
168  if (fileName!="") \
169  { \
170  itk::ReadImage<ImageType>(fileName, imageName); \
171  std::vector<int> sizeTmp = itk::GetVectorImage3DVolumeSize(imageName); \
172  utlGlobalException(!utl::IsSameVector(_size,sizeTmp), "images have different size"); \
173  imageItName = itk::VectorImageRegionIteratorWithIndex<ImageType> (imageName, imageName->GetLargestPossibleRegion(),axis); \
174  }
175 
176 
177 #define __SetImageOrScalar(ImageType, fileName, imageName, imageItName, valueName) \
178  typename ImageType::Pointer imageName = ImageType::New(); \
179  itk::VectorImageRegionIteratorWithIndex<ImageType> imageItName; \
180  double valueName=0; \
181  if (fileName!="") \
182  { \
183  if (utl::IsNumber(fileName)) \
184  { \
185  valueName = utl::ConvertStringToNumber<double>(fileName); \
186  } \
187  else \
188  { \
189  itk::ReadImage<ImageType>(fileName, imageName); \
190  std::vector<int> sizeTmp = itk::GetVectorImage3DVolumeSize(imageName); \
191  utlGlobalException(!utl::IsSameVector(_size,sizeTmp), "images have different size"); \
192  imageItName = itk::VectorImageRegionIteratorWithIndex<ImageType> (imageName, imageName->GetLargestPossibleRegion()); \
193  } \
194  }
195 
196 
197 // #define __XYZVectorToScalarFunctor(argName, opName, axisStr, funcName) \
198 // do { \
199 // if (argName##Arg.isSet()) \
200 // { \
201 // typename NDImageType::Pointer outImagePermute = NDImageType::New(); \
202 // typename NDImageType::Pointer outImagePermute2 = NDImageType::New(); \
203 // SetOperationWithChecking(_Operation, opName); \
204 // PermuteImage(image, axisStr, imagePermute); \
205 // UnaryVectorOPImage<NDImageType, NDImageType>(imagePermute, outImagePermute, _MaskImageFile, funcName); \
206 // PermuteImage(outImagePermute, axisStr, outImagePermute2); \
207 // ConvertImage(outImagePermute2, outImage); \
208 // } \
209 // } while(0)
210 
211 
212 
213 #define __FunctionFromStringOpImage(funcStr) \
214 do \
215  { \
216  itk::FunctorFromStringOPImage(imageVec, outImage, funcStr, mask, _numberOfThreads); \
217  } while ( 0 );
218 
219 
220 #define __UnaryScalarFunctor(argNoAxis, opNoAxis, funcName) \
221  if (argNoAxis##Arg.isSet()) \
222  { \
223  _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()==1); \
224  utl::Functor::ScalarFunctorWrapper<funcName > func; \
225  itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, func, mask, _numberOfThreads); \
226  }
227 
228 
229 #define __XYZTUnaryVectorFunctor(axisNum, argNoAxis, opNoAxis, funcName) \
230  if (argNoAxis##Arg.isSet()) \
231  { \
232  _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()==1); \
233  funcName func; \
234  itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, func, mask, _numberOfThreads, axisNum); \
235  }
236 
237 
238 #define __XYZTUnaryVectorFunctorWithArguments(axisNum, argNoAxis, opNoAxis, funcName) \
239  if (argNoAxis##Arg.isSet()) \
240  { \
241  _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()==1); \
242  funcName func; \
243  func.SetArguments(argNoAxis); \
244  itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, func, mask, _numberOfThreads, axisNum); \
245  }
246 
247 
248 #define __XYZTMultiVectorFunctor(axisNum, argNoAxis, opNoAxis, funcName) \
249  if (argNoAxis##Arg.isSet()) \
250  { \
251  _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()>=2); \
252  funcName func; \
253  itk::MultiVariableVectorOPImage<ImageType, ImageOutType>(imageVec, outImage, func, mask, _numberOfThreads, axisNum); \
254  }
255 
256 
257 #define __XYZTMultiVectorFunctor_FixedSize(axisNum, argNoAxis, opNoAxis, funcName, inSize) \
258  if (argNoAxis##Arg.isSet()) \
259  { \
260  _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()==inSize); \
261  funcName func; \
262  itk::MultiVariableVectorOPImage<ImageType, ImageOutType>(imageVec, outImage, func, mask, _numberOfThreads, axisNum); \
263  }
264 
265 
266 // #define __XYZTBinaryVectorFunctor(axisNum, argNoAxis, opNoAxis, funcName) \
267 // if (argNoAxis##Arg.isSet()) \
268 // { \
269 // _SetOperationWithChecking(_Operation, OP_##opNoAxis, imageVec.size()==2); \
270 // utlGlobalException(imageVec.size()!=2, "need 2 input images"); \
271 // BinaryVectorOPImage<ImageType, ImageOutType>(image, imageVec[1], outImage, _MaskImageFile, funcName, axisNum); \
272 // }
273 
274 // [>* Unary operations. Output image has the same size as the input image. <]
275 // template <class ImageType, class ImageOutType, class OpFunctor>
276 // void
277 // UnaryOPImage(const itk::SmartPointer<ImageType>& image, itk::SmartPointer<ImageOutType>& outImage, std::string _MaskImageFile, const OpFunctor& func)
278 // {
279 // itk::VectorImageRegionIteratorWithIndex<ImageType> it(image, image->GetLargestPossibleRegion());
280 
281 // typedef itk::Image<double, 4> ScalarImageType;
282 // __SetImage(ScalarImageType, _MaskImageFile, mask, maskIt, -1);
283 
284 // outImage = ImageOutType::New();
285 // itk::CopyImageInformation(image, outImage);
286 // outImage->Allocate();
287 // itk::VectorImageRegionIteratorWithIndex<ImageOutType> outIt(outImage, outImage->GetLargestPossibleRegion());
288 
289 // typename ImageType::IndexType index;
290 // itk::VariableLengthVector<double> inPixel, outPixel;
291 // outPixel.SetSize(itk::GetVectorImageVectorSize(image));
292 // for (it.GoToBegin(), maskIt.GoToBegin(), outIt.GoToBegin();
293 // !it.IsAtEnd();
294 // ++it, ++maskIt, ++outIt)
295 // {
296 // if (_MaskImageFile!="" && maskIt.Get()<=0)
297 // { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
298 
299 // index = it.GetIndex();
300 // if (IsOutsideBox(index, _realBox))
301 // { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
302 
303 // it.GetVector(inPixel);
304 
305 // for ( int i = 0; i < inPixel.GetSize(); ++i )
306 // outPixel[i] = func(inPixel[i]);
307 
308 // if (utl::IsLogDebug())
309 // {
310 // std::cout << "index = " << index << std::endl << std::flush;
311 // itk::PrintVariableLengthVector(inPixel, "inPixel");
312 // itk::PrintVariableLengthVector(outPixel, "outPixel");
313 // }
314 // outIt.SetVector(outPixel);
315 // }
316 // }
317 
319 template <class ImageType, class ImageOutType, class OpFunctor>
320 void
321 BinaryOPImage(const itk::SmartPointer<ImageType>& image, itk::SmartPointer<ImageOutType>& outImage, std::string _MaskImageFile, std::string _opImageFile, const OpFunctor& func)
322 {
323  utlGlobalException(_opImageFile=="", "cannot be empty");
324  bool isValue = utl::IsNumber(_opImageFile);
325 
326  itk::VectorImageRegionIteratorWithIndex<ImageType> it(image, image->GetLargestPossibleRegion());
327 
328  typedef itk::Image<double, 4> ScalarImageType;
329  __SetImage(ScalarImageType, _MaskImageFile, mask, maskIt, -1);
330 
331  __SetImageOrScalar(ImageType, _opImageFile, opImage, opIt, opValue);
332 
333  int imageVecSize = itk::GetVectorImageVectorSize(image);
334  if (!isValue)
335  {
336  int opImageVecSize = itk::GetVectorImageVectorSize(opImage);
337  utlSAGlobalException(opImageVecSize!=1 && opImageVecSize!=imageVecSize)
338  (opImageVecSize)(imageVecSize).msg("wrong size of inputImage and opImage");
339  }
340 
341  outImage = ImageOutType::New();
342  itk::CopyImageInformation(image, outImage);
343  outImage->Allocate();
344  itk::VectorImageRegionIteratorWithIndex<ImageOutType> outIt(outImage, outImage->GetLargestPossibleRegion());
345 
346  typename ImageType::IndexType index;
347  itk::VariableLengthVector<double> inPixel, outPixel, inTmp;
348  outPixel.SetSize(imageVecSize);
349  for (it.GoToBegin(), maskIt.GoToBegin(), opIt.GoToBegin(), outIt.GoToBegin();
350  !it.IsAtEnd();
351  ++it, ++maskIt, ++opIt, ++outIt)
352  {
353  if (_MaskImageFile!="" && maskIt.Get()<=0)
354  { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
355 
356  index = it.GetIndex();
357  if (IsOutsideBox(index, _realBox))
358  { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
359 
360  it.GetVector(inPixel);
361 
362  if (isValue)
363  {
364  for ( int i = 0; i < inPixel.GetSize(); ++i )
365  outPixel[i] = func(inPixel[i], opValue);
366  }
367  else
368  {
369  opIt.GetVector(inTmp);
370  for ( int i = 0; i < inPixel.GetSize(); ++i )
371  outPixel[i] = inTmp.Size()==1 ? func(inPixel[i], inTmp[0]) : func(inPixel[i], inTmp[i]);
372  }
373 
374  if (utl::IsLogDebug())
375  {
376  std::cout << "index = " << index << std::endl << std::flush;
377  itk::PrintVariableLengthVector(inPixel, "inPixel");
378  itk::PrintVariableLengthVector(outPixel, "outPixel");
379  }
380  outIt.SetVector(outPixel);
381  }
382 
383 }
384 
385 // template <class ImageType, class ImageOutType, class OpFunctor>
386 // void
387 // BinaryVectorOPImage(const itk::SmartPointer<ImageType>& image, const itk::SmartPointer<ImageType>& image2, itk::SmartPointer<ImageOutType>& outImage, std::string _MaskImageFile, const OpFunctor& func, int vectorAxis=3)
388 // {
389 // utlGlobalException(!itk::VerifyImageSize(image, image2, false), "the two images should have the same shape.");
390 // utlSAAssert(vectorAxis>=0 && vectorAxis<=3)(vectorAxis).msg("wrong vectorAxis");
391 // itk::VectorImageRegionIteratorWithIndex<ImageType> it(image, image->GetLargestPossibleRegion(), vectorAxis);
392 // itk::VectorImageRegionIteratorWithIndex<ImageType> it2(image2, image2->GetLargestPossibleRegion(), vectorAxis);
393 
394 // typedef itk::Image<double, 4> ScalarImageType;
395 // __SetImage(ScalarImageType, _MaskImageFile, mask, maskIt, vectorAxis);
396 
397 // int vecInputSize = itk::GetVectorImageVectorSize(image);
398 
399 // outImage = ImageOutType::New();
400 // itk::CopyImageInformation(image, outImage);
401 // int outVectorSize = func.GetOutputDimension(vecInputSize);
402 // int outDim= (vectorAxis==3) ? func.GetOutputDimension(vecInputSize) : vecInputSize;
403 // if (vectorAxis!=3)
404 // {
405 // typename ImageOutType::RegionType regionTmp = outImage->GetLargestPossibleRegion();
406 // typename ImageOutType::SizeType sizeTmp = regionTmp.GetSize();
407 // sizeTmp[vectorAxis]=func.GetOutputDimension(sizeTmp[vectorAxis]);
408 // outVectorSize = func.GetOutputDimension(sizeTmp[vectorAxis]);
409 // regionTmp.SetSize(sizeTmp);
410 // outImage->SetRegions(regionTmp);
411 // }
412 // utlPrintVar(utl::IsLogDebug(), vectorAxis, outDim, vecInputSize, func.GetOutputDimension(vecInputSize));
413 // utlSAException(outDim<=0)(outDim).msg("wrong outDim");
414 // SetVectorImageVectorSize(outImage, outDim);
415 // outImage->Allocate();
416 // // outImage->Print(std::cout<<"outImage =\n");
417 // itk::VectorImageRegionIteratorWithIndex<ImageOutType> outIt(outImage, outImage->GetLargestPossibleRegion(), vectorAxis);
418 
419 // typename ImageType::IndexType index;
420 // itk::VariableLengthVector<double> inPixel, inPixel2, outPixel;
421 // utl::Vector<double> inVec, inVec2, outVec(outVectorSize);
422 // outPixel.SetSize(outVectorSize);
423 // outPixel.Fill(0.0);
424 // for (it.GoToBegin(), it2.GoToBegin(), maskIt.GoToBegin(), outIt.GoToBegin();
425 // !it.IsAtEnd();
426 // ++it, ++it2, ++maskIt, ++outIt)
427 // {
428 // if (_MaskImageFile!="" && maskIt.Get()<=0)
429 // { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
430 
431 // index = it.GetIndex();
432 // if (IsOutsideBox(index, _realBox))
433 // { outPixel.Fill(0.0); outIt.SetVector(outPixel); continue; }
434 
435 // for ( int i = 0; i < (vectorAxis==3?1:vecInputSize); ++i )
436 // {
437 // it.GetVector(inPixel, i);
438 // it2.GetVector(inPixel2, i);
439 
440 // inVec = utl::VariableLengthVectorToUtlVector(inPixel);
441 // inVec2 = utl::VariableLengthVectorToUtlVector(inPixel2);
442 // outVec = func(inVec, inVec2);
443 // outPixel = utl::UtlVectorToVariableLengthVector(outVec);
444 
445 // if (utl::IsLogDebug())
446 // {
447 // std::cout << "index = " << index << ", i = " << i << std::endl << std::flush;
448 // itk::PrintVariableLengthVector(inPixel, "inPixel");
449 // itk::PrintVariableLengthVector(inPixel2, "inPixel2");
450 // itk::PrintVariableLengthVector(outPixel, "outPixel");
451 // }
452 // outIt.SetVector(outPixel,i);
453 // }
454 // }
455 // }
456 
457 template <class ImageType, class ImageOutType>
458 int
459 ImageMath(int argc, char const* argv[])
460 {
461  PARSE_ARGS;
462  utl::LogLevel = _Verbose;
463  _numberOfThreads = _NumberOfThreads;
464 
465  if (_AxisArg.isSet())
466  {
467  _axis = utl::StringToUpperCase(_Axis);
469  }
470 
471  std::string _InputImageFile, _OutputImageFile;
472  GetImageFiles(_ImageFiles, _InputImageFile, _OutputImageFile);
473 
474  typedef itk::Image<double,4> Scalar4DImageType;
475  typename Scalar4DImageType::Pointer mask = Scalar4DImageType::New();
476  if (_MaskImageFileArg.isSet())
477  {
478  itk::ReadImage(_MaskImageFile, mask);
479  }
480 
481  typedef typename ImageType::Pointer ImagePointer;
482  typename ImageType::Pointer image = ImageType::New();
483  typename NDImageType::Pointer imagePermute = NDImageType::New();
484  typename ImageOutType::Pointer outImage = ImageOutType::New();
485 
486  typedef utl::Functor::VectorUnaryFunctionWrapper<> VectorUnaryFunctionWrapper;
487  typedef utl::Functor::VectorMultiVariableFunctionWrapper<> VectorMultiVariableFunctionWrapper;
488 
489  std::vector<ImagePointer> imageVec(_ImageFiles.size()-1);
490  for ( int i = 0; i < _ImageFiles.size()-1; ++i )
491  {
492  itk::ReadImage<ImageType>(_ImageFiles[i], imageVec[i]);
493  }
494  image = imageVec[0];
495  _size = itk::GetVectorImage3DVolumeSize(image);
496  int dimT = itk::GetVectorImageVectorSize(image);
497 
498  utlGlobalException(_Box.size()!=6, "wrong size of _Box");
499  _realBox = _Box;
500  _realBox[0] = _realBox[0]<0 ? 0 : (_realBox[0]>_size[0]-1 ? _size[0]-1 : _realBox[0]);
501  _realBox[1] = _realBox[1]==-1 ? _size[0]-1 : (_realBox[1]<0 ? 0 : (_realBox[1]>_size[0]-1 ? _size[0]-1 : _realBox[1]) );
502  _realBox[2] = _realBox[2]<0 ? 0 : (_realBox[2]>_size[1]-1 ? _size[1]-1 : _realBox[2]);
503  _realBox[3] = _realBox[3]==-1 ? _size[1]-1 : (_realBox[3]<0 ? 0 : (_realBox[3]>_size[1]-1 ? _size[1]-1 : _realBox[3]) );
504  _realBox[4] = _realBox[4]<0 ? 0 : (_realBox[4]>_size[2]-1 ? _size[2]-1 : _realBox[4]);
505  _realBox[5] = _realBox[5]==-1 ? _size[2]-1 : (_realBox[5]<0 ? 0 : (_realBox[5]>_size[2]-1 ? _size[2]-1 : _realBox[5]) );
506 
507  // binary
508  if (_AddImageFileArg.isSet())
509  {
510  _SetOperationWithChecking(_Operation, OP_ADD, imageVec.size()==1);
511  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _AddImageFile, std::plus<double>());
512  }
513  if (_MinusImageFileArg.isSet())
514  {
515  _SetOperationWithChecking(_Operation, OP_MINUS, imageVec.size()==1);
516  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _MinusImageFile, std::minus<double>());
517  }
518  if (_MultiplyImageFileArg.isSet())
519  {
520  _SetOperationWithChecking(_Operation, OP_MULTIPLY, imageVec.size()==1);
521  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _MultiplyImageFile, std::multiplies<double>());
522  }
523  if (_DivideImageFileArg.isSet())
524  {
525  _SetOperationWithChecking(_Operation, OP_DIVIDE, imageVec.size()==1);
526  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _DivideImageFile, std::divides<double>());
527  }
528  if (_MaxImageFileArg.isSet())
529  {
530  _SetOperationWithChecking(_Operation, OP_MAX, imageVec.size()==1);
531  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _MaxImageFile, utl::Functor::Max<double>());
532  }
533  if (_MinImageFileArg.isSet())
534  {
535  _SetOperationWithChecking(_Operation, OP_MIN, imageVec.size()==1);
536  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _MinImageFile, utl::Functor::Min<double>());
537  }
538  if (_PowImageFileArg.isSet())
539  {
540  _SetOperationWithChecking(_Operation, OP_POW, imageVec.size()==1);
541  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, _PowImageFile, utl::Functor::Pow<double>());
542  }
543 
544 
545  // math expression
546  if (_FunctorArg.isSet())
547  {
549  utlPrintVar(true, _Functor);
550  __FunctionFromStringOpImage(_Functor);
551  }
552 
553  // unary
559 
560  // unary vector
561  if (_NormArg.isSet())
562  {
563  _SetOperationWithChecking(_Operation, OP_NORM, imageVec.size()==1);
564  if (_Norm=="L2") itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, utl::Functor::TwoNorm<utl::Vector<double> >(), mask, _numberOfThreads, _axisNumber);
565  if (_Norm=="L1") itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, utl::Functor::OneNorm<utl::Vector<double> >(), mask, _numberOfThreads, _axisNumber);
566  if (_Norm=="L0") itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, utl::Functor::ZeroNorm<utl::Vector<double> >(), mask, _numberOfThreads, _axisNumber);
567  if (_Norm=="INF") itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, outImage, utl::Functor::InfNorm<utl::Vector<double> >(), mask, _numberOfThreads, _axisNumber);
568  }
574 
576 
578 
579  // __XYZTBinaryVectorFunctor(_axisNumber, _DotProduct, DOTPRODUCT, utl::Functor::DotProduct<utl::Vector<double> >());
581 
582  if (_CropArg.isSet())
583  {
584  _SetOperationWithChecking(_Operation, OP_CROP, imageVec.size()==1);
585  std::vector<int> cropBox(8,-1);
586  utlGlobalException(_Crop.size()>8, "wrong size of _Crop, should be no more than 8");
587  for ( int i = 0; i < _Crop.size(); ++i )
588  cropBox[i] = _Crop[i];
589  cropBox[0] = cropBox[0]<0 ? 0 : (cropBox[0]>_size[0]-1 ? _size[0]-1 : cropBox[0]);
590  cropBox[1] = cropBox[1]==-1 ? _size[0]-1 : (cropBox[1]<0 ? 0 : (cropBox[1]>_size[0]-1 ? _size[0]-1 : cropBox[1]) );
591  cropBox[2] = cropBox[2]<0 ? 0 : (cropBox[2]>_size[1]-1 ? _size[1]-1 : cropBox[2]);
592  cropBox[3] = cropBox[3]==-1 ? _size[1]-1 : (cropBox[3]<0 ? 0 : (cropBox[3]>_size[1]-1 ? _size[1]-1 : cropBox[3]) );
593  cropBox[4] = cropBox[4]<0 ? 0 : (cropBox[4]>_size[2]-1 ? _size[2]-1 : cropBox[4]);
594  cropBox[5] = cropBox[5]==-1 ? _size[2]-1 : (cropBox[5]<0 ? 0 : (cropBox[5]>_size[2]-1 ? _size[2]-1 : cropBox[5]) );
595  cropBox[6] = cropBox[6]<0 ? 0 : (cropBox[6]>dimT-1 ? dimT-1 : cropBox[6]);
596  cropBox[7] = cropBox[7]==-1 ? dimT-1 : (cropBox[7]<0 ? 0 : (cropBox[7]>dimT-1 ? dimT-1 : cropBox[7]) );
597 
598  utlPrintVar(true, dimT);
599  utl::PrintVector(_Crop, "_Crop");
600  utl::PrintVector(cropBox, "cropBox");
601  if (cropBox[6]==0 && cropBox[7]==dimT-1)
602  {
603  // no crop in t-axis
604  typename ImageType::Pointer imageTmp = ImageType::New();
605  typename ImageType::IndexType desiredStart;
606  typename ImageType::SizeType desiredSize;
607 
608  for ( int i = 0; i < ImageType::ImageDimension; ++i )
609  desiredStart[i] = cropBox[2*i];
610  for ( int i = 0; i < ImageType::ImageDimension; ++i )
611  desiredSize[i] = cropBox[2*i+1]-cropBox[2*i]+1;
612 
613  typename ImageType::RegionType desiredRegion(desiredStart, desiredSize);
614 
615  typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > FilterType;
616  typename FilterType::Pointer filter = FilterType::New();
617  filter->SetRegionOfInterest(desiredRegion);
618  filter->SetInput(image);
619  filter->Update();
620  imageTmp = filter->GetOutput();
621  ConvertImage(imageTmp, outImage);
622  }
623  else
624  {
625  // crop in t-axis
627  op.SetOffset(cropBox[6]);
628  op.SetChunkSize(cropBox[7]-cropBox[6]+1);
629  op.SetSpace(dimT);
630  typename ImageOutType::Pointer imageTmp = ImageOutType::New();
631  itk::UnaryVectorOPImage<ImageType, ImageOutType>(image, imageTmp, op, mask, _numberOfThreads);
632 
633  if (cropBox[0]==0 && cropBox[1]==_size[0]-1 && cropBox[2]==0 && cropBox[3]==_size[1]-1 && cropBox[4]==0 && cropBox[5]==_size[2]-1)
634  {
635  // no crop in x/y/z axis
636  outImage = imageTmp;
637  }
638  else
639  {
640  // crop in x/y/z axis
641  typename ImageOutType::IndexType desiredStart;
642  typename ImageOutType::SizeType desiredSize;
643 
644  for ( int i = 0; i < ImageOutType::ImageDimension; ++i )
645  desiredStart[i] = cropBox[2*i];
646  for ( int i = 0; i < ImageOutType::ImageDimension; ++i )
647  desiredSize[i] = cropBox[2*i+1]-cropBox[2*i]+1;
648 
649  // correct t-axis, since crop in t-axis has been done already.
650  if (ImageOutType::ImageDimension==4)
651  {
652  desiredStart[3] = 0;
653  desiredSize[3] = cropBox[7]-cropBox[6]+1;
654  }
655 
656  typename ImageOutType::RegionType desiredRegion(desiredStart, desiredSize);
657 
658  typedef itk::RegionOfInterestImageFilter< ImageOutType, ImageOutType > FilterType;
659  typename FilterType::Pointer filter = FilterType::New();
660  filter->SetRegionOfInterest(desiredRegion);
661  filter->SetInput(imageTmp);
662  filter->Update();
663  outImage = filter->GetOutput();
664  }
665  }
666  }
667 
668  if (_Operation!=OP_NULL)
669  itk::SaveImage(outImage, _OutputImageFile);
670  else if (_ImageFiles.size()==2)
671  {
672  // if no operation is set, then conversion of image formats (with mask if provided)
673  BinaryOPImage<ImageType, ImageOutType>(image, outImage, _MaskImageFile, utl::ConvertNumberToString(1.0), std::multiplies<double>());
674  // __FunctionFromStringOpImage("x");
675  itk::SaveImage(outImage, _OutputImageFile);
676  }
677  else
678  utlGlobalException(true, "wrong operator or inputs");
679 
680  return 0;
681 }
682 
686 int
687 main (int argc, char const* argv[])
688 {
689  // GenerateCLP
690  PARSE_ARGS;
691 
692  // utlGlobalException(!_OutputImageFileArg.isSet(), "need to set output -o");
693  std::string _InputImageFile, _OutputImageFile;
694  GetImageFiles(_ImageFiles, _InputImageFile, _OutputImageFile);
695 
696  if (!_OutputFormatArg.isSet() || (_OutputFormatArg.isSet() && _OutputFormat=="NONE"))
697  {
698  if (itk::IsVectorImage(_InputImageFile) || itk::Is3DImage(_InputImageFile))
699  return ImageMath<VectorImageType, VectorImageType>(argc, argv);
700  else
701  return ImageMath<NDImageType, NDImageType>(argc, argv);
702  }
703  else
704  {
705  if (itk::IsVectorImage(_InputImageFile) && _OutputFormat=="VECTOR")
706  return ImageMath<VectorImageType, VectorImageType>(argc, argv);
707  else if (itk::IsVectorImage(_InputImageFile) && _OutputFormat=="4D")
708  return ImageMath<VectorImageType, NDImageType>(argc, argv);
709  else if (!itk::IsVectorImage(_InputImageFile) && _OutputFormat=="VECTOR")
710  return ImageMath<NDImageType, VectorImageType>(argc, argv);
711  else if (!itk::IsVectorImage(_InputImageFile) && _OutputFormat=="4D")
712  return ImageMath<NDImageType, NDImageType>(argc, argv);
713  }
714 
715  return 0;
716 }
717 
NDArray<T,1> is a vector class which uses blas mkl.
Definition: utlVector.h:36
bool Is3DImage(const std::string &filename)
Definition: utlITK.h:259
#define __UnaryScalarFunctor(argNoAxis, opNoAxis, funcName)
Created "06-10-2016.
helper functions specifically used in dmritool
bool IsNumber(const std::string &input)
Definition: utlCore.h:726
#define __SetImage(ImageType, fileName, imageName, imageItName, axis)
void MultiVolumeToVectorImage(const SmartPointer< Image< ScalarType, Dim+1 > > &image1, SmartPointer< VectorImage< ScalarType, Dim > > &image2)
#define __XYZTMultiVectorFunctor_FixedSize(axisNum, argNoAxis, opNoAxis, funcName, inSize)
std::string ConvertNumberToString(const T value, const int precision=6)
Definition: utlCore.h:740
#define utlSAGlobalException(expr)
Definition: utlCoreMacro.h:362
bool SaveImage(const SmartPointer< ImageType > &image, const std::string &filename, const std::string &printInfo="Writing Image:")
Definition: utlITK.h:157
#define __XYZTUnaryVectorFunctor(axisNum, argNoAxis, opNoAxis, funcName)
void ConvertImage< VectorImageType, NDImageType >(const itk::SmartPointer< VectorImageType > &image, itk::SmartPointer< NDImageType > &imageOut)
#define _SetOperationWithChecking(op, value, numberOfInputsCond)
Definition: 4DImageMath.cxx:93
void PrintVariableLengthVector(const VariableLengthVector< T >vec, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout)
Definition: utlITK.h:45
void SetOperationWithChecking(int &op, int value)
Definition: 4DImageMath.cxx:83
int ImageMath(int argc, char const *argv[])
std::string StringToUpperCase(const std::string &str)
Definition: utlCore.h:301
int _axisNumber
Definition: 4DImageMath.cxx:80
bool IsLogDebug(const int level=utl::LogLevel)
Definition: utlCoreMacro.h:213
#define __FunctionFromStringOpImage(funcStr)
virtual void SetSpace(const int _arg)
Definition: utlFunctors.h:478
bool ReadImage(const std::string &filename, SmartPointer< ImageType > &image, const std::string &printInfo="Reading Image:")
Definition: utlITK.h:110
void CopyImageInformation(const SmartPointer< ImageWithInfoType > &imageFrom, SmartPointer< ImageType > &imageTo)
Definition: utlITK.h:552
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
virtual void SetOffset(const int _arg)
Definition: utlFunctors.h:476
int GetNumberFromAxis(std::string axis)
int _numberOfThreads
Definition: 4DImageMath.cxx:78
#define utlPrintVar(cond,...)
Definition: utlCoreMacro.h:309
std::vector< int > GetVectorImage3DVolumeSize(const SmartPointer< ImageType > &image)
Definition: utlITK.h:331
void BinaryOPImage(const itk::SmartPointer< ImageType > &image, itk::SmartPointer< ImageOutType > &outImage, std::string _MaskImageFile, std::string _opImageFile, const OpFunctor &func)
static int LogLevel
Definition: utlCoreMacro.h:203
void GetImageFiles(const std::vector< std::string > &imageVec, std::string &inImage0, std::string &outImage)
int main(int argc, char const *argv[])
4DImage math. It works for both itk::Image<double,4> and itk::VectorImage<double,3> ...
int GetVectorImageVectorSize(const SmartPointer< ImageType > &image, const int axis=-1)
Definition: utlITK.h:278
void PrintVector(const std::vector< T > &vec, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout, bool showStats=true)
Definition: utlCore.h:1002
#define __XYZTUnaryVectorFunctorWithArguments(axisNum, argNoAxis, opNoAxis, funcName)
#define __SetImageOrScalar(ImageType, fileName, imageName, imageItName, valueName)
#define __XYZTMultiVectorFunctor(axisNum, argNoAxis, opNoAxis, funcName)
void ConvertImage< NDImageType, VectorImageType >(const itk::SmartPointer< NDImageType > &image, itk::SmartPointer< VectorImageType > &imageOut)
std::string _axis
Definition: 4DImageMath.cxx:79
itk::Image< ScalarType, 4 > NDImageType
Definition: 4DImageMath.cxx:29
void VectorToMultiVolumeImage(const SmartPointer< VectorImage< ScalarType, Dim > > &image1, SmartPointer< Image< ScalarType, Dim+1 > > &image2)
obtain a part of a vector.
Definition: utlFunctors.h:412
bool IsVectorImage(const std::string &filename)
Definition: utlITK.h:253
itk::VectorImage< ScalarType, 3 > VectorImageType
Definition: 4DImageMath.cxx:30
int _Operation
Definition: 4DImageMath.cxx:75
void ConvertImage(const itk::SmartPointer< ImageType > &image, itk::SmartPointer< Image2Type > &imageOut)
bool IsOutsideBox(const IndexType &index, const std::vector< int > &box)
Definition: 4DImageMath.cxx:34
std::vector< int > _size
Definition: 4DImageMath.cxx:77
std::vector< int > _realBox
Definition: 4DImageMath.cxx:76
A multi-dimensional iterator templated over image type. It provides the same interfaces for both itk:...
virtual void SetChunkSize(const int _arg)
Definition: utlFunctors.h:477
#define ABS(a)
Definition: utils.h:51
void SetVector(const PixelVectorType &value, const int offIndex=0) const
double ScalarType
Definition: 4DImageMath.cxx:28