ENH: bandmathx, added new operators

Modules/Filtering/MathParserX/include/otbParserXPlugins.h

@@ -385,6 +385,27 @@ public:
   };
 
 
+class vect2scal : public mup::ICallback
+  {
+public:
+    vect2scal():ICallback(mup::cmFUNC, "vect2scal", 1)
+    {}
+
+    virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int itkNotUsed(a_iArgc));
+
+    const mup::char_type* GetDesc() const
+    {
+      return "vect2scal - Convert one dimensional vector to scalar";
+    }
+
+    mup::IToken* Clone() const
+    {
+      return new vect2scal(*this);
+    }
+  };
+  
+  
+
 class vcos : public mup::ICallback
   {
 public:
@@ -405,6 +426,25 @@ public:
   };
 
 
+class vacos : public mup::ICallback
+  {
+public:
+    vacos():ICallback(mup::cmFUNC, "vacos", 1)
+    {}
+
+    virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int itkNotUsed(a_iArgc));
+
+    const mup::char_type* GetDesc() const
+    {
+      return "vacos - Arccosinus for noncomplex vectors & matrices";
+    }
+
+    mup::IToken* Clone() const
+    {
+      return new vacos(*this);
+    }
+  };
+
 class vsin : public mup::ICallback
   {
 public:
@@ -424,6 +464,25 @@ public:
     }
   };
 
+class vasin : public mup::ICallback
+  {
+public:
+    vasin():ICallback(mup::cmFUNC, "vasin", 1)
+    {}
+
+    virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int itkNotUsed(a_iArgc));
+
+    const mup::char_type* GetDesc() const
+    {
+      return "vasin - Arcsinus for noncomplex vectors & matrices";
+    }
+
+    mup::IToken* Clone() const
+    {
+      return new vasin(*this);
+    }
+  };
+
 
 class vtan : public mup::ICallback
   {
@@ -445,6 +504,26 @@ public:
   };
 
 
+class vatan : public mup::ICallback
+  {
+public:
+    vatan():ICallback(mup::cmFUNC, "vatan", 1)
+    {}
+
+    virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int itkNotUsed(a_iArgc));
+
+    const mup::char_type* GetDesc() const
+    {
+      return "vatan - Arctangent for noncomplex vectors & matrices";
+    }
+
+    mup::IToken* Clone() const
+    {
+      return new vatan(*this);
+    }
+  };
+  
+
 class vtanh : public mup::ICallback
   {
 public:

Modules/Filtering/MathParserX/src/otbParserX.cxx

@@ -77,7 +77,11 @@ public:
     m_MuParserX.DefineFun(new vnorm);
     m_MuParserX.DefineFun(new vmin);
     m_MuParserX.DefineFun(new vmax);
+    m_MuParserX.DefineFun(new vect2scal);
     m_MuParserX.DefineFun(new vcos);
+    m_MuParserX.DefineFun(new vacos);
+    m_MuParserX.DefineFun(new vasin);
+    m_MuParserX.DefineFun(new vatan);
     m_MuParserX.DefineFun(new vsin);
     m_MuParserX.DefineFun(new vtan);
     m_MuParserX.DefineFun(new vtanh);

Modules/Filtering/MathParserX/src/otbParserXPlugins.cxx

@@ -734,6 +734,33 @@ void vmax::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
     }
 
 
+void vect2scal::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+    {
+      if(a_iArgc != 1)
+        return;
+
+      assert(a_pArg[0]->GetType()=='m');
+
+      // Get the argument from the argument input vector
+      const mup::matrix_type a = a_pArg[0]->GetArray();
+
+
+      int nbrows = a.GetRows();
+      int nbcols = a.GetCols();
+      
+      assert(nbrows == 1);
+      assert(nbcols == 1);
+
+      mup::float_type res;
+
+      res = a.At(0,0).GetFloat();
+          
+
+      // The return value is passed by writing it to the reference ret
+      *ret = res;
+    }
+
+
 void vcos::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
     {
       if(a_iArgc != 1)
@@ -758,6 +785,32 @@ void vcos::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
       // The return value is passed by writing it to the reference ret
       *ret = res;
     }
+    
+    
+void vacos::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+    {
+      if(a_iArgc != 1)
+        return;
+
+      assert(a_pArg[0]->GetType()=='m');
+
+      // Get the argument from the argument input vector
+      const mup::matrix_type a = a_pArg[0]->GetArray();
+
+
+      int nbrows = a.GetRows();
+      int nbcols = a.GetCols();
+
+      mup::matrix_type res(nbrows,nbcols,0.);
+
+      for (int k=0; k<nbcols; ++k)
+        for (int p=0; p<nbrows; ++p)
+          res.At(p,k) = vcl_acos(a.At(p,k).GetFloat());
+          
+
+      // The return value is passed by writing it to the reference ret
+      *ret = res;
+    }
 
 
 void vsin::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
@@ -784,6 +837,32 @@ void vsin::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
     }
 
 
+void vasin::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+    {
+      if(a_iArgc != 1)
+        return;
+
+      assert(a_pArg[0]->GetType()=='m');
+
+      // Get the argument from the argument input vector
+      const mup::matrix_type a = a_pArg[0]->GetArray();
+
+
+      int nbrows = a.GetRows();
+      int nbcols = a.GetCols();
+
+      mup::matrix_type res(nbrows,nbcols,0.);
+
+      for (int k=0; k<nbcols; ++k)
+        for (int p=0; p<nbrows; ++p)
+          res.At(p,k) = vcl_asin(a.At(p,k).GetFloat());
+          
+
+      // The return value is passed by writing it to the reference ret
+      *ret = res;
+    }
+
+
 void vtan::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
     {
       if(a_iArgc != 1)
@@ -808,6 +887,30 @@ void vtan::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
     }
 
 
+void vatan::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+    {
+      if(a_iArgc != 1)
+        return;
+
+      assert(a_pArg[0]->GetType()=='m');
+
+      // Get the argument from the argument input vector
+      const mup::matrix_type a = a_pArg[0]->GetArray();
+
+      int nbrows = a.GetRows();
+      int nbcols = a.GetCols();
+
+      mup::matrix_type res(nbrows,nbcols,0.);
+
+      for (int k=0; k<nbcols; ++k)
+        for (int p=0; p<nbrows; ++p)
+          res.At(p,k) = vcl_atan(a.At(p,k).GetFloat());
+
+      // The return value is passed by writing it to the reference ret
+      *ret = res;
+    }
+
+
 void vtanh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
     {
       if(a_iArgc != 1)

Modules/Filtering/MathParserX/test/otbBandMathXImageFilter.cxx

@@ -299,7 +299,7 @@ int otbBandMathXImageFilterConv( int itkNotUsed(argc), char* argv [])
   //filter->SetConstant("expo",expo);
   //filter->SetExpression("conv(kernel1,imageAb1N3x5,imageAb2N3x5); im2b1^1.1; vcos(canal3); mean(imageAb2N3x3); var(imageAb2N3x3); median(imageAb2N3x3)");
   filter->ImportContext(inputFilename); //Equivalent to three commands above
-  filter->SetExpression("(vmax(canal3b1N3x5)+vmin(canal3b1N3x5)) div {2.0} + {imageAb3Var} dv 2.0 + {imageAb2Mini / im2b1Maxi} mlt 3.4 + {imageAb3Mean / imageAb1Sum * imageAb3Var} pw 1.2");
+  filter->SetExpression("(vmax(canal3b1N3x5)+vmin(canal3b1N3x5)) div {2.0} + {imageAb3Var} dv 2.0 + {imageAb2Mini / im2b1Maxi} mlt 3.4 + {imageAb3Mean / imageAb1Sum * imageAb3Var} pw 1.2 ; vect2scal(vacos({0.5}) + vasin({0.5}) + vatan({0.5})) > 2.0 ?{1}:{0}");
   filter->Update();
 
   if (filter->GetNumberOfOutputs() != 2)
@@ -312,7 +312,7 @@ int otbBandMathXImageFilterConv( int itkNotUsed(argc), char* argv [])
   if (output1->GetNumberOfComponentsPerPixel() != 7)
     itkGenericExceptionMacro(<< "Wrong number of components per pixel (input 1).");
 
-  if (output2->GetNumberOfComponentsPerPixel() != 1)
+  if (output2->GetNumberOfComponentsPerPixel() != 2)
     itkGenericExceptionMacro(<< "Wrong number of components per pixel (input 2).");
 
   std::cout << "\n---  Standard Use\n";
@@ -370,6 +370,10 @@ int otbBandMathXImageFilterConv( int itkNotUsed(argc), char* argv [])
       vect2.push_back(it3.GetPixel(i)[0]);
     std::sort(vect2.begin(),vect2.end());
     px2[0] = (vect2.back() + vect2.front())/2.0 + imageAb3Var / 2.0 + (imageAb2Mini / im2b1Maxi)*3.4 + vcl_pow(imageAb3Mean / imageAb1Sum * imageAb3Var,1.2);
+    if ( vcl_acos(0.5)+vcl_asin(0.5)+vcl_atan(0.5) > 2.0 )
+		px2[1]=1.0;
+	else
+		px2[1]=0.0;
 
 
 
@@ -390,10 +394,13 @@ int otbBandMathXImageFilterConv( int itkNotUsed(argc), char* argv [])
 
     //expression 2
     double result8 = itoutput2.GetCenterPixel()[0];
+    double result9 = itoutput2.GetCenterPixel()[1];
     double expected8 = px2[0];
+    double expected9 = px2[1];
     double error8 = (result8 - expected8) * (result8 - expected8) / (result8 + expected8);
+    double error9 = (result9 - expected9) * (result9 - expected9) / (result9 + expected9);
 
-    if ( ( error1 > 1E-9 ) || ( error2 > 1E-9 ) || ( error3 > 1E-9 ) || ( error4 > 1E-9 ) || ( error5 > 1E-9 ) || ( error6 > 1E-9 ) || ( error7 > 1E-9 ) || (error8 > 1E-9) )
+    if ( ( error1 > 1E-9 ) || ( error2 > 1E-9 ) || ( error3 > 1E-9 ) || ( error4 > 1E-9 ) || ( error5 > 1E-9 ) || ( error6 > 1E-9 ) || ( error7 > 1E-9 ) || (error8 > 1E-9) || (error9 > 1E-9))
       {
       itkGenericExceptionMacro( << "TEST FAILLED" << std::endl
          << "Error1 = " << error1  << std::endl
@@ -419,7 +426,10 @@ int otbBandMathXImageFilterConv( int itkNotUsed(argc), char* argv [])
          << "     Expected7 =  " << expected7 << std::endl
          << "Error8 = " << error8 << std::endl
          << "     Result8 =  "   << result8
-         << "     Expected8 =  " << expected8 << std::endl);
+         << "     Expected8 =  " << expected8 << std::endl
+         << "Error9 = " << error9 << std::endl
+         << "     Result9 =  "   << result9
+         << "     Expected9 =  " << expected9 << std::endl);
       }
   }