Commit e0e7e11d authored by Christophe Palmann's avatar Christophe Palmann

TEST: bandmathx, use of global stats OK

parent daf58361
......@@ -25,6 +25,7 @@
#include "itkImageToImageFilter.h"
#include "itkArray.h"
#include "otbStreamingStatisticsVectorImageFilter.h"
#include "otbParserX.h"
#include <vector>
......@@ -33,43 +34,21 @@
namespace otb
{
/** \class BandMathImageFilterX
* \brief Performs a mathematical operation on the input images
* \brief Performs mathematical operations on the input images
* according to the formula specified by the user.
*
* This filter is based on the mathematical parser library muParser.
* This filter is based on the mathematical parser library muParserX.
* 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
* \url{http:*articles.beltoforion.de/article.php?a=muparserx}.
*
* 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.
*
* the corresponding SetNthInput method. For the jth (j=1..T) input image, if
* no associated variable name has been specified, a default variable
* name is given by concatenating the prefix "im" with the
* corresponding input index plus one (for instance, im1 is related to the first input).
* If the jth input image is multidimensional, then the variable imj represents a vector whose components are related to its bands.
* In order to access the kth band, the variable observes the following pattern : imjbk.
*
* \sa Parser
*
......@@ -107,6 +86,11 @@ public:
typedef ParserX ParserType;
typedef typename ParserType::ValueType ValueType;
/** Typedef for statistic computing. */
typedef StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType;
typedef typename StreamingStatisticsVectorImageFilterType::Pointer StreamingStatisticsVectorImageFilterPointerType;
typedef typename StreamingStatisticsVectorImageFilterType::MatrixType MatrixType;
/** 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);
......
......@@ -29,8 +29,6 @@
#include "itkProgressReporter.h"
#include "otbMacro.h"
#include <otbStreamingStatisticsVectorImageFilter.h>
#include <iostream>
#include <fstream>
#include <string>
......@@ -97,7 +95,9 @@ void BandMathImageFilterX<TImage>
{
Superclass::PrintSelf(os, indent);
os << indent << "Expression: " << m_Expression[0] << std::endl; //TODO
os << indent << "Expressions: " << std::endl;
for (unsigned int i=0; i<m_Expression.size(); i++)
os << indent << m_Expression[i] << std::endl;
os << indent << "Computed values follow:" << std::endl;
os << indent << "UnderflowCount: " << m_UnderflowCount << std::endl;
os << indent << "OverflowCount: " << m_OverflowCount << std::endl;
......@@ -154,7 +154,7 @@ void BandMathImageFilterX<TImage>
//Mandatory before call of GetNumberOfComponentsPerPixel
//Really important not to call the filter's UpdateOutputInformation method here:
//this method is not ready until all inputs, variables and expressions are set.
imagebis->UpdateOutputInformation();
imagebis->UpdateOutputInformation();
//imibj
for (int j=0; j<imagebis->GetNumberOfComponentsPerPixel(); j++)
......@@ -191,6 +191,8 @@ void BandMathImageFilterX<TImage>
statsNames.push_back("Mini");
statsNames.push_back("Maxi");
statsNames.push_back("Mean");
statsNames.push_back("Sum");
statsNames.push_back("Var");
for (int j=0; j<imagebis->GetNumberOfComponentsPerPixel(); j++)
for (int t=0; t<statsNames.size(); t++)
......@@ -245,6 +247,7 @@ void BandMathImageFilterX<TImage>
}
template< typename TImage >
void BandMathImageFilterX<TImage>
::SetMatrix(const std::string& name, const std::string& definition)
......@@ -357,7 +360,7 @@ void BandMathImageFilterX<TImage>
iss << " " << m_VAllowedVarNameAddedByUser[i].value.At(k,0);
for(int p=1; p<m_VAllowedVarNameAddedByUser[i].value.GetCols(); p++)
iss << " , " << m_VAllowedVarNameAddedByUser[i].value.At(k,p);
iss << ";";
iss << " ;";
}
str=iss.str();
str.erase(str.size()-1);
......@@ -548,7 +551,7 @@ void BandMathImageFilterX<TImage>
m_VFinalAllowedVarName.clear();
// m_VFinalAllowedVarName = m_VAllowedVarNameAuto + m_VAllowedVarNameAddedByUser
// m_VFinalAllowedVarName = variable names dictionary
// m_VFinalAllowedVarName = variable names dictionary
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++)
......@@ -607,7 +610,7 @@ void BandMathImageFilterX<TImage>
// Important to remember that variables of m_VVarName come from a call of GetExprVar method
// Only useful variables are allocated in this filter
int nbVar = m_VVarName.size();
int nbVar = m_VVarName.size();
m_StatsVarDetected.clear();
......@@ -678,14 +681,14 @@ void BandMathImageFilterX<TImage>
if (m_AImage[i][j].type == 8 ) // global stats
{
m_AImage[i][j].value = ValueType(initValue);
//m_AImage[i][j].info[0] = Image ID : useful to know which images must have their regions set to largest possible region (see GenerateInputRequestedRegion)
m_AImage[i][j].value = ValueType(initValue);
//m_AImage[i][j].info[0] = Image ID : useful to know which images must have their regions set to largest possible region (see GenerateInputRequestedRegion)
bool found = false;
for (int r=0; r<m_StatsVarDetected.size() && !found; r++)
for (int r=0; r<m_StatsVarDetected.size() && !found; r++)
if (m_StatsVarDetected[r] == m_AImage[i][j].info[0])
found = true;
if (!found)
m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
if (!found)
m_StatsVarDetected.push_back(m_AImage[i][j].info[0]);
}
......@@ -789,13 +792,14 @@ void BandMathImageFilterX< TImage >
ImagePointer inputPtr = const_cast<TImage *>(this->GetInput(m_StatsVarDetected[i]));
inputPtr->SetRequestedRegionToLargestPossibleRegion();
}
else
itkExceptionMacro(<< "Requested input #" << m_StatsVarDetected[i] << ", but only " << this->GetNumberOfInputs() << " inputs are available." << std::endl);
else
itkExceptionMacro(<< "Requested input #" << m_StatsVarDetected[i] << ", but only " << this->GetNumberOfInputs() << " inputs are available." << std::endl);
}
}
template< typename TImage >
void BandMathImageFilterX<TImage>
::BeforeThreadedGenerateData()
......@@ -828,13 +832,13 @@ void BandMathImageFilterX<TImage>
for (unsigned int i=0; i<m_StatsVarDetected.size(); i++)
{
typedef otb::StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType;
typename StreamingStatisticsVectorImageFilterType::Pointer filter = StreamingStatisticsVectorImageFilterType::New();
StreamingStatisticsVectorImageFilterPointerType filter = StreamingStatisticsVectorImageFilterType::New();
filter->SetInput( this->GetNthInput(m_StatsVarDetected[i]) );
filter->Update();
PixelType pix; //Variable length vector
MatrixType mat;
for(int t=0; t<m_AImage.size(); t++) // for each thread
for(int v=0; v<m_AImage[t].size(); v++) // for each variable
......@@ -868,6 +872,26 @@ void BandMathImageFilterX<TImage>
m_AImage[t][v].value = pix[b];
break;
break;
case 3: // sum
pix = filter->GetSum();
for (unsigned int b=0; b<pix.GetSize(); b++) // for each band
if (m_AImage[t][v].info[1] == b) // info[1] : band ID
m_AImage[t][v].value = pix[b];
break;
case 4: // stddev
mat = filter->GetCovariance();
for (unsigned int b=0; b<mat.Cols(); b++) // for each band
if (m_AImage[t][v].info[1] == b) // info[1] : band ID
m_AImage[t][v].value = mat(b,b);
break;
}
}
}
......@@ -899,7 +923,7 @@ void BandMathImageFilterX<TImage>
m_OverflowCount += m_ThreadOverflow[i];
}
if((m_UnderflowCount != 0) || (m_OverflowCount!=0)) //TODO
if((m_UnderflowCount != 0) || (m_OverflowCount!=0))
{
std::stringstream sstm;
sstm << std::endl
......@@ -949,7 +973,7 @@ void BandMathImageFilterX<TImage>
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 : better to use ImageBoundaryFacesCalculator
VNit.back().NeedToUseBoundaryConditionOn();
}
......
......@@ -253,6 +253,8 @@ int otbBandMathImageFilterXConv( int itkNotUsed(argc), char* argv [])
typename PixelType::ValueType imageAb2Mini = itk::NumericTraits<typename PixelType::ValueType>::max();
typename PixelType::ValueType imageAb3Mean = 0.0, n=0.0;
typename PixelType::ValueType imageAb3Var = 0.0;
typename PixelType::ValueType imageAb1Sum = 0.0;
typename PixelType::ValueType im2b1Maxi = itk::NumericTraits<typename PixelType::ValueType>::min();
for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3)
......@@ -273,13 +275,20 @@ int otbBandMathImageFilterXConv( int itkNotUsed(argc), char* argv [])
imageAb3Mean += it1.GetCenterPixel()[2];
n++;
// Var of im1 band 3
imageAb3Var += vcl_pow(it1.GetCenterPixel()[2],2.0);
// Maximum of im2 band 1
if (im2b1Maxi < it2.GetCenterPixel()[0])
im2b1Maxi = it2.GetCenterPixel()[0];
//Sum of im1 band1
imageAb1Sum += it1.GetCenterPixel()[0];
}
imageAb3Mean = imageAb3Mean /n;
imageAb3Var = (n/(n-1)) * (imageAb3Var /n - imageAb3Mean*imageAb3Mean); //unbiased
FilterType::Pointer filter = FilterType::New();
std::cout << "Number Of Threads : " << filter->GetNumberOfThreads() << std::endl;
......@@ -293,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} ");
filter->SetExpression("(vmax(canal3b1N3x5)+vmin(canal3b1N3x5)) div {2.0} + {imageAb2Mini / im2b1Maxi} + {imageAb3Mean / imageAb1Sum * imageAb3Var} ");
filter->Update();
if (filter->GetNumberOfOutputs() != 2)
......@@ -367,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;
px2[0] = (vect2.back() + vect2.front())/2.0 + (imageAb2Mini / im2b1Maxi) + imageAb3Mean / imageAb1Sum * imageAb3Var;
......@@ -436,7 +445,7 @@ int otbBandMathImageFilterXTxt( int itkNotUsed(argc), char* argv [])
FilterType::Pointer filter = FilterType::New();
filter->ImportContext(inputFilename);
filter->ImportContext(inputFilename);
filter->ExportContext(outputFilename);
return EXIT_SUCCESS;
......
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