From 7e74d95583d5612f17071b38d22444d89a652d40 Mon Sep 17 00:00:00 2001
From: Christophe Palmann <christophe.palmann@c-s.fr>
Date: Wed, 29 Oct 2014 11:34:05 +0100
Subject: [PATCH] WIP : otbBandMathX --> neighborhood access
---
Code/BasicFilters/otbBandMathImageFilterX.h | 7 +-
Code/BasicFilters/otbBandMathImageFilterX.txx | 224 ++++++++++++++----
Code/Common/otbParserX.cxx | 2 +-
Code/Common/otbParserXPlugins.h | 208 ++++++++++------
Testing/Code/BasicFilters/CMakeLists.txt | 16 +-
.../BasicFilters/otbBandMathImageFilterX.cxx | 148 ++++++++++++
.../BasicFilters/otbBasicFiltersTests13.cxx | 3 +-
7 files changed, 476 insertions(+), 132 deletions(-)
diff --git a/Code/BasicFilters/otbBandMathImageFilterX.h b/Code/BasicFilters/otbBandMathImageFilterX.h
index cab9e032ff..92feeb57bb 100644
--- a/Code/BasicFilters/otbBandMathImageFilterX.h
+++ b/Code/BasicFilters/otbBandMathImageFilterX.h
@@ -112,6 +112,9 @@ public:
/** Set the expression to be parsed */
void SetExpression(const std::string& expression);
+ /** Set a matrix (or a vector) */
+ void SetMatrix(const std::string& name, const std::string& definition);
+
/** Return the expression to be parsed */
std::string GetExpression(int) const;
@@ -158,7 +161,9 @@ private :
std::vector<ParserType::Pointer> m_VParser;
std::vector< std::vector<adhocStruct2> > m_AImage;
std::vector< adhocStruct > m_VVarName;
- std::vector< adhocStruct > m_VAllowedVarName;
+ std::vector< adhocStruct > m_VAllowedVarNameAuto;
+ std::vector< adhocStruct > m_VAllowedVarNameAddedByUser;
+ std::vector< adhocStruct > m_VFinalAllowedVarName; // m_VAllowedVarNameAuto + m_VAllowedVarNameAddedByUser
std::vector< adhocStruct > m_VNotAllowedVarName;
unsigned int m_NbVar;
std::vector< int > m_outputsDimensions;
diff --git a/Code/BasicFilters/otbBandMathImageFilterX.txx b/Code/BasicFilters/otbBandMathImageFilterX.txx
index 63a2df8505..a1b7ea6ffe 100644
--- a/Code/BasicFilters/otbBandMathImageFilterX.txx
+++ b/Code/BasicFilters/otbBandMathImageFilterX.txx
@@ -24,6 +24,7 @@
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
+#include "itkConstNeighborhoodIterator.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "otbMacro.h"
@@ -55,12 +56,12 @@ BandMathImageFilterX<TImage>
adhocStruct ahcX;
ahcX.name = "idxX";
ahcX.type = 0;
- m_VAllowedVarName.push_back(ahcX);
+ m_VAllowedVarNameAuto.push_back(ahcX);
adhocStruct ahcY;
ahcY.name = "idxY";
ahcY.type = 1;
- m_VAllowedVarName.push_back(ahcY);
+ m_VAllowedVarNameAuto.push_back(ahcY);
//this->SetNumberOfThreads(1);
@@ -113,7 +114,7 @@ void BandMathImageFilterX<TImage>
ahcPhyX.name = sstmPhyX.str();
ahcPhyX.type = 2;
ahcPhyX.info[0] = idx; // Input image #ID
- m_VAllowedVarName.push_back(ahcPhyX);
+ m_VAllowedVarNameAuto.push_back(ahcPhyX);
std::stringstream sstmPhyY;
adhocStruct ahcPhyY;
@@ -121,7 +122,7 @@ void BandMathImageFilterX<TImage>
ahcPhyY.name = sstmPhyY.str();
ahcPhyY.type = 3;
ahcPhyY.info[0] = idx; // Input image #ID
- m_VAllowedVarName.push_back(ahcPhyY);
+ m_VAllowedVarNameAuto.push_back(ahcPhyY);
//imi
std::stringstream sstm;
@@ -130,7 +131,7 @@ void BandMathImageFilterX<TImage>
ahc.name = sstm.str();
ahc.type = 4;
ahc.info[0] = idx; // Input image #ID
- m_VAllowedVarName.push_back(ahc);
+ m_VAllowedVarNameAuto.push_back(ahc);
//imibj
for (int j=0; j<image->GetNumberOfComponentsPerPixel(); j++)
@@ -142,25 +143,25 @@ void BandMathImageFilterX<TImage>
ahc.type = 5;
ahc.info[0] = idx; // Input image #ID
ahc.info[1] = j; // Band #ID
- m_VAllowedVarName.push_back(ahc);
+ m_VAllowedVarNameAuto.push_back(ahc);
}
//imibjNkxp
- int size=0;
+ int size=10;
for (int j=0; j<image->GetNumberOfComponentsPerPixel(); j++)
- for(int k=1; k<=size; k++)
- for(int p=1; p<=size; p++)
+ for(int x=0; x<=size; x++)
+ for(int y=0; y<=size; y++)
{
std::stringstream sstm;
adhocStruct ahc;
- sstm << varName << "b" << (j+1) << "N" << k << "x" << p;
+ sstm << varName << "b" << (j+1) << "N" << 2*x+1 << "x" << 2*y+1;
ahc.name = sstm.str();
ahc.type = 6;
ahc.info[0] = idx; // Input image #ID
- ahc.info[1] = j; // Band #ID
- ahc.info[2] = k;
- ahc.info[3] = p;
- m_VAllowedVarName.push_back(ahc);
+ 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);
}
}
@@ -184,6 +185,68 @@ void BandMathImageFilterX<TImage>
this->Modified();
}
+template< typename TImage >
+void BandMathImageFilterX<TImage>
+::SetMatrix(const std::string& name, const std::string& definition)
+{
+ if (name.empty())
+ itkExceptionMacro(<< "Please, set a name for your matrix/vector." << std::endl);
+
+ if (definition.empty())
+ itkExceptionMacro(<< "Please, set the definition of your matrix/vector." << std::endl);
+
+ 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;
+ std::istringstream iss( def );
+ 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;
+ //std::cout << val << std::endl;
+ mat.back().push_back(val);
+ }
+ }
+
+ for(int i=0; i<mat.size(); i++)
+ {
+ std::cout << std::endl;
+ for(int j=0; j<mat[i].size(); j++)
+ std::cout << mat[i][j] << " ";
+ }
+ std::cout << std::endl;
+
+ //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);
+
+
+ std::stringstream sstm;
+ adhocStruct ahc;
+ sstm << name;
+ ahc.name = sstm.str();
+ 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)
+ m_VAllowedVarNameAddedByUser.push_back(ahc);
+
+}
+
template< typename TImage >
std::string BandMathImageFilterX<TImage>
::GetExpression(int IDExpression) const
@@ -191,6 +254,7 @@ std::string BandMathImageFilterX<TImage>
return m_Expression.at(IDExpression);
}
+
template< typename TImage >
std::vector<std::string>& BandMathImageFilterX<TImage>
::GetVarNames() const
@@ -200,7 +264,6 @@ std::vector<std::string>& BandMathImageFilterX<TImage>
res.push_back(m_VVarName[y].name);
return res;
-
}
@@ -226,27 +289,34 @@ void BandMathImageFilterX<TImage>
// Generate variables names
m_VVarName.clear();
m_VNotAllowedVarName.clear();
+ m_VFinalAllowedVarName.clear();
+
+ //TODO
+ for(int i=0; i<m_VAllowedVarNameAddedByUser.size(); i++)
+ m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAddedByUser[i]);
+ for(int i=0; i<m_VAllowedVarNameAuto.size(); i++)
+ m_VFinalAllowedVarName.push_back(m_VAllowedVarNameAuto[i]);
this->SetNumberOfRequiredOutputs((int) m_Expression.size());
- for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression)
+ for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression) // For each expression
{
ParserType::Pointer dummyParser = ParserType::New();
dummyParser->SetExpr(this->GetExpression(IDExpression));
mup::var_maptype vmap = dummyParser->GetExprVar();
for (mup::var_maptype::iterator item = vmap.begin(); item!=vmap.end(); ++item)
- {
+ {//std::cout << "item->first : " << item->first << std::endl;
bool OK=false; int i=0;
- while( ( !OK ) && (i<m_VAllowedVarName.size()) )
+ while( ( !OK ) && (i<m_VFinalAllowedVarName.size()) )
{
- if (item->first == m_VAllowedVarName[i].name)
+ if (item->first == m_VFinalAllowedVarName[i].name)
OK=true;
else
i++;
}
- if (OK) {AddVariable(m_VAllowedVarName[i]);}
+ if (OK) {AddVariable(m_VFinalAllowedVarName[i]); }
else {
adhocStruct ahc;
ahc.name = item->first;
@@ -254,9 +324,10 @@ void BandMathImageFilterX<TImage>
}
}
- }
+ }// At this step, m_VVarName has been built
+ //Checking formulas consistency
if (m_VNotAllowedVarName.size()>0)
{
std::stringstream sstm;
@@ -267,7 +338,7 @@ void BandMathImageFilterX<TImage>
itkExceptionMacro(<< sstm.str());
}
-
+
// Register variables for each parser (important : one parser per thread)
m_VParser.clear();
unsigned int nbThreads = this->GetNumberOfThreads();
@@ -282,10 +353,10 @@ void BandMathImageFilterX<TImage>
m_AImage.resize(nbThreads);
- double initValue = 1.0;
- for(int i = 0; i < nbThreads; ++i)
+ double initValue = 0.1;
+ for(int i = 0; i < nbThreads; ++i) // For each thread
{
- m_AImage[i].resize(m_NbVar);
+ m_AImage[i].resize(m_NbVar); // For each variable
for(int j=0; j < m_NbVar; ++j)
{
@@ -325,14 +396,23 @@ void BandMathImageFilterX<TImage>
if (m_AImage[i][j].type == 6 ) // neighborhood
{
- //TODO
+ m_AImage[i][j].value = ValueType(m_AImage[i][j].info[3],m_AImage[i][j].info[2],initValue);
+ }
+
+ if (m_AImage[i][j].type == 7 ) // TODO
+ {
+ m_AImage[i][j].value = ValueType(m_AImage[i][j].info[1],m_AImage[i][j].info[0],0.0);
+ m_AImage[i][j].value.At(0,0)=0.1;
+ m_AImage[i][j].value.At(0,1)=0.5;
+ m_AImage[i][j].value.At(0,2)=0.1;
+
}
m_VParser.at(i)->DefineVar(m_AImage[i][j].name, &(m_AImage[i][j].value));
- initValue += 0.01;
- if (initValue>10.0)
- initValue=1.0;
+ initValue += 0.001;
+ if (initValue>1.0)
+ initValue=0.1;
}
}
@@ -366,10 +446,15 @@ void BandMathImageFilterX< TImage >
case 'm':
mup::matrix_type vect = value.GetArray();
- m_outputsDimensions.push_back(vect.GetCols());
+ if ( vect.GetRows() == 1 ) //Vector
+ m_outputsDimensions.push_back(vect.GetCols());
+ else //Matrix
+ itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
break;
}
+
+ std::cout << "Type = " << value.GetType() << " dimension = " << m_outputsDimensions.back() << std::endl;
}
}
@@ -383,7 +468,7 @@ void BandMathImageFilterX< TImage >
typedef itk::ImageBase< TImage::ImageDimension > ImageBaseType;
typename ImageBaseType::Pointer outputPtr;
- PrepareParsers(); // addition
+ PrepareParsers(); // addition
OutputsDimensions(); // addition
// Allocate the output memory
@@ -400,7 +485,7 @@ void BandMathImageFilterX< TImage >
if ( outputPtr )
{
outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
- outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]);
+ outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]); // addition
outputPtr->Allocate();
}
@@ -485,7 +570,9 @@ void BandMathImageFilterX<TImage>
ValueType value;
unsigned int nbInputImages = this->GetNumberOfInputs();
- // Iterators
+ //----------------- --------- -----------------//
+ //----------------- Iterators -----------------//
+ //----------------- --------- -----------------//
typedef itk::ImageRegionConstIterator<TImage> ImageRegionConstIteratorType;
std::vector< ImageRegionConstIteratorType > Vit;
Vit.resize(nbInputImages);
@@ -499,17 +586,36 @@ void BandMathImageFilterX<TImage>
VoutIt[j] = ImageRegionConstIteratorType (this->GetOutput(j), outputRegionForThread);
+ //Special case : neighborhoods
+
+ //std::vector< adhocStruct > ngbhName;
+ std::vector< itk::ConstNeighborhoodIterator<TImage> > VNit;
+ for(int j=0; j<m_VVarName.size(); ++j)
+ if (m_VVarName[j].type == 6)
+ {
+ 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
+ VNit.back().NeedToUseBoundaryConditionOn(); //TODO
+ }
+
+
// Support progress methods/callbacks
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
- // Pixel affectation
- for(int j=0; j < nbInputImages; ++j) {Vit[j].GoToBegin();}
- for(int j=0; j < m_Expression.size(); ++j) {VoutIt[j].GoToBegin();}
- //outIt.GoToBegin();
+
+ //----------------- --------------------- -----------------//
+ //----------------- Variable affectations -----------------//
+ //----------------- --------------------- -----------------//
+ for(int j=0; j < nbInputImages; ++j) { Vit[j].GoToBegin(); }
+ for(int j=0; j < m_Expression.size(); ++j) { VoutIt[j].GoToBegin(); }
+ for(int j=0; j < VNit.size(); ++j) { VNit[j].GoToBegin(); }
while(!Vit.at(0).IsAtEnd()) // For each pixel
{
+ int ngbhNameIndex=0; int index=0;
for(int j=0; j < m_AImage[threadId].size(); ++j) // For each variable, perform a copy
{
@@ -546,14 +652,35 @@ void BandMathImageFilterX<TImage>
case 6 : //neighborhood
// TODO
+
+ // 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")
+
+ for(int cols=0; cols<m_AImage[threadId][j].info[2]; ++cols)
+ for(int rows=0; rows<m_AImage[threadId][j].info[3]; ++rows)
+ {
+ m_AImage[threadId][j].value.At(rows,cols) = VNit[ngbhNameIndex].GetPixel(index)[m_AImage[threadId][j].info[1]];
+ index++;
+ }
+
+ ngbhNameIndex++;
+ break;
+
+ case 7 :
+ //Nothing to do
break;
default :
itkExceptionMacro(<< "Type of the variable is unknown");
break;
}
- }
+ }//End while
+
+ //----------------- ----------- -----------------//
+ //----------------- Evaluations -----------------//
+ //----------------- ----------- -----------------//
for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression)
{
@@ -579,18 +706,22 @@ void BandMathImageFilterX<TImage>
break;
case 'c':
- itkExceptionMacro(<< "Complex numbers not supported." << std::endl);
+ itkExceptionMacro(<< "Complex numbers are not supported." << std::endl);
break;
case 'm':
mup::matrix_type vect = value.GetArray();
- for(int p=0; p<vect.GetCols(); ++p)
- VoutIt[IDExpression].Get()[p] = vect.At(0,p);
- break;
+ if ( vect.GetRows() == 1 ) //Vector
+ for(int p=0; p<vect.GetCols(); ++p)
+ VoutIt[IDExpression].Get()[p] = vect.At(0,p).GetFloat();
+ else //Matrix
+ itkExceptionMacro(<< "Result of the evaluation can't be a matrix." << std::endl);
+ break;
}
-
+
+ //----------------- Pixel affectations -----------------//
for(int p=0; p<VoutIt[IDExpression].Get().GetSize(); ++p)
{
// Case value is equal to -inf or inferior to the minimum value
@@ -611,9 +742,10 @@ void BandMathImageFilterX<TImage>
}
- for(int j=0; j < nbInputImages; ++j) {++Vit[j];}
- for(int j=0; j < m_Expression.size(); ++j) {++VoutIt[j];}
-
+ for(int j=0; j < nbInputImages; ++j) { ++Vit[j]; }
+ for(int j=0; j < m_Expression.size(); ++j) { ++VoutIt[j]; }
+ for(int j=0; j < VNit.size(); ++j) { ++VNit[j]; }
+
progress.CompletedPixel();
}
diff --git a/Code/Common/otbParserX.cxx b/Code/Common/otbParserX.cxx
index b3082a34bf..9b4fb4d314 100644
--- a/Code/Common/otbParserX.cxx
+++ b/Code/Common/otbParserX.cxx
@@ -58,7 +58,7 @@ public:
virtual void InitFun()
{
m_MuParserX.DefineFun(new ndvi);
- //m_MuParserX.DefineFun(new conv);
+ m_MuParserX.DefineFun(new conv);
m_MuParserX.DefineFun(new bands);
m_MuParserX.DefineFun(new vcos);
m_MuParserX.DefineFun(new vsin);
diff --git a/Code/Common/otbParserXPlugins.h b/Code/Common/otbParserXPlugins.h
index 698365eee8..d2495f73d5 100644
--- a/Code/Common/otbParserXPlugins.h
+++ b/Code/Common/otbParserXPlugins.h
@@ -65,53 +65,75 @@ public:
};
-
-class ndvi : public mup::ICallback
+class conv : public mup::ICallback
{
public:
- ndvi():ICallback(mup::cmFUNC, "ndvi", 2)
+ conv():ICallback(mup::cmFUNC, "conv", -1)
{}
virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
// Get the argument from the argument input vector
- mup::float_type r = a_pArg[0]->GetFloat();
+ mup::matrix_type m1 = a_pArg[0]->GetArray();
- mup::float_type niri = a_pArg[1]->GetFloat();
+ int nbrows = m1.GetRows();
+ int nbcols = m1.GetCols();
- // The return value is passed by writing it to the reference ret
- if ( vcl_abs(r + niri) < 1E-6 )
- *ret = 0.;
- else
- *ret = (niri-r)/(niri+r);
+ //matrix_type m3(nbrows,nbcols,0);
+
+ mup::matrix_type res(a_iArgc-1,1,0);
+
+ for (int k=1; k<a_iArgc; ++k)
+ {
+
+ float sum=0.0;
+
+ mup::matrix_type m2 = a_pArg[k]->GetArray();
+
+ for (int i=0; i<nbrows; i++)
+ for (int j=0; j<nbcols; j++)
+ sum += m1.At(i,j).GetFloat() * m2.At(i,j).GetFloat();
+
+
+ res.At(k-1,0)=sum;
+ }
+
+
+ // The return value is passed by writing it to the reference ret
+ *ret = res;
}
const mup::char_type* GetDesc() const
{
- return "NDVI - Normalized Difference Vegetation Index";
+ return "conv(m1,m2) - A matrix convolution";
}
mup::IToken* Clone() const
{
- return new ndvi(*this);
+ return new conv(*this);
}
};
-//--------------------------------------------------------------------------------------------------------//
-class vcos : public mup::ICallback
+
+class MatDiv : public mup::IOprtBin
{
-public:
- vcos():ICallback(mup::cmFUNC, "vcos", 1)
+ public:
+ MatDiv():IOprtBin(_T("div"), (int)(mup::prMUL_DIV), mup::oaLEFT)
{}
- virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+ virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
- // Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
+ const mup::matrix_type b = a_pArg[1]->GetArray();
+ /*if (!arg1->IsNonComplexScalar())
+ throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg1->GetType(), 'f', 1));
+ if (!arg2->IsNonComplexScalar())
+ throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg2->GetType(), 'f', 2));*/
+
int nbrows = a.GetRows();
int nbcols = a.GetCols();
@@ -119,37 +141,43 @@ public:
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
- res.At(p,k) = vcl_cos(a.At(p,k).GetFloat());
+ res.At(p,k) = a.At(p,k).GetFloat() / b.At(p,k).GetFloat();
// The return value is passed by writing it to the reference ret
- *ret = res;
+ *ret = res;
}
- const mup::char_type* GetDesc() const
- {
- return "vcos - Cosinus for noncomplex vectors & matrices";
+ virtual const mup::char_type* GetDesc() const
+ {
+ return _T("x ./ y - Division for noncomplex vectors & matrices");
}
-
- mup::IToken* Clone() const
- {
- return new vcos(*this);
+
+ virtual mup::IToken* Clone() const
+ {
+ return new MatDiv(*this);
}
};
-class vsin : public mup::ICallback
+
+class MatConv : public mup::IOprtBin
{
-public:
- vsin():ICallback(mup::cmFUNC, "vsin", 1)
+ public:
+ MatConv():IOprtBin(_T("mult"), (int)(mup::prMUL_DIV), mup::oaLEFT)
{}
- virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+ virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
{
- // Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
+ const mup::matrix_type b = a_pArg[1]->GetArray();
+ /*if (!arg1->IsNonComplexScalar())
+ throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg1->GetType(), 'f', 1));
+ if (!arg2->IsNonComplexScalar())
+ throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg2->GetType(), 'f', 2));*/
+
int nbrows = a.GetRows();
int nbcols = a.GetCols();
@@ -157,41 +185,71 @@ public:
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
- res.At(p,k) = vcl_sin(a.At(p,k).GetFloat());
-
+ res.At(p,k) = a.At(p,k).GetFloat() * b.At(p,k).GetFloat();
+
// The return value is passed by writing it to the reference ret
- *ret = res;
+ *ret = res;
+ }
+
+ virtual const mup::char_type* GetDesc() const
+ {
+ return _T("x mult y - Convolution for noncomplex vectors & matrices");
+ }
+
+ virtual mup::IToken* Clone() const
+ {
+ return new MatConv(*this);
+ }
+ };
+
+
+
+class ndvi : public mup::ICallback
+ {
+public:
+ ndvi():ICallback(mup::cmFUNC, "ndvi", 2)
+ {}
+
+ virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
+ {
+ // Get the argument from the argument input vector
+ mup::float_type r = a_pArg[0]->GetFloat();
+
+ mup::float_type niri = a_pArg[1]->GetFloat();
+
+
+ // The return value is passed by writing it to the reference ret
+ if ( vcl_abs(r + niri) < 1E-6 )
+ *ret = 0.;
+ else
+ *ret = (niri-r)/(niri+r);
}
const mup::char_type* GetDesc() const
{
- return "vsin - Sinus for noncomplex vectors & matrices";
+ return "NDVI - Normalized Difference Vegetation Index";
}
mup::IToken* Clone() const
{
- return new vsin(*this);
+ return new ndvi(*this);
}
};
-class MatDiv : public mup::IOprtBin
+//--------------------------------------------------------------------------------------------------------//
+class vcos : public mup::ICallback
{
- public:
- MatDiv():IOprtBin(_T("div"), (int)(mup::prMUL_DIV), mup::oaLEFT)
+public:
+ vcos():ICallback(mup::cmFUNC, "vcos", 1)
{}
- virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
+ virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
+ // Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
- const mup::matrix_type b = a_pArg[1]->GetArray();
- /*if (!arg1->IsNonComplexScalar())
- throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg1->GetType(), 'f', 1));
- if (!arg2->IsNonComplexScalar())
- throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg2->GetType(), 'f', 2));*/
-
int nbrows = a.GetRows();
int nbcols = a.GetCols();
@@ -199,41 +257,37 @@ class MatDiv : public mup::IOprtBin
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
- res.At(p,k) = a.At(p,k).GetFloat() / b.At(p,k).GetFloat();
+ res.At(p,k) = vcl_cos(a.At(p,k).GetFloat());
// The return value is passed by writing it to the reference ret
- *ret = res;
+ *ret = res;
}
- virtual const mup::char_type* GetDesc() const
- {
- return _T("x ./ y - Division for noncomplex vectors & matrices");
+ const mup::char_type* GetDesc() const
+ {
+ return "vcos - Cosinus for noncomplex vectors & matrices";
}
-
- virtual mup::IToken* Clone() const
- {
- return new MatDiv(*this);
+
+ mup::IToken* Clone() const
+ {
+ return new vcos(*this);
}
};
-class MatConv : public mup::IOprtBin
+
+class vsin : public mup::ICallback
{
- public:
- MatConv():IOprtBin(_T("mult"), (int)(mup::prMUL_DIV), mup::oaLEFT)
+public:
+ vsin():ICallback(mup::cmFUNC, "vsin", 1)
{}
- virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int)
+ virtual void Eval(mup::ptr_val_type &ret, const mup::ptr_val_type *a_pArg, int a_iArgc)
{
+ // Get the argument from the argument input vector
const mup::matrix_type a = a_pArg[0]->GetArray();
- const mup::matrix_type b = a_pArg[1]->GetArray();
- /*if (!arg1->IsNonComplexScalar())
- throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg1->GetType(), 'f', 1));
- if (!arg2->IsNonComplexScalar())
- throw ParserError( ErrorContext(ecTYPE_CONFLICT_FUN, -1, GetIdent(), arg2->GetType(), 'f', 2));*/
-
int nbrows = a.GetRows();
int nbcols = a.GetCols();
@@ -241,24 +295,26 @@ class MatConv : public mup::IOprtBin
for (int k=0; k<nbcols; ++k)
for (int p=0; p<nbrows; ++p)
- res.At(p,k) = a.At(p,k).GetFloat() * b.At(p,k).GetFloat();
-
+ res.At(p,k) = vcl_sin(a.At(p,k).GetFloat());
+
// The return value is passed by writing it to the reference ret
- *ret = res;
+ *ret = res;
}
- virtual const mup::char_type* GetDesc() const
- {
- return _T("x mult y - Convolution for noncomplex vectors & matrices");
+ const mup::char_type* GetDesc() const
+ {
+ return "vsin - Sinus for noncomplex vectors & matrices";
}
-
- virtual mup::IToken* Clone() const
- {
- return new MatConv(*this);
+
+ mup::IToken* Clone() const
+ {
+ return new vsin(*this);
}
};
+
+
}//end namespace otb
#endif
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index d147b73c4a..1096ef9b93 100644
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -2117,22 +2117,24 @@ add_test(bfTvImageAndVectorImageOperationFilterTest ${BASICFILTERS_TESTS12}
add_test(bfTuBandMathImageFilter ${BASICFILTERS_TESTS13}
otbBandMathImageFilterNew)
-add_test(bfTuBandMathImageFilterX ${BASICFILTERS_TESTS13}
- otbBandMathImageFilterXNew)
-
add_test(bfTvBandMathImageFilter ${BASICFILTERS_TESTS13}
otbBandMathImageFilter)
-add_test(bfTvBandMathImageFilterX ${BASICFILTERS_TESTS13}
- otbBandMathImageFilterX)
-
-
add_test(bfTvBandMathImageFilterWithIdx ${BASICFILTERS_TESTS13}
otbBandMathImageFilterWithIdx
${TEMP}/bfTvBandMathImageFilterWithIdx1.tif
${TEMP}/bfTvBandMathImageFilterWithIdx2.tif
)
+add_test(bfTuBandMathImageFilterXNew ${BASICFILTERS_TESTS13}
+ otbBandMathImageFilterXNew)
+
+add_test(bfTvBandMathImageFilterX ${BASICFILTERS_TESTS13}
+ otbBandMathImageFilterX)
+
+add_test(bfTvotbBandMathImageFilterXConv ${BASICFILTERS_TESTS13}
+ otbBandMathImageFilterXConv)
+
add_test(bfTvBandMathImageFilterXWithIdx ${BASICFILTERS_TESTS13}
otbBandMathImageFilterXWithIdx
${TEMP}/bfTvBandMathImageFilterWithIdx1.tif
diff --git a/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx b/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx
index 6ab50093e4..9930766e55 100644
--- a/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx
+++ b/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx
@@ -199,6 +199,154 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) [])
}
+int otbBandMathImageFilterXConv( int itkNotUsed(argc), char* itkNotUsed(argv) [])
+{
+
+ typedef otb::VectorImage<double, 2> ImageType;
+ typedef ImageType::PixelType PixelType;
+ typedef otb::BandMathImageFilterX<ImageType> FilterType;
+
+ unsigned int i;
+ const unsigned int N = 100, D1=3, D2=1, D3=1;
+ unsigned int FAIL_FLAG = 0;
+
+ ImageType::SizeType size;
+ size.Fill(N);
+ ImageType::IndexType index;
+ index.Fill(0);
+ ImageType::RegionType region;
+ region.SetSize(size);
+ region.SetIndex(index);
+
+ ImageType::Pointer image1 = ImageType::New();
+ ImageType::Pointer image2 = ImageType::New();
+ ImageType::Pointer image3 = ImageType::New();
+
+ image1->SetLargestPossibleRegion( region );
+ image1->SetBufferedRegion( region );
+ image1->SetRequestedRegion( region );
+ image1->SetNumberOfComponentsPerPixel(D1);
+ image1->Allocate();
+
+ image2->SetLargestPossibleRegion( region );
+ image2->SetBufferedRegion( region );
+ image2->SetRequestedRegion( region );
+ image2->SetNumberOfComponentsPerPixel(D2);
+ image2->Allocate();
+
+ image3->SetLargestPossibleRegion( region );
+ image3->SetBufferedRegion( region );
+ image3->SetRequestedRegion( region );
+ image3->SetNumberOfComponentsPerPixel(D3);
+ image3->Allocate();
+
+ typedef itk::ConstNeighborhoodIterator<ImageType> IteratorType;
+ typename IteratorType::RadiusType radius;
+ radius[0]=1; // Size x direction
+ radius[1]=0; // Size y direction
+
+ IteratorType it1(radius, image1, region);
+ IteratorType it2(radius, image2, region);
+ IteratorType it3(radius, image3, region);
+
+ for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3)
+ {
+ ImageType::IndexType i1 = it1.GetIndex();
+ ImageType::IndexType i2 = it2.GetIndex();
+ ImageType::IndexType i3 = it3.GetIndex();
+
+ it1.GetCenterPixel()[0] = i1[0] + i1[1] -50; it1.GetCenterPixel()[1] = i1[0] * i1[1] -50; it1.GetCenterPixel()[2] = i1[0] / (i1[1]+1)+5;
+ it2.GetCenterPixel()[0] = i2[0] * i2[1];
+ it3.GetCenterPixel()[0] = i3[0] + i3[1] * i3[1];
+
+ }
+
+
+ FilterType::Pointer filter = FilterType::New();
+ std::cout << "Number Of Threads : " << filter->GetNumberOfThreads() << std::endl;
+
+
+ filter->SetNthInput(0, image1);
+ filter->SetNthInput(1, image2);
+ filter->SetNthInput(2, image3, "canal3");
+ filter->SetMatrix("kernel1N3x1","{0.1 , 0.5 , 0.1}");
+ filter->SetExpression("conv(kernel1N3x1,im1b1N3x1)");
+ filter->Update();
+
+ if (filter->GetNumberOfOutputs() != 1)
+ itkGenericExceptionMacro( <<std::endl << "Error = wrong number of outputs ");
+
+
+ ImageType::Pointer output1 = filter->GetOutput(0);
+
+ if (output1->GetNumberOfComponentsPerPixel() != 1)
+ itkGenericExceptionMacro( <<std::endl << "Error = wrong number of components per pixel ");
+
+ std::cout << "\n--- Standard Use\n";
+ std::cout << "Parsed Expression : " << filter->GetExpression(0) << std::endl;
+
+
+
+ //Sub-test 1
+ IteratorType itoutput1(radius, output1, region);
+
+ for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(), itoutput1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3, ++itoutput1)
+ {
+ ImageType::IndexType i1 = it1.GetIndex();
+ ImageType::IndexType i2 = it2.GetIndex();
+ ImageType::IndexType i3 = it3.GetIndex();
+ PixelType px1(D1),px2(D2),px3(D3);
+
+ px1[0]=0; float coefs[3]={0.1,0.5,0.1};
+ for(int i=0; i<it1.Size(); ++i)
+ px1[0] += coefs[i]*it1.GetPixel(i)[0];
+ //px1[1] = ( i1[0] * i1[1] -50 ); px1[2] = ( i1[0] / (i1[1]+1)+5 );
+
+
+ px2[0] = ( i2[0] * i2[1] );
+ px3[0] = ( i3[0] + i3[1] * i3[1] );
+
+ double result1 = itoutput1.GetCenterPixel()[0], result2 = itoutput1.GetCenterPixel()[1], result3 = itoutput1.GetCenterPixel()[2];
+ double error1,error2,error3;
+
+ double expected1 = px1[0];
+ /*double expected2 = px1[1];
+ double expected3 = px1[2]; */
+
+ /*std::cout << "Pixel_1 = " << it1.Get()[0] << " Pixel_2 = " << it2.Get()[0] << " Pixel_3 = " << it3.Get()[0]
+ << " Result = " << itoutput1.Get()[0] << " Expected = " << expected1 << std::endl;*/
+
+
+ error1 = (result1 - expected1) * (result1 - expected1) / (result1 + expected1);
+ /*error2 = (result2 - expected2) * (result2 - expected2) / (result2 + expected2);
+ error3 = (result3 - expected3) * (result3 - expected3) / (result3 + expected3);*/
+
+ if ( ( error1 > 1E-9 ) /*|| ( error2 > 1E-9 ) || ( error3 > 1E-9 )*/)
+ {
+ itkGenericExceptionMacro( <<std::endl
+<< "Error = " << error1 << " > 1E-9 -> TEST FAILLED" << std::endl
+ << "Pixel_[-1]_band1 = " << it1.GetPixel(0)[0]
+ << " Pixel_[0]_band1 = " << it1.GetPixel(1)[0]
+ << " Pixel_[1]_band1 = " << it1.GetPixel(2)[0]
+ << " Result = " << result1
+ << " Expected = " << expected1 << std::endl);
+/*<< "Error = " << error2 << " > 1E-9 -> TEST FAILLED" << std::endl
+ << "Pixel_[-1]_band2 = " << it1.GetPixel(0)[1]
+ << " Pixel_[0]_band2 = " << it1.GetPixel(1)[1]
+ << " Pixel_[1]_band2 = " << it1.GetPixel(2)[1]
+ << " Result = " << result2
+ << " Expected = " << expected2 << std::endl
+<< "Error = " << error3 << " > 1E-9 -> TEST FAILLED" << std::endl
+ << "Pixel_[-1]_band3 = " << it1.GetPixel(0)[2]
+ << " Pixel_[0]_band3 = " << it1.GetPixel(1)[2]
+ << " Pixel_[1]_band3 = " << it1.GetPixel(2)[2]
+ << " Result = " << result3
+ << " Expected = " << expected3 << std::endl);*/
+ }
+ }
+
+ return EXIT_SUCCESS;
+}
int otbBandMathImageFilterXWithIdx( int itkNotUsed(argc), char* argv[])
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests13.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests13.cxx
index 3784e866ee..2d95f95d03 100644
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests13.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests13.cxx
@@ -25,9 +25,10 @@
void RegisterTests()
{
REGISTER_TEST(otbBandMathImageFilterNew);
- REGISTER_TEST(otbBandMathImageFilterXNew);
REGISTER_TEST(otbBandMathImageFilter);
+ REGISTER_TEST(otbBandMathImageFilterXNew);
REGISTER_TEST(otbBandMathImageFilterX);
+ REGISTER_TEST(otbBandMathImageFilterXConv);
REGISTER_TEST(otbBandMathImageFilterWithIdx);
REGISTER_TEST(otbBandMathImageFilterXWithIdx);
REGISTER_TEST(otbComplexToIntensityFilterTest);
--
GitLab