diff --git a/Code/BasicFilters/otbBandMathImageFilterX.h b/Code/BasicFilters/otbBandMathImageFilterX.h index b5603f26b2d305b0b03bdcfcef80624ce9aef9a3..62a66d8f6b4dc7f08ba4fe8d15a018cc86a97810 100644 --- a/Code/BasicFilters/otbBandMathImageFilterX.h +++ b/Code/BasicFilters/otbBandMathImageFilterX.h @@ -116,7 +116,7 @@ public: void SetExpression(const std::string& expression); /** Return the expression to be parsed */ - std::string GetExpression() const; + std::string GetExpression(int) const; /** Return the nth filter input associated variable name */ std::string GetNthInputName(unsigned int idx) const; @@ -129,6 +129,8 @@ protected : virtual ~BandMathImageFilterX(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; + void AllocateOutputs(); + void BeforeThreadedGenerateData(); void ThreadedGenerateData(const ImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); void AfterThreadedGenerateData(); @@ -151,9 +153,9 @@ private : BandMathImageFilterX(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented - void addVariable(adhocStruct&); - void generateVariables(); - + void AddVariable(adhocStruct&); + void GenerateVariables(); + void OutputsDimensions(); std::vector<std::string> m_Expression; std::vector<ParserType::Pointer> m_VParser; @@ -162,6 +164,7 @@ private : std::vector< adhocStruct > m_VAllowedVarName; std::vector< adhocStruct > m_VNotAllowedVarName; unsigned int m_NbVar; + std::vector< int > m_outputsDimensions; SpacingType m_Spacing; OrigineType m_Origin; diff --git a/Code/BasicFilters/otbBandMathImageFilterX.txx b/Code/BasicFilters/otbBandMathImageFilterX.txx index 08e58521d4419017d85fbb8fc00da44c7f2c2925..c6df9b262756a96bad4906bb8fa0906a18a04b03 100644 --- a/Code/BasicFilters/otbBandMathImageFilterX.txx +++ b/Code/BasicFilters/otbBandMathImageFilterX.txx @@ -62,7 +62,7 @@ BandMathImageFilterX<TImage> ahcY.type = 1; m_VAllowedVarName.push_back(ahcY); - //this->SetNumberOfThreads(1); + this->SetNumberOfThreads(1); } @@ -73,6 +73,7 @@ BandMathImageFilterX<TImage> { } + template <class TImage> void BandMathImageFilterX<TImage> ::PrintSelf(std::ostream& os, itk::Indent indent) const @@ -182,16 +183,19 @@ template< typename TImage > void BandMathImageFilterX<TImage> ::SetExpression(const std::string& expression) { - //if (m_Expression != expression) - m_Expression.push_back(expression); + m_Expression.push_back(expression); //TODO + + if (m_Expression.size()>1) + this->SetNthOutput( (int) (m_Expression.size()) -1, ( TImage::New() ).GetPointer() ); + this->Modified(); } template< typename TImage > std::string BandMathImageFilterX<TImage> -::GetExpression() const +::GetExpression(int IDExpression) const { - return m_Expression[0]; + return m_Expression.at(IDExpression); //TODO } template< typename TImage > @@ -204,7 +208,7 @@ std::string BandMathImageFilterX<TImage> template< typename TImage > void BandMathImageFilterX<TImage> -::addVariable(adhocStruct &ahc) +::AddVariable(adhocStruct &ahc) { bool found=false; for(int i=0; i<m_VVarName.size(); ++i) @@ -218,16 +222,19 @@ void BandMathImageFilterX<TImage> template< typename TImage > void BandMathImageFilterX<TImage> -::generateVariables() +::GenerateVariables() { + // Generate variables names m_VVarName.clear(); m_VNotAllowedVarName.clear(); + this->SetNumberOfRequiredOutputs((int) m_Expression.size()); + for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression) { ParserType::Pointer dummyParser = ParserType::New(); - dummyParser->SetExpr(this->GetExpression()); + dummyParser->SetExpr(this->GetExpression(IDExpression)); mup::var_maptype vmap = dummyParser->GetExprVar(); for (mup::var_maptype::iterator item = vmap.begin(); item!=vmap.end(); ++item) @@ -241,18 +248,21 @@ void BandMathImageFilterX<TImage> i++; } - if (OK) {addVariable(m_VAllowedVarName[i]);} //<<<<<<------------ + if (OK) {AddVariable(m_VAllowedVarName[i]);} else { adhocStruct ahc; ahc.name = item->first; m_VNotAllowedVarName.push_back(ahc); } } + } +/*for(int y=0; y<m_VAllowedVarName.size(); y++) + std::cout << "--> " << m_VAllowedVarName[y].name << " " << m_VAllowedVarName[y].type << std::endl; -/*for(int y=0; y<m_VVarName.size(); y++) - std::cout << "--> " << m_VVarName[y].name << " " << m_VVarName[y].type << std::endl;*/ +for(int y=0; y<m_VVarName.size(); y++) + std::cout << "---------> " << m_VVarName[y].name << " " << m_VVarName[y].type << std::endl;*/ if (m_VNotAllowedVarName.size()>0) { @@ -264,42 +274,10 @@ void BandMathImageFilterX<TImage> itkExceptionMacro(<< sstm.str()); } -} - -template< typename TImage > -void BandMathImageFilterX<TImage> -::BeforeThreadedGenerateData() -{ - // Some useful variable + + // Register variables for each parser (important : one parser per thread) + m_VParser.clear(); unsigned int nbThreads = this->GetNumberOfThreads(); - unsigned int nbInputImages = this->GetNumberOfInputs(); - - // Check if input image dimensions matches - unsigned int inputSize[2]; - inputSize[0] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(0); - inputSize[1] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(1); - - for(unsigned int p = 1; p < nbInputImages; p++) - { - if((inputSize[0] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0)) - || (inputSize[1] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1))) - { - itkExceptionMacro(<< "Input images must have the same dimensions." << std::endl - << "band #1 is [" << inputSize[0] << ";" << inputSize[1] << "]" << std::endl - << "band #" << p+1 << " is [" - << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0) << ";" - << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1) << "]"); - } - } - - // Allocate and initialize the thread temporaries - m_ThreadUnderflow.SetSize(nbThreads); - m_ThreadUnderflow.Fill(0); - m_ThreadOverflow.SetSize(nbThreads); - m_ThreadOverflow.Fill(0); - - - // Important : one parser for each thread typename std::vector<ParserType::Pointer>::iterator itParser; m_VParser.resize(nbThreads); for(itParser = m_VParser.begin(); itParser < m_VParser.end(); itParser++) @@ -307,16 +285,15 @@ void BandMathImageFilterX<TImage> *itParser = ParserType::New(); } - // Generate variables, and register them for each parser (reminder : one per thread) - generateVariables(); m_NbVar = m_VVarName.size(); m_AImage.resize(nbThreads); + double initValue = 1.0; for(int i = 0; i < nbThreads; ++i) { m_AImage[i].resize(m_NbVar); - m_VParser[i]->SetExpr(m_Expression[0]); + //m_VParser[i]->SetExpr(m_Expression[0]); //To be vired for(int j=0; j < m_NbVar; ++j) { @@ -328,7 +305,7 @@ void BandMathImageFilterX<TImage> if ( (m_AImage[i][j].type == 0 ) || (m_AImage[i][j].type == 1) ) // indices (idxX & idxY) { - m_AImage[i][j].value = ValueType(0.0); + m_AImage[i][j].value = ValueType(initValue); } if (m_AImage[i][j].type == 2) //imiPhyX @@ -346,12 +323,12 @@ void BandMathImageFilterX<TImage> if (m_AImage[i][j].type == 4 ) // vector { unsigned int nbBands = this->GetNthInput(m_AImage[i][j].info[0])->GetNumberOfComponentsPerPixel(); - m_AImage[i][j].value = ValueType(1,nbBands,0.0); + m_AImage[i][j].value = ValueType(1,nbBands,initValue); } if (m_AImage[i][j].type == 5 ) // pixel { - m_AImage[i][j].value = ValueType(0.0); + m_AImage[i][j].value = ValueType(initValue); } if (m_AImage[i][j].type == 6 ) // neighborhood @@ -360,9 +337,121 @@ void BandMathImageFilterX<TImage> } 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; } } - + +} + + +template< typename TImage > +void BandMathImageFilterX< TImage > +::OutputsDimensions() +{ + m_outputsDimensions.clear(); + + for(int i=0; i<m_Expression.size(); ++i) + { + m_VParser.at(0)->SetExpr(m_Expression[i]); + ValueType value = m_VParser.at(0)->Eval(); + + switch (value.GetType()) + { //ValueType + case 'i': + m_outputsDimensions.push_back(1); + break; + + case 'f': + m_outputsDimensions.push_back(1); + break; + + case 'c': + itkExceptionMacro(<< "Complex numbers not supported." << std::endl); + break; + + case 'm': + mup::matrix_type vect = value.GetArray(); + m_outputsDimensions.push_back(vect.GetCols()); + break; + + } + } + + /*for(int i=0; i<m_Expression.size(); ++i) + std::cout << "m_outputsDimensions[i] = " << m_outputsDimensions[i] << std::endl;*/ + +} + + + +template< typename TImage > +void BandMathImageFilterX< TImage > +::AllocateOutputs() +{ + typedef itk::ImageBase< TImage::ImageDimension > ImageBaseType; + typename ImageBaseType::Pointer outputPtr; + + GenerateVariables(); + OutputsDimensions(); + + // Allocate the output memory + int i=0; + for ( itk::OutputDataObjectIterator it(this); !it.IsAtEnd(); it++ ) + { + // Check whether the output is an image of the appropriate + // dimension (use ProcessObject's version of the GetInput() + // method since it returns the input as a pointer to a + // DataObject as opposed to the subclass version which + // static_casts the input to an TInputImage). + outputPtr = dynamic_cast< ImageBaseType * >( it.GetOutput() ); + + if ( outputPtr ) + { + outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() ); + outputPtr->SetNumberOfComponentsPerPixel(m_outputsDimensions[i]); + outputPtr->Allocate(); + } + + i++; + } +} + +template< typename TImage > +void BandMathImageFilterX<TImage> +::BeforeThreadedGenerateData() +{ + // Some useful variable + unsigned int nbThreads = this->GetNumberOfThreads(); + unsigned int nbInputImages = this->GetNumberOfInputs(); + + // Check if input image dimensions matches + unsigned int inputSize[2]; + inputSize[0] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(0); + inputSize[1] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(1); + + for(unsigned int p = 1; p < nbInputImages; p++) + { + if((inputSize[0] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0)) + || (inputSize[1] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1))) + { + itkExceptionMacro(<< "Input images must have the same dimensions." << std::endl + << "band #1 is [" << inputSize[0] << ";" << inputSize[1] << "]" << std::endl + << "band #" << p+1 << " is [" + << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0) << ";" + << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1) << "]"); + } + } + + // Allocate and initialize the thread temporaries + m_ThreadUnderflow.SetSize(nbThreads); + m_ThreadUnderflow.Fill(0); + m_ThreadOverflow.SetSize(nbThreads); + m_ThreadOverflow.Fill(0); + + //GenerateVariables(); } @@ -386,7 +475,7 @@ void BandMathImageFilterX<TImage> if((m_UnderflowCount != 0) || (m_OverflowCount!=0)) otbWarningMacro(<< std::endl << "The Following Parsed Expression : " - << this->GetExpression() << std::endl + << this->GetExpression(0) << std::endl //TODO << "Generated " << m_UnderflowCount << " Underflow(s) " << "And " << m_OverflowCount << " Overflow(s) " << std::endl << "The Parsed Expression, The Inputs And The Output " @@ -411,14 +500,24 @@ void BandMathImageFilterX<TImage> Vit[j] = ImageRegionConstIteratorType (this->GetNthInput(j), outputRegionForThread); } - itk::ImageRegionIterator<TImage> outIt (this->GetOutput(), outputRegionForThread); + std::vector< ImageRegionConstIteratorType > VoutIt; + VoutIt.resize(m_Expression.size()); + for(int j=0; j < VoutIt.size(); ++j) + { + VoutIt[j] = ImageRegionConstIteratorType (this->GetOutput(j), outputRegionForThread); +//std::cout << " this->GetOutput(j)->GetNumberOfComponentsPerPixel() = " << this->GetOutput(j)->GetNumberOfComponentsPerPixel() << std::endl; + } + + //itk::ImageRegionIterator<TImage> outIt (this->GetOutput(0), outputRegionForThread); // Support progress methods/callbacks itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); // Pixel affectation for(int j=0; j < nbInputImages; ++j) {Vit[j].GoToBegin();} - outIt.GoToBegin(); + for(int j=0; j < m_Expression.size(); ++j) {VoutIt[j].GoToBegin();} + //outIt.GoToBegin(); + while(!Vit.at(0).IsAtEnd()) // For each pixel { @@ -466,61 +565,66 @@ void BandMathImageFilterX<TImage> } } + for(int IDExpression=0; IDExpression<m_Expression.size(); ++IDExpression) + { - try - { - value = m_VParser[threadId]->Eval(); - } - catch(itk::ExceptionObject& err) - { - itkExceptionMacro(<< err); - } + m_VParser[threadId]->SetExpr(m_Expression[IDExpression]); - switch (value.GetType()) - { //ValueType - case 'i': - outIt.Get()[0] = value.GetInteger(); - break; + try + { + value = m_VParser[threadId]->Eval(); + } + catch(itk::ExceptionObject& err) + { + itkExceptionMacro(<< err); + } - case 'f': - outIt.Get()[0] = value.GetFloat(); - break; + switch (value.GetType()) + { //ValueType + case 'i': + VoutIt[IDExpression].Get()[0] = value.GetInteger(); + break; - case 'c': - itkExceptionMacro(<< "Complex numbers not supported." << std::endl); - break; + case 'f': + VoutIt[IDExpression].Get()[0] = value.GetFloat(); + break; - case 'm': - mup::matrix_type vect = value.GetArray(); - for(int p=0; p<vect.GetCols(); ++p) - outIt.Get()[p] = vect.At(0,p); - break; + case 'c': + itkExceptionMacro(<< "Complex numbers 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; + } - for(int p=0; p<outIt.Get().GetSize(); ++p) - { - // Case value is equal to -inf or inferior to the minimum value - // allowed by the pixelType cast - if (outIt.Get()[p] < double(itk::NumericTraits<typename PixelType::ValueType>::NonpositiveMin())) + for(int p=0; p<VoutIt[IDExpression].Get().GetSize(); ++p) { - outIt.Get()[p] = itk::NumericTraits<typename PixelType::ValueType>::NonpositiveMin(); - m_ThreadUnderflow[threadId]++; + // Case value is equal to -inf or inferior to the minimum value + // allowed by the pixelType cast + if (VoutIt[IDExpression].Get()[p] < double(itk::NumericTraits<typename PixelType::ValueType>::NonpositiveMin())) + { + VoutIt[IDExpression].Get()[p] = itk::NumericTraits<typename PixelType::ValueType>::NonpositiveMin(); + m_ThreadUnderflow[threadId]++; + } + // Case value is equal to inf or superior to the maximum value + // allowed by the pixelType cast + else if (VoutIt[IDExpression].Get()[p] > double(itk::NumericTraits<typename PixelType::ValueType>::max())) + { + VoutIt[IDExpression].Get()[p] = itk::NumericTraits<typename PixelType::ValueType>::max(); + m_ThreadOverflow[threadId]++; + } } - // Case value is equal to inf or superior to the maximum value - // allowed by the pixelType cast - else if (outIt.Get()[p] > double(itk::NumericTraits<typename PixelType::ValueType>::max())) - { - outIt.Get()[p] = itk::NumericTraits<typename PixelType::ValueType>::max(); - m_ThreadOverflow[threadId]++; - } - } + } for(int j=0; j < nbInputImages; ++j) {++Vit[j];} - ++outIt; + for(int j=0; j < m_Expression.size(); ++j) {++VoutIt[j];} + //++outIt; progress.CompletedPixel(); } diff --git a/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx b/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx index 28f588da94575040af9effc97caa8762635d0dd8..336658a2457efe952c820fa273c8e98144bbf47e 100644 --- a/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx +++ b/Testing/Code/BasicFilters/otbBandMathImageFilterX.cxx @@ -104,16 +104,23 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) []) filter->SetNthInput(1, image2); filter->SetNthInput(2, image3, "canal3"); - filter->SetExpression("vcos(2 * pi * im1) div (2 * pi * im2 + {1E-3}) mult vsin(pi * canal3) + ndvi(im1b1, im2b1) * sqrt(2) * canal3"); + filter->SetExpression("vcos(2 * pi * im1) div (2 * pi * im2 + {1E-3}) mult vsin(pi * canal3) + ndvi(im1b1, im2b1) * sqrt(2) * canal3"); //Sub-test 1 + filter->SetExpression("im1b1 / im2b1"); //Sub-test 2 (Edge Effect Handling) filter->Update(); + //if (filter->GetNumberOfOutputs() != 2) + + + ImageType::Pointer output1 = filter->GetOutput(0); + ImageType::Pointer output2 = filter->GetOutput(1); + std::cout << "\n--- Standard Use\n"; - std::cout << "Parsed Expression : " << filter->GetExpression() << std::endl; + std::cout << "Parsed Expression : " << filter->GetExpression(0) << std::endl; - ImageType::Pointer output = filter->GetOutput(); - IteratorType it(output, region); + //Sub-test 1 + IteratorType itoutput1(output1, region); - for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(), it.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3, ++it) + 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(); @@ -124,7 +131,7 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) []) px2[0] = ( i2[0] * i2[1] ); px3[0] = ( i3[0] + i3[1] * i3[1] ); - double result = it.Get()[0]; + double result = itoutput1.Get()[0]; double ndvi_expected; double error; @@ -138,7 +145,7 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) []) /*std::cout << "Pixel_1 = " << it1.Get()[0] << " Pixel_2 = " << it2.Get()[0] << " Pixel_3 = " << it3.Get()[0] - << " Result = " << it.Get()[0] << " Expected = " << expected << std::endl;*/ + << " Result = " << itoutput1.Get()[0] << " Expected = " << expected << std::endl;*/ error = (result - expected) * (result - expected) / (result + expected); @@ -149,7 +156,7 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) []) << "Pixel_1 = " << it1.Get()[0] << " Pixel_2 = " << it2.Get()[0] << " Pixel_3 = " << it3.Get()[0] - << " Result = " << it.Get()[0] + << " Result = " << itoutput1.Get()[0] << " Expected = " << expected << std::endl ); FAIL_FLAG++; } @@ -161,23 +168,25 @@ int otbBandMathImageFilterX( int itkNotUsed(argc), char* itkNotUsed(argv) []) FAIL_FLAG = 0; - + //Sub-test 2 /** Edge Effect Handling */ + + IteratorType itoutput2(output2, region); std::cout << "\n--- Edge Effect Handling\n"; std::cout << "- +/-inf section\n"; - filter->SetExpression("im1b1 / im2b1"); - filter->Update(); + //filter->SetExpression("im1b1 / im2b1"); + //filter->Update(); std::cout << "- nan section\n"; - it1.GoToBegin(); it2.GoToBegin(); it.GoToBegin(); - for(i=1; i<=50; ++i , ++it1, ++it2, ++it){} - if(vnl_math_isnan(it.Get()[0])) - std::cout << "Pixel_1 = " << it1.Get() << " Pixel_2 = " << it2.Get() << " Result = " << it.Get() << " Expected = nan\n"; + it1.GoToBegin(); it2.GoToBegin(); itoutput2.GoToBegin(); + for(i=1; i<=50; ++i , ++it1, ++it2, ++itoutput2){} + if(vnl_math_isnan(itoutput2.Get()[0])) + std::cout << "Pixel_1 = " << it1.Get() << " Pixel_2 = " << it2.Get() << " Result = " << itoutput2.Get() << " Expected = nan\n"; else itkGenericExceptionMacro( << "\nError > Bad Edge Effect Handling -> Test Failled\n" << "Pixel_1 = " << it1.Get() << " Pixel_2 = " << it2.Get() - << " Result = " << it.Get() << " Expected = nan\n" ); + << " Result = " << itoutput2.Get() << " Expected = nan\n" ); std::cout << std::endl; return EXIT_SUCCESS; @@ -255,7 +264,7 @@ int otbBandMathImageFilterXBatchMode( int itkNotUsed(argc), char* itkNotUsed(arg filter->Update(); std::cout << "\n--- Standard Use\n"; - std::cout << "Parsed Expression : " << filter->GetExpression() << std::endl; + std::cout << "Parsed Expression : " << filter->GetExpression(0) << std::endl; ImageType::Pointer output = filter->GetOutput(); IteratorType it(output, region); @@ -400,25 +409,17 @@ int otbBandMathImageFilterXWithIdx( int itkNotUsed(argc), char* argv[]) filter->SetNthInput(1, image2); filter->SetNthInput(2, image3); - // #ifdef OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS filter->SetExpression("(sqrt(idxX*idxX+idxY*idxY) < 50) ? im1b1 : im2b1"); - //#else - //filter->SetExpression("if(sqrt(idxX*idxX+idxY*idxY) < 50, b2, b3)"); - //#endif WriterType::Pointer writer = WriterType::New(); writer->SetInput(filter->GetOutput()); writer->SetFileName(outfname1); writer->Update(); - //#ifdef OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS - filter->SetExpression("(sqrt(im2PhyX*im2PhyX+im2PhyY*im2PhyY) < 25) ? im2b1 : im3b1"); - /*#else - filter->SetExpression("if(sqrt(idxPhyX*idxPhyX+idxPhyY*idxPhyY) < 25, b2, b3)"); - #endif*/ + /*filter->SetExpression("(sqrt(im2PhyX*im2PhyX+im2PhyY*im2PhyY) < 25) ? im2b1 : im3b1"); WriterType::Pointer writer2 = WriterType::New(); writer2->SetInput(filter->GetOutput()); writer2->SetFileName(outfname2); - writer2->Update(); + writer2->Update();*/ return EXIT_SUCCESS; }