OTB  7.0.0
Orfeo Toolbox
otbBandMathXImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbBandMathXImageFilter_hxx
23 #define otbBandMathXImageFilter_hxx
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkImageRegionConstIterator.h"
28 #include "itkImageScanlineConstIterator.h"
29 #include "itkImageScanlineIterator.h"
30 #include "itkImageRegionConstIteratorWithOnlyIndex.h"
31 #include "itkConstNeighborhoodIterator.h"
32 #include "itkNumericTraits.h"
33 #include "itkProgressReporter.h"
34 #include "otbMacro.h"
35 
36 #include <iostream>
37 #include <fstream>
38 #include <string>
39 
40 namespace otb
41 {
42 
44 template <class TImage>
46 {
47  // This number will be incremented each time an image
48  // is added over the one minimumrequired
49  this->SetNumberOfRequiredInputs(1);
50 
51  m_UnderflowCount = 0;
52  m_OverflowCount = 0;
53  m_ThreadUnderflow.SetSize(1);
54  m_ThreadOverflow.SetSize(1);
55 
56 
57  // idxX and idxY
58  adhocStruct ahcX;
59  ahcX.name = "idxX";
60  ahcX.type = 0;
61  m_VAllowedVarNameAuto.push_back(ahcX);
62 
63  adhocStruct ahcY;
64  ahcY.name = "idxY";
65  ahcY.type = 1;
66  m_VAllowedVarNameAuto.push_back(ahcY);
67 
68  m_SizeNeighbourhood = 10;
69 
70  m_ManyExpressions = true;
71 }
72 
74 template <class TImage>
76 {
77  m_Expression.clear();
78  m_VParser.clear();
80 
81  for (unsigned int i = 0; i < m_AImage.size(); ++i)
82  m_AImage[i].clear();
83  m_AImage.clear();
84 
85  m_VVarName.clear();
86  m_VAllowedVarNameAuto.clear();
87  m_VAllowedVarNameAddedByUser.clear();
88  m_VFinalAllowedVarName.clear();
89  m_VNotAllowedVarName.clear();
90  m_outputsDimensions.clear();
91 }
92 
93 
94 template <class TImage>
95 void BandMathXImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
96 {
97  Superclass::PrintSelf(os, indent);
98 
99  os << indent << "Expressions: " << std::endl;
100  for (unsigned int i = 0; i < m_Expression.size(); i++)
101  os << indent << m_Expression[i] << std::endl;
102  os << indent << "Computed values follow:" << std::endl;
103  os << indent << "UnderflowCount: " << m_UnderflowCount << std::endl;
104  os << indent << "OverflowCount: " << m_OverflowCount << std::endl;
105  os << indent << "itk::NumericTraits<typename PixelValueType>::NonpositiveMin() : " << itk::NumericTraits<PixelValueType>::NonpositiveMin() << std::endl;
106  os << indent << "itk::NumericTraits<typename PixelValueType>::max() : " << itk::NumericTraits<PixelValueType>::max() << std::endl;
107 }
108 
109 template <class TImage>
111 {
112 
113  std::stringstream sstm;
114  sstm << "im" << (idx + 1);
115  this->SetNthInput(idx, image, sstm.str());
116 }
117 
118 template <class TImage>
119 void BandMathXImageFilter<TImage>::SetNthInput(DataObjectPointerArraySizeType idx, const ImageType* image, const std::string& varName)
120 {
121 
122  ImageType* imagebis = const_cast<ImageType*>(image); // Useful for call of UpdateOutputInformation() (see below)
123  this->SetInput(idx, imagebis);
124 
125  // imiPhyX and imiPhyY
126  std::stringstream sstmPhyX;
127  adhocStruct ahcPhyX;
128  sstmPhyX << varName << "PhyX";
129  ahcPhyX.name = sstmPhyX.str();
130  ahcPhyX.type = 2;
131  ahcPhyX.info[0] = idx; // Input image #ID
132  m_VAllowedVarNameAuto.push_back(ahcPhyX);
133 
134  std::stringstream sstmPhyY;
135  adhocStruct ahcPhyY;
136  sstmPhyY << varName << "PhyY";
137  ahcPhyY.name = sstmPhyY.str();
138  ahcPhyY.type = 3;
139  ahcPhyY.info[0] = idx; // Input image #ID
140  m_VAllowedVarNameAuto.push_back(ahcPhyY);
141 
142  // imi
143  std::stringstream sstm_glob;
144  adhocStruct ahc_glob;
145  sstm_glob << varName;
146  ahc_glob.name = sstm_glob.str();
147  ahc_glob.type = 4;
148  ahc_glob.info[0] = idx; // Input image #ID
149  m_VAllowedVarNameAuto.push_back(ahc_glob);
150 
151  // Mandatory before call of GetNumberOfComponentsPerPixel
152  // Really important not to call the filter's UpdateOutputInformation method here:
153  // this method is not ready until all inputs, variables and expressions are set.
154  imagebis->UpdateOutputInformation();
155 
156  // imibj
157  for (unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
158  {
159  std::stringstream sstm;
160  adhocStruct ahc;
161  sstm << varName << "b" << (j + 1);
162  ahc.name = sstm.str();
163  ahc.type = 5;
164  ahc.info[0] = idx; // Input image #ID
165  ahc.info[1] = j; // Band #ID
166  m_VAllowedVarNameAuto.push_back(ahc);
167  }
168 
169  // imibjNkxp
170  for (unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
171  for (unsigned int x = 0; x <= m_SizeNeighbourhood; x++)
172  for (unsigned int y = 0; y <= m_SizeNeighbourhood; y++)
173  {
174  std::stringstream sstm;
175  adhocStruct ahc;
176  sstm << varName << "b" << (j + 1) << "N" << 2 * x + 1 << "x" << 2 * y + 1;
177  ahc.name = sstm.str();
178  ahc.type = 6;
179  ahc.info[0] = idx; // Input image #ID
180  ahc.info[1] = j; // Band #ID
181  ahc.info[2] = 2 * x + 1; // Size x direction (matrix convention = cols)
182  ahc.info[3] = 2 * y + 1; // Size y direction (matrix convention = rows)
183  m_VAllowedVarNameAuto.push_back(ahc);
184  }
185 
186  // imibjSTATS
187  std::vector<std::string> statsNames;
188  statsNames.push_back("Mini");
189  statsNames.push_back("Maxi");
190  statsNames.push_back("Mean");
191  statsNames.push_back("Sum");
192  statsNames.push_back("Var");
193 
194  for (unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
195  for (unsigned int t = 0; t < statsNames.size(); t++)
196  {
197  std::stringstream sstm;
198  adhocStruct ahc;
199  sstm << varName << "b" << (j + 1) << statsNames[t];
200  ahc.name = sstm.str();
201  ahc.type = 8;
202  ahc.info[0] = idx; // Input image #ID
203  ahc.info[1] = j; // Band #ID
204  ahc.info[2] = t; // Sub-type : 0 Mini, 1 Maxi, 2 Mean ...
205  m_VAllowedVarNameAuto.push_back(ahc);
206  }
207 }
208 
209 template <typename TImage>
211 {
212  return const_cast<TImage*>(this->GetInput(idx));
213 }
214 
215 template <typename TImage>
217 {
218  m_ManyExpressions = flag;
219 }
220 
221 template <typename TImage>
222 void BandMathXImageFilter<TImage>::SetExpression(const std::string& expression)
223 {
224  std::string expressionToBePushed = expression;
225 
226  if (expression.find(";") != std::string::npos)
227  {
228  std::ostringstream oss;
229  oss << "cat(";
230  for (unsigned int i = 0; i < expression.size(); ++i)
231  if (expression[i] == ';')
232  oss << ",";
233  else
234  oss << expression[i];
235 
236  oss << ")";
237  expressionToBePushed = oss.str();
238  }
239 
240  if (m_ManyExpressions)
241  m_Expression.push_back(expressionToBePushed);
242  else if (m_Expression.size() == 0)
243  m_Expression.push_back(expressionToBePushed);
244 
245  if (m_Expression.size() > 1)
246  this->SetNthOutput((DataObjectPointerArraySizeType)(m_Expression.size()) - 1, (TImage::New()).GetPointer());
247 
248  this->Modified();
249 }
250 
251 template <typename TImage>
253 {
254  m_Expression.clear();
255  this->Modified();
256 }
257 template <typename TImage>
258 void BandMathXImageFilter<TImage>::SetMatrix(const std::string& name, const std::string& definition)
259 {
260 
261  for (unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
262  if (name.compare(m_VAllowedVarNameAddedByUser[i].name) == 0)
263  itkExceptionMacro(<< "Variable name '" << name << "' already used." << std::endl);
264 
265 
266  if ((definition.find("{") != 0) || (definition.find("}")) != definition.size() - 1)
267  itkExceptionMacro(<< "Definition of a matrix must begin with { and end with } characters." << std::endl);
268 
269  // Get rid of { and } characters
270  std::string def;
271  for (unsigned int i = 1; i < definition.size() - 1; ++i)
272  def.push_back(definition[i]);
273 
274 
275  std::vector<std::vector<double>> mat;
276  std::istringstream iss(def);
277  std::string rows;
278  while (std::getline(iss, rows, ';'))
279  {
280  mat.push_back(std::vector<double>(0));
281  std::istringstream iss2(rows);
282  std::string elmt;
283  while (std::getline(iss2, elmt, ','))
284  {
285  std::istringstream iss3(elmt);
286  double val;
287  iss3 >> val;
288  mat.back().push_back(val);
289  }
290  }
291 
292 
293  // Check dimensions of the matrix
294  for (unsigned int i = 0; i < mat.size() - 1; i++)
295  if (mat[i].size() != mat[i + 1].size())
296  itkExceptionMacro(<< "Each row must have the same number of cols : " << definition << std::endl);
297 
298 
299  // Registration
300  adhocStruct ahc;
301  ahc.name = name;
302  ahc.type = 7;
303  ahc.info[0] = mat[0].size(); // Size x direction (matrix convention = cols)
304  ahc.info[1] = mat.size(); // Size y direction (matrix convention = rows)
305 
306  ahc.value = ValueType(ahc.info[1], ahc.info[0], 0.0);
307  for (unsigned int i = 0; i < mat.size(); i++)
308  for (unsigned int j = 0; j < mat[i].size(); j++)
309  ahc.value.At(i, j) = mat[i][j];
310 
311  m_VAllowedVarNameAddedByUser.push_back(ahc);
312 }
313 
314 
315 template <typename TImage>
316 void BandMathXImageFilter<TImage>::SetConstant(const std::string& name, double value)
317 {
318  for (unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
319  if (name.compare(m_VAllowedVarNameAddedByUser[i].name) == 0)
320  itkExceptionMacro(<< "Variable name '" << name << "' already used." << std::endl);
321 
322  adhocStruct ahc;
323  ahc.name = name;
324  ahc.type = 7;
325  ahc.value = value;
326 
327  m_VAllowedVarNameAddedByUser.push_back(ahc);
328 }
329 
330 
331 template <typename TImage>
332 void BandMathXImageFilter<TImage>::ExportContext(const std::string& filename)
333 {
334 
335  std::vector<std::string> vectI, vectF, vectM, vectFinal;
336 
337  for (unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
338  {
339  std::ostringstream iss;
340  std::string str;
341 
342  switch (m_VAllowedVarNameAddedByUser[i].value.GetType())
343  {
344  case 'i':
345  iss << "#I " << m_VAllowedVarNameAddedByUser[i].name << " " << m_VAllowedVarNameAddedByUser[i].value.GetInteger();
346  str = iss.str();
347  vectI.push_back(str);
348  break;
349  case 'f':
350  iss << "#F " << m_VAllowedVarNameAddedByUser[i].name << " " << m_VAllowedVarNameAddedByUser[i].value.GetFloat();
351  str = iss.str();
352  vectF.push_back(str);
353  break;
354  case 'c':
355  itkExceptionMacro(<< "Complex numbers not supported." << std::endl);
356  break;
357  case 'm':
358  iss << "#M " << m_VAllowedVarNameAddedByUser[i].name << " "
359  << "{";
360  for (int k = 0; k < m_VAllowedVarNameAddedByUser[i].value.GetRows(); k++)
361  {
362  iss << " " << m_VAllowedVarNameAddedByUser[i].value.At(k, 0);
363  for (int p = 1; p < m_VAllowedVarNameAddedByUser[i].value.GetCols(); p++)
364  iss << " , " << m_VAllowedVarNameAddedByUser[i].value.At(k, p);
365  iss << ";";
366  }
367  str = iss.str();
368  str.erase(str.size() - 1);
369  str.push_back('}');
370  vectM.push_back(str);
371  break;
372  }
373  }
374 
375  // Sorting : I F M and E
376  for (unsigned int i = 0; i < vectI.size(); ++i)
377  vectFinal.push_back(vectI[i]);
378  for (unsigned int i = 0; i < vectF.size(); ++i)
379  vectFinal.push_back(vectF[i]);
380  for (unsigned int i = 0; i < vectM.size(); ++i)
381  vectFinal.push_back(vectM[i]);
382  for (unsigned int i = 0; i < m_Expression.size(); ++i)
383  {
384  std::ostringstream iss;
385  iss << "#E " << m_Expression[i] << std::endl;
386  std::string str = iss.str();
387  vectFinal.push_back(str);
388  }
389 
390  std::ofstream exportFile(filename, std::ios::out | std::ios::trunc);
391  if (exportFile)
392  {
393  for (unsigned int i = 0; i < vectFinal.size(); ++i)
394  exportFile << vectFinal[i] << std::endl;
395 
396  exportFile.close();
397  }
398  else
399  itkExceptionMacro(<< "Could not open " << filename << "." << std::endl);
400 }
401 
402 template <typename TImage>
403 void BandMathXImageFilter<TImage>::ImportContext(const std::string& filename)
404 {
405  std::ifstream importFile(filename, std::ios::in);
406 
407  std::string wholeline, line, name, matrixdef;
408  int pos, pos2, lineID = 0, nbSuccesses = 0;
409  double value;
410 
411  if (importFile)
412  {
413 
414  while (std::getline(importFile, wholeline))
415  {
416  lineID++;
417 
418  pos = wholeline.find_first_not_of(' ');
419 
420  if (pos != (int)std::string::npos)
421  {
422  line = wholeline.substr(pos);
423 
424  if ((line[0] == '#') && ((line[1] == 'I') || (line[1] == 'i') || (line[1] == 'F') || (line[1] == 'f')))
425  {
426 
427  pos = line.find_first_not_of(' ', 2);
428 
429  if (pos == (int)std::string::npos)
430  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID << " : please, set the name and the value of the constant." << std::endl);
431 
432  std::string sub = line.substr(pos);
433 
434  pos = sub.find_first_of(' ');
435  name = sub.substr(0, pos);
436 
437  if (sub.find_first_of('{', pos) != std::string::npos)
438  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID
439  << " : symbol #F found, but find vector/matrix definition. Please, set an integer or a float number." << std::endl);
440 
441  if (sub.find_first_not_of(' ', pos) == std::string::npos)
442  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID << " : please, set the value of the constant." << std::endl)
443 
444  std::istringstream iss(sub.substr(pos));
445  iss >> value;
446 
447  SetConstant(name, value);
448  nbSuccesses++;
449  }
450  else if ((line[0] == '#') && ((line[1] == 'M') || (line[1] == 'm')))
451  {
452 
453  pos = line.find_first_not_of(' ', 2);
454 
455  if (pos == (int)std::string::npos)
456  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID << " : please, set the name and the definition of the vector/matrix."
457  << std::endl);
458 
459  std::string sub = line.substr(pos);
460 
461  pos = sub.find_first_of(' ');
462  name = sub.substr(0, pos);
463  pos2 = sub.find_first_of('{');
464  if (pos2 != (int)std::string::npos)
465  matrixdef = sub.substr(pos2);
466  else
467  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID << " : symbol #M found, but couldn't not find vector/matrix definition."
468  << std::endl);
469 
470  SetMatrix(name, matrixdef);
471  nbSuccesses++;
472  }
473  else if ((line[0] == '#') && ((line[1] == 'E') || (line[1] == 'e')))
474  {
475  pos = line.find_first_not_of(' ', 2);
476 
477  if (pos == (int)std::string::npos)
478  itkExceptionMacro(<< "In file '" << filename << "', line " << lineID << " : symbol #E found, but couldn't not find any expression." << std::endl);
479 
480  std::string sub = line.substr(pos);
481 
482  SetExpression(sub);
483  nbSuccesses++;
484  }
485  }
486  } // while
487 
488  importFile.close();
489 
490 
491  if (nbSuccesses == 0)
492  itkExceptionMacro(<< "No constant or expression could be set; please, ensure that the file '" << filename << "' is correct." << std::endl);
493  }
494  else
495  itkExceptionMacro(<< "Could not open " << filename << "." << std::endl);
496 }
497 
498 
499 template <typename TImage>
500 std::string BandMathXImageFilter<TImage>::GetExpression(unsigned int IDExpression) const
501 {
502  if (IDExpression < m_Expression.size())
503  return m_Expression[IDExpression];
504  return "";
505 }
506 
507 
508 template <typename TImage>
509 std::vector<std::string> BandMathXImageFilter<TImage>::GetVarNames() const
510 {
511  std::vector<std::string> res;
512  for (int y = 0; y < m_VVarName.size(); y++)
513  res.push_back(m_VVarName[y].name);
514 
515  return res;
516 }
517 
518 
519 template <typename TImage>
521 {
522  bool found = false;
523  for (unsigned int i = 0; i < m_VVarName.size(); ++i)
524  if (m_VVarName[i].name == ahc.name)
525  found = true;
526 
527  if (!found)
528  m_VVarName.push_back(ahc);
529 }
530 
531 
532 template <typename TImage>
534 {
535 
536  if (m_Expression.size() == 0)
537  itkExceptionMacro(<< "No expression set; please set at least one expression." << std::endl);
538 
539 
540  // Generate variables names
541  m_VVarName.clear();
542  m_VNotAllowedVarName.clear();
543  m_VFinalAllowedVarName.clear();
544 
545  // m_VFinalAllowedVarName = m_VAllowedVarNameAuto + m_VAllowedVarNameAddedByUser
546  // m_VFinalAllowedVarName = variable names dictionary
547  for (unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
548  m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAddedByUser[i]);
549  for (unsigned int i = 0; i < m_VAllowedVarNameAuto.size(); i++)
550  m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAuto[i]);
551 
552  unsigned int nbExpr = m_Expression.size();
553  for (unsigned int IDExpression = 0; IDExpression < nbExpr; ++IDExpression) // For each expression
554  {
555  ParserType::Pointer dummyParser = ParserType::New();
556  dummyParser->SetExpr(this->GetExpression(IDExpression));
557 
558  mup::var_maptype vmap = dummyParser->GetExprVar();
559  for (mup::var_maptype::iterator item = vmap.begin(); item != vmap.end(); ++item)
560  {
561  bool OK = false;
562  int i = 0;
563  while ((!OK) && (i < (int)m_VFinalAllowedVarName.size()))
564  {
565  if (item->first == m_VFinalAllowedVarName[i].name)
566  OK = true;
567  else
568  i++;
569  }
570 
571  if (OK)
572  {
573  AddVariable(m_VFinalAllowedVarName[i]);
574  }
575  else
576  {
577  adhocStruct ahc;
578  ahc.name = item->first;
579  m_VNotAllowedVarName.push_back(ahc);
580  }
581  }
582 
583  } // At this step, m_VVarName has been built
584 
585 
586  // Checking formulas consistency
587  if (m_VNotAllowedVarName.size() > 0)
588  {
589  std::stringstream sstm;
590  sstm << "Following variables not allowed : ";
591  for (unsigned int i = 0; i < m_VNotAllowedVarName.size(); ++i)
592  sstm << m_VNotAllowedVarName[i].name << " ";
593  sstm << std::endl;
594  itkExceptionMacro(<< sstm.str());
595  }
596 
597 
598  // Register variables for each parser (important : one parser per thread and per expression)
599  m_VParser.clear();
600  unsigned int nbThreads = this->GetNumberOfThreads();
601  for (unsigned int k = 0; k < nbThreads; k++)
602  {
603  std::vector<ParserType::Pointer> parserList;
604  for (unsigned int i = 0; i < nbExpr; i++)
605  {
606  parserList.push_back(ParserType::New());
607  }
608  m_VParser.push_back(parserList);
609  }
610 
611  // Important to remember that variables of m_VVarName come from a call of GetExprVar method
612  // Only useful variables are allocated in this filter
613  int nbVar = m_VVarName.size();
614 
615  m_StatsVarDetected.clear();
616  m_NeighDetected.clear();
617  m_NeighExtremaSizes.clear();
618  unsigned int nbInputImages = this->GetNumberOfInputs();
619  RadiusType dummyRadius;
620  dummyRadius[0] = 1;
621  dummyRadius[1] = 1;
622  m_NeighExtremaSizes.resize(nbInputImages, dummyRadius);
623 
624  // Reset
625  for (unsigned int i = 0; i < m_AImage.size(); ++i)
626  m_AImage[i].clear();
627  m_AImage.clear();
628 
629  m_AImage.resize(nbThreads);
630 
631  double initValue = 0.1;
632  for (unsigned int i = 0; i < nbThreads; ++i) // For each thread
633  {
634 
635  m_AImage[i].resize(nbVar); // For each variable
636 
637  for (int j = 0; j < nbVar; ++j)
638  {
639  m_AImage[i][j].name = m_VVarName[j].name;
640  m_AImage[i][j].type = m_VVarName[j].type;
641  for (int t = 0; t < 5; ++t)
642  m_AImage[i][j].info[t] = m_VVarName[j].info[t];
643 
644 
645  if ((m_AImage[i][j].type == 0) || (m_AImage[i][j].type == 1)) // indices (idxX & idxY)
646  {
647  m_AImage[i][j].value = ValueType(initValue);
648  }
649 
650  if (m_AImage[i][j].type == 2) // imiPhyX
651  {
652  SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
653  m_AImage[i][j].value = ValueType(static_cast<double>(spacing[0]));
654  }
655 
656  if (m_AImage[i][j].type == 3) // imiPhyY
657  {
658  SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
659  m_AImage[i][j].value = ValueType(static_cast<double>(spacing[1]));
660  }
661 
662  if (m_AImage[i][j].type == 4) // vector
663  {
664  unsigned int nbBands = this->GetNthInput(m_AImage[i][j].info[0])->GetNumberOfComponentsPerPixel();
665  m_AImage[i][j].value = ValueType(1, nbBands, initValue);
666  }
667 
668  if (m_AImage[i][j].type == 5) // pixel
669  {
670  m_AImage[i][j].value = ValueType(initValue);
671  }
672 
673  if (m_AImage[i][j].type == 6) // neighborhood
674  {
675  m_AImage[i][j].value = ValueType(m_AImage[i][j].info[3], m_AImage[i][j].info[2], initValue);
676 
677  // m_AImage[i][j].info[0] = Image ID
678  bool found = false;
679  for (unsigned int r = 0; r < m_NeighDetected.size() && !found; r++)
680  if (m_NeighDetected[r] == (unsigned int)m_AImage[i][j].info[0])
681  found = true;
682  if (!found)
683  m_NeighDetected.push_back(m_AImage[i][j].info[0]);
684 
685  // find biggest radius for a given input image (idis given by info[0])
686  if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] < (unsigned int)((m_VVarName[j].info[2] - 1) / 2)) // Size x direction (otb convention)
687  m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] = (unsigned int)((m_VVarName[j].info[2] - 1) / 2);
688 
689  if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] < (unsigned int)((m_VVarName[j].info[3] - 1) / 2)) // Size y direction (otb convention)
690  m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] = (unsigned int)((m_VVarName[j].info[3] - 1) / 2);
691  }
692 
693  if (m_AImage[i][j].type == 7) // user defined variables
694  {
695 
696  for (int t = 0; t < (int)m_VAllowedVarNameAddedByUser.size(); t++)
697  if (m_VAllowedVarNameAddedByUser[t].name.compare(m_AImage[i][j].name) == 0)
698  m_AImage[i][j].value = m_VAllowedVarNameAddedByUser[t].value;
699  }
700 
701  if (m_AImage[i][j].type == 8) // global stats
702  {
703  m_AImage[i][j].value = ValueType(initValue);
704  // m_AImage[i][j].info[0] = Image ID
705  bool found = false;
706  for (unsigned int r = 0; r < m_StatsVarDetected.size() && !found; r++)
707  if (m_StatsVarDetected[r] == m_AImage[i][j].info[0])
708  found = true;
709  if (!found)
710  m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
711  }
712 
713 
714  // Register variable
715  for (unsigned int k = 0; k < nbExpr; k++)
716  {
717  m_VParser[i][k]->DefineVar(m_AImage[i][j].name, &(m_AImage[i][j].value));
718  }
719 
720 
721  initValue += 0.001;
722  if (initValue > 1.0)
723  initValue = 0.1;
724  }
725  }
726 
727  // Set expressions
728  for (unsigned int k = 0; k < nbThreads; k++)
729  {
730  for (unsigned int i = 0; i < nbExpr; i++)
731  {
732  m_VParser[k][i]->SetExpr(m_Expression[i]);
733  }
734  }
735 }
736 
737 
738 template <typename TImage>
740 {
741  // Must instantiate stats variables of the parsers
742  // Note : at this stage, inputs have already been set to largest possible regions.
743  for (unsigned int i = 0; i < m_StatsVarDetected.size(); i++)
744  {
745 
746  StreamingStatisticsVectorImageFilterPointerType filter = StreamingStatisticsVectorImageFilterType::New();
747 
748  filter->SetInput(this->GetNthInput(m_StatsVarDetected[i]));
749  filter->Update();
750 
751  PixelType pix; // Variable length vector
752  MatrixType mat;
753 
754  for (unsigned int t = 0; t < m_AImage.size(); t++) // for each thread
755  for (unsigned int v = 0; v < m_AImage[t].size(); v++) // for each variable
756  if ((m_AImage[t][v].type == 8) && (m_AImage[t][v].info[0] == m_StatsVarDetected[i])) // type 8 : flag identifying a glob stat; info[0] : input ID
757  {
758  switch (m_AImage[t][v].info[2]) // info[2] sub-type (see also SetNthInput method above)
759  {
760  case 0: // mini
761 
762  pix = filter->GetMinimum();
763  for (int b = 0; b < (int)pix.GetSize(); b++) // for each band
764  if (m_AImage[t][v].info[1] == b) // info[1] : band ID
765  m_AImage[t][v].value = pix[b];
766 
767  break;
768 
769  case 1: // maxi
770 
771  pix = filter->GetMaximum();
772  for (int b = 0; b < (int)pix.GetSize(); b++) // for each band
773  if (m_AImage[t][v].info[1] == b) // info[1] : band ID
774  m_AImage[t][v].value = pix[b];
775 
776  break;
777 
778  case 2: // mean
779 
780  pix = filter->GetMean();
781  for (int b = 0; b < (int)pix.GetSize(); b++) // for each band
782  if (m_AImage[t][v].info[1] == b) // info[1] : band ID
783  m_AImage[t][v].value = pix[b];
784 
785  break;
786 
787  break;
788 
789  case 3: // sum
790 
791  pix = filter->GetSum();
792  for (int b = 0; b < (int)pix.GetSize(); b++) // for each band
793  if (m_AImage[t][v].info[1] == b) // info[1] : band ID
794  m_AImage[t][v].value = pix[b];
795 
796  break;
797 
798  case 4: // stddev
799 
800  mat = filter->GetCovariance();
801  for (int b = 0; b < (int)mat.Cols(); b++) // for each band
802  if (m_AImage[t][v].info[1] == b) // info[1] : band ID
803  m_AImage[t][v].value = mat(b, b);
804 
805  break;
806  }
807  }
808  }
809 }
810 
811 template <typename TImage>
813 {
814 
815  this->SetNumberOfRequiredOutputs((int)m_Expression.size());
816 
817  m_outputsDimensions.clear();
818 
819  for (int i = 0; i < (int)m_Expression.size(); ++i)
820  {
821  ValueType value = m_VParser[0][i]->EvalRef();
822 
823  switch (value.GetType())
824  { // ValueType
825  case 'b':
826  itkExceptionMacro(<< "Booleans not supported." << std::endl);
827  break;
828 
829  case 'i':
830  m_outputsDimensions.push_back(1);
831  break;
832 
833  case 'f':
834  m_outputsDimensions.push_back(1);
835  break;
836 
837  case 'c':
838  itkExceptionMacro(<< "Complex numbers not supported." << std::endl);
839  break;
840 
841  case 'm':
842  {
843  const mup::matrix_type& vect = value.GetArray();
844  if (vect.GetRows() == 1) // Vector
845  m_outputsDimensions.push_back(vect.GetCols());
846  else // Matrix
847  itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
848  }
849  break;
850 
851  default:
852  itkExceptionMacro(<< "Unknown output type : " << value.GetType() << std::endl);
853  break;
854  }
855 
856  // std::cout << "Type = " << value.GetType() << " dimension = " << m_outputsDimensions.back() << std::endl;
857  }
858 }
859 
860 template <typename TImage>
862 {
863  // Check if input image dimensions match
864  unsigned int nbInputImages = this->GetNumberOfInputs();
865  unsigned int inputSize[2];
866  inputSize[0] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(0);
867  inputSize[1] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(1);
868 
869  for (unsigned int p = 1; p < nbInputImages; p++)
870  if ((inputSize[0] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0)) ||
871  (inputSize[1] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1)))
872  {
873 
874  itkExceptionMacro(<< "Input images must have the same dimensions." << std::endl
875  << "band #1 is [" << inputSize[0] << ";" << inputSize[1] << "]" << std::endl
876  << "band #" << p + 1 << " is [" << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0) << ";"
877  << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1) << "]");
878  }
879 }
880 
881 template <typename TImage>
883 {
884  Superclass::GenerateOutputInformation();
885 
886 
887  CheckImageDimensions();
888  PrepareParsers();
889  if (GlobalStatsDetected())
890  PrepareParsersGlobStats();
891  OutputsDimensions();
892 
893 
894  typedef itk::ImageBase<TImage::ImageDimension> ImageBaseType;
895  typename ImageBaseType::Pointer outputPtr;
896 
897  int i = 0;
898  for (itk::OutputDataObjectIterator it(this); !it.IsAtEnd(); i++, it++)
899  {
900  // Check whether the output is an image of the appropriate
901  // dimension (use ProcessObject's version of the GetInput()
902  // method since it returns the input as a pointer to a
903  // DataObject as opposed to the subclass version which
904  // static_casts the input to an TImage).
905  outputPtr = dynamic_cast<ImageBaseType*>(it.GetOutput());
906 
907  if (outputPtr)
908  outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]);
909  }
910 }
911 
912 
913 template <typename TImage>
915 {
916  // call the superclass' implementation of this method
917  Superclass::GenerateInputRequestedRegion();
918 
919  for (unsigned int i = 0; i < m_NeighDetected.size(); i++)
920  if (m_NeighDetected[i] < this->GetNumberOfInputs())
921  {
922 
923  // get pointers to the input and output
924  typename Superclass::InputImagePointer inputPtr = const_cast<TImage*>(this->GetNthInput(m_NeighDetected[i]));
925 
926 
927  ImageRegionType inputRequestedRegion;
928  inputRequestedRegion = inputPtr->GetRequestedRegion();
929 
930  // pad the input requested region by the operator radius
931  inputRequestedRegion.PadByRadius(m_NeighExtremaSizes[m_NeighDetected[i]]);
932 
933  // crop the input requested region at the input's largest possible region
934  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
935  {
936  inputPtr->SetRequestedRegion(inputRequestedRegion);
937  return;
938  }
939  else
940  {
941  // Couldn't crop the region (requested region is outside the largest
942  // possible region). Throw an exception.
943 
944  // store what we tried to request (prior to trying to crop)
945  inputPtr->SetRequestedRegion(inputRequestedRegion);
946 
947  // build an exception
948  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
949  std::ostringstream msg, msg2;
950  msg << static_cast<const char*>(this->GetNameOfClass()) << "::GenerateInputRequestedRegion()";
951  e.SetLocation(msg.str());
952  msg2 << "Requested region is (at least partially) outside the largest possible region (input #" << m_NeighDetected[i] << ").";
953  e.SetDescription(msg2.str());
954  e.SetDataObject(inputPtr);
955  throw e;
956  }
957  }
958  else
959  itkExceptionMacro(<< "Requested input #" << m_NeighDetected[i] << ", but only " << this->GetNumberOfInputs() << " inputs are available." << std::endl);
960 }
961 
962 
963 template <typename TImage>
965 {
966 
967  unsigned int nbThreads = this->GetNumberOfThreads();
968  // Allocate and initialize the thread temporaries
969  m_ThreadUnderflow.SetSize(nbThreads);
970  m_ThreadUnderflow.Fill(0);
971  m_ThreadOverflow.SetSize(nbThreads);
972  m_ThreadOverflow.Fill(0);
973 }
974 
975 
976 template <typename TImage>
978 {
979  unsigned int nbThreads = this->GetNumberOfThreads();
980  unsigned int i;
981 
982  m_UnderflowCount = 0;
983  m_OverflowCount = 0;
984 
985  // Accumulate counts for each thread
986  for (i = 0; i < nbThreads; ++i)
987  {
988  m_UnderflowCount += m_ThreadUnderflow[i];
989  m_OverflowCount += m_ThreadOverflow[i];
990  }
991 
992  if ((m_UnderflowCount != 0) || (m_OverflowCount != 0))
993  {
994  std::stringstream sstm;
995  sstm << std::endl << "The Following Parsed Expression : ";
996  for (unsigned int t = 0; t < m_Expression.size(); ++t)
997  sstm << this->GetExpression(t) << std::endl;
998  sstm << "Generated " << m_UnderflowCount << " Underflow(s) "
999  << "And " << m_OverflowCount << " Overflow(s) " << std::endl
1000  << "The Parsed Expression, The Inputs And The Output "
1001  << "Type May Be Incompatible !";
1002 
1003  otbWarningMacro(<< sstm.str());
1004  }
1005 }
1006 
1007 template <typename TImage>
1008 void BandMathXImageFilter<TImage>::ThreadedGenerateData(const ImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
1009 {
1010 
1011  ValueType value;
1012  unsigned int nbInputImages = this->GetNumberOfInputs();
1013 
1014  //----------------- --------- -----------------//
1015  //----------------- Iterators -----------------//
1016  //----------------- --------- -----------------//
1017  typedef itk::ImageScanlineConstIterator<TImage> ImageScanlineConstIteratorType;
1018  typedef itk::ImageScanlineIterator<TImage> ImageScanlineIteratorType;
1019  typedef itk::ImageRegionConstIteratorWithOnlyIndex<TImage> IndexIteratorType;
1020  std::vector<ImageScanlineConstIteratorType> Vit;
1021  Vit.resize(nbInputImages);
1022  for (unsigned int j = 0; j < nbInputImages; ++j)
1023  Vit[j] = ImageScanlineConstIteratorType(this->GetNthInput(j), outputRegionForThread);
1024 
1025 
1026  std::vector<ImageScanlineIteratorType> VoutIt;
1027  VoutIt.resize(m_Expression.size());
1028  for (unsigned int j = 0; j < VoutIt.size(); ++j)
1029  VoutIt[j] = ImageScanlineIteratorType(this->GetOutput(j), outputRegionForThread);
1030 
1031 
1032  // Special case : neighborhoods
1033  std::vector<itk::ConstNeighborhoodIterator<TImage>> VNit;
1034  for (unsigned int j = 0; j < m_VVarName.size(); ++j)
1035  if (m_VVarName[j].type == 6)
1036  {
1037  RadiusType radius;
1038  radius[0] = (int)((m_VVarName[j].info[2] - 1) / 2); // Size x direction (otb convention)
1039  radius[1] = (int)((m_VVarName[j].info[3] - 1) / 2); // Size y direction (otb convention)
1040  VNit.push_back(
1041  itk::ConstNeighborhoodIterator<TImage>(radius, this->GetNthInput(m_VVarName[j].info[0]), outputRegionForThread)); // info[0] = Input image ID
1042  VNit.back().NeedToUseBoundaryConditionOn();
1043  }
1044 
1045  // Index only iterator
1046  IndexIteratorType indexIterator(this->GetNthInput(0), outputRegionForThread);
1047 
1048  // Support progress methods/callbacks
1049  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
1050 
1051  // iterator on variables
1052  typename std::vector<adhocStruct>::iterator iterVarStart = m_AImage[threadId].begin();
1053  typename std::vector<adhocStruct>::iterator iterVarEnd = m_AImage[threadId].end();
1054  typename std::vector<adhocStruct>::iterator iterVar = iterVarStart;
1055 
1056  // temporary output vectors
1057  std::vector<PixelType> tmpOutputs(m_Expression.size());
1058  for (unsigned int k = 0; k < m_Expression.size(); ++k)
1059  tmpOutputs[k].SetSize(m_outputsDimensions[k]);
1060 
1061  //----------------- --------------------- -----------------//
1062  //----------------- Variable affectations -----------------//
1063  //----------------- --------------------- -----------------//
1064  for (unsigned int j = 0; j < nbInputImages; ++j)
1065  {
1066  Vit[j].GoToBegin();
1067  }
1068  for (unsigned int j = 0; j < m_Expression.size(); ++j)
1069  {
1070  VoutIt[j].GoToBegin();
1071  }
1072  for (unsigned int j = 0; j < VNit.size(); ++j)
1073  {
1074  VNit[j].GoToBegin();
1075  }
1076  indexIterator.GoToBegin();
1077 
1078  while (!Vit[0].IsAtEnd()) // For each pixel
1079  {
1080  while (!Vit[0].IsAtEndOfLine()) // For each line
1081  {
1082  int ngbhNameIndex = 0;
1083  int index;
1084 
1085  iterVar = iterVarStart;
1086  while (iterVar != iterVarEnd)
1087  {
1088  switch (iterVar->type)
1089  {
1090  case 0: // idxX
1091  iterVar->value = static_cast<double>(indexIterator.GetIndex()[0]);
1092  break;
1093 
1094  case 1: // idxY
1095  iterVar->value = static_cast<double>(indexIterator.GetIndex()[1]);
1096  break;
1097 
1098  case 2: // Spacing X (imiPhyX)
1099  // Nothing to do (already set inside BeforeThreadedGenerateData)"
1100  break;
1101 
1102  case 3: // Spacing Y (imiPhyY)
1103  // Nothing to do (already set inside BeforeThreadedGenerateData)"
1104  break;
1105 
1106  case 4: // vector
1107  // iterVar->info[0] : Input image #ID
1108  for (int p = 0; p < iterVar->value.GetCols(); ++p)
1109  iterVar->value.At(0, p) = Vit[iterVar->info[0]].Get()[p];
1110  break;
1111 
1112  case 5: // pixel
1113  // iterVar->info[0] : Input image #ID
1114  // iterVar->info[1] : Band #ID
1115  iterVar->value = Vit[iterVar->info[0]].Get()[iterVar->info[1]];
1116  break;
1117 
1118  case 6: // neighborhood
1119 
1120  // iterVar->info[1] : Band #ID
1121  if (iterVar->info[2] * iterVar->info[3] != (int)VNit[ngbhNameIndex].Size())
1122  itkExceptionMacro(<< "Size of muparserx variable is different from its related otb neighborhood iterator")
1123 
1124  index = 0;
1125  for (int rows = 0; rows < iterVar->info[3]; ++rows)
1126  for (int cols = 0; cols < iterVar->info[2]; ++cols)
1127  {
1128  iterVar->value.At(rows, cols) = VNit[ngbhNameIndex].GetPixel(index)[iterVar->info[1]];
1129  index++;
1130  }
1131 
1132  ngbhNameIndex++;
1133  break;
1134 
1135  case 7:
1136  // Nothing to do : user defined variable or constant, which have already been set inside PrepareParsers (see above)
1137  break;
1138 
1139  case 8:
1140  // Nothing to do : variable has already been set inside PrepareParsersGlobStats method (see above)
1141  break;
1142 
1143  default:
1144  itkExceptionMacro(<< "Type of the variable is unknown");
1145  break;
1146  }
1147 
1148  iterVar++;
1149  } // End while on vars
1150 
1151  //----------------- ----------- -----------------//
1152  //----------------- Evaluations -----------------//
1153  //----------------- ----------- -----------------//
1154  for (unsigned int IDExpression = 0; IDExpression < m_Expression.size(); ++IDExpression)
1155  {
1156  value = m_VParser[threadId][IDExpression]->EvalRef();
1157 
1158  switch (value.GetType())
1159  { // ValueType
1160  case 'i':
1161  tmpOutputs[IDExpression][0] = value.GetInteger();
1162  break;
1163 
1164  case 'f':
1165  tmpOutputs[IDExpression][0] = value.GetFloat();
1166  break;
1167 
1168  case 'c':
1169  itkExceptionMacro(<< "Complex numbers are not supported." << std::endl);
1170  break;
1171 
1172  case 'm':
1173  {
1174  const mup::matrix_type& vect = value.GetArray();
1175 
1176  if (vect.GetRows() == 1) // Vector
1177  for (int p = 0; p < vect.GetCols(); ++p)
1178  tmpOutputs[IDExpression][p] = vect.At(0, p).GetFloat();
1179  else // Matrix
1180  itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
1181  }
1182  break;
1183  }
1184 
1185  //----------------- Pixel affectations -----------------//
1186  for (unsigned int p = 0; p < m_outputsDimensions[IDExpression]; ++p)
1187  {
1188  // Case value is equal to -inf or inferior to the minimum value
1189  // allowed by the PixelValueType cast
1190  if (tmpOutputs[IDExpression][p] < double(itk::NumericTraits<PixelValueType>::NonpositiveMin()))
1191  {
1192  tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::NonpositiveMin();
1193  m_ThreadUnderflow[threadId]++;
1194  }
1195  // Case value is equal to inf or superior to the maximum value
1196  // allowed by the PixelValueType cast
1197  else if (tmpOutputs[IDExpression][p] > double(itk::NumericTraits<PixelValueType>::max()))
1198  {
1199  tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::max();
1200  m_ThreadOverflow[threadId]++;
1201  }
1202  }
1203  VoutIt[IDExpression].Set(tmpOutputs[IDExpression]);
1204  }
1205 
1206  for (unsigned int j = 0; j < nbInputImages; ++j)
1207  {
1208  ++Vit[j];
1209  }
1210  for (unsigned int j = 0; j < m_Expression.size(); ++j)
1211  {
1212  ++VoutIt[j];
1213  }
1214  for (unsigned int j = 0; j < VNit.size(); ++j)
1215  {
1216  ++VNit[j];
1217  }
1218  ++indexIterator;
1219 
1220  progress.CompletedPixel();
1221  }
1222  for (unsigned int j = 0; j < nbInputImages; ++j)
1223  {
1224  Vit[j].NextLine();
1225  }
1226  for (unsigned int j = 0; j < m_Expression.size(); ++j)
1227  {
1228  VoutIt[j].NextLine();
1229  }
1230  }
1231 }
1232 
1233 } // end namespace otb
1234 
1235 #endif
StreamingStatisticsVectorImageFilterType::MatrixType MatrixType
void SetExpression(const std::string &expression)
static constexpr char const * name() noexcept
itk::SmartPointer< Self > Pointer
Definition: otbParserX.h:65
void SetConstant(const std::string &name, double value)
void PrintSelf(std::ostream &os, itk::Indent indent) const override
void SetMatrix(const std::string &name, const std::string &definition)
void SetNthInput(DataObjectPointerArraySizeType idx, const ImageType *image)
ParserType::ValueType ValueType
std::string GetExpression(unsigned int IDExpression) const
void ExportContext(const std::string &filename)
StreamingStatisticsVectorImageFilterType::Pointer StreamingStatisticsVectorImageFilterPointerType
itk::ImageBase< 2 > ImageBaseType
void ThreadedGenerateData(const ImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
itk::ConstNeighborhoodIterator< TImage >::RadiusType RadiusType
ImageType::RegionType ImageRegionType
ImageType::SpacingType SpacingType
ImageType * GetNthInput(DataObjectPointerArraySizeType idx)
#define otbWarningMacro(x)
Definition: otbMacro.h:65
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
std::vector< std::string > GetVarNames() const
static constexpr GLenum value() noexcept
void ImportContext(const std::string &filename)
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType