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