otbBandMathXImageFilter.txx 34.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.

  Some parts of this code are derived from ITK. See ITKCopyright.txt
  for details.


     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
21 22 23
#ifndef __otbBandMathXImageFilter_txx
#define __otbBandMathXImageFilter_txx
#include "otbBandMathXImageFilter.h"
24 25 26

#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
27
#include "itkConstNeighborhoodIterator.h"
28 29 30 31 32
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "otbMacro.h"

#include <iostream>
33
#include <fstream>
34 35 36 37 38 39 40
#include <string>

namespace otb
{

/** Constructor */
template <class TImage>
41 42
BandMathXImageFilter<TImage>
::BandMathXImageFilter()
43 44 45 46 47 48 49 50 51
{
  //This number will be incremented each time an image
  //is added over the one minimumrequired
  this->SetNumberOfRequiredInputs( 1 );

  m_UnderflowCount = 0;
  m_OverflowCount = 0;
  m_ThreadUnderflow.SetSize(1);
  m_ThreadOverflow.SetSize(1);
52 53 54 55 56 57


  //idxX and idxY
  adhocStruct ahcX;
  ahcX.name = "idxX";
  ahcX.type = 0;
58
  m_VAllowedVarNameAuto.push_back(ahcX);
59 60 61 62

  adhocStruct ahcY;
  ahcY.name = "idxY";
  ahcY.type = 1;
63
  m_VAllowedVarNameAuto.push_back(ahcY);
64

65
  m_SizeNeighbourhood=10;
66

67 68 69 70
}

/** Destructor */
template <class TImage>
71 72
BandMathXImageFilter<TImage>
::~BandMathXImageFilter()
73
{
74 75 76 77 78 79 80 81 82 83
  m_Expression.clear();
  m_VParser.clear();

  for(int i=0; i<m_AImage.size(); ++i)
    m_AImage[i].clear();
  m_AImage.clear();

  m_VVarName.clear();
  m_VAllowedVarNameAuto.clear();
  m_VAllowedVarNameAddedByUser.clear();
OTB Bot's avatar
STYLE  
OTB Bot committed
84
  m_VFinalAllowedVarName.clear();
85 86 87 88
  m_VNotAllowedVarName.clear();
  m_outputsDimensions.clear();


89 90
}

91

92
template <class TImage>
93
void BandMathXImageFilter<TImage>
94 95 96 97
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);

98 99
  os << indent << "Expressions: " << std::endl;
  for (unsigned int i=0; i<m_Expression.size(); i++)
OTB Bot's avatar
STYLE  
OTB Bot committed
100
    os << indent << m_Expression[i] << std::endl;
101 102 103
  os << indent << "Computed values follow:"                            << std::endl;
  os << indent << "UnderflowCount: "  << m_UnderflowCount              << std::endl;
  os << indent << "OverflowCount: "   << m_OverflowCount               << std::endl;
104 105 106 107
  os << indent << "itk::NumericTraits<typename PixelValueType>::NonpositiveMin()  :  "
               << itk::NumericTraits<PixelValueType>::NonpositiveMin()      << std::endl;
  os << indent << "itk::NumericTraits<typename PixelValueType>::max()  :             "
         << itk::NumericTraits<PixelValueType>::max()                 << std::endl;
108 109 110
}

template <class TImage>
111
void BandMathXImageFilter<TImage>
112 113
::SetNthInput(unsigned int idx, const ImageType * image)
{
114 115 116 117

  std::stringstream sstm;
  sstm << "im" << (idx+1);
  this->SetNthInput(idx, image, sstm.str());
118 119 120
}

template <class TImage>
121
void BandMathXImageFilter<TImage>
122 123
::SetNthInput(unsigned int idx, const ImageType * image, const std::string& varName)
{
124 125 126

  ImageType * imagebis = const_cast<ImageType *>( image ); // Useful for call of UpdateOutputInformation() (see below)
  this->SetInput(idx, imagebis );
127

128 129 130 131 132 133 134
  //imiPhyX and imiPhyY
  std::stringstream sstmPhyX;
  adhocStruct ahcPhyX;
  sstmPhyX << varName << "PhyX";
  ahcPhyX.name = sstmPhyX.str();
  ahcPhyX.type = 2;
  ahcPhyX.info[0] = idx; // Input image #ID
135
  m_VAllowedVarNameAuto.push_back(ahcPhyX);
136 137 138 139 140 141 142

  std::stringstream sstmPhyY;
  adhocStruct ahcPhyY;
  sstmPhyY << varName << "PhyY";
  ahcPhyY.name = sstmPhyY.str();
  ahcPhyY.type = 3;
  ahcPhyY.info[0] = idx; // Input image #ID
143
  m_VAllowedVarNameAuto.push_back(ahcPhyY);
144 145 146 147 148 149 150 151

  //imi
  std::stringstream sstm;
  adhocStruct ahc;
  sstm << varName;
  ahc.name = sstm.str();
  ahc.type = 4;
  ahc.info[0] = idx; // Input image #ID
152
  m_VAllowedVarNameAuto.push_back(ahc);
153

154
  //Mandatory before call of GetNumberOfComponentsPerPixel
155
  //Really important not to call the filter's UpdateOutputInformation method here:
156
  //this method is not ready until all inputs, variables and expressions are set.
OTB Bot's avatar
STYLE  
OTB Bot committed
157
  imagebis->UpdateOutputInformation();
158

159
  //imibj
160
  for (int j=0; j<imagebis->GetNumberOfComponentsPerPixel(); j++)
161 162 163 164 165 166 167 168
  {
    std::stringstream sstm;
    adhocStruct ahc;
    sstm << varName << "b" << (j+1);
    ahc.name = sstm.str();
    ahc.type = 5;
    ahc.info[0] = idx; // Input image #ID
    ahc.info[1] = j; // Band #ID
169
    m_VAllowedVarNameAuto.push_back(ahc);
170 171 172
  }

  //imibjNkxp
173
  for (int j=0; j<imagebis->GetNumberOfComponentsPerPixel(); j++)
174 175
    for(int x=0; x<=m_SizeNeighbourhood; x++)
      for(int y=0; y<=m_SizeNeighbourhood; y++)
176 177 178
      {
        std::stringstream sstm;
        adhocStruct ahc;
179
        sstm << varName << "b" << (j+1) << "N" << 2*x+1 << "x" << 2*y+1;
180 181 182
        ahc.name = sstm.str();
        ahc.type = 6;
        ahc.info[0] = idx; // Input image #ID
183 184 185 186
        ahc.info[1] = j;   // Band #ID
        ahc.info[2] = 2*x+1;  // Size x direction (matrix convention = cols)
        ahc.info[3] = 2*y+1;  // Size y direction (matrix convention = rows)
        m_VAllowedVarNameAuto.push_back(ahc);
187 188
      }

189 190 191 192 193
  //imibjSTATS
  std::vector< std::string > statsNames;
  statsNames.push_back("Mini");
  statsNames.push_back("Maxi");
  statsNames.push_back("Mean");
194 195
  statsNames.push_back("Sum");
  statsNames.push_back("Var");
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210

  for (int j=0; j<imagebis->GetNumberOfComponentsPerPixel(); j++)
    for (int t=0; t<statsNames.size(); t++)
    {
      std::stringstream sstm;
      adhocStruct ahc;
      sstm << varName << "b" << (j+1) << statsNames[t];
      ahc.name = sstm.str();
      ahc.type = 8;
      ahc.info[0] = idx; // Input image #ID
      ahc.info[1] = j;   // Band #ID
      ahc.info[2] = t;   // Sub-type : 0 Mini, 1 Maxi, 2 Mean ...
      m_VAllowedVarNameAuto.push_back(ahc);
    }

211 212 213
}

template <typename TImage>
214
TImage * BandMathXImageFilter<TImage>
215 216 217 218 219
::GetNthInput(unsigned int idx)
{
  return const_cast<TImage *>(this->GetInput(idx));
}

220

221
template< typename TImage >
222
void BandMathXImageFilter<TImage>
OTB Bot's avatar
STYLE  
OTB Bot committed
223
::SetExpression(const std::string& expression)
224
{
225

226
  if (expression.find(";") != std::string::npos)
227 228 229 230
  {
    std::ostringstream oss;
    oss << "cat(";
    for(int i=0; i < expression.size(); ++i)
231
      if (expression[i] == ';')
232
        oss << ",";
OTB Bot's avatar
STYLE  
OTB Bot committed
233
      else
234 235 236 237 238 239 240 241
        oss << expression[i];

    oss << ")";
    m_Expression.push_back(oss.str());

  }
  else
    m_Expression.push_back(expression);
242 243 244 245

  if (m_Expression.size()>1)
    this->SetNthOutput( (int) (m_Expression.size()) -1, ( TImage::New() ).GetPointer() );

246 247 248
  this->Modified();
}

249

250
template< typename TImage >
251
void BandMathXImageFilter<TImage>
252 253 254
::SetMatrix(const std::string& name, const std::string& definition)
{

255 256 257 258
  for(int i=0; i<m_VAllowedVarNameAddedByUser.size(); i++)
    if (name.compare(m_VAllowedVarNameAddedByUser[i].name) == 0)
      itkExceptionMacro(<< "Variable name '"<< name << "' already used." << std::endl);

259 260 261 262 263 264 265 266 267 268 269

  if ( (definition.find("{") != 0) || (definition.find("}")) != definition.size()-1 )
    itkExceptionMacro(<< "Definition of a matrix must begin with { and end with } characters." << std::endl);

  //Get rid of { and } characters
  std::string def;
  for(int i=1; i<definition.size()-1; ++i)
    def.push_back(definition[i]);


  std::vector< std::vector<double> > mat;
OTB Bot's avatar
STYLE  
OTB Bot committed
270
  std::istringstream iss( def );
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
  std::string rows;
  while (std::getline( iss, rows, ';' ) )
  {
      mat.push_back(std::vector<double>(0));
      std::istringstream iss2( rows );
      std::string elmt;
        while (std::getline( iss2, elmt, ',' ) )
        {
            std::istringstream iss3( elmt );
            double val;
            iss3 >> val;
            mat.back().push_back(val);
        }
  }
  
286

287 288 289 290 291 292
  //Check dimensions of the matrix
  for (int i=0; i<mat.size()-1; i++)
    if (mat[i].size() != mat[i+1].size())
      itkExceptionMacro(<< "Each row must have the same number of cols : " << definition << std::endl);
  

293
  //Registration
294
  adhocStruct ahc;
295
  ahc.name = name;
296 297 298
  ahc.type = 7;
  ahc.info[0] = mat[0].size();  // Size x direction (matrix convention = cols)
  ahc.info[1] = mat.size();     // Size y direction (matrix convention = rows)
299 300 301 302 303 304

  ahc.value = ValueType(ahc.info[1],ahc.info[0],0.0);
  for(int i=0; i<mat.size(); i++)
    for(int j=0; j<mat[i].size(); j++)
      ahc.value.At(i,j)=mat[i][j];
 
305 306 307 308
  m_VAllowedVarNameAddedByUser.push_back(ahc);

}

309 310

template< typename TImage >
311
void BandMathXImageFilter<TImage>
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
::SetConstant(const std::string& name, double value)
{
  for(int i=0; i<m_VAllowedVarNameAddedByUser.size(); i++)
    if (name.compare(m_VAllowedVarNameAddedByUser[i].name) == 0)
      itkExceptionMacro(<< "Variable name '"<< name << "' already used." << std::endl);

  adhocStruct ahc;
  ahc.name = name;
  ahc.type = 7;
  ahc.value = value;
 
  m_VAllowedVarNameAddedByUser.push_back(ahc);

}


template< typename TImage >
329
void BandMathXImageFilter<TImage>
OTB Bot's avatar
STYLE  
OTB Bot committed
330
::ExportContext(const std::string& filename)
331 332 333 334 335
{

  std::vector< std::string > vectI,vectF,vectM, vectFinal;

  for(int i=0; i<m_VAllowedVarNameAddedByUser.size(); i++)
OTB Bot's avatar
STYLE  
OTB Bot committed
336
           {
337 338 339 340 341 342 343
        std::ostringstream iss;
        std::string        str;

        switch (m_VAllowedVarNameAddedByUser[i].value.GetType())
        {
          case 'i':
            iss << "#I " << m_VAllowedVarNameAddedByUser[i].name << " " << m_VAllowedVarNameAddedByUser[i].value.GetInteger();
OTB Bot's avatar
STYLE  
OTB Bot committed
344
            str=iss.str();
345 346 347
            vectI.push_back(str);
          break;
          case 'f':
OTB Bot's avatar
STYLE  
OTB Bot committed
348 349
            iss << "#F " << m_VAllowedVarNameAddedByUser[i].name << " " << m_VAllowedVarNameAddedByUser[i].value.GetFloat();
            str=iss.str();
350 351 352 353 354 355 356 357 358 359 360 361
            vectF.push_back(str);
          break;
          case 'c':
            itkExceptionMacro(<< "Complex numbers not supported." << std::endl);
          break;
          case 'm':
            iss << "#M " << m_VAllowedVarNameAddedByUser[i].name << " " << "{";
            for(int k=0; k<m_VAllowedVarNameAddedByUser[i].value.GetRows(); k++)
            {
              iss << " " << m_VAllowedVarNameAddedByUser[i].value.At(k,0);
              for(int p=1; p<m_VAllowedVarNameAddedByUser[i].value.GetCols(); p++)
                iss << " , " <<  m_VAllowedVarNameAddedByUser[i].value.At(k,p);
362
              iss << " ;";
363
            }
OTB Bot's avatar
STYLE  
OTB Bot committed
364
            str=iss.str();
365 366 367 368
            str.erase(str.size()-1);
            str.push_back('}');
            vectM.push_back(str);
          break;
OTB Bot's avatar
STYLE  
OTB Bot committed
369
        }
370 371 372

    }

OTB Bot's avatar
STYLE  
OTB Bot committed
373
  // Sorting : I F M and E
OTB Bot's avatar
STYLE  
OTB Bot committed
374
  for(int i=0; i<vectI.size(); ++i)
375
    vectFinal.push_back(vectI[i]);
OTB Bot's avatar
STYLE  
OTB Bot committed
376
  for(int i=0; i<vectF.size(); ++i)
377
    vectFinal.push_back(vectF[i]);
OTB Bot's avatar
STYLE  
OTB Bot committed
378
  for(int i=0; i<vectM.size(); ++i)
379 380 381 382 383 384 385 386 387 388
    vectFinal.push_back(vectM[i]);
  for(int i=0; i < m_Expression.size(); ++i)
    {
      std::ostringstream iss;
      iss << "#E " << m_Expression[i] << std::endl;
      std::string str=iss.str();
      vectFinal.push_back(str);
    }

  std::ofstream exportFile(filename.c_str(), std::ios::out | std::ios::trunc);
OTB Bot's avatar
STYLE  
OTB Bot committed
389
  if(exportFile)
390 391 392 393 394 395
    {
      for(int i=0; i<vectFinal.size(); ++i)
        exportFile << vectFinal[i] << std::endl;

      exportFile.close();
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
396
    else
397 398 399 400
      itkExceptionMacro(<< "Could not open " << filename << "." << std::endl);
}

template< typename TImage >
401
void BandMathXImageFilter<TImage>
OTB Bot's avatar
STYLE  
OTB Bot committed
402
::ImportContext(const std::string& filename)
403 404 405
{
  std::ifstream importFile(filename.c_str(), std::ios::in);

406 407
  std::string wholeline,line,sub,name,matrixdef;
  int pos,pos2,lineID=0,nbSuccesses=0;
408 409
  double value;

OTB Bot's avatar
STYLE  
OTB Bot committed
410
  if(importFile)
411 412
    {
 
413
      while(std::getline(importFile,wholeline))
414 415 416
      {
        lineID++;

417 418 419 420 421
        pos = wholeline.find_first_not_of(' ');
        
        if (pos != std::string::npos)
        {
          line = wholeline.substr(pos);
422

423 424 425 426
          if ( (line[0] == '#') && ((line[1] == 'I') || (line[1] == 'i') || (line[1] == 'F') || (line[1] == 'f')) )
            {
              
              pos = line.find_first_not_of(' ',2);
427

428 429
              if (pos == std::string::npos)
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID << " : please, set the name and the value of the constant." << std::endl);
430

431
              sub = line.substr(pos);
432

433
              pos = sub.find_first_of(' ');
OTB Bot's avatar
STYLE  
OTB Bot committed
434
              name = sub.substr(0,pos);
435 436
          
              if (sub.find_first_of('{',pos) != std::string::npos)
OTB Bot's avatar
STYLE  
OTB Bot committed
437
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID
438
                << " : symbol #F found, but find vector/matrix definition. Please, set an integer or a float number." << std::endl);
439

440 441
              if (sub.find_first_not_of(' ',pos) == std::string::npos )
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID << " : please, set the value of the constant." << std::endl)
442

443
              std::istringstream iss( sub.substr(pos) );
OTB Bot's avatar
STYLE  
OTB Bot committed
444
              iss >> value;
445 446 447
                
              SetConstant(name,value);
              nbSuccesses++;
448

449 450 451
            }
          else if ( (line[0] == '#') && ((line[1] == 'M') || (line[1] == 'm')) )
            {
452

453
              pos = line.find_first_not_of(' ',2);
454

455 456
              if (pos == std::string::npos)
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID << " : please, set the name and the definition of the vector/matrix." << std::endl);
457

458
              std::string sub = line.substr(pos);
459

460 461 462 463 464 465 466
              pos = sub.find_first_of(' ');
              name = sub.substr(0,pos);
              pos2 = sub.find_first_of('{');
              if (pos2 != std::string::npos)
                matrixdef = sub.substr(pos2);
              else
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID << " : symbol #M found, but couldn't not find vector/matrix definition." << std::endl);
467

468 469
              SetMatrix(name,matrixdef);
              nbSuccesses++;
470

471 472 473 474
            }
          else if ( (line[0] == '#') && ((line[1] == 'E') || (line[1] == 'e')) )
            {
              pos = line.find_first_not_of(' ',2);
475

476 477
              if (pos == std::string::npos)
                itkExceptionMacro(<< "In file '"<< filename << "', line " << lineID << " : symbol #E found, but couldn't not find any expression." << std::endl);
478

479
              sub = line.substr(pos);
480

481 482 483 484 485 486 487
              SetExpression(sub);
              nbSuccesses++;

            }
 
        }
      }//while
488 489

      importFile.close();
490 491 492 493 494


    if (nbSuccesses == 0)
      itkExceptionMacro(<< "No constant or expression could be set; please, ensure that the file '" << filename << "' is correct." << std::endl);

495
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
496
    else
497 498
      itkExceptionMacro(<< "Could not open " << filename << "." << std::endl);

499

500 501 502
}


503
template< typename TImage >
504
std::string BandMathXImageFilter<TImage>
505
::GetExpression(int IDExpression) const
506
{
OTB Bot's avatar
STYLE  
OTB Bot committed
507
  return m_Expression.at(IDExpression);
508 509
}

510

511
template< typename TImage >
512
std::vector<std::string> BandMathXImageFilter<TImage>
513
::GetVarNames() const
514
{
515 516 517 518 519
  std::vector<std::string> res;
  for(int y=0; y<m_VVarName.size(); y++)
    res.push_back(m_VVarName[y].name);

  return res;
520 521 522
}


523
template< typename TImage >
524
void BandMathXImageFilter<TImage>
525
::AddVariable(adhocStruct &ahc)
526 527 528 529 530 531 532 533
{
    bool found=false;
    for(int i=0; i<m_VVarName.size(); ++i)
      if (m_VVarName[i].name == ahc.name)
        found=true;

    if (!found)
      m_VVarName.push_back(ahc);
534
    
535 536 537
}


538
template< typename TImage >
539
void BandMathXImageFilter<TImage>
540
::PrepareParsers()
541 542
{

543
  if (m_Expression.size() == 0)
544 545
    itkExceptionMacro(<< "No expression set; please set at least one expression." << std::endl);

546

547
  // Generate variables names
548 549
  m_VVarName.clear();
  m_VNotAllowedVarName.clear();
550 551
  m_VFinalAllowedVarName.clear();

552
  // m_VFinalAllowedVarName = m_VAllowedVarNameAuto + m_VAllowedVarNameAddedByUser
OTB Bot's avatar
STYLE  
OTB Bot committed
553
  // m_VFinalAllowedVarName = variable names dictionary
554
  for(int i=0; i<m_VAllowedVarNameAddedByUser.size(); i++)
OTB Bot's avatar
STYLE  
OTB Bot committed
555
    m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAddedByUser[i]);
556 557
  for(int i=0; i<m_VAllowedVarNameAuto.size(); i++)
    m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAuto[i]);
558

559

560
  for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression) // For each expression
561
  {
562
      ParserType::Pointer dummyParser = ParserType::New();
563
      dummyParser->SetExpr(this->GetExpression(IDExpression));
564 565 566

      mup::var_maptype vmap = dummyParser->GetExprVar();
      for (mup::var_maptype::iterator item = vmap.begin(); item!=vmap.end(); ++item)
567
      {
568
        bool OK=false; int i=0;
569
        while( ( !OK ) && (i<m_VFinalAllowedVarName.size()) )
570
        {
571
          if (item->first == m_VFinalAllowedVarName[i].name)
572 573 574 575 576
            OK=true;
          else
            i++;
        }
        
OTB Bot's avatar
STYLE  
OTB Bot committed
577
        if (OK) {AddVariable(m_VFinalAllowedVarName[i]); }
578 579
        else {
                adhocStruct ahc;
OTB Bot's avatar
STYLE  
OTB Bot committed
580
                ahc.name = item->first;
581 582 583
                m_VNotAllowedVarName.push_back(ahc);
              }
      }
584

585
  }// At this step, m_VVarName has been built
586 587


588
  //Checking formulas consistency
589 590 591 592 593 594 595 596 597 598
  if (m_VNotAllowedVarName.size()>0)
  {
    std::stringstream sstm;
    sstm << "Following variables not allowed : ";
    for(int i=0; i<m_VNotAllowedVarName.size(); ++i)
      sstm << m_VNotAllowedVarName[i].name << " ";
    sstm << std::endl;
    itkExceptionMacro(<< sstm.str());
  }

599

600 601
  // Register variables for each parser (important : one parser per thread)
  m_VParser.clear();
602
  unsigned int nbThreads = this->GetNumberOfThreads();
603 604
  typename std::vector<ParserType::Pointer>::iterator        itParser;
  m_VParser.resize(nbThreads);
605 606
  for(itParser = m_VParser.begin(); itParser < m_VParser.end(); itParser++)
    {
607
      *itParser = ParserType::New();
608 609
    }

610 611
  // Important to remember that variables of m_VVarName come from a call of GetExprVar method
  // Only useful variables are allocated in this filter
OTB Bot's avatar
STYLE  
OTB Bot committed
612
  int nbVar = m_VVarName.size();
613 614

  m_StatsVarDetected.clear();
615

616 617 618
  //Reset
  for(int i=0; i<m_AImage.size(); ++i)
    m_AImage[i].clear();
619
  m_AImage.clear();
620

621
  m_AImage.resize(nbThreads);
622

623 624
  double initValue = 0.1;
  for(int i = 0; i < nbThreads; ++i) // For each thread
625
  {
626
 
627
    m_AImage[i].resize(nbVar); // For each variable
628

629
    for(int j=0; j < nbVar; ++j)
630
      {
631 632
        m_AImage[i][j].name = m_VVarName[j].name;
        m_AImage[i][j].type  = m_VVarName[j].type;
OTB Bot's avatar
STYLE  
OTB Bot committed
633
        for (int t=0; t<5; ++t)
634 635 636
          m_AImage[i][j].info[t]=m_VVarName[j].info[t];

       
637

638 639
        if ( (m_AImage[i][j].type == 0 ) || (m_AImage[i][j].type == 1) ) // indices (idxX & idxY)
        {
640
            m_AImage[i][j].value = ValueType(initValue);
641 642
        }

OTB Bot's avatar
STYLE  
OTB Bot committed
643
        if (m_AImage[i][j].type == 2) //imiPhyX
644 645 646 647 648 649 650 651 652 653 654 655 656 657
        {
          SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSpacing();
          m_AImage[i][j].value = ValueType(static_cast<double>(spacing[0]));
        }

        if (m_AImage[i][j].type == 3) //imiPhyY
        {
          SpacingType spacing = this->GetNthInput(m_AImage[i][j].info[0])->GetSpacing();
          m_AImage[i][j].value = ValueType(static_cast<double>(spacing[1]));
        }

        if (m_AImage[i][j].type == 4 ) // vector
        {
            unsigned int nbBands = this->GetNthInput(m_AImage[i][j].info[0])->GetNumberOfComponentsPerPixel();
658
            m_AImage[i][j].value = ValueType(1,nbBands,initValue);
659 660 661 662
        }

        if (m_AImage[i][j].type == 5 ) // pixel
        {
663
            m_AImage[i][j].value = ValueType(initValue);
664 665 666 667
        }

        if (m_AImage[i][j].type == 6 ) // neighborhood
        {
668 669 670
            m_AImage[i][j].value = ValueType(m_AImage[i][j].info[3],m_AImage[i][j].info[2],initValue);
        }

671
        if (m_AImage[i][j].type == 7 ) // user defined variables
672 673
        {

674 675 676
          for(int t=0; t<m_VAllowedVarNameAddedByUser.size(); t++)
            if (m_VAllowedVarNameAddedByUser[t].name.compare(m_AImage[i][j].name) == 0)
              m_AImage[i][j].value = m_VAllowedVarNameAddedByUser[t].value;
677

678
              
679 680
        }

681 682
        if (m_AImage[i][j].type == 8 ) // global stats
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
683 684
            m_AImage[i][j].value = ValueType(initValue);
            //m_AImage[i][j].info[0] = Image ID : useful to know which images must have their regions set to largest possible region (see GenerateInputRequestedRegion)
685
            bool found = false;
OTB Bot's avatar
STYLE  
OTB Bot committed
686
            for (int r=0; r<m_StatsVarDetected.size() && !found; r++)
687 688
                if (m_StatsVarDetected[r] == m_AImage[i][j].info[0])
                  found = true;
OTB Bot's avatar
STYLE  
OTB Bot committed
689 690
            if (!found)
              m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
691 692 693 694
        }


        //Register variable
OTB Bot's avatar
STYLE  
OTB Bot committed
695
        m_VParser.at(i)->DefineVar(m_AImage[i][j].name, &(m_AImage[i][j].value));
696

697

698 699 700
        initValue += 0.001;
        if (initValue>1.0)
          initValue=0.1;
701
      }
702
  }
703 704 705 706 707

}


template< typename TImage >
708
void BandMathXImageFilter< TImage >
709 710
::OutputsDimensions()
{
711 712 713

  this->SetNumberOfRequiredOutputs((int) m_Expression.size());

714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736
  m_outputsDimensions.clear();

  for(int i=0; i<m_Expression.size(); ++i)
  {
    m_VParser.at(0)->SetExpr(m_Expression[i]);
    ValueType value = m_VParser.at(0)->Eval();

    switch (value.GetType())
    {   //ValueType
        case 'i':
        m_outputsDimensions.push_back(1);
        break;

        case 'f':
        m_outputsDimensions.push_back(1);
        break;

        case 'c':
        itkExceptionMacro(<< "Complex numbers not supported." << std::endl);
        break;

        case 'm':
        mup::matrix_type vect = value.GetArray();
737
        if ( vect.GetRows() == 1 ) //Vector
OTB Bot's avatar
STYLE  
OTB Bot committed
738
          m_outputsDimensions.push_back(vect.GetCols());
739 740
        else //Matrix
          itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
741 742 743
        break;

    }
744

745
    //std::cout << "Type = " << value.GetType() << " dimension = " << m_outputsDimensions.back() << std::endl;
746 747 748 749 750 751
  }

}


template< typename TImage >
752
void BandMathXImageFilter< TImage >
753
::GenerateOutputInformation(void)
754
{
755 756
  Superclass::GenerateOutputInformation();

757 758 759
  typedef itk::ImageBase< TImage::ImageDimension > ImageBaseType;
  typename ImageBaseType::Pointer outputPtr;

760
  PrepareParsers();    // addition
761
  OutputsDimensions(); // addition
762 763

  int i=0;
764
  for ( itk::OutputDataObjectIterator it(this); !it.IsAtEnd(); i++, it++ )
765 766 767 768 769
    {
      // Check whether the output is an image of the appropriate
      // dimension (use ProcessObject's version of the GetInput()
      // method since it returns the input as a pointer to a
      // DataObject as opposed to the subclass version which
770
      // static_casts the input to an TImage).
771 772 773
      outputPtr = dynamic_cast< ImageBaseType * >( it.GetOutput() );

        if ( outputPtr )
774
          outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]);
775
    }
776

777 778
}

779

780
template< typename TImage >
781
void BandMathXImageFilter< TImage >
782 783 784 785 786
::GenerateInputRequestedRegion()
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

787
  for (int i=0; i<m_StatsVarDetected.size(); i++) //Must request largest possible regions (only for concerned images)
788 789 790 791 792 793
  {
      if ( m_StatsVarDetected[i] < this->GetNumberOfInputs() )
      {
          ImagePointer  inputPtr = const_cast<TImage *>(this->GetInput(m_StatsVarDetected[i]));
          inputPtr->SetRequestedRegionToLargestPossibleRegion();
      }
OTB Bot's avatar
STYLE  
OTB Bot committed
794 795
      else
        itkExceptionMacro(<< "Requested input #" << m_StatsVarDetected[i] << ", but only " << this->GetNumberOfInputs() << " inputs are available." << std::endl);
796 797 798 799 800
  }

}


801
template< typename TImage >
802
void BandMathXImageFilter<TImage>
803 804
::BeforeThreadedGenerateData()
{
805
  // Some useful variables
806 807 808
  unsigned int nbThreads = this->GetNumberOfThreads();
  unsigned int nbInputImages = this->GetNumberOfInputs();

809
  // Check if input image dimensions match
810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
  unsigned int inputSize[2];
  inputSize[0] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(0);
  inputSize[1] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(1);

  for(unsigned int p = 1; p < nbInputImages; p++)
    {
    if((inputSize[0] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0))
       || (inputSize[1] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1)))
      {
      itkExceptionMacro(<< "Input images must have the same dimensions." << std::endl
                        << "band #1 is [" << inputSize[0] << ";" << inputSize[1] << "]" << std::endl
                        << "band #" << p+1 << " is ["
                        << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0) << ";"
                        << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1) << "]");
      }
    }

OTB Bot's avatar
STYLE  
OTB Bot committed
827
  if (globalStatsDetected())
828
  // Must instantiate stats variables of the parsers
OTB Bot's avatar
STYLE  
OTB Bot committed
829
  // Note : at this stage, inputs have already been set to largest possible regions.
830 831
    for (unsigned int i=0; i<m_StatsVarDetected.size(); i++)
    {
832

833
      StreamingStatisticsVectorImageFilterPointerType filter = StreamingStatisticsVectorImageFilterType::New();
834 835 836 837 838

      filter->SetInput( this->GetNthInput(m_StatsVarDetected[i]) );
      filter->Update();

      PixelType pix; //Variable length vector
839
      MatrixType mat;
840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872
      
      for(int t=0; t<m_AImage.size(); t++) // for each thread
        for(int v=0; v<m_AImage[t].size(); v++) // for each variable
          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
          {
              switch ( m_AImage[t][v].info[2] ) // info[2] sub-type (see also SetNthInput method above)
              {
                  case 0: // mini

                    pix = filter->GetMinimum();
                    for (unsigned int b=0; b<pix.GetSize(); b++) // for each band
                      if (m_AImage[t][v].info[1] == b) // info[1] : band ID
                        m_AImage[t][v].value = pix[b];

                  break;

                  case 1: // maxi

                    pix = filter->GetMaximum();
                    for (unsigned int b=0; b<pix.GetSize(); b++) // for each band
                      if (m_AImage[t][v].info[1] == b) // info[1] : band ID
                        m_AImage[t][v].value = pix[b];

                  break;

                  case 2: // mean

                    pix = filter->GetMean();
                    for (unsigned int b=0; b<pix.GetSize(); b++) // for each band
                      if (m_AImage[t][v].info[1] == b) // info[1] : band ID
                        m_AImage[t][v].value = pix[b];

                  break;
873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892

                  break;

                  case 3: // sum

                    pix = filter->GetSum();
                    for (unsigned int b=0; b<pix.GetSize(); b++) // for each band
                      if (m_AImage[t][v].info[1] == b) // info[1] : band ID
                        m_AImage[t][v].value = pix[b];

                  break;

                  case 4: // stddev

                   mat = filter->GetCovariance();
                    for (unsigned int b=0; b<mat.Cols(); b++) // for each band
                      if (m_AImage[t][v].info[1] == b) // info[1] : band ID
                        m_AImage[t][v].value = mat(b,b);

                  break;
893 894
              }
          }
895 896 897
    }


898 899 900 901 902 903
  // Allocate and initialize the thread temporaries
  m_ThreadUnderflow.SetSize(nbThreads);
  m_ThreadUnderflow.Fill(0);
  m_ThreadOverflow.SetSize(nbThreads);
  m_ThreadOverflow.Fill(0);

904 905
}

906

907
template< typename TImage >
908
void BandMathXImageFilter<TImage>
909 910 911 912 913 914 915 916 917 918 919 920 921 922 923
::AfterThreadedGenerateData()
{
  unsigned int nbThreads = this->GetNumberOfThreads();
  unsigned int i;

  m_UnderflowCount = 0;
  m_OverflowCount = 0;

  // Accumulate counts for each thread
  for(i = 0; i < nbThreads; ++i)
    {
    m_UnderflowCount += m_ThreadUnderflow[i];
    m_OverflowCount += m_ThreadOverflow[i];
    }

OTB Bot's avatar
STYLE  
OTB Bot committed
924
  if((m_UnderflowCount != 0) || (m_OverflowCount!=0))
925 926 927 928 929 930 931
  {
    std::stringstream sstm;
    sstm << std::endl
        << "The Following Parsed Expression  :  ";
    for(int t=0; t<m_Expression.size(); ++t)
        sstm << this->GetExpression(t) << std::endl;
    sstm << "Generated " << m_UnderflowCount << " Underflow(s) "
932 933
        << "And " << m_OverflowCount        << " Overflow(s) "   << std::endl
        << "The Parsed Expression, The Inputs And The Output "
934 935 936 937
        << "Type May Be Incompatible !";

    otbWarningMacro(<< sstm.str());
  }
938 939 940
}

template< typename TImage >
941
void BandMathXImageFilter<TImage>
942 943 944
::ThreadedGenerateData(const ImageRegionType& outputRegionForThread,
           itk::ThreadIdType threadId)
{
945 946

  ValueType value;
947 948
  unsigned int nbInputImages = this->GetNumberOfInputs();

OTB Bot's avatar
STYLE  
OTB Bot committed
949 950 951
  //----------------- --------- -----------------//
  //----------------- Iterators -----------------//
  //----------------- --------- -----------------//
952 953 954
  typedef itk::ImageRegionConstIterator<TImage> ImageRegionConstIteratorType;
  std::vector< ImageRegionConstIteratorType > Vit;
  Vit.resize(nbInputImages);
955
  for(int j=0; j < nbInputImages; ++j)
956 957
    Vit[j] = ImageRegionConstIteratorType (this->GetNthInput(j), outputRegionForThread);
    
958

959 960 961
  std::vector< ImageRegionConstIteratorType > VoutIt;
  VoutIt.resize(m_Expression.size());
  for(int j=0; j < VoutIt.size(); ++j)
962 963
    VoutIt[j] = ImageRegionConstIteratorType (this->GetOutput(j), outputRegionForThread);
    
964

965 966 967
  //Special case : neighborhoods
  std::vector< itk::ConstNeighborhoodIterator<TImage> > VNit;
  for(int j=0; j<m_VVarName.size(); ++j)
OTB Bot's avatar
STYLE  
OTB Bot committed
968
    if (m_VVarName[j].type == 6)
969 970 971 972 973
     {
        typename itk::ConstNeighborhoodIterator<TImage>::RadiusType radius;
        radius[0]=(int) ((m_VVarName[j].info[2]-1)/2); // Size x direction (otb convention)
        radius[1]=(int) ((m_VVarName[j].info[3]-1)/2); // Size y direction (otb convention)
        VNit.push_back( itk::ConstNeighborhoodIterator<TImage>(radius, this->GetNthInput(m_VVarName[j].info[0]),outputRegionForThread)); // info[0] = Input image ID
OTB Bot's avatar
STYLE  
OTB Bot committed
974
        VNit.back().NeedToUseBoundaryConditionOn();
975 976 977
     }


978
  // Support progress methods/callbacks
979 980
  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());

981

OTB Bot's avatar
STYLE  
OTB Bot committed
982 983 984
  //----------------- --------------------- -----------------//
  //----------------- Variable affectations -----------------//
  //----------------- --------------------- -----------------//
985 986
  for(int j=0; j < nbInputImages; ++j)       {  Vit[j].GoToBegin();     }
  for(int j=0; j < m_Expression.size(); ++j) {  VoutIt[j].GoToBegin();  }
OTB Bot's avatar
STYLE  
OTB Bot committed
987
  for(int j=0; j < VNit.size(); ++j)         {  VNit[j].GoToBegin();    }
988

989 990 991
  while(!Vit.at(0).IsAtEnd()) // For each pixel
  {

OTB Bot's avatar
STYLE  
OTB Bot committed
992
    int ngbhNameIndex=0; int index;
993
    for(int j=0; j < m_AImage[threadId].size(); ++j) // For each variable, perform a copy
994 995
    {

OTB Bot's avatar
STYLE  
OTB Bot committed
996
       switch (m_AImage[threadId][j].type)
997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014
        {

          case 0 : //idxX
            m_AImage[threadId][j].value = static_cast<double>(Vit[0].GetIndex()[0]);
          break;

          case 1 : //idxY
            m_AImage[threadId][j].value = static_cast<double>(Vit[0].GetIndex()[1]);
          break;

          case 2 : //Spacing X (imiPhyX)
            //Nothing to do (already set inside BeforeThreadedGenerateData)"
          break;

          case 3 : //Spacing Y (imiPhyY)
            //Nothing to do (already set inside BeforeThreadedGenerateData)"
          break;

OTB Bot's avatar
STYLE  
OTB Bot committed
1015
          case 4 : //vector
1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027
            // m_AImage[threadId][j].info[0] : Input image #ID
            for(int p=0; p < m_AImage[threadId][j].value.GetCols(); ++p)
              m_AImage[threadId][j].value.At(0,p) = Vit[m_AImage[threadId][j].info[0]].Get()[p];
          break;

          case 5 : //pixel
            // m_AImage[threadId][j].info[0] : Input image #ID
            // m_AImage[threadId][j].info[1] : Band #ID
            m_AImage[threadId][j].value = Vit[m_AImage[threadId][j].info[0]].Get()[m_AImage[threadId][j].info[1]];
          break;

          case 6 : //neighborhood
1028
    
1029 1030 1031 1032
          // m_AImage[threadId][j].info[1] : Band #ID
          if (m_AImage[threadId][j].info[2]*m_AImage[threadId][j].info[3] != VNit[ngbhNameIndex].Size() )
            itkExceptionMacro(<< "Size of muparserx variable is different from its related otb neighborhood iterator")

1033 1034 1035
          index=0;
          for(int rows=0; rows<m_AImage[threadId][j].info[3]; ++rows)
            for(int cols=0; cols<m_AImage[threadId][j].info[2]; ++cols)
1036 1037 1038 1039 1040 1041 1042 1043 1044
              {
                m_AImage[threadId][j].value.At(rows,cols) = VNit[ngbhNameIndex].GetPixel(index)[m_AImage[threadId][j].info[1]];
                index++;
              }

          ngbhNameIndex++;
          break;

          case 7 :
1045
          //Nothing to do : user defined variable or constant, which have already been set
1046 1047
          break;

1048 1049 1050 1051
          case 8 :
          //Nothing to do : variable has already been set inside BeforeThreadedGenerateData method (see above)
          break;

1052 1053 1054 1055
          default :
            itkExceptionMacro(<< "Type of the variable is unknown");
          break;
        }
1056
    }//End while
1057

1058

OTB Bot's avatar
STYLE  
OTB Bot committed
1059 1060 1061
  //----------------- ----------- -----------------//
  //----------------- Evaluations -----------------//
  //----------------- ----------- -----------------//
1062 1063
  for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression)
  {
1064

1065
        m_VParser[threadId]->SetExpr(m_Expression[IDExpression]);
1066

1067
        value = m_VParser[threadId]->Eval();
1068

1069 1070 1071 1072 1073
        switch (value.GetType())
        {   //ValueType
            case 'i':
            VoutIt[IDExpression].Get()[0] = value.GetInteger();
            break;
1074

1075 1076 1077
            case 'f':
            VoutIt[IDExpression].Get()[0] = value.GetFloat();
            break;
1078

1079
            case 'c':
1080
            itkExceptionMacro(<< "Complex numbers are not supported." << std::endl);
1081
            break;
1082

1083 1084
            case 'm':
            mup::matrix_type vect = value.GetArray();
1085

1086 1087
            if ( vect.GetRows() == 1 ) //Vector
              for(int p=0; p<vect.GetCols(); ++p)
OTB Bot's avatar
STYLE  
OTB Bot committed
1088
                VoutIt[IDExpression].Get()[p] = vect.At(0,p).GetFloat();
1089 1090 1091
            else //Matrix
              itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
            break;
1092
        }
1093

1094
 
OTB Bot's avatar
STYLE  
OTB Bot committed
1095
        //----------------- Pixel affectations -----------------//
1096
        for(int p=0; p<VoutIt[IDExpression].Get().GetSize(); ++p)
1097
        {
1098
            // Case value is equal to -inf or inferior to the minimum value
1099 1100
            // allowed by the PixelValueType cast
            if (VoutIt[IDExpression].Get()[p] < double(itk::NumericTraits<PixelValueType>::NonpositiveMin()))
1101
            {
1102
                VoutIt[IDExpression].Get()[p] = itk::NumericTraits<PixelValueType>::NonpositiveMin();
1103 1104 1105
                m_ThreadUnderflow[threadId]++;
            }
            // Case value is equal to inf or superior to the maximum value
1106 1107
            // allowed by the PixelValueType cast
            else if (VoutIt[IDExpression].Get()[p] > double(itk::NumericTraits<PixelValueType>::max()))
1108
            {
1109
               VoutIt[IDExpression].Get()[p] = itk::NumericTraits<PixelValueType>::max();
1110
               m_ThreadOverflow[threadId]++;
OTB Bot's avatar
STYLE  
OTB Bot committed
1111
            }
1112
        }
1113
    }
1114

1115 1116
    for(int j=0; j < nbInputImages; ++j)        {   ++Vit[j];    }
    for(int j=0; j < m_Expression.size(); ++j)  {   ++VoutIt[j]; }
OTB Bot's avatar
STYLE  
OTB Bot committed
1117
    for(int j=0; j < VNit.size(); ++j)          {   ++VNit[j];   }
1118

1119 1120 1121
    progress.CompletedPixel();
  }

1122 1123 1124 1125 1126
}

}// end namespace otb

#endif