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);
50 this->DynamicMultiThreadingOff();
53 m_ThreadUnderflow.SetSize(1);
54 m_ThreadOverflow.SetSize(1);
62 m_VAllowedVarNameAuto.push_back(ahcX);
67 m_VAllowedVarNameAuto.push_back(ahcY);
69 m_SizeNeighbourhood = 10;
71 m_ManyExpressions =
true;
75 template <
class TImage>
82 for (
unsigned int i = 0; i < m_AImage.size(); ++i)
87 m_VAllowedVarNameAuto.clear();
88 m_VAllowedVarNameAddedByUser.clear();
89 m_VFinalAllowedVarName.clear();
90 m_VNotAllowedVarName.clear();
91 m_outputsDimensions.clear();
95 template <
class TImage>
98 Superclass::PrintSelf(os, indent);
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;
110 template <
class TImage>
114 std::stringstream sstm;
115 sstm <<
"im" << (idx + 1);
116 this->SetNthInput(idx, image, sstm.str());
119 template <
class TImage>
124 this->SetInput(idx, imagebis);
127 std::stringstream sstmPhyX;
129 sstmPhyX << varName <<
"PhyX";
130 ahcPhyX.
name = sstmPhyX.str();
132 ahcPhyX.
info[0] = idx;
133 m_VAllowedVarNameAuto.push_back(ahcPhyX);
135 std::stringstream sstmPhyY;
137 sstmPhyY << varName <<
"PhyY";
138 ahcPhyY.
name = sstmPhyY.str();
140 ahcPhyY.
info[0] = idx;
141 m_VAllowedVarNameAuto.push_back(ahcPhyY);
144 std::stringstream sstm_glob;
146 sstm_glob << varName;
147 ahc_glob.
name = sstm_glob.str();
149 ahc_glob.
info[0] = idx;
150 m_VAllowedVarNameAuto.push_back(ahc_glob);
155 imagebis->UpdateOutputInformation();
158 for (
unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
160 std::stringstream sstm;
162 sstm << varName <<
"b" << (j + 1);
163 ahc.
name = sstm.str();
167 m_VAllowedVarNameAuto.push_back(ahc);
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++)
175 std::stringstream sstm;
177 sstm << varName <<
"b" << (j + 1) <<
"N" << 2 * x + 1 <<
"x" << 2 * y + 1;
178 ahc.
name = sstm.str();
182 ahc.
info[2] = 2 * x + 1;
183 ahc.
info[3] = 2 * y + 1;
184 m_VAllowedVarNameAuto.push_back(ahc);
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");
195 for (
unsigned int j = 0; j < imagebis->GetNumberOfComponentsPerPixel(); j++)
196 for (
unsigned int t = 0; t < statsNames.size(); t++)
198 std::stringstream sstm;
200 sstm << varName <<
"b" << (j + 1) << statsNames[t];
201 ahc.
name = sstm.str();
206 m_VAllowedVarNameAuto.push_back(ahc);
210 template <
typename TImage>
213 return const_cast<TImage*
>(this->GetInput(idx));
216 template <
typename TImage>
219 m_ManyExpressions = flag;
222 template <
typename TImage>
225 std::string expressionToBePushed = expression;
227 if (expression.find(
";") != std::string::npos)
229 std::ostringstream oss;
231 for (
unsigned int i = 0; i < expression.size(); ++i)
232 if (expression[i] ==
';')
235 oss << expression[i];
238 expressionToBePushed = oss.str();
241 if (m_ManyExpressions)
242 m_Expression.push_back(expressionToBePushed);
243 else if (m_Expression.size() == 0)
244 m_Expression.push_back(expressionToBePushed);
246 if (m_Expression.size() > 1)
252 template <
typename TImage>
255 m_Expression.clear();
258 template <
typename TImage>
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);
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);
272 for (
unsigned int i = 1; i < definition.size() - 1; ++i)
273 def.push_back(definition[i]);
276 std::vector<std::vector<double>> mat;
277 std::istringstream iss(def);
279 while (std::getline(iss, rows,
';'))
281 mat.push_back(std::vector<double>(0));
282 std::istringstream iss2(rows);
284 while (std::getline(iss2, elmt,
','))
286 std::istringstream iss3(elmt);
289 mat.back().push_back(val);
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);
304 ahc.
info[0] = mat[0].size();
305 ahc.
info[1] = mat.size();
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];
312 m_VAllowedVarNameAddedByUser.push_back(ahc);
316 template <
typename TImage>
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);
328 m_VAllowedVarNameAddedByUser.push_back(ahc);
332 template <
typename TImage>
336 std::vector<std::string> vectI, vectF, vectM, vectFinal;
338 for (
unsigned int i = 0; i < m_VAllowedVarNameAddedByUser.size(); i++)
340 std::ostringstream iss;
343 switch (m_VAllowedVarNameAddedByUser[i].value.GetType())
346 iss <<
"#I " << m_VAllowedVarNameAddedByUser[i].name <<
" " << m_VAllowedVarNameAddedByUser[i].value.GetInteger();
348 vectI.push_back(str);
351 iss <<
"#F " << m_VAllowedVarNameAddedByUser[i].name <<
" " << m_VAllowedVarNameAddedByUser[i].value.GetFloat();
353 vectF.push_back(str);
356 itkExceptionMacro(<<
"Complex numbers not supported." << std::endl);
359 iss <<
"#M " << m_VAllowedVarNameAddedByUser[i].name <<
" "
361 for (
int k = 0; k < m_VAllowedVarNameAddedByUser[i].value.GetRows(); k++)
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);
369 str.erase(str.size() - 1);
371 vectM.push_back(str);
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)
385 std::ostringstream iss;
386 iss <<
"#E " << m_Expression[i] << std::endl;
387 std::string str = iss.str();
388 vectFinal.push_back(str);
391 std::ofstream exportFile(filename, std::ios::out | std::ios::trunc);
394 for (
unsigned int i = 0; i < vectFinal.size(); ++i)
395 exportFile << vectFinal[i] << std::endl;
400 itkExceptionMacro(<<
"Could not open " << filename <<
"." << std::endl);
403 template <
typename TImage>
406 std::ifstream importFile(filename, std::ios::in);
408 std::string wholeline, line, name, matrixdef;
409 int pos, pos2, lineID = 0, nbSuccesses = 0;
415 while (std::getline(importFile, wholeline))
419 pos = wholeline.find_first_not_of(
' ');
421 if (pos != (
int)std::string::npos)
423 line = wholeline.substr(pos);
425 if ((line[0] ==
'#') && ((line[1] ==
'I') || (line[1] ==
'i') || (line[1] ==
'F') || (line[1] ==
'f')))
428 pos = line.find_first_not_of(
' ', 2);
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);
433 std::string sub = line.substr(pos);
435 pos = sub.find_first_of(
' ');
436 name = sub.substr(0, pos);
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);
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)
445 std::istringstream iss(sub.substr(pos));
448 SetConstant(name, value);
451 else if ((line[0] ==
'#') && ((line[1] ==
'M') || (line[1] ==
'm')))
454 pos = line.find_first_not_of(
' ', 2);
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."
460 std::string sub = line.substr(pos);
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);
468 itkExceptionMacro(<<
"In file '" << filename <<
"', line " << lineID <<
" : symbol #M found, but couldn't not find vector/matrix definition."
471 SetMatrix(name, matrixdef);
474 else if ((line[0] ==
'#') && ((line[1] ==
'E') || (line[1] ==
'e')))
476 pos = line.find_first_not_of(
' ', 2);
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);
481 std::string sub = line.substr(pos);
492 if (nbSuccesses == 0)
493 itkExceptionMacro(<<
"No constant or expression could be set; please, ensure that the file '" << filename <<
"' is correct." << std::endl);
496 itkExceptionMacro(<<
"Could not open " << filename <<
"." << std::endl);
500 template <
typename TImage>
503 if (IDExpression < m_Expression.size())
504 return m_Expression[IDExpression];
509 template <
typename TImage>
512 std::vector<std::string> res;
513 for (
int y = 0; y < m_VVarName.size(); y++)
514 res.push_back(m_VVarName[y].name);
520 template <
typename TImage>
524 for (
unsigned int i = 0; i < m_VVarName.size(); ++i)
525 if (m_VVarName[i].name == ahc.
name)
529 m_VVarName.push_back(ahc);
533 template <
typename TImage>
537 if (m_Expression.size() == 0)
538 itkExceptionMacro(<<
"No expression set; please set at least one expression." << std::endl);
543 m_VNotAllowedVarName.clear();
544 m_VFinalAllowedVarName.clear();
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]);
553 unsigned int nbExpr = m_Expression.size();
554 for (
unsigned int IDExpression = 0; IDExpression < nbExpr; ++IDExpression)
557 dummyParser->SetExpr(this->GetExpression(IDExpression));
559 mup::var_maptype vmap = dummyParser->GetExprVar();
560 for (mup::var_maptype::iterator item = vmap.begin(); item != vmap.end(); ++item)
564 while ((!OK) && (i < (
int)m_VFinalAllowedVarName.size()))
566 if (item->first == m_VFinalAllowedVarName[i].name)
574 AddVariable(m_VFinalAllowedVarName[i]);
579 ahc.
name = item->first;
580 m_VNotAllowedVarName.push_back(ahc);
588 if (m_VNotAllowedVarName.size() > 0)
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 <<
" ";
595 itkExceptionMacro(<< sstm.str());
601 unsigned int nbThreads = this->GetNumberOfWorkUnits();
602 for (
unsigned int k = 0; k < nbThreads; k++)
604 std::vector<ParserType::Pointer> parserList;
605 for (
unsigned int i = 0; i < nbExpr; i++)
607 parserList.push_back(ParserType::New());
609 m_VParser.push_back(parserList);
614 int nbVar = m_VVarName.size();
616 m_StatsVarDetected.clear();
617 m_NeighDetected.clear();
618 m_NeighExtremaSizes.clear();
619 unsigned int nbInputImages = this->GetNumberOfInputs();
623 m_NeighExtremaSizes.resize(nbInputImages, dummyRadius);
626 for (
unsigned int i = 0; i < m_AImage.size(); ++i)
630 m_AImage.resize(nbThreads);
632 double initValue = 0.1;
633 for (
unsigned int i = 0; i < nbThreads; ++i)
636 m_AImage[i].resize(nbVar);
638 for (
int j = 0; j < nbVar; ++j)
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];
646 if ((m_AImage[i][j].type == 0) || (m_AImage[i][j].type == 1))
648 m_AImage[i][j].value =
ValueType(initValue);
651 if (m_AImage[i][j].type == 2)
653 SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
654 m_AImage[i][j].value =
ValueType(
static_cast<double>(spacing[0]));
657 if (m_AImage[i][j].type == 3)
659 SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSignedSpacing();
660 m_AImage[i][j].value =
ValueType(
static_cast<double>(spacing[1]));
663 if (m_AImage[i][j].type == 4)
665 unsigned int nbBands = this->GetNthInput(m_AImage[i][j].info[0])->GetNumberOfComponentsPerPixel();
666 m_AImage[i][j].value =
ValueType(1, nbBands, initValue);
669 if (m_AImage[i][j].type == 5)
671 m_AImage[i][j].value =
ValueType(initValue);
674 if (m_AImage[i][j].type == 6)
676 m_AImage[i][j].value =
ValueType(m_AImage[i][j].info[3], m_AImage[i][j].info[2], initValue);
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])
684 m_NeighDetected.push_back(m_AImage[i][j].info[0]);
687 if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] < (
unsigned int)((m_VVarName[j].info[2] - 1) / 2))
688 m_NeighExtremaSizes[m_AImage[i][j].info[0]][0] = (
unsigned int)((m_VVarName[j].info[2] - 1) / 2);
690 if (m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] < (
unsigned int)((m_VVarName[j].info[3] - 1) / 2))
691 m_NeighExtremaSizes[m_AImage[i][j].info[0]][1] = (
unsigned int)((m_VVarName[j].info[3] - 1) / 2);
694 if (m_AImage[i][j].type == 7)
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;
702 if (m_AImage[i][j].type == 8)
704 m_AImage[i][j].value =
ValueType(initValue);
707 for (
unsigned int r = 0; r < m_StatsVarDetected.size() && !found; r++)
708 if (m_StatsVarDetected[r] == m_AImage[i][j].info[0])
711 m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
716 for (
unsigned int k = 0; k < nbExpr; k++)
718 m_VParser[i][k]->DefineVar(m_AImage[i][j].name, &(m_AImage[i][j].value));
729 for (
unsigned int k = 0; k < nbThreads; k++)
731 for (
unsigned int i = 0; i < nbExpr; i++)
733 m_VParser[k][i]->SetExpr(m_Expression[i]);
739 template <
typename TImage>
744 for (
unsigned int i = 0; i < m_StatsVarDetected.size(); i++)
749 filter->SetInput(this->GetNthInput(m_StatsVarDetected[i]));
755 for (
unsigned int t = 0; t < m_AImage.size(); t++)
756 for (
unsigned int v = 0; v < m_AImage[t].size(); v++)
757 if ((m_AImage[t][v].type == 8) && (m_AImage[t][v].info[0] == m_StatsVarDetected[i]))
759 switch (m_AImage[t][v].info[2])
763 pix = filter->GetMinimum();
764 for (
int b = 0; b < (int)pix.GetSize(); b++)
765 if (m_AImage[t][v].info[1] == b)
766 m_AImage[t][v].value = pix[b];
772 pix = filter->GetMaximum();
773 for (
int b = 0; b < (int)pix.GetSize(); b++)
774 if (m_AImage[t][v].info[1] == b)
775 m_AImage[t][v].value = pix[b];
781 pix = filter->GetMean();
782 for (
int b = 0; b < (int)pix.GetSize(); b++)
783 if (m_AImage[t][v].info[1] == b)
784 m_AImage[t][v].value = pix[b];
792 pix = filter->GetSum();
793 for (
int b = 0; b < (int)pix.GetSize(); b++)
794 if (m_AImage[t][v].info[1] == b)
795 m_AImage[t][v].value = pix[b];
801 mat = filter->GetCovariance();
802 for (
int b = 0; b < (int)mat.Cols(); b++)
803 if (m_AImage[t][v].info[1] == b)
804 m_AImage[t][v].value = mat(b, b);
812 template <
typename TImage>
816 this->SetNumberOfRequiredOutputs((
int)m_Expression.size());
818 m_outputsDimensions.clear();
820 for (
int i = 0; i < (int)m_Expression.size(); ++i)
822 ValueType value = m_VParser[0][i]->EvalRef();
824 switch (value.GetType())
827 itkExceptionMacro(<<
"Booleans not supported." << std::endl);
831 m_outputsDimensions.push_back(1);
835 m_outputsDimensions.push_back(1);
839 itkExceptionMacro(<<
"Complex numbers not supported." << std::endl);
844 const mup::matrix_type& vect = value.GetArray();
845 if (vect.GetRows() == 1)
846 m_outputsDimensions.push_back(vect.GetCols());
848 itkExceptionMacro(<<
"Result of the evaluation can't be a matrix." << std::endl);
853 itkExceptionMacro(<<
"Unknown output type : " << value.GetType() << std::endl);
861 template <
typename TImage>
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);
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)))
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) <<
"]");
882 template <
typename TImage>
885 Superclass::GenerateOutputInformation();
888 CheckImageDimensions();
890 if (GlobalStatsDetected())
891 PrepareParsersGlobStats();
895 typedef itk::ImageBase<TImage::ImageDimension>
ImageBaseType;
896 typename ImageBaseType::Pointer outputPtr;
899 for (itk::OutputDataObjectIterator it(
this); !it.IsAtEnd(); i++, it++)
909 outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]);
914 template <
typename TImage>
918 Superclass::GenerateInputRequestedRegion();
920 for (
unsigned int i = 0; i < m_NeighDetected.size(); i++)
921 if (m_NeighDetected[i] < this->GetNumberOfInputs())
925 typename Superclass::InputImagePointer inputPtr =
const_cast<TImage*
>(this->GetNthInput(m_NeighDetected[i]));
929 inputRequestedRegion = inputPtr->GetRequestedRegion();
932 inputRequestedRegion.PadByRadius(m_NeighExtremaSizes[m_NeighDetected[i]]);
935 if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
937 inputPtr->SetRequestedRegion(inputRequestedRegion);
946 inputPtr->SetRequestedRegion(inputRequestedRegion);
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);
960 itkExceptionMacro(<<
"Requested input #" << m_NeighDetected[i] <<
", but only " << this->GetNumberOfInputs() <<
" inputs are available." << std::endl);
964 template <
typename TImage>
968 unsigned int nbThreads = this->GetNumberOfWorkUnits();
970 m_ThreadUnderflow.SetSize(nbThreads);
971 m_ThreadUnderflow.Fill(0);
972 m_ThreadOverflow.SetSize(nbThreads);
973 m_ThreadOverflow.Fill(0);
977 template <
typename TImage>
980 unsigned int nbThreads = this->GetNumberOfWorkUnits();
983 m_UnderflowCount = 0;
987 for (i = 0; i < nbThreads; ++i)
989 m_UnderflowCount += m_ThreadUnderflow[i];
990 m_OverflowCount += m_ThreadOverflow[i];
993 if ((m_UnderflowCount != 0) || (m_OverflowCount != 0))
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 !";
1008 template <
typename TImage>
1013 unsigned int nbInputImages = this->GetNumberOfInputs();
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);
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);
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)
1039 radius[0] = (int)((m_VVarName[j].info[2] - 1) / 2);
1040 radius[1] = (int)((m_VVarName[j].info[3] - 1) / 2);
1042 itk::ConstNeighborhoodIterator<TImage>(radius, this->GetNthInput(m_VVarName[j].info[0]), outputRegionForThread));
1043 VNit.back().NeedToUseBoundaryConditionOn();
1047 IndexIteratorType indexIterator(this->GetNthInput(0), outputRegionForThread);
1050 itk::ProgressReporter progress(
this, threadId, outputRegionForThread.GetNumberOfPixels());
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;
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]);
1065 for (
unsigned int j = 0; j < nbInputImages; ++j)
1069 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1071 VoutIt[j].GoToBegin();
1073 for (
unsigned int j = 0; j < VNit.size(); ++j)
1075 VNit[j].GoToBegin();
1077 indexIterator.GoToBegin();
1079 while (!Vit[0].IsAtEnd())
1081 while (!Vit[0].IsAtEndOfLine())
1083 int ngbhNameIndex = 0;
1086 iterVar = iterVarStart;
1087 while (iterVar != iterVarEnd)
1089 switch (iterVar->type)
1092 iterVar->value =
static_cast<double>(indexIterator.GetIndex()[0]);
1096 iterVar->value =
static_cast<double>(indexIterator.GetIndex()[1]);
1109 for (
int p = 0; p < iterVar->value.GetCols(); ++p)
1110 iterVar->value.At(0, p) = Vit[iterVar->info[0]].Get()[p];
1116 iterVar->value = Vit[iterVar->info[0]].Get()[iterVar->info[1]];
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")
1126 for (
int rows = 0; rows < iterVar->info[3]; ++rows)
1127 for (
int cols = 0; cols < iterVar->info[2]; ++cols)
1129 iterVar->value.At(rows, cols) = VNit[ngbhNameIndex].GetPixel(index)[iterVar->info[1]];
1145 itkExceptionMacro(<<
"Type of the variable is unknown");
1155 for (
unsigned int IDExpression = 0; IDExpression < m_Expression.size(); ++IDExpression)
1157 value = m_VParser[threadId][IDExpression]->EvalRef();
1159 switch (value.GetType())
1162 tmpOutputs[IDExpression][0] = value.GetInteger();
1166 tmpOutputs[IDExpression][0] = value.GetFloat();
1170 itkExceptionMacro(<<
"Complex numbers are not supported." << std::endl);
1175 const mup::matrix_type& vect = value.GetArray();
1177 if (vect.GetRows() == 1)
1178 for (
int p = 0; p < vect.GetCols(); ++p)
1179 tmpOutputs[IDExpression][p] = vect.At(0, p).GetFloat();
1181 itkExceptionMacro(<<
"Result of the evaluation can't be a matrix." << std::endl);
1187 for (
unsigned int p = 0; p < m_outputsDimensions[IDExpression]; ++p)
1191 if (tmpOutputs[IDExpression][p] <
double(itk::NumericTraits<PixelValueType>::NonpositiveMin()))
1193 tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::NonpositiveMin();
1194 m_ThreadUnderflow[threadId]++;
1198 else if (tmpOutputs[IDExpression][p] >
double(itk::NumericTraits<PixelValueType>::max()))
1200 tmpOutputs[IDExpression][p] = itk::NumericTraits<PixelValueType>::max();
1201 m_ThreadOverflow[threadId]++;
1204 VoutIt[IDExpression].Set(tmpOutputs[IDExpression]);
1207 for (
unsigned int j = 0; j < nbInputImages; ++j)
1211 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1215 for (
unsigned int j = 0; j < VNit.size(); ++j)
1221 progress.CompletedPixel();
1223 for (
unsigned int j = 0; j < nbInputImages; ++j)
1227 for (
unsigned int j = 0; j < m_Expression.size(); ++j)
1229 VoutIt[j].NextLine();
ImageType::SpacingType SpacingType
void AfterThreadedGenerateData() override
void SetExpression(const std::string &expression)
void ThreadedGenerateData(const ImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
void GenerateInputRequestedRegion() override
ImageType::RegionType ImageRegionType
StreamingStatisticsVectorImageFilterType::Pointer StreamingStatisticsVectorImageFilterPointerType
ImageType * GetNthInput(DataObjectPointerArraySizeType idx)
void ImportContext(const std::string &filename)
itk::ConstNeighborhoodIterator< TImage >::RadiusType RadiusType
~BandMathXImageFilter() override
void AddVariable(adhocStruct &)
void CheckImageDimensions()
void SetManyExpressions(bool flag)
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 GenerateOutputInformation() override
void PrepareParsersGlobStats()
void BeforeThreadedGenerateData() override
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
ImageType::PixelType PixelType
std::string GetExpression(unsigned int IDExpression) const
ParserType::ValueType ValueType
itk::SmartPointer< Self > Pointer
itk::ImageBase< 2 > ImageBaseType
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
#define otbWarningMacro(x)