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