From 5d3bfa864b1250cea2e179c1775ce3413a810344 Mon Sep 17 00:00:00 2001
From: Jonathan Guinet <jonathan.guinet@c-s.fr>
Date: Thu, 7 Apr 2011 10:01:15 +0200
Subject: [PATCH] ENH MuPaser based Functor for binary mask generation and
connected component segmentation application
---
.../otbConnectedComponentMuParserFunctor.h | 194 ++++++++++++++++++
Code/BasicFilters/otbMaskMuParserFilter.h | 142 +++++++++++++
Code/BasicFilters/otbMaskMuParserFilter.txx | 169 +++++++++++++++
Code/BasicFilters/otbMaskMuParserFunctor.h | 185 +++++++++++++++++
4 files changed, 690 insertions(+)
create mode 100644 Code/BasicFilters/otbConnectedComponentMuParserFunctor.h
create mode 100644 Code/BasicFilters/otbMaskMuParserFilter.h
create mode 100644 Code/BasicFilters/otbMaskMuParserFilter.txx
create mode 100644 Code/BasicFilters/otbMaskMuParserFunctor.h
diff --git a/Code/BasicFilters/otbConnectedComponentMuParserFunctor.h b/Code/BasicFilters/otbConnectedComponentMuParserFunctor.h
new file mode 100644
index 0000000000..98080cf6fa
--- /dev/null
+++ b/Code/BasicFilters/otbConnectedComponentMuParserFunctor.h
@@ -0,0 +1,194 @@
+/*=========================================================================
+
+ 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 __otbConnectedComponentMuParserFunctor_h
+#define __otbConnectedComponentMuParserFunctor_h
+
+#include "otbParser.h"
+#include "otbMacro.h"
+#include "itkMacro.h"
+
+
+#include <vnl/algo/vnl_lsqr.h>
+#include <vnl/vnl_sparse_matrix_linear_system.h>
+#include <vnl/vnl_least_squares_function.h>
+
+
+namespace otb
+{
+/** \class ConnectedComponentMuParserFunctor
+ * \brief functor used as input to itk connected component segmentation module.
+ *
+ * This functor 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
+ *
+ * image pixels :
+ * p1bX, p2bX, where X denote band index.
+ *
+ * Specific Functor additional input :
+ * distance : euclidian distance
+ * spectralAngle : NOT IMPLEMENTED
+ *
+ * OTB additional functions:
+ * ndvi(r, niri)
+ *
+ * OTB additional constants:
+ * e - log2e - log10e - ln2 - ln10 - pi - euler
+ *
+ *
+ * \sa Parser
+ *
+ */
+
+/** \class ConnectedComponentMuParserFunctor
+ * \brief This functor use MuParser as criteria for itk connected component module
+ */
+ namespace Functor
+ {
+
+ template<class TInput>
+ class ITK_EXPORT ConnectedComponentMuParserFunctor
+ {
+
+ public:
+ typedef Parser ParserType;
+ typedef ConnectedComponentMuParserFunctor Self;
+
+ std::string GetNameOfClass()
+ {
+ return "ConnectedComponentMuParserFunctor";
+ }
+
+ inline bool operator()(const TInput &p1, const TInput &p2)
+ {
+
+
+ double value;
+
+ if(p1.GetSize()!=m_NbOfBands)
+ {
+ this->SetNumberOfBands(p1.GetSize());
+ }
+
+ // we fill the buffer
+ for(unsigned int i=0; i<m_NbOfBands; i++)
+ {
+ m_AImageP1[i]= static_cast<double> (p1[i]);
+ m_AImageP2[i]= static_cast<double> (p2[i]);
+ }
+
+ m_Distance = 0.0;
+ for(unsigned int i = 0; i<m_NbOfBands;++i)
+ {
+ m_Distance +=(p1[i]-p2[i])*(p1[i]-p2[i]);
+ }
+ m_Distance = vcl_sqrt(m_Distance);
+
+ // cast
+ try
+ {
+ value = m_Parser->Eval();
+ }
+ catch(itk::ExceptionObject& err)
+ {
+ itkExceptionMacro(<< err);
+ }
+
+ return static_cast<bool> (value);
+ }
+
+
+ void SetExpression(const std::string expression)
+ {
+ m_Expression=expression;
+ m_Parser->SetExpr(m_Expression);
+ }
+
+ /** Return the expression to be parsed */
+ std::string GetExpression() const
+ {
+ return m_Expression;
+ }
+
+ void SetNumberOfBands(unsigned int NbOfBands)
+ {
+
+ m_NbOfBands=NbOfBands;
+ std::ostringstream varName;
+
+
+ m_AImageP1.resize(NbOfBands,0.0);
+ m_AImageP2.resize(NbOfBands,0.0);
+
+ for(unsigned int i=0;i<NbOfBands;i++)
+ {
+ varName << "p1b" << i+1 ;
+ m_Parser->DefineVar(varName.str(), &(m_AImageP1[i]));
+ varName.str("");
+ varName << "p2b" << i+1 ;
+ m_Parser->DefineVar(varName.str(), &(m_AImageP2[i]));
+ varName.str("");
+ }
+ // customized data
+ //m_NbVar++;
+ //this->SetDataSize(m_NbVar);
+ m_Parser->DefineVar("distance", &m_Distance);
+ m_Parser->DefineVar("spectralAngle", &m_SpectralAngle);
+ //this->SetVarName(m_NbVar-1,"spectralDistance");
+
+ }
+
+ ConnectedComponentMuParserFunctor()
+ {
+ m_Parser = ParserType::New();
+ m_NbOfBands =0;
+ };
+
+ ~ConnectedComponentMuParserFunctor()
+ {};
+
+ protected:
+
+
+ private:
+
+ ConnectedComponentMuParserFunctor(const Self &);//purposely not implemented
+ void operator =(const Self &); //purposely not implemented
+
+ std::string m_Expression;
+ ParserType::Pointer m_Parser;
+ std::vector<double> m_AImageP1;
+ std::vector<double> m_AImageP2;
+ double m_Distance;
+ double m_SpectralAngle;
+ std::vector<std::string > m_VarName;
+ unsigned int m_NbOfBands;
+ double m_ParserResult;
+
+ };
+ } // end of Functor namespace
+
+
+}//end namespace otb
+
+
+#endif
diff --git a/Code/BasicFilters/otbMaskMuParserFilter.h b/Code/BasicFilters/otbMaskMuParserFilter.h
new file mode 100644
index 0000000000..26d5044362
--- /dev/null
+++ b/Code/BasicFilters/otbMaskMuParserFilter.h
@@ -0,0 +1,142 @@
+/*=========================================================================
+
+ 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 __otbMaskMuParserFilter_h
+#define __otbMaskMuParserFilter_h
+
+
+#include "itkProgressReporter.h"
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "itkArray.h"
+
+#include "otbParser.h"
+#include "otbMacro.h"
+#include "itkMacro.h"
+#include "otbMaskMuParserFunctor.h"
+
+#include "itkImageToImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+
+namespace otb
+{
+/** \class MaskMuParserImageFilter
+ * \brief Performs a mathematical operation on the input images
+ * according to the formula specified by the user. values different from 0 are set to 1
+ *
+ * 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
+ *
+ * an input vector image and a Mu Parser compliant fomula is needed
+ * each band of vector image is name im1bX, where X is the band index
+ * for example im1b2 correspond to the second band of the input image.
+ * 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 :
+ * "im1b1<140"
+ *
+ *
+ * \sa Parser
+ *
+ * \ingroup Streamed
+ * \ingroup Threaded
+ */
+
+
+
+template< class TInputImage,class TOutputImage,class TFunction=Functor::MaskMuParserFunctor<
+ typename TInputImage::PixelType > >
+class ITK_EXPORT MaskMuParserFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef MaskMuParserFilter Self;
+ typedef itk::ImageToImageFilter< TInputImage, TOutputImage> 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(MaskMuParserFilter, itk::ImageToImageFilter);
+
+ /** Some convenient typedefs. */
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename InputImageType::PixelType PixelType;
+ typedef typename InputImageType::Pointer InputImagePointer;
+ typedef typename InputImageType::ConstPointer InputImageConstPointer;
+ typedef TOutputImage OutputImageType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+ typedef typename OutputImageType::Pointer OutputImagePointer;
+
+ typedef TFunction FunctorType;
+ typedef typename FunctorType::Pointer FunctorPointer;
+
+ typedef MaskMuParserFilter<InputImageType,OutputImageType,FunctorType> MaskMuParserFilterType;
+
+
+ /** Set the expression to be parsed */
+ void SetExpression(const std::string expression);
+
+ /** Return the expression to be parsed */
+ std::string GetExpression() const;
+
+protected :
+ MaskMuParserFilter();
+ virtual ~MaskMuParserFilter();
+ virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ void BeforeThreadedGenerateData();
+ void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread,int threadId);
+ void AfterThreadedGenerateData();
+
+
+private :
+ MaskMuParserFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ std::vector<FunctorPointer> m_VFunctor;
+ std::string m_Expression;
+ long m_UnderflowCount;
+ long m_OverflowCount;
+ itk::Array<long> m_ThreadUnderflow;
+ itk::Array<long> m_ThreadOverflow;
+};
+
+}//end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMaskMuParserFilter.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbMaskMuParserFilter.txx b/Code/BasicFilters/otbMaskMuParserFilter.txx
new file mode 100644
index 0000000000..1d4d01d72c
--- /dev/null
+++ b/Code/BasicFilters/otbMaskMuParserFilter.txx
@@ -0,0 +1,169 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.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 __otbMaskMuParserFilter_txx
+#define __otbMaskMuParserFilter_txx
+
+#include "otbMaskMuParserFilter.h"
+#include <iostream>
+#include <string>
+
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkNumericTraits.h"
+#include "itkProgressReporter.h"
+
+namespace otb
+{
+
+// constructor
+template < class TInputImage,class TOutputImage,class TFunction>
+MaskMuParserFilter<TInputImage,TOutputImage,TFunction>::MaskMuParserFilter()
+{
+ m_UnderflowCount = 0;
+ m_OverflowCount = 0;
+ m_ThreadUnderflow.SetSize(1);
+ m_ThreadOverflow.SetSize(1);
+}
+
+// Destructor
+template < class TInputImage,class TOutputImage,class TFunction>
+MaskMuParserFilter<TInputImage,TOutputImage,TFunction>::~MaskMuParserFilter()
+{
+
+}
+
+template < class TInputImage,class TOutputImage,class TFunction>
+void MaskMuParserFilter<TInputImage,TOutputImage,TFunction>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+
+ os << indent << "Expression: " << m_Expression << std::endl;
+}
+
+
+
+template< class TInputImage,class TOutputImage,class TFunction>
+void MaskMuParserFilter<TInputImage,TOutputImage,TFunction>
+::SetExpression(const std::string expression)
+ {
+ if (m_Expression != expression)
+ m_Expression = expression;
+ this->Modified();
+ }
+
+template< class TInputImage,class TOutputImage,class TFunction>
+std::string MaskMuParserFilter<TInputImage,TOutputImage,TFunction>
+::GetExpression() const
+ {
+ return m_Expression;
+ }
+
+
+/**
+ * BeforeThreadedGenerateData
+ */
+template < class TInputImage,class TOutputImage,class TFunction>
+void MaskMuParserFilter<TInputImage,TOutputImage,TFunction>
+::BeforeThreadedGenerateData()
+ {
+
+ typename std::vector<FunctorPointer>::iterator itFunctor;
+ unsigned int nbThreads = this->GetNumberOfThreads();
+ unsigned int nbInputImages = this->GetNumberOfInputs();
+ unsigned int nbOfBands =this->GetInput()->GetNumberOfComponentsPerPixel();
+ unsigned int thread_index;
+ std::ostringstream varName;
+
+ // Allocate and initialize the thread temporaries
+ m_ThreadUnderflow.SetSize(nbThreads);
+ m_ThreadUnderflow.Fill(0);
+ m_ThreadOverflow.SetSize(nbThreads);
+ m_ThreadOverflow.Fill(0);
+ m_VFunctor.resize(nbThreads);
+
+
+ for(itFunctor = m_VFunctor.begin(); itFunctor < m_VFunctor.end(); itFunctor++)
+ {
+ *itFunctor = FunctorType::New();
+ }
+
+ for(thread_index = 0; thread_index < nbThreads; thread_index++)
+ {
+
+ m_VFunctor.at(thread_index)->SetExpression(m_Expression);
+ m_VFunctor.at(thread_index)->SetNumberOfBands(nbOfBands);
+ // we associate image band with its name im1b1 for band 1
+ for(unsigned int i=0;i<nbOfBands;i++)
+ {
+ varName << "b" << i+1; // +1 var[0]=im1b1
+ m_VFunctor.at(thread_index)->SetVarName(i,varName.str());
+ varName.str("");
+ }
+
+ if(!m_VFunctor.at(thread_index)->CheckExpression())
+ {
+ std::cerr<<"parser can't evaluate formula "<<this->GetExpression()<<"in thread "<<thread_index<<std::endl;
+ }
+ }
+
+ }
+
+template < class TInputImage,class TOutputImage,class TFunction>
+void MaskMuParserFilter<TInputImage,TOutputImage,TFunction>
+::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread,int threadId)
+ {
+
+
+ InputImageConstPointer inputPtr = this->GetInput();
+ OutputImagePointer outputPtr = this->GetOutput(0);
+
+ // Define the portion of the input to walk for this thread, using
+ // the CallCopyOutputRegionToInputRegion method allows for the input
+ // and output images to be different dimensions
+ InputImageRegionType inputRegionForThread;
+ this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
+
+ // Define the iterators
+ itk::ImageRegionConstIterator<TInputImage> inputIt(inputPtr, inputRegionForThread);
+ itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
+
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ inputIt.GoToBegin();
+ outputIt.GoToBegin();
+
+ FunctorPointer functorP = m_VFunctor.at(threadId);
+ FunctorType& functor = *functorP;
+ while( !inputIt.IsAtEnd() )
+ {
+ outputIt.Set( functor( inputIt.Get() ) );
+ ++inputIt;
+ ++outputIt;
+ progress.CompletedPixel(); // potential exception thrown here
+ }
+ }
+
+template < class TInputImage,class TOutputImage,class TFunction>
+void MaskMuParserFilter<TInputImage,TOutputImage,TFunction>
+::AfterThreadedGenerateData()
+ {
+ }
+
+} // end namespace otb
+
+#endif
diff --git a/Code/BasicFilters/otbMaskMuParserFunctor.h b/Code/BasicFilters/otbMaskMuParserFunctor.h
new file mode 100644
index 0000000000..63632da2d1
--- /dev/null
+++ b/Code/BasicFilters/otbMaskMuParserFunctor.h
@@ -0,0 +1,185 @@
+/*=========================================================================
+
+ 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 __otbMaskMuParserFunctor_h
+#define __otbMaskMuParserFunctor_h
+
+#include "otbParser.h"
+#include "otbMacro.h"
+#include "itkMacro.h"
+
+
+#include "itkLightObject.h"
+#include "itkObjectFactory.h"
+
+
+namespace otb
+{
+/** \class MaskMuParserFunctor
+ * \brief functor used to create binary mask for input of connected component segmentation module.
+ *
+ * This functor 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
+ *
+ * \sa Parser
+ *
+ */
+
+/** \class MaskMuParserFunctor
+ * \brief This functor use MuParser to evaluate and process mathematical expression
+ */
+namespace Functor
+{
+
+template<class TInput>
+class ITK_EXPORT MaskMuParserFunctor : public itk::LightObject
+{
+public:
+typedef Parser ParserType;
+typedef MaskMuParserFunctor Self;
+typedef itk::LightObject 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(MaskMuParserFunctor,itk::LightObject);
+
+inline bool operator()(const TInput &p)
+{
+
+ double value;
+
+ if(p.Size() !=m_NbVar)
+ {
+ itkExceptionMacro(<< "wrong number of bands" << std::endl);
+ }
+
+ // we fill the buffer
+ for(unsigned int i=0; i<m_NbVar; i++)
+ {
+
+ m_AImage[i]= static_cast<double> (p[i]);
+ }
+
+ // cast
+ try
+ {
+ value = m_Parser->Eval();
+ }
+ catch(itk::ExceptionObject& err)
+ {
+ itkExceptionMacro(<< err);
+ }
+ return static_cast<bool> (value);
+
+}
+
+void SetExpression(const std::string expression)
+{
+ m_Expression=expression;
+ m_Parser->SetExpr(m_Expression);
+}
+
+/** Return the expression to be parsed */
+std::string GetExpression() const
+{
+ return m_Expression;
+}
+
+void SetVarName(unsigned int idx, const std::string varName)
+{
+ m_VarName[idx] = varName;
+ m_Parser->DefineVar(m_VarName[idx], &(m_AImage[idx]));
+}
+
+std::string GetVarName(unsigned int idx)
+{
+ return m_VarName[idx];
+}
+
+void SetNumberOfBands(unsigned int NumberOfBands)
+{
+ m_NbVar=NumberOfBands;
+ m_VarName.resize(NumberOfBands,"");
+ m_AImage.resize(NumberOfBands,0.0);
+}
+
+/** Check the expression */
+bool CheckExpression()
+{
+
+ for(unsigned int i=0; i<m_NbVar; i++)
+ {
+ m_AImage[i]=0;
+ }
+
+ try
+ {
+ m_Parser->Eval();
+
+ }
+ catch(itk::ExceptionObject& err)
+ {
+ return false;
+ }
+
+ return true;
+}
+
+protected:
+MaskMuParserFunctor()
+{
+m_Parser = ParserType::New();
+};
+
+~MaskMuParserFunctor()
+{};
+
+
+private:
+
+MaskMuParserFunctor(const Self &);//purposely not implemented
+void operator =(const Self &); //purposely not implemented
+
+std::string m_Expression;
+ParserType::Pointer m_Parser;
+std::vector<double> m_AImage;
+std::vector<std::string > m_VarName;
+unsigned int m_NbVar;
+double m_ParserResult;
+
+};
+} // end of Functor namespace
+
+
+}//end namespace otb
+
+
+#endif
--
GitLab