22 #ifndef otbBandMathXImageFilter_hxx
23 #define otbBandMathXImageFilter_hxx
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"
44 template <
class TImage>
49 this->SetNumberOfRequiredInputs(1);
53 m_ThreadUnderflow.SetSize(1);
54 m_ThreadOverflow.SetSize(1);
61 m_VAllowedVarNameAuto.push_back(ahcX);
66 m_VAllowedVarNameAuto.push_back(ahcY);
68 m_SizeNeighbourhood = 10;
70 m_ManyExpressions =
true;
74 template <
class TImage>
81 for (
unsigned int i = 0; i < m_AImage.size(); ++i)
86 m_VAllowedVarNameAuto.clear();
87 m_VAllowedVarNameAddedByUser.clear();
88 m_VFinalAllowedVarName.clear();
89 m_VNotAllowedVarName.clear();
90 m_outputsDimensions.clear();
94 template <
class TImage>
97 Superclass::PrintSelf(os, indent);
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;
109 template <
class TImage>
113 std::stringstream sstm;
114 sstm <<
"im" << (idx + 1);
115 this->SetNthInput(idx, image, sstm.str());
118 template <
class TImage>
123 this->SetInput(idx, imagebis);
126 std::stringstream sstmPhyX;
128 sstmPhyX << varName <<
"PhyX";
129 ahcPhyX.
name = sstmPhyX.str();
131 ahcPhyX.
info[0] = idx;
132 m_VAllowedVarNameAuto.push_back(ahcPhyX);
134 std::stringstream sstmPhyY;
136 sstmPhyY << varName <<
"PhyY";
137 ahcPhyY.
name = sstmPhyY.str();
139 ahcPhyY.
info[0] = idx;
140 m_VAllowedVarNameAuto.push_back(ahcPhyY);
143 std::stringstream sstm_glob;
145 sstm_glob << varName;
146 ahc_glob.
name = sstm_glob.str();
148 ahc_glob.
info[0] = idx;
149 m_VAllowedVarNameAuto.push_back(ahc_glob);
154 imagebis->UpdateOutputInformation();
157 for (
unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
159 std::stringstream sstm;
161 sstm << varName <<
"b" << (j + 1);
162 ahc.
name = sstm.str();
166 m_VAllowedVarNameAuto.push_back(ahc);
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++)
174 std::stringstream sstm;
176 sstm << varName <<
"b" << (j + 1) <<
"N" << 2 * x + 1 <<
"x" << 2 * y + 1;
177 ahc.
name = sstm.str();
181 ahc.
info[2] = 2 * x + 1;
182 ahc.
info[3] = 2 * y + 1;
183 m_VAllowedVarNameAuto.push_back(ahc);
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");
194 for (
unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
195 for (
unsigned int t = 0; t < statsNames.size(); t++)
197 std::stringstream sstm;
199 sstm << varName <<
"b" << (j + 1) << statsNames[t];
200 ahc.
name = sstm.str();
205 m_VAllowedVarNameAuto.push_back(ahc);
209 template <
typename TImage>
212 return const_cast<TImage*
>(this->GetInput(idx));
215 template <
typename TImage>
218 m_ManyExpressions = flag;
221 template <
typename TImage>
224 std::string expressionToBePushed = expression;
226 if (expression.find(
";") != std::string::npos)
228 std::ostringstream oss;
230 for (
unsigned int i = 0; i < expression.size(); ++i)
231 if (expression[i] ==
';')
234 oss << expression[i];
237 expressionToBePushed = oss.str();
240 if (m_ManyExpressions)
241 m_Expression.push_back(expressionToBePushed);
242 else if (m_Expression.size() == 0)
243 m_Expression.push_back(expressionToBePushed);
245 if (m_Expression.size() > 1)
251 template <
typename TImage>
254 m_Expression.clear();
257 template <
typename TImage>
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);
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);
271 for (
unsigned int i = 1; i < definition.size() - 1; ++i)
272 def.push_back(definition[i]);
275 std::vector<std::vector<double>> mat;
276 std::istringstream iss(def);
278 while (std::getline(iss, rows,
';'))
280 mat.push_back(std::vector<double>(0));
281 std::istringstream iss2(rows);
283 while (std::getline(iss2, elmt,
','))
285 std::istringstream iss3(elmt);
288 mat.back().push_back(val);
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);
303 ahc.
info[0] = mat[0].size();
304 ahc.
info[1] = mat.size();
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];
311 m_VAllowedVarNameAddedByUser.push_back(ahc);
315 template <
typename TImage>
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);
327 m_VAllowedVarNameAddedByUser.push_back(ahc);
331 template <
typename TImage>
335 std::vector<std::string> vectI, vectF, vectM, vectFinal;
337 for (
unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
339 std::ostringstream iss;
342 switch (m_VAllowedVarNameAddedByUser[i].value.GetType())
345 iss <<
"#I " << m_VAllowedVarNameAddedByUser[i].name <<
" " << m_VAllowedVarNameAddedByUser[i].value.GetInteger();
347 vectI.push_back(str);
350 iss <<
"#F " << m_VAllowedVarNameAddedByUser[i].name <<
" " << m_VAllowedVarNameAddedByUser[i].value.GetFloat();
352 vectF.push_back(str);
355 itkExceptionMacro(<<
"Complex numbers not supported." << std::endl);
358 iss <<
"#M " << m_VAllowedVarNameAddedByUser[i].name <<
" "
360 for (
int k = 0; k < m_VAllowedVarNameAddedByUser[i].value.GetRows(); k++)
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);
368 str.erase(str.size() - 1);
370 vectM.push_back(str);
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)
384 std::ostringstream iss;
385 iss <<
"#E " << m_Expression[i] << std::endl;
386 std::string str = iss.str();
387 vectFinal.push_back(str);
390 std::ofstream exportFile(filename, std::ios::out | std::ios::trunc);
393 for (
unsigned int i = 0; i < vectFinal.size(); ++i)
394 exportFile << vectFinal[i] << std::endl;
399 itkExceptionMacro(<<
"Could not open " << filename <<
"." << std::endl);
402 template <
typename TImage>
405 std::ifstream importFile(filename, std::ios::in);
407 std::string wholeline, line, name, matrixdef;
408 int pos, pos2, lineID = 0, nbSuccesses = 0;
414 while (std::getline(importFile, wholeline))
418 pos = wholeline.find_first_not_of(
' ');
420 if (pos != (
int)std::string::npos)
422 line = wholeline.substr(pos);
424 if ((line[0] ==
'#') && ((line[1] ==
'I') || (line[1] ==
'i') || (line[1] ==
'F') || (line[1] ==
'f')))
427 pos = line.find_first_not_of(
' ', 2);
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);
432 std::string sub = line.substr(pos);
434 pos = sub.find_first_of(
' ');
435 name = sub.substr(0, pos);
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);
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)
444 std::istringstream iss(sub.substr(pos));
447 SetConstant(name, value);
450 else if ((line[0] ==
'#') && ((line[1] ==
'M') || (line[1] ==
'm')))
453 pos = line.find_first_not_of(
' ', 2);
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."
459 std::string sub = line.substr(pos);
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);
467 itkExceptionMacro(<<
"In file '" << filename <<
"', line " << lineID <<
" : symbol #M found, but couldn't not find vector/matrix definition."
470 SetMatrix(name, matrixdef);
473 else if ((line[0] ==
'#') && ((line[1] ==
'E') || (line[1] ==
'e')))
475 pos = line.find_first_not_of(
' ', 2);
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);
480 std::string sub = line.substr(pos);
491 if (nbSuccesses == 0)
492 itkExceptionMacro(<<
"No constant or expression could be set; please, ensure that the file '" << filename <<
"' is correct." << std::endl);
495 itkExceptionMacro(<<
"Could not open " << filename <<
"." << std::endl);
499 template <
typename TImage>
502 if (IDExpression < m_Expression.size())
503 return m_Expression[IDExpression];
508 template <
typename TImage>
511 std::vector<std::string> res;
512 for (
int y = 0; y < m_VVarName.size(); y++)
513 res.push_back(m_VVarName[y].name);
519 template <
typename TImage>
523 for (
unsigned int i = 0; i < m_VVarName.size(); ++i)
524 if (m_VVarName[i].name == ahc.
name)
528 m_VVarName.push_back(ahc);
532 template <
typename TImage>
536 if (m_Expression.size() == 0)
537 itkExceptionMacro(<<
"No expression set; please set at least one expression." << std::endl);
542 m_VNotAllowedVarName.clear();
543 m_VFinalAllowedVarName.clear();
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]);
552 unsigned int nbExpr = m_Expression.size();
553 for (
unsigned int IDExpression = 0; IDExpression < nbExpr; ++IDExpression)
556 dummyParser->SetExpr(this->GetExpression(IDExpression));
558 mup::var_maptype vmap = dummyParser->GetExprVar();
559 for (mup::var_maptype::iterator item = vmap.begin(); item != vmap.end(); ++item)
563 while ((!OK) && (i < (
int)m_VFinalAllowedVarName.size()))
565 if (item->first == m_VFinalAllowedVarName[i].name)
573 AddVariable(m_VFinalAllowedVarName[i]);
578 ahc.
name = item->first;
579 m_VNotAllowedVarName.push_back(ahc);
587 if (m_VNotAllowedVarName.size() > 0)
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 <<
" ";
594 itkExceptionMacro(<< sstm.str());
600 unsigned int nbThreads = this->GetNumberOfThreads();
601 for (
unsigned int k = 0; k < nbThreads; k++)
603 std::vector<ParserType::Pointer> parserList;
604 for (
unsigned int i = 0; i < nbExpr; i++)
606 parserList.push_back(ParserType::New());
608 m_VParser.push_back(parserList);
613 int nbVar = m_VVarName.size();
615 m_StatsVarDetected.clear();
616 m_NeighDetected.clear();
617 m_NeighExtremaSizes.clear();
618 unsigned int nbInputImages = this->GetNumberOfInputs();
622 m_NeighExtremaSizes.resize(nbInputImages, dummyRadius);
625 for (
unsigned int i = 0; i < m_AImage.size(); ++i)
629 m_AImage.resize(nbThreads);
631 double initValue = 0.1;
632 for (
unsigned int i = 0; i < nbThreads; ++i)
635 m_AImage[i].resize(nbVar);
637 for (
int j = 0; j < nbVar; ++j)
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];
645 if ((m_AImage[i][j].type == 0) || (m_AImage[i][j].type == 1))
647 m_AImage[i][j].value =
ValueType(initValue);
650 if (m_AImage[i][j].type == 2)
652 SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
653 m_AImage[i][j].value =
ValueType(
static_cast<double>(spacing[0]));
656 if (m_AImage[i][j].type == 3)
658 SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
659 m_AImage[i][j].value =
ValueType(
static_cast<double>(spacing[1]));
662 if (m_AImage[i][j].type == 4)
664 unsigned int nbBands = this->GetNthInput(m_AImage[i][j].info[0])->GetNumberOfComponentsPerPixel();
665 m_AImage[i][j].value =
ValueType(1, nbBands, initValue);
668 if (m_AImage[i][j].type == 5)
670 m_AImage[i][j].value =
ValueType(initValue);
673 if (m_AImage[i][j].type == 6)
675 m_AImage[i][j].value =
ValueType(m_AImage[i][j].info[3], m_AImage[i][j].info[2], initValue);
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])
683 m_NeighDetected.push_back(m_AImage[i][j].info[0]);
686 if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] < (
unsigned int)((m_VVarName[j].info[2] - 1) / 2))
687 m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] = (
unsigned int)((m_VVarName[j].info[2] - 1) / 2);
689 if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] < (
unsigned int)((m_VVarName[j].info[3] - 1) / 2))
690 m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] = (
unsigned int)((m_VVarName[j].info[3] - 1) / 2);
693 if (m_AImage[i][j].type == 7)
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;
701 if (m_AImage[i][j].type == 8)
703 m_AImage[i][j].value =
ValueType(initValue);
706 for (
unsigned int r = 0; r < m_StatsVarDetected.size() && !found; r++)
707 if (m_StatsVarDetected[r] == m_AImage[i][j].info[0])
710 m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
715 for (
unsigned int k = 0; k < nbExpr; k++)
717 m_VParser[i][k]->DefineVar(m_AImage[i][j].name, &(m_AImage[i][j].value));
728 for (
unsigned int k = 0; k < nbThreads; k++)
730 for (
unsigned int i = 0; i < nbExpr; i++)
732 m_VParser[k][i]->SetExpr(m_Expression[i]);
738 template <
typename TImage>
743 for (
unsigned int i = 0; i < m_StatsVarDetected.size(); i++)
748 filter->SetInput(this->GetNthInput(m_StatsVarDetected[i]));
754 for (
unsigned int t = 0; t < m_AImage.size(); t++)
755 for (
unsigned int v = 0; v < m_AImage[t].size(); v++)
756 if ((m_AImage[t][v].type == 8) && (m_AImage[t][v].info[0] == m_StatsVarDetected[i]))
758 switch (m_AImage[t][v].info[2])
762 pix = filter->GetMinimum();
763 for (
int b = 0; b < (int)pix.GetSize(); b++)
764 if (m_AImage[t][v].info[1] == b)
765 m_AImage[t][v].value = pix[b];
771 pix = filter->GetMaximum();
772 for (
int b = 0; b < (int)pix.GetSize(); b++)
773 if (m_AImage[t][v].info[1] == b)
774 m_AImage[t][v].value = pix[b];
780 pix = filter->GetMean();
781 for (
int b = 0; b < (int)pix.GetSize(); b++)
782 if (m_AImage[t][v].info[1] == b)
783 m_AImage[t][v].value = pix[b];
791 pix = filter->GetSum();
792 for (
int b = 0; b < (int)pix.GetSize(); b++)
793 if (m_AImage[t][v].info[1] == b)
794 m_AImage[t][v].value = pix[b];
800 mat = filter->GetCovariance();
801 for (
int b = 0; b < (int)mat.Cols(); b++)
802 if (m_AImage[t][v].info[1] == b)
803 m_AImage[t][v].value = mat(b, b);
811 template <
typename TImage>
815 this->SetNumberOfRequiredOutputs((
int)m_Expression.size());
817 m_outputsDimensions.clear();
819 for (
int i = 0; i < (int)m_Expression.size(); ++i)
821 ValueType value = m_VParser[0][i]->EvalRef();
823 switch (value.GetType())
826 itkExceptionMacro(<<
"Booleans not supported." << std::endl);
830 m_outputsDimensions.push_back(1);
834 m_outputsDimensions.push_back(1);
838 itkExceptionMacro(<<
"Complex numbers not supported." << std::endl);
843 const mup::matrix_type& vect = value.GetArray();
844 if (vect.GetRows() == 1)
845 m_outputsDimensions.push_back(vect.GetCols());
847 itkExceptionMacro(<<
"Result of the evaluation can't be a matrix." << std::endl);
852 itkExceptionMacro(<<
"Unknown output type : " << value.GetType() << std::endl);
860 template <
typename TImage>
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);
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)))
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) <<
"]");
881 template <
typename TImage>
884 Superclass::GenerateOutputInformation();
887 CheckImageDimensions();
889 if (GlobalStatsDetected())
890 PrepareParsersGlobStats();
894 typedef itk::ImageBase<TImage::ImageDimension>
ImageBaseType;
895 typename ImageBaseType::Pointer outputPtr;
898 for (itk::OutputDataObjectIterator it(
this); !it.IsAtEnd(); i++, it++)
908 outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]);
913 template <
typename TImage>
917 Superclass::GenerateInputRequestedRegion();
919 for (
unsigned int i = 0; i < m_NeighDetected.size(); i++)
920 if (m_NeighDetected[i] < this->GetNumberOfInputs())
924 typename Superclass::InputImagePointer inputPtr =
const_cast<TImage*
>(this->GetNthInput(m_NeighDetected[i]));
928 inputRequestedRegion = inputPtr->GetRequestedRegion();
931 inputRequestedRegion.PadByRadius(m_NeighExtremaSizes[m_NeighDetected[i]]);
934 if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
936 inputPtr->SetRequestedRegion(inputRequestedRegion);
945 inputPtr->SetRequestedRegion(inputRequestedRegion);
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);
959 itkExceptionMacro(<<
"Requested input #" << m_NeighDetected[i] <<
", but only " << this->GetNumberOfInputs() <<
" inputs are available." << std::endl);
963 template <
typename TImage>
967 unsigned int nbThreads = this->GetNumberOfThreads();
969 m_ThreadUnderflow.SetSize(nbThreads);
970 m_ThreadUnderflow.Fill(0);
971 m_ThreadOverflow.SetSize(nbThreads);
972 m_ThreadOverflow.Fill(0);
976 template <
typename TImage>
979 unsigned int nbThreads = this->GetNumberOfThreads();
982 m_UnderflowCount = 0;
986 for (i = 0; i < nbThreads; ++i)
988 m_UnderflowCount += m_ThreadUnderflow[i];
989 m_OverflowCount += m_ThreadOverflow[i];
992 if ((m_UnderflowCount != 0) || (m_OverflowCount != 0))
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 !";
1007 template <
typename TImage>
1012 unsigned int nbInputImages = this->GetNumberOfInputs();
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);
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);
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)
1038 radius[0] = (int)((m_VVarName[j].info[2] - 1) / 2);
1039 radius[1] = (int)((m_VVarName[j].info[3] - 1) / 2);
1041 itk::ConstNeighborhoodIterator<TImage>(radius, this->GetNthInput(m_VVarName[j].info[0]), outputRegionForThread));
1042 VNit.back().NeedToUseBoundaryConditionOn();
1046 IndexIteratorType indexIterator(this->GetNthInput(0), outputRegionForThread);
1049 itk::ProgressReporter progress(
this, threadId, outputRegionForThread.GetNumberOfPixels());
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;
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]);
1064 for (
unsigned int j = 0; j < nbInputImages; ++j)
1068 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1070 VoutIt[j].GoToBegin();
1072 for (
unsigned int j = 0; j < VNit.size(); ++j)
1074 VNit[j].GoToBegin();
1076 indexIterator.GoToBegin();
1078 while (!Vit[0].IsAtEnd())
1080 while (!Vit[0].IsAtEndOfLine())
1082 int ngbhNameIndex = 0;
1085 iterVar = iterVarStart;
1086 while (iterVar != iterVarEnd)
1088 switch (iterVar->type)
1091 iterVar->value =
static_cast<double>(indexIterator.GetIndex()[0]);
1095 iterVar->value =
static_cast<double>(indexIterator.GetIndex()[1]);
1108 for (
int p = 0; p < iterVar->value.GetCols(); ++p)
1109 iterVar->value.At(0, p) = Vit[iterVar->info[0]].Get()[p];
1115 iterVar->value = Vit[iterVar->info[0]].Get()[iterVar->info[1]];
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")
1125 for (
int rows = 0; rows < iterVar->info[3]; ++rows)
1126 for (
int cols = 0; cols < iterVar->info[2]; ++cols)
1128 iterVar->value.At(rows, cols) = VNit[ngbhNameIndex].GetPixel(index)[iterVar->info[1]];
1144 itkExceptionMacro(<<
"Type of the variable is unknown");
1154 for (
unsigned int IDExpression = 0; IDExpression < m_Expression.size(); ++IDExpression)
1156 value = m_VParser[threadId][IDExpression]->EvalRef();
1158 switch (value.GetType())
1161 tmpOutputs[IDExpression][0] = value.GetInteger();
1165 tmpOutputs[IDExpression][0] = value.GetFloat();
1169 itkExceptionMacro(<<
"Complex numbers are not supported." << std::endl);
1174 const mup::matrix_type& vect = value.GetArray();
1176 if (vect.GetRows() == 1)
1177 for (
int p = 0; p < vect.GetCols(); ++p)
1178 tmpOutputs[IDExpression][p] = vect.At(0, p).GetFloat();
1180 itkExceptionMacro(<<
"Result of the evaluation can't be a matrix." << std::endl);
1186 for (
unsigned int p = 0; p < m_outputsDimensions[IDExpression]; ++p)
1190 if (tmpOutputs[IDExpression][p] <
double(itk::NumericTraits<PixelValueType>::NonpositiveMin()))
1192 tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::NonpositiveMin();
1193 m_ThreadUnderflow[threadId]++;
1197 else if (tmpOutputs[IDExpression][p] >
double(itk::NumericTraits<PixelValueType>::max()))
1199 tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::max();
1200 m_ThreadOverflow[threadId]++;
1203 VoutIt[IDExpression].Set(tmpOutputs[IDExpression]);
1206 for (
unsigned int j = 0; j < nbInputImages; ++j)
1210 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1214 for (
unsigned int j = 0; j < VNit.size(); ++j)
1220 progress.CompletedPixel();
1222 for (
unsigned int j = 0; j < nbInputImages; ++j)
1226 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1228 VoutIt[j].NextLine();