Commit 3c6fc044 authored by Christophe Palmann's avatar Christophe Palmann

ENH: bandmathx, uniformization of operators, also added new ones

parent e0e7e11d
......@@ -61,15 +61,19 @@ public:
m_MuParserX.DefineFun(new conv);
m_MuParserX.DefineFun(new bands);
m_MuParserX.DefineFun(new cat);
m_MuParserX.DefineOprt(new ElementWiseDivision);
m_MuParserX.DefineOprt(new ElementWiseMultiplication);
m_MuParserX.DefineOprt(new ElementWisePower);
m_MuParserX.DefineFun(new mean);
m_MuParserX.DefineFun(new var);
m_MuParserX.DefineFun(new median);
m_MuParserX.DefineOprt(new ElementWiseDivision);
m_MuParserX.DefineOprt(new ElementWiseMultiplication);
m_MuParserX.DefineOprt(new ElementWisePower);
m_MuParserX.DefineOprt(new DivisionByScalar);
m_MuParserX.DefineOprt(new MultiplicationByScalar);
m_MuParserX.DefineOprt(new PowerByScalar);
m_MuParserX.DefineFun(new vmin);
m_MuParserX.DefineFun(new vmax);
m_MuParserX.DefineFun(new vcos);
m_MuParserX.DefineFun(new vsin);
m_MuParserX.DefineFun(new vtan);
......
......@@ -38,9 +38,12 @@ void bands::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_
mup::matrix_type b = a_pArg[1]->GetArray();
int nbrows = b.GetRows(); assert(nbrows==1); // Bands selection is done by a row vector
int nbrows = b.GetRows();
int nbcols = b.GetCols();
assert(a.GetRows()==1); // Bands selection is done on a row vector
assert(nbrows==1); // Bands selection is done by a row vector
mup::matrix_type res(1,nbcols,0.);
for (int k=0; k<nbcols; ++k)
......@@ -70,17 +73,11 @@ void conv::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
float sum=0.0;
assert(a_pArg[k]->GetType()=='m');
mup::matrix_type m2 = a_pArg[k]->GetArray();
if ( (m2.GetRows() != nbrows) || (m2.GetCols() != nbcols) )
{
mup::ErrorContext err;
err.Errc = mup::ecMATRIX_DIMENSION_MISMATCH;
err.Arg = a_iArgc;
err.Ident = GetIdent();
throw mup::ParserError(err);
}
assert(m2.GetRows() == nbrows);
assert(m2.GetCols() == nbcols);
for (int i=0; i<nbrows; i++)
for (int j=0; j<nbcols; j++)
......@@ -98,6 +95,7 @@ void conv::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void ElementWiseDivision::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
assert(a_pArg[0]->GetType()=='m');
assert(a_pArg[1]->GetType()=='m');
......@@ -111,20 +109,15 @@ void ElementWiseDivision::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *
int nbrows2 = b.GetRows();
int nbcols2 = b.GetCols();
if ( (nbrows != nbrows2) || (nbcols != nbcols2) )
{
mup::ErrorContext err;
err.Errc = mup::ecMATRIX_DIMENSION_MISMATCH;
err.Ident = GetIdent();
throw mup::ParserError(err);
}
assert(nbrows == nbrows2);
assert(nbcols == nbcols2);
mup::matrix_type res(nbrows,nbcols,0.);
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
{
assert(vcl_abs(b.At(p,k).GetFloat()) > 1e-6);
assert(vcl_abs(b.At(p,k).GetFloat()) > 1e-9);
res.At(p,k) = a.At(p,k).GetFloat() / b.At(p,k).GetFloat();
}
......@@ -134,6 +127,49 @@ void ElementWiseDivision::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *
}
void DivisionByScalar::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
assert(a_pArg[0]->GetType()=='m');
const mup::matrix_type a = a_pArg[0]->GetArray();
mup::matrix_type b;
double scalar;
switch (a_pArg[1]->GetType())
{
case 'i':
scalar = (double) a_pArg[1]->GetInteger();
break;
case 'f':
scalar = a_pArg[1]->GetFloat();
break;
case 'm':
b = a_pArg[1]->GetArray();
assert(b.GetRows() == 1);
assert(b.GetCols() == 1);
scalar = b.At(0,0);
break;
}
assert(vcl_abs(scalar) > 1e-9);
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) = a.At(p,k).GetFloat() / scalar;
// The return value is passed by writing it to the reference ret
*ret = res;
}
void ElementWiseMultiplication::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
......@@ -149,13 +185,8 @@ void ElementWiseMultiplication::Eval(mup::ptr_val_type &ret, const mup::ptr_val_
int nbrows2 = b.GetRows();
int nbcols2 = b.GetCols();
if ( (nbrows != nbrows2) || (nbcols != nbcols2) )
{
mup::ErrorContext err;
err.Errc = mup::ecMATRIX_DIMENSION_MISMATCH;
err.Ident = GetIdent();
throw mup::ParserError(err);
}
assert(nbrows == nbrows2);
assert(nbcols == nbcols2);
mup::matrix_type res(nbrows,nbcols,0.);
......@@ -167,27 +198,73 @@ void ElementWiseMultiplication::Eval(mup::ptr_val_type &ret, const mup::ptr_val_
*ret = res;
}
void MultiplicationByScalar::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
assert(a_pArg[0]->GetType()=='m');
const mup::matrix_type a = a_pArg[0]->GetArray();
mup::matrix_type b;
double scalar;
switch (a_pArg[1]->GetType())
{
case 'i':
scalar = (double) a_pArg[1]->GetInteger();
break;
case 'f':
scalar = a_pArg[1]->GetFloat();
break;
case 'm':
b = a_pArg[1]->GetArray();
assert(b.GetRows() == 1);
assert(b.GetCols() == 1);
scalar = b.At(0,0);
break;
}
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) = a.At(p,k).GetFloat() * scalar;
// The return value is passed by writing it to the reference ret
*ret = res;
}
void ElementWisePower::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==2);
assert(a_pArg[0]->GetType()=='m');
assert(a_pArg[1]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
const double pw = a_pArg[1]->GetFloat();
const mup::matrix_type b = a_pArg[1]->GetArray();
int nbrows = a.GetRows();
int nbcols = a.GetCols();
int nbrows2 = b.GetRows();
int nbcols2 = b.GetCols();
assert(nbrows == nbrows2);
assert(nbcols == nbcols2);
mup::matrix_type res(nbrows,nbcols,0.);
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
{
assert(a.At(p,k).GetFloat() >= 0 );
res.At(p,k) = vcl_pow(a.At(p,k).GetFloat(),pw);
res.At(p,k) = vcl_pow(a.At(p,k).GetFloat(),b.At(p,k).GetFloat());
}
// The return value is passed by writing it to the reference ret
......@@ -195,8 +272,52 @@ void ElementWisePower::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_p
}
void PowerByScalar::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
assert(a_pArg[0]->GetType()=='m');
const mup::matrix_type a = a_pArg[0]->GetArray();
mup::matrix_type b;
double scalar;
switch (a_pArg[1]->GetType())
{
case 'i':
scalar = (double) a_pArg[1]->GetInteger();
break;
case 'f':
scalar = a_pArg[1]->GetFloat();
break;
case 'm':
b = a_pArg[1]->GetArray();
assert(b.GetRows() == 1);
assert(b.GetCols() == 1);
scalar = b.At(0,0);
break;
}
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_pow(a.At(p,k).GetFloat() , scalar);
// The return value is passed by writing it to the reference ret
*ret = res;
}
void ndvi::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==2);
// Get the argument from the argument input vector
mup::float_type r = a_pArg[0]->GetFloat();
mup::float_type niri = a_pArg[1]->GetFloat();
......@@ -536,6 +657,9 @@ void vmax::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
//--------------------------------------------------------------------------------------------------------//
void vcos::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -557,10 +681,12 @@ void vcos::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vsin::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
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();
......@@ -577,10 +703,13 @@ void vsin::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vtan::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
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();
......@@ -597,6 +726,9 @@ void vtan::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vtanh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -617,6 +749,9 @@ void vtanh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_
void vsinh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -637,6 +772,10 @@ void vsinh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_
void vcosh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -657,10 +796,12 @@ void vcosh::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_
void vlog::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
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();
......@@ -677,6 +818,9 @@ void vlog::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vlog10::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -697,10 +841,12 @@ void vlog10::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a
void vabs::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
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();
......@@ -717,6 +863,9 @@ void vabs::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vexp::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
assert(a_pArg[0]->GetType()=='m');
// Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
......@@ -737,10 +886,13 @@ void vexp::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_i
void vsqrt::Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
assert(a_iArgc==1);
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();
......
......@@ -87,6 +87,26 @@ class ElementWiseDivision : public mup::IOprtBin
};
class DivisionByScalar : public mup::IOprtBin
{
public:
DivisionByScalar():IOprtBin(_T("dv"), (int)(mup::prMUL_DIV), mup::oaLEFT)
{}
virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int);
virtual const mup::char_type* GetDesc() const
{
return _T("x dv y - division of vectors / matrices by a scalar");
}
virtual mup::IToken* Clone() const
{
return new DivisionByScalar(*this);
}
};
class ElementWiseMultiplication : public mup::IOprtBin
{
public:
......@@ -107,6 +127,26 @@ class ElementWiseMultiplication : public mup::IOprtBin
};
class MultiplicationByScalar : public mup::IOprtBin
{
public:
MultiplicationByScalar():IOprtBin(_T("mlt"), (int)(mup::prMUL_DIV), mup::oaLEFT)
{}
virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int);
virtual const mup::char_type* GetDesc() const
{
return _T("x mlt y - multiplication of vectors / matrices by a scalar");
}
virtual mup::IToken* Clone() const
{
return new MultiplicationByScalar(*this);
}
};
class ElementWisePower : public mup::IOprtBin
{
public:
......@@ -127,6 +167,26 @@ public:
};
class PowerByScalar : public mup::IOprtBin
{
public:
PowerByScalar():IOprtBin(_T("pw"), (int)(mup::prMUL_DIV), mup::oaLEFT)
{}
virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int);
virtual const mup::char_type* GetDesc() const
{
return _T("x pw y - power of vectors / matrices by a scalar");
}
virtual mup::IToken* Clone() const
{
return new PowerByScalar(*this);
}
};
class ndvi : public mup::ICallback
{
public:
......
......@@ -302,7 +302,7 @@ int otbBandMathImageFilterXConv( 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} + {imageAb2Mini / im2b1Maxi} + {imageAb3Mean / imageAb1Sum * imageAb3Var} ");
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->Update();
if (filter->GetNumberOfOutputs() != 2)
......@@ -376,7 +376,7 @@ int otbBandMathImageFilterXConv( int itkNotUsed(argc), char* argv [])
for (int i=0; i<it3.Size(); i++)
vect2.push_back(it3.GetPixel(i)[0]);
std::sort(vect2.begin(),vect2.end());
px2[0] = (vect2.back() + vect2.front())/2.0 + (imageAb2Mini / im2b1Maxi) + imageAb3Mean / imageAb1Sum * imageAb3Var;
px2[0] = (vect2.back() + vect2.front())/2.0 + imageAb3Var / 2.0 + (imageAb2Mini / im2b1Maxi)*3.4 + vcl_pow(imageAb3Mean / imageAb1Sum * imageAb3Var,1.2);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment