WIP: otbBandMathImageFilterX / otbParserX / otbParserXPlugins (all based on muParserX)

CMake/FindMuParserX.cmake

+# - Find MuParserX
+# Find the native MuParser includes and library
+#
+#   MUPARSER_FOUND       - True if MuParser found.
+#   MUPARSER_INCLUDE_DIR - where to find tinyxml.h, etc.
+#   MUPARSER_LIBRARIES   - List of libraries when using MuParser.
+#
+
+if( MUPARSERX_INCLUDE_DIR )
+    # Already in cache, be silent
+    set( MuParserX_FIND_QUIETLY TRUE )
+endif()
+
+find_path( MUPARSERX_INCLUDE_DIR mpParser.h
+           PATH_SUFFIXES mpParser )
+
+find_library( MUPARSERX_LIBRARIES
+              NAMES muparserx
+              PATH_SUFFIXES muparserx )
+
+# handle the QUIETLY and REQUIRED arguments and set MUPARSER_FOUND to TRUE if
+# all listed variables are TRUE
+include( FindPackageHandleStandardArgs )
+FIND_PACKAGE_HANDLE_STANDARD_ARGS( MuParserX DEFAULT_MSG MUPARSERX_LIBRARIES MUPARSERX_INCLUDE_DIR )
+
+mark_as_advanced( MUPARSERX_INCLUDE_DIR MUPARSERX_LIBRARIES )

CMake/ImportMuParserX.cmake

+message(STATUS "Importing MuParserX...")
+find_package(MuParserX)
+
+if(MUPARSERX_FOUND)
+  option(OTB_USE_EXTERNAL_MUPARSERX "Use external MuParserX library." ON)
+else()
+  option(OTB_USE_EXTERNAL_MUPARSERX "Use external MuParserX library." OFF)
+endif()
+mark_as_advanced(OTB_USE_EXTERNAL_MUPARSERX)
+
+if(OTB_USE_EXTERNAL_MUPARSERX)
+  if(MUPARSERX_FOUND)
+    # Starting with muparser 2.0.0,
+    #  intrinsic operators "and", "or", "xor" have been removed
+    #  and intrinsic operators "&&" and "||" have been introduced as replacements
+    set(CMAKE_REQUIRED_INCLUDES ${MUPARSERX_INCLUDE_DIR})
+    set(CMAKE_REQUIRED_LIBRARIES ${MUPARSERX_LIBRARIES})
+    unset(CMAKE_REQUIRED_FLAGS)
+    unset(CMAKE_REQUIRED_DEFINES)
+    #file(READ ${CMAKE_CURRENT_SOURCE_DIR}/CMake/otbTestMuParserHasCxxLogicalOperators.cxx
+    #    OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS_SOURCEFILE)
+    #check_cxx_source_runs(
+    #  "${OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS_SOURCEFILE}"
+    #  OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS
+    #  )
+    unset(CMAKE_REQUIRED_INCLUDES)
+    unset(CMAKE_REQUIRED_LIBRARIES)
+    unset(CMAKE_REQUIRED_FLAGS)
+    unset(CMAKE_REQUIRED_DEFINES)
+
+    #if(OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS)
+    #  message(STATUS "  MuParser version is >= 2.0.0 : "
+    #                 "uses '&&' and '||' logical operators, and C++ like "
+    #                 "ternary if-then-else operator")
+    #else()
+    #  message(STATUS "  MuParser version is < 2.0.0  : "
+    #                 "uses 'and' and 'or' logical operators, and ternary "
+    #                 "operator 'if( ; ; )'")
+    #endif()
+  else()
+    message(FATAL_ERROR "Can't build OTB without MuParserX. Instal it "
+                        "on your system, or disable the option "
+                        "OTB_USE_EXTERNAL_MUPARSERX to use internal one")
+  endif()
+else()
+  set(MUPARSERX_LIBRARIES otbmuparserx)
+  #unset(OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS) # internal version is < 2.0.0
+  message(STATUS "  Using MuParserX internal version ")
+endif()
+

CMake/otbIncludeDirectories.cmake

@@ -95,6 +95,14 @@ else()
       ${OTB_SOURCE_DIR}/Utilities/otbmuparser)
 endif()
 
+if(OTB_USE_EXTERNAL_MUPARSERX)
+  set(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
+      ${MUPARSERX_INCLUDE_DIR})
+else()
+  set(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
+      ${OTB_SOURCE_DIR}/Utilities/otbmuparserx)
+endif()
+
 if(OTB_USE_EXTERNAL_LIBKML)
   set(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
       ${LIBKML_INCLUDE_DIR})

CMakeLists.txt

@@ -216,6 +216,7 @@ include(ImportQt4)
 include(ImportSiftFast)
 include(ImportTinyXML)
 include(ImportMuParser)
+include(ImportMuParserX)
 include(ImportLibKML)
 include(ImportOpenCV)
 #-----------------------------------------------------------------------------

Code/BasicFilters/CMakeLists.txt

@@ -3,8 +3,9 @@
 file(GLOB OTBBasicFilters_SRCS "*.cxx" )
 
 add_library(OTBBasicFilters ${OTBBasicFilters_SRCS})
-target_link_libraries(OTBBasicFilters OTBCommon otbedison ${MUPARSER_LIBRARIES})
+target_link_libraries(OTBBasicFilters OTBCommon otbedison ${MUPARSER_LIBRARIES} ${MUPARSERX_LIBRARIES})
 
+ 
 # Explicit link to fftw is required
 # The link of ITKFFT to fftw does not propagate well :
 # in an external project using OTB with FFT features,

Code/BasicFilters/otbBandMathImageFilterX.h

+/*=========================================================================
+
+  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.
+
+=========================================================================*/
+
+#ifndef __otbBandMathImageFilterX_h
+#define __otbBandMathImageFilterX_h
+
+#include "itkInPlaceImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkArray.h"
+
+#include "otbParserX.h"
+
+namespace otb
+{
+/** \class BandMathImageFilterX
+ * \brief Performs a mathematical operation on the input images
+ * according to the formula specified by the user.
+ *
+ * This filter is based on the mathematical parser library muParser.
+ * The built in functions and operators list is available at:
+ * http://muparser.sourceforge.net/mup_features.html#idDef2
+ *
+ * OTB additional functions:
+ * ndvi(r, niri)
+ *
+ * OTB additional constants:
+ * e - log2e - log10e - ln2 - ln10 - pi - euler
+ *
+ * In order to use this filter, at least one input image is to be
+ * set. An associated variable name can be specified or not by using
+ * the corresponding SetNthInput method. For the nth input image, if
+ * no associated variable name has been spefified, a default variable
+ * name is given by concatenating the letter "b" (for band) and the
+ * corresponding input index.
+ * Next step is to set the expression according to the variable
+ * names. For example, in the default case with three input images the
+ * following expression is valid :
+ * "ndvi(b1, b2)*b3"
+ *
+ * As an additional functionality, the filter also granted access to
+ * indexes information under special virtual bands named idxX, idxY
+ * for the images indexes and idxPhyX, idxPhyY for the physical
+ * indexes.
+ * It allows the user to perform, for example a spatial processing
+ * aiming to suppress a determined area :
+ * "if(sqrt((idxPhyX-105.3)*(idxPhyX-105.3)+
+ *          (idxPhyY-207.1)*(idxPhyY-207.1))>100, b1, 0)"
+ * This expression replace the physical zone around the point of
+ * physical index (105.3; 207.1) by a black area
+ * This functionality assumes that all the band involved have the same
+ * spacing and origin.
+ *
+ *
+ * \sa Parser
+ *
+ * \ingroup Streamed
+ * \ingroup Threaded
+ */
+
+template< class TImage >
+class ITK_EXPORT BandMathImageFilterX
+  : public itk::InPlaceImageFilter< TImage >
+{
+public:
+  /** Standard class typedefs. */
+  typedef BandMathImageFilterX< TImage >                 Self;
+  typedef itk::InPlaceImageFilter< TImage >             Superclass;
+  typedef itk::SmartPointer< Self >                     Pointer;
+  typedef itk::SmartPointer< const Self >               ConstPointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(BandMathImageFilterX, InPlaceImageFilter);
+
+  /** Some convenient typedefs. */
+  typedef TImage                                  ImageType;
+  typedef typename ImageType::ConstPointer        ImagePointer;
+  typedef typename ImageType::RegionType          ImageRegionType;
+  typedef typename ImageType::PixelType           PixelType;
+  typedef typename ImageType::IndexType           IndexType;
+  typedef typename ImageType::PointType           OrigineType;
+  typedef typename ImageType::SpacingType         SpacingType;
+  typedef ParserX                                  ParserType;
+  typedef typename ParserType::ValueType          ValueType;
+
+  /** Set the nth filter input with or without a specified associated variable name */
+  void SetNthInput( unsigned int idx, const ImageType * image);
+  void SetNthInput( unsigned int idx, const ImageType * image, const std::string& varName);
+
+  /** Change the nth filter input associated variable name */
+  void SetNthInputName(unsigned int idx, const std::string& expression);
+
+  /** Set the expression to be parsed */
+  void SetExpression(const std::string& expression);
+
+  /** Return the expression to be parsed */
+  std::string GetExpression() const;
+
+  /** Return the nth filter input associated variable name */
+  std::string GetNthInputName(unsigned int idx) const;
+
+  /** Return a pointer on the nth filter input */
+  ImageType * GetNthInput(unsigned int idx);
+
+protected :
+  BandMathImageFilterX();
+  virtual ~BandMathImageFilterX();
+  virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  void BeforeThreadedGenerateData();
+  void ThreadedGenerateData(const ImageRegionType& outputRegionForThread, itk::ThreadIdType threadId );
+  void AfterThreadedGenerateData();
+
+private :
+  BandMathImageFilterX(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+  std::string                           m_Expression;
+  std::vector<ParserType::Pointer>      m_VParser;
+  std::vector< std::vector<ValueType> >    m_AImage;
+  std::vector< std::string >            m_VVarName;
+  unsigned int                          m_NbVar;
+
+  SpacingType                           m_Spacing;
+  OrigineType                           m_Origin;
+
+  long                                  m_UnderflowCount;
+  long                                  m_OverflowCount;
+  itk::Array<long>                      m_ThreadUnderflow;
+  itk::Array<long>                      m_ThreadOverflow;
+};
+
+}//end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbBandMathImageFilterX.txx"
+#endif
+
+#endif

Code/BasicFilters/otbBandMathImageFilterX.txx

+/*=========================================================================
+
+  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.
+
+=========================================================================*/
+#ifndef __otbBandMathImageFilterX_txx
+#define __otbBandMathImageFilterX_txx
+#include "otbBandMathImageFilterX.h"
+
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkNumericTraits.h"
+#include "itkProgressReporter.h"
+#include "otbMacro.h"
+
+
+#include <iostream>
+#include <string>
+
+namespace otb
+{
+
+/** Constructor */
+template <class TImage>
+BandMathImageFilterX<TImage>
+::BandMathImageFilterX()
+{
+  //This number will be incremented each time an image
+  //is added over the one minimumrequired
+  this->SetNumberOfRequiredInputs( 1 );
+  this->InPlaceOff();
+
+  m_UnderflowCount = 0;
+  m_OverflowCount = 0;
+  m_ThreadUnderflow.SetSize(1);
+  m_ThreadOverflow.SetSize(1);
+}
+
+/** Destructor */
+template <class TImage>
+BandMathImageFilterX<TImage>
+::~BandMathImageFilterX()
+{
+}
+
+template <class TImage>
+void BandMathImageFilterX<TImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+
+  os << indent << "Expression: "      << m_Expression                  << std::endl;
+  os << indent << "Computed values follow:"                            << std::endl;
+  os << indent << "UnderflowCount: "  << m_UnderflowCount              << std::endl;
+  os << indent << "OverflowCount: "   << m_OverflowCount               << std::endl;
+  os << indent << "itk::NumericTraits<PixelType>::NonpositiveMin()  :  "
+               << itk::NumericTraits<PixelType>::NonpositiveMin()      << std::endl;
+  os << indent << "itk::NumericTraits<PixelType>::max()  :             "
+         << itk::NumericTraits<PixelType>::max()                 << std::endl;
+}
+
+template <class TImage>
+void BandMathImageFilterX<TImage>
+::SetNthInput(unsigned int idx, const ImageType * image)
+{
+  this->SetInput(idx, const_cast<TImage *>( image ));
+  unsigned int nbInput = this->GetNumberOfInputs();
+  m_VVarName.resize(nbInput+4);
+  std::ostringstream varName;
+  varName << "b" << nbInput;
+  m_VVarName[idx] = varName.str();
+
+  m_VVarName[idx+1] = "idxX";
+  m_VVarName[idx+2] = "idxY";
+  m_VVarName[idx+3] = "idxPhyX";
+  m_VVarName[idx+4] = "idxPhyY";
+}
+
+template <class TImage>
+void BandMathImageFilterX<TImage>
+::SetNthInput(unsigned int idx, const ImageType * image, const std::string& varName)
+{
+  this->SetInput(idx, const_cast<TImage *>( image ));
+  m_VVarName.resize(this->GetNumberOfInputs()+4);
+  m_VVarName[idx] = varName;
+
+  m_VVarName[idx+1] = "idxX";
+  m_VVarName[idx+2] = "idxY";
+  m_VVarName[idx+3] = "idxPhyX";
+  m_VVarName[idx+4] = "idxPhyY";
+}
+
+template <class TImage>
+void BandMathImageFilterX<TImage>
+::SetNthInputName(unsigned int idx, const std::string& varName)
+{
+  m_VVarName[idx] = varName;
+}
+
+template <typename TImage>
+TImage * BandMathImageFilterX<TImage>
+::GetNthInput(unsigned int idx)
+{
+  return const_cast<TImage *>(this->GetInput(idx));
+}
+
+template< typename TImage >
+void BandMathImageFilterX<TImage>
+::SetExpression(const std::string& expression)
+{
+  if (m_Expression != expression)
+    m_Expression = expression;
+  this->Modified();
+}
+
+template< typename TImage >
+std::string BandMathImageFilterX<TImage>
+::GetExpression() const
+{
+  return m_Expression;
+}
+
+template< typename TImage >
+std::string BandMathImageFilterX<TImage>
+::GetNthInputName(unsigned int idx) const
+{
+  return m_VVarName.at(idx);
+}
+
+template< typename TImage >
+void BandMathImageFilterX<TImage>
+::BeforeThreadedGenerateData()
+{
+  typename std::vector<ParserType::Pointer>::iterator        itParser;
+  typename std::vector< std::vector<PixelType> >::iterator   itVImage;
+  unsigned int nbThreads = this->GetNumberOfThreads();
+  unsigned int nbInputImages = this->GetNumberOfInputs();
+  unsigned int nbAccessIndex = 4; //to give access to image and physical index
+  unsigned int i, j;
+  unsigned int inputSize[2];
+  std::vector< std::string > tmpIdxVarNames;
+
+  tmpIdxVarNames.resize(nbAccessIndex);
+  tmpIdxVarNames.at(0) = "idxX";
+  tmpIdxVarNames.at(1) = "idxY";
+  tmpIdxVarNames.at(2) = "idxPhyX";
+  tmpIdxVarNames.at(3) = "idxPhyY";
+
+  // Check if input image dimensions matches
+  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) << "]");
+      }
+    }
+
+  // Store images specs
+  m_Spacing = this->GetNthInput(0)->GetSpacing();
+  m_Origin = this->GetNthInput(0)->GetOrigin();
+
+  // Allocate and initialize the thread temporaries
+  m_ThreadUnderflow.SetSize(nbThreads);
+  m_ThreadUnderflow.Fill(0);
+  m_ThreadOverflow.SetSize(nbThreads);
+  m_ThreadOverflow.Fill(0);
+  m_VParser.resize(nbThreads);
+  m_AImage.resize(nbThreads);
+  m_NbVar = nbInputImages+nbAccessIndex;
+  m_VVarName.resize(m_NbVar);
+
+  for(itParser = m_VParser.begin(); itParser < m_VParser.end(); itParser++)
+    {
+    *itParser = ParserType::New();
+    }
+
+  for(i = 0; i < nbThreads; ++i)
+    {
+    m_AImage.at(i).resize(m_NbVar);
+    m_VParser.at(i)->SetExpr(m_Expression);
+
+    for(j=0; j < nbInputImages; ++j)
+      {
+      m_VParser.at(i)->DefineVar(m_VVarName.at(j), &m_AImage.at(i).at(j)); 
+      }
+
+    for(j=nbInputImages; j < nbInputImages+nbAccessIndex; ++j)
+      {
+      m_VVarName.at(j) = tmpIdxVarNames.at(j-nbInputImages);
+      m_VParser.at(i)->DefineVar(m_VVarName.at(j), &m_AImage.at(i).at(j));
+      }
+
+    }
+
+ 
+}
+
+template< typename TImage >
+void BandMathImageFilterX<TImage>
+::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];
+    }
+
+  if((m_UnderflowCount != 0) || (m_OverflowCount!=0))
+    otbWarningMacro(<< std::endl
+        << "The Following Parsed Expression  :  "
+        << this->GetExpression()                                 << std::endl
+        << "Generated " << m_UnderflowCount << " Underflow(s) "
+        << "And " << m_OverflowCount        << " Overflow(s) "   << std::endl
+        << "The Parsed Expression, The Inputs And The Output "
+        << "Type May Be Incompatible !");
+}
+
+template< typename TImage >
+void BandMathImageFilterX<TImage>
+::ThreadedGenerateData(const ImageRegionType& outputRegionForThread,
+           itk::ThreadIdType threadId)
+{
+  double value;
+  unsigned int j;
+  unsigned int nbInputImages = this->GetNumberOfInputs();
+
+  typedef itk::ImageRegionConstIterator<TImage> ImageRegionConstIteratorType;
+
+  std::vector< ImageRegionConstIteratorType > Vit;
+  Vit.resize(nbInputImages);
+  for(j=0; j < nbInputImages; ++j)
+    {
+    Vit[j] = ImageRegionConstIteratorType (this->GetNthInput(j), outputRegionForThread);
+    }
+
+  itk::ImageRegionIterator<TImage> ot (this->GetOutput(), outputRegionForThread);
+
+  // support progress methods/callbacks
+  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+  while(!Vit.at(0).IsAtEnd())
+    {
+    for(j=0; j < nbInputImages; ++j)
+      {
+      m_AImage.at(threadId).at(j) = static_cast<double>(Vit.at(j).Get());
+      }
+
+    // Image Indexes
+    for(j=0; j < 2; ++j)
+      {
+      m_AImage.at(threadId).at(nbInputImages+j)   = static_cast<double>(Vit.at(0).GetIndex()[j]);
+      }
+    for(j=0; j < 2; ++j)
+      {
+      m_AImage.at(threadId).at(nbInputImages+2+j) = static_cast<double>(m_Origin[j])
+        +static_cast<double>(Vit.at(0).GetIndex()[j]) * static_cast<double>(m_Spacing[j]);
+      }
+
+    try
+      {
+      value = m_VParser.at(threadId)->Eval().GetFloat();
+      }
+      catch(itk::ExceptionObject& err)
+      {
+      itkExceptionMacro(<< err);
+      }
+
+    // Case value is equal to -inf or inferior to the minimum value
+    // allowed by the pixelType cast
+    if (value < double(itk::NumericTraits<PixelType>::NonpositiveMin()))
+      {
+      ot.Set(itk::NumericTraits<PixelType>::NonpositiveMin());
+      m_ThreadUnderflow[threadId]++;
+      }
+    // Case value is equal to inf or superior to the maximum value
+    // allowed by the pixelType cast
+    else if (value > double(itk::NumericTraits<PixelType>::max()))
+      {
+      ot.Set(itk::NumericTraits<PixelType>::max());
+      m_ThreadOverflow[threadId]++;
+      }
+    else
+      {
+      ot.Set(static_cast<PixelType>(value));
+      }
+
+    for(j=0; j < nbInputImages; ++j)
+      {
+      ++(Vit.at(j));
+      }
+
+    ++ot;
+
+    progress.CompletedPixel();
+    }
+}
+
+}// end namespace otb
+
+#endif

Code/Common/CMakeLists.txt

@@ -12,7 +12,7 @@ add_library(OTBCommon ${OTBCommon_SRCS})
 #    PROPERTIES
 #    LINK_INTERFACE_LIBRARIES ""
 #)