diff --git a/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.h b/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..569ce70428e1326a1875c8b03d85c3943406cdaf --- /dev/null +++ b/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.h @@ -0,0 +1,218 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + Some parts of this code are derived from ITK. See ITKCopyright.txt + for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANT2ABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __otbStreamingInnerProductVectorImageFilter_h +#define __otbStreamingInnerProductVectorImageFilter_h + + +#include "otbPersistentImageFilter.h" +#include "otbPersistentFilterStreamingDecorator.h" +#include "itkNumericTraits.h" +#include "itkArray.h" +#include "itkSimpleDataObjectDecorator.h" +#include "itkImageRegionSplitter.h" +#include "itkVariableSizeMatrix.h" +#include "itkVariableLengthVector.h" +#include "vnl/vnl_matrix.h" + + +namespace otb +{ + +/** \class StreamingInnerProductVectorImageFilter + * \brief Compute the inner product of a large image using streaming + * + * This filter persists its temporary data. It means that if you Update it n times on n different + * requested regions, the output statistics will be the statitics of the whole set of n regions. + * + * To reset the temporary data, one should call the Reset() function. + * + * To get the statistics once the regions have been processed via the pipeline, use the Synthetize() method. + * + * \sa PersistentImageFilter + * \ingroup Streamed + * \ingroup Multithreaded + * \ingroup MathematicalStatisticsImageFilters + * + */ +template<class TInputImage> +class ITK_EXPORT PersistentInnerProductVectorImageFilter : + public PersistentImageFilter<TInputImage, TInputImage> +{ +public: + /** Standard Self typedef */ + typedef PersistentInnerProductVectorImageFilter Self; + typedef PersistentImageFilter<TInputImage,TInputImage> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(PersistentInnerProductVectorImageFilter,PersistentImageFilter); + + /** Image related typedefs. */ + typedef TInputImage ImageType; + typedef typename TInputImage::Pointer InputImagePointer; + typedef typename TInputImage::RegionType RegionType; + typedef typename TInputImage::PixelType PixelType; + + itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); + + /** Image related typedefs. */ + itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension ); + + /** Smart Pointer type to a DataObject. */ + typedef typename itk::DataObject::Pointer DataObjectPointer; + + /** Type definition for a double matrix. */ + typedef vnl_matrix<double> MatrixType; + typedef typename std::vector<MatrixType> ArrayMatrixType; + + /** Type of DataObjects used for scalar outputs */ + typedef itk::SimpleDataObjectDecorator<MatrixType> MatrixObjectType; + + /** Return the computed inner product matrix. */ + + MatrixType GetInnerProduct() const + { + return this->GetInnerProductOutput()->Get(); + }; + MatrixObjectType* GetInnerProductOutput(); + const MatrixObjectType* GetInnerProductOutput() const; + + + /** Make a DataObject of the correct type to be used as the specified + * output. + */ + virtual DataObjectPointer MakeOutput(unsigned int idx); + + /** Pass the input through unmodified. Do this by Grafting in the + * AllocateOutputs method. + */ + virtual void AllocateOutputs(); + virtual void GenerateOutputInformation(); + virtual void Synthetize(void); + virtual void Reset(void); + +protected: + PersistentInnerProductVectorImageFilter(); + virtual ~PersistentInnerProductVectorImageFilter() {}; + virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; + /** Multi-thread version GenerateData. */ + void ThreadedGenerateData (const RegionType& outputRegionForThread,int threadId); + +private: + PersistentInnerProductVectorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + ArrayMatrixType m_ThreadInnerProduct; + +}; // end of class PersistentStatisticsVectorImageFilter + +/**===========================================================================*/ + +/** \class StreamingInnerProductVectorImageFilter + * \brief This class streams the whole input image through the PersistentStatisticsImageFilter. + * + * This way, it allows to compute the inner product of this image. It calls the + * Reset() method of the PersistentStatisticsImageFilter before streaming the image and the + * Synthetize() method of the PersistentStatisticsImageFilter after having streamed the image + * to compute the statistics. The accessor on the results are wrapping the accessors of the + * internal PersistentStatisticsImageFilter. + * + * \sa PersistentStatisticsVectorImageFilter + * \sa PersistentImageFilter + * \sa PersistentFilterStreamingDecorator + * \sa StreamingImageVirtualWriter + * \ingroup Streamed + * \ingroup Multithreaded + * \ingroup MathematicalStatisticsImageFilters + */ + +template<class TInputImage> +class ITK_EXPORT StreamingInnerProductVectorImageFilter : + public PersistentFilterStreamingDecorator< PersistentInnerProductVectorImageFilter<TInputImage> > +{ +public: + /** Standard Self typedef */ + typedef StreamingInnerProductVectorImageFilter Self; + typedef PersistentFilterStreamingDecorator + < PersistentInnerProductVectorImageFilter<TInputImage> > Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + + /** Type macro */ + itkNewMacro(Self); + + /** Creation through object factory macro */ + itkTypeMacro(StreamingInnerProductVectorImageFilter,PersistentFilterStreamingDecorator); + + typedef TInputImage InputImageType; + typedef typename Superclass::FilterType StatFilterType; + typedef typename StatFilterType::MatrixType MatrixType; + + + /** Type of DataObjects used for scalar outputs */ + typedef typename StatFilterType::MatrixObjectType MatrixObjectType; + + void SetInput(TInputImage * input) + { + this->GetFilter()->SetInput(input); + } + TInputImage * GetInput() + { + return this->GetFilter()->GetInput(); + } + + /** Return the computed inner product. */ + MatrixType GetInnerProduct() const + { + return this->GetFilter()->GetInnerProductOutput()->Get(); + }; + MatrixObjectType* GetInnerProductOutput() + { + return this->GetFilter()->GetInnerProductOutput(); + }; + const MatrixObjectType* GetInnerProductOutput() const + { + return this->GetFilter()->GetInnerProductOutput(); + }; + +protected: + /** Constructor */ + StreamingInnerProductVectorImageFilter() {}; + /** Destructor */ + virtual ~StreamingInnerProductVectorImageFilter() {}; + +private: + StreamingInnerProductVectorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented +}; + +} // end namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbStreamingInnerProductVectorImageFilter.txx" +#endif + +#endif diff --git a/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.txx b/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..83918751951b6f386f3a2abf50edef1cd503c5dc --- /dev/null +++ b/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.txx @@ -0,0 +1,243 @@ +/*========================================================================= + +Program: ORFEO Toolbox +Language: C++ +Date: $Date$ +Version: $Revision$ + + +Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. +See OTBCopyright.txt for details. + +Some parts of this code are derived from ITK. See ITKCopyright.txt +for details. + + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __otbStreamingInnerProductVectorImageFilter_txx +#define __otbStreamingInnerProductVectorImageFilter_txx +#include "otbStreamingInnerProductVectorImageFilter.h" + +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkNumericTraits.h" +#include "itkProgressReporter.h" +#include "otbMacro.h" + +namespace otb +{ + +template<class TInputImage> +PersistentInnerProductVectorImageFilter<TInputImage> +::PersistentInnerProductVectorImageFilter() +{ + // first output is a copy of the image, DataObject created by + // superclass + // + // allocate the data objects for the outputs which are + // just decorators around pixel types + + // allocate the data objects for the outputs which are + // just decorators around matrix types + typename ImageType::Pointer output1 = static_cast<ImageType*>(this->MakeOutput(0).GetPointer()); + this->itk::ProcessObject::SetNthOutput(0, output1.GetPointer()); + typename MatrixObjectType::Pointer output2 = static_cast<MatrixObjectType*>(this->MakeOutput(1).GetPointer()); + this->itk::ProcessObject::SetNthOutput(1, output2.GetPointer()); +} + +template<class TInputImage> +itk::DataObject::Pointer +PersistentInnerProductVectorImageFilter<TInputImage> +::MakeOutput(unsigned int output) +{ + switch (output) + { + case 0: + return static_cast<itk::DataObject*>(TInputImage::New().GetPointer()); + break; + case 1: + return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer()); + break; + default: + // might as well make an image + return static_cast<itk::DataObject*>(TInputImage::New().GetPointer()); + break; + } +} + +template<class TInputImage> +typename PersistentInnerProductVectorImageFilter<TInputImage>::MatrixObjectType* +PersistentInnerProductVectorImageFilter<TInputImage> +::GetInnerProductOutput() +{ + return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1)); +} + +template<class TInputImage> +const typename PersistentInnerProductVectorImageFilter<TInputImage>::MatrixObjectType* +PersistentInnerProductVectorImageFilter<TInputImage> +::GetInnerProductOutput() const +{ + return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1)); +} + + +template<class TInputImage> +void +PersistentInnerProductVectorImageFilter<TInputImage> +::GenerateOutputInformation() +{ + Superclass::GenerateOutputInformation(); + if ( this->GetInput() ) + { + this->GetOutput()->CopyInformation(this->GetInput()); + this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion()); + + if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()==0) + { + this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion()); + } + } +} + +template<class TInputImage> +void +PersistentInnerProductVectorImageFilter<TInputImage> +::AllocateOutputs() +{ + // This is commented to prevent the streaming of the whole image for the first stream strip + // It shall not cause any problem because the output image of this filter is not intended to be used. + //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() ); + //this->GraftOutput( image ); + // Nothing that needs to be allocated for the remaining outputs +} + +template<class TInputImage> +void +PersistentInnerProductVectorImageFilter<TInputImage> +::Reset() +{ + TInputImage * inputPtr = const_cast<TInputImage * >(this->GetInput()); + inputPtr->UpdateOutputInformation(); + + if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()==0) + { + this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion()); + } + + unsigned int numberOfThreads = this->GetNumberOfThreads(); + unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel(); + // Set the number of training image + MatrixType tempMatrix; + tempMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages); + tempMatrix.fill(0); + m_ThreadInnerProduct = ArrayMatrixType(numberOfThreads, tempMatrix); + + MatrixType initMatrix; + initMatrix.set_size(numberOfTrainingImages, numberOfTrainingImages); + initMatrix.fill(0); + this->GetInnerProductOutput()->Set( initMatrix ); + +} + +template<class TInputImage> +void +PersistentInnerProductVectorImageFilter<TInputImage> +::Synthetize() +{ + // Compute Inner product Matrix + TInputImage * inputPtr = const_cast<TInputImage * >(this->GetInput()); + unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel(); + unsigned int numberOfThreads = this->GetNumberOfThreads(); + MatrixType innerProduct; + innerProduct.set_size( numberOfTrainingImages, numberOfTrainingImages ); + innerProduct.fill( 0 ); + + // Concatenate threaded matrix + for ( unsigned int thread = 0; thread < numberOfThreads; thread++) + { + innerProduct += m_ThreadInnerProduct[thread]; + } + + //--------------------------------------------------------------------- + // Fill the rest of the inner product matrix and make it symmetric + //--------------------------------------------------------------------- + for(unsigned int band_x = 0; band_x < (numberOfTrainingImages - 1); band_x++) + { + for(unsigned int band_y = band_x+1; band_y < numberOfTrainingImages; band_y++) + { + innerProduct[band_x][band_y] = innerProduct[band_y][band_x]; + }// end band_y loop + }// end band_x loop + + if( ( numberOfTrainingImages - 1 ) != 0 ) + { + innerProduct /= ( numberOfTrainingImages - 1 ); + } + else + { + innerProduct.fill(0); + } + + // Set the output + this->GetInnerProductOutput()->Set( innerProduct ); +} + +template<class TInputImage> +void +PersistentInnerProductVectorImageFilter<TInputImage> +::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId) +{ + /** + * Grab the input + */ + InputImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput() ); + // support progress methods/callbacks + itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); + unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel(); + + itk::ImageRegionConstIterator<TInputImage> it (inputPtr, outputRegionForThread); + it.GoToBegin(); + // do the work + while (!it.IsAtEnd()) + { + PixelType vectorValue = it.Get(); + double mean(0.); + for (unsigned int i=0; i<vectorValue.GetSize(); i++) + { + mean += static_cast<double>(vectorValue[i]); + } + mean /= static_cast<double>(vectorValue.GetSize()); + + // Matrix iteration + for(unsigned int band_x = 0; band_x < numberOfTrainingImages; band_x++) + { + for(unsigned int band_y = 0; band_y <= band_x; band_y++ ) + { + m_ThreadInnerProduct[threadId][band_x][band_y] += + (static_cast<double>(vectorValue[band_x]) - mean) * (static_cast<double>(vectorValue[band_y]) - mean); + } // end: band_y loop + } // end: band_x loop + ++it; + progress.CompletedPixel(); + }// end: looping through the image +} + +template <class TImage> +void +PersistentInnerProductVectorImageFilter<TImage> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os,indent); + + os << indent << "InnerProduct: " << this->GetInnerProductOutput()->Get() << std::endl; + +} + + +}// end namespace otb +#endif diff --git a/Code/Radiometry/otbVegetationIndex.h b/Code/Radiometry/otbVegetationIndex.h index ba2763e05680eac2a46fa74939c6a604112faee5..68d268a4a7a4dac65fe384b028c284deb4fb07d5 100644 --- a/Code/Radiometry/otbVegetationIndex.h +++ b/Code/Radiometry/otbVegetationIndex.h @@ -295,6 +295,84 @@ private: double m_Gamma; }; + + + +/** \class TSARVI + * \brief This functor calculate the Transformed Soil Atmospherical Resistant Vegetation Index (TSARVI) + * + * [Yoram J. Kaufman and Didier Tanré, 1992] + * + * \ingroup Functor + */ +template <class TInput1, class TInput2, class TInput3, class TOutput> +class TSARVI +{ +public: + TSARVI() : m_X(0.08), m_Gamma(0.5) {}; + ~TSARVI() {}; + inline TOutput operator()(const TInput1 &r, const TInput2 &b, const TInput3 &nir) + { + double dr = static_cast<double>(r); + double db = static_cast<double>(b); + double dnir = static_cast<double>(nir); + double dRB = dr - m_Gamma*(db - dr); + double denominator = dRB + m_A*dnir - m_A*m_B + m_X*(1.+m_A*m_A); + if ( denominator == 0. ) + { + return static_cast<TOutput>(0.); + } + return ( static_cast<TOutput>( (m_A*(dnir - m_A*dRB - m_B))/denominator ) ); + } + /** Set/Get A and B parameters */ + void SetA(const double A) + { + m_A = A; + } + double GetA(void)const + { + return (m_A); + } + void SetB(const double B) + { + m_B = B; + } + double GetB(void)const + { + return (m_B); + } + /** Set/Get X parameter */ + void SetX(const double X) + { + m_X = X; + } + double GetX(void)const + { + return (m_X); + } + /** Set/Get the gamma parameter */ + void SetGamma(const double gamma) + { + m_Gamma = gamma; + } + double GetGamma(void)const + { + return (m_Gamma); + } + +private: + + /** A and B parameters */ + double m_A; + double m_B; + /** X parameter */ + double m_X; + /** Gamma parameter */ + double m_Gamma; + +}; + + /** \class EVI * \brief This functor calculate the Enhanced Vegetation Index (EVI) * @@ -321,7 +399,6 @@ public: return static_cast<TOutput>(0.); } return ( static_cast<TOutput>( m_G * (dnir - dr)/denominator ) ); -//return ( static_cast<TOutput>( dnir ) ); } /** Set/Get G parameter */ void SetG(const double g) diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt index 96c31626dddff3ede128eedc3a0cf91fa60aa267..3031a6dee25ce71170baa868d9765f34dccc020c 100644 --- a/Testing/Code/BasicFilters/CMakeLists.txt +++ b/Testing/Code/BasicFilters/CMakeLists.txt @@ -1184,6 +1184,17 @@ ADD_TEST(bfTvFunctionWithNeighborhoodToImageFilterNew ${BASICFILTERS_TESTS11} otbFunctionWithNeighborhoodToImageFilterNew ) +ADD_TEST(bfTuStreamingInnerProductVectorImageFilterNew ${BASICFILTERS_TESTS11} + otbStreamingInnerProductVectorImageFilterNew +) + +ADD_TEST(bfTvStreamingInnerProductVectorImageFilter ${BASICFILTERS_TESTS11} + otbStreamingInnerProductVectorImageFilter + ${INPUTDATA}/poupees_sub.png + ${TEMP}/bfTvStreamingInnerProductVectorImageFilterResults.txt +) + + # A enrichir SET(BasicFilters_SRCS1 @@ -1354,6 +1365,8 @@ otbKeyPointDensityImageFilterNew.cxx otbKeyPointDensityImageFilterTest.cxx otbImagePCAShapeModelEstimatorTest.cxx otbFunctionWithNeighborhoodToImageFilterNew.cxx +otbStreamingInnerProductVectorImageFilterNew.cxx +otbStreamingInnerProductVectorImageFilter.cxx ) diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests11.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests11.cxx index 94a81b14a92118c134641c662b5c78fb0f2146cf..a5670072ea8cacad484b7f75ad310ba32ef5dea3 100644 --- a/Testing/Code/BasicFilters/otbBasicFiltersTests11.cxx +++ b/Testing/Code/BasicFilters/otbBasicFiltersTests11.cxx @@ -37,4 +37,6 @@ void RegisterTests() REGISTER_TEST(otbKeyPointDensityImageFilterTest); REGISTER_TEST(otbImagePCAShapeModelEstimatorTest); REGISTER_TEST(otbFunctionWithNeighborhoodToImageFilterNew); + REGISTER_TEST(otbStreamingInnerProductVectorImageFilterNew); + REGISTER_TEST(otbStreamingInnerProductVectorImageFilter); } diff --git a/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..ff94bec5753f58cfb16693e9b4cf8ea231ec7eae --- /dev/null +++ b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx @@ -0,0 +1,56 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbStreamingInnerProductVectorImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include <fstream> + +int otbStreamingInnerProductVectorImageFilter( int argc, char* argv[] ) +{ + const char * inputFileName = argv[1]; + const char * outfname = argv[2]; + + typedef double PixelType; + const unsigned int Dimension = 2; + + // Typedef + typedef otb::VectorImage<PixelType,Dimension> ImageType; + typedef otb::ImageFileReader< ImageType > ReaderType; + typedef otb::StreamingInnerProductVectorImageFilter<ImageType> FilterType; + + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(inputFileName); + + // Instanciation object + FilterType::Pointer filter = FilterType::New(); + + filter->GetStreamer()->SetNumberOfStreamDivisions(10); + filter->SetInput(reader->GetOutput()); + filter->Update(); + + + std::ofstream file; + file.open(outfname); + file.precision(5); + file<<std::fixed; + file<<"Inner Product: Dim ["<<filter->GetInnerProduct().size()<<"]:"<<std::endl; + file<<filter->GetInnerProduct()<<std::endl; + file.close(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilterNew.cxx b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0b7c15c002d3d6f785ca2c79b3a392993030df6b --- /dev/null +++ b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilterNew.cxx @@ -0,0 +1,34 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbStreamingInnerProductVectorImageFilter.h" +#include "otbVectorImage.h" + +int otbStreamingInnerProductVectorImageFilterNew( int argc, char* argv[] ) +{ + typedef double PixelType; + const unsigned int Dimension = 2; + + // Typedef + typedef otb::VectorImage<PixelType,Dimension> ImageType; + typedef otb::StreamingInnerProductVectorImageFilter<ImageType> FilterType; + + // Instanciation object + FilterType::Pointer filter = FilterType::New(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/CMakeLists.txt b/Testing/Code/Radiometry/CMakeLists.txt index c05f18b6ae2c2870a69f3a4dbab5cef3136fc1c5..82d1ff611b1c8573073ca8f515b3afa04b58c148 100644 --- a/Testing/Code/Radiometry/CMakeLists.txt +++ b/Testing/Code/Radiometry/CMakeLists.txt @@ -181,7 +181,7 @@ ADD_TEST(raTvARVI_MultiChannelRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY ${INPUTDATA}/poupees_sub.png ${TEMP}/raMultiChannelRAndBAndNIRVegetationIndex_ARVI_poupees_subc1c2c3.tif 1 2 3 - 0.6 # Gamma parameter + 0.6 # Gamma parameter ) # ------- otb::ImageToLuminanceImageFilter ------------------------------ @@ -487,7 +487,7 @@ ADD_TEST(raTvAtmosphericCorrectionSequencementTest ${RADIOMETRY_TESTS4} ) -# ------- otb::RAndBAndNIRVegetationIndexImageFilter ------------------------------ +# ------- EVI RAndBAndNIRVegetationIndexImageFilter ------------------------------ ADD_TEST(raTvEVIRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_TESTS4} --compare-image ${EPSILON} ${BASELINE}/raRAndBAndNIRVegetationIndex_EVI_poupees_subc1c2c3.tif ${TEMP}/raRAndBAndNIRVegetationIndex_EVI_poupees_subc1c2c3.tif @@ -503,7 +503,7 @@ ADD_TEST(raTvEVIRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_TESTS4} 1.0 ) -# ------- otb::MultiChannelRAndBAndNIRVegetationIndexImageFilter ------------------------------ +# ------- EVI MultiChannelRAndBAndNIRVegetationIndexImageFilter ------------------------------ ADD_TEST(raTvEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_TESTS4} --compare-image ${EPSILON} ${BASELINE}/raRAndBAndNIRVegetationIndex_EVI_qb_RoadExtract.tif ${TEMP}/raRAndBAndNIRVegetationIndex_EVI_qb_RoadExtract.tif @@ -520,6 +520,37 @@ ADD_TEST(raTvEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_T 400.0 # canopy background adjustment ) +# ------- TSARVI RAndBAndNIRVegetationIndexImageFilter ------------------------------ +ADD_TEST(raTvTSARVIRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_TESTS4} + --compare-image ${EPSILON} ${BASELINE}/raRAndBAndNIRVegetationIndex_TSARVI_poupees_subc1c2c3.tif + ${TEMP}/raRAndBAndNIRVegetationIndex_TSARVI_poupees_subc1c2c3.tif + otbTSARVIRAndBAndNIRVegetationIndexImageFilter + ${INPUTDATA}/poupees_sub_c1.png + ${INPUTDATA}/poupees_sub_c2.png + ${INPUTDATA}/poupees_sub_c3.png + ${TEMP}/raRAndBAndNIRVegetationIndex_TSARVI_poupees_subc1c2c3.tif + 0.7 # a ( pente de la droite des sols nus dans l'espace RB/PIR ) + 0.9 # b ( ordonnée à l'origine de la droite des sols nus dans l'espace RB/PIR ) + 0.08 # x coeff a priori constant + 0.5 # gamma +) + +# ------- TSARVI MultiChannelRAndBAndNIRVegetationIndexImageFilter ------------------------------ +ADD_TEST(raTvTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter ${RADIOMETRY_TESTS4} + --compare-image ${EPSILON} ${BASELINE}/raRAndBAndNIRVegetationIndex_TSARVI_qb_RoadExtract.tif + ${TEMP}/raRAndBAndNIRVegetationIndex_TSARVI_qb_RoadExtract.tif + otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter + ${INPUTDATA}/qb_RoadExtract.img.hdr + ${TEMP}/raRAndBAndNIRVegetationIndex_TSARVI_qb_RoadExtract.tif + 3 # red + 1 # blue + 4 # nir + 0.7 # a ( pente de la droite des sols nus dans l'espace RB/PIR ) + 0.9 # b ( ordonnée à l'origine de la droite des sols nus dans l'espace RB/PIR ) + 0.08 # x coeff a priori constant + 0.5 # gamma +) + @@ -563,6 +594,8 @@ otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx otbAtmosphericCorrectionSequencement.cxx otbEVIRAndBAndNIRVegetationIndexImageFilter.cxx otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx +otbTSARVIRAndBAndNIRVegetationIndexImageFilter.cxx +otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx ) INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}") diff --git a/Testing/Code/Radiometry/otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx b/Testing/Code/Radiometry/otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx index c8d18acff476275c6591a14acb1d44c01aa6de67..bd3d6503c45b9fa0d68c939a7d0857ce6c5593f9 100644 --- a/Testing/Code/Radiometry/otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx +++ b/Testing/Code/Radiometry/otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx @@ -45,7 +45,6 @@ int generic_EVIMultiChannelRAndBAndNIRVegetationIndexImageFilter(int argc, char unsigned int redChannel(::atoi(argv[3])); unsigned int blueChannel(::atoi(argv[4])); unsigned int nirChannel(::atoi(argv[5])); -std::cout << "ORDER : "<<redChannel<<" "<<blueChannel<<" "<<nirChannel<<std::endl; double g(::atof(argv[6])); double c1(::atof(argv[7])); diff --git a/Testing/Code/Radiometry/otbRadiometryTests4.cxx b/Testing/Code/Radiometry/otbRadiometryTests4.cxx index c51fb1b5a29a17e11848460b77f322886111d710..406d589c40f86d222a6efb829b852805a3b85235 100644 --- a/Testing/Code/Radiometry/otbRadiometryTests4.cxx +++ b/Testing/Code/Radiometry/otbRadiometryTests4.cxx @@ -32,5 +32,7 @@ void RegisterTests() REGISTER_TEST(otbAtmosphericCorrectionSequencementTest); REGISTER_TEST(otbEVIRAndBAndNIRVegetationIndexImageFilter); REGISTER_TEST(otbEVIMultiChannelRAndBAndNIRVegetationIndexImageFilter); + REGISTER_TEST(otbTSARVIRAndBAndNIRVegetationIndexImageFilter); + REGISTER_TEST(otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter); } diff --git a/Testing/Code/Radiometry/otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx b/Testing/Code/Radiometry/otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7785d4f83890ba4e4bb5a97b9f635051687ffc49 --- /dev/null +++ b/Testing/Code/Radiometry/otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx @@ -0,0 +1,79 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "itkExceptionObject.h" + +#include "otbMultiChannelRAndBAndNIRVegetationIndexImageFilter.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbVegetationIndex.h" + + +int otbTSARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef otb::VectorImage<double ,Dimension> InputImageType; + typedef otb::Image<double,Dimension> OutputImageType; + typedef otb::ImageFileReader<InputImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + typedef otb::Functor::TSARVI< InputImageType::InternalPixelType, + InputImageType::InternalPixelType, + InputImageType::InternalPixelType, + OutputImageType::PixelType > FunctorType; + typedef otb::MultiChannelRAndBAndNIRVegetationIndexImageFilter<InputImageType,OutputImageType,FunctorType> + MultiChannelRAndBAndNIRVegetationIndexImageFilterType; + + // Instantiating object + MultiChannelRAndBAndNIRVegetationIndexImageFilterType::Pointer filter = MultiChannelRAndBAndNIRVegetationIndexImageFilterType::New(); + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + + const char * inputFilename = argv[1]; + const char * outputFilename = argv[2]; + + unsigned int redChannel(::atoi(argv[3])); + unsigned int blueChannel(::atoi(argv[4])); + unsigned int nirChannel(::atoi(argv[5])); + + double a(::atof(argv[6])); + double b(::atof(argv[7])); + double x(::atof(argv[8])); + double gamma(::atof(argv[9])); + + reader->SetFileName( inputFilename ); + writer->SetFileName( outputFilename ); + filter->SetRedIndex(redChannel); + filter->SetBlueIndex(blueChannel); + filter->SetNIRIndex(nirChannel); + filter->SetInput( reader->GetOutput() ); + + filter->GetFunctor().SetA(a); + filter->GetFunctor().SetB(b); + filter->GetFunctor().SetX(x); + filter->GetFunctor().SetGamma(gamma); + + writer->SetInput( filter->GetOutput() ); + writer->Update(); + + return EXIT_SUCCESS; +} + + + + diff --git a/Testing/Code/Radiometry/otbTSARVIRAndBAndNIRVegetationIndexImageFilter.cxx b/Testing/Code/Radiometry/otbTSARVIRAndBAndNIRVegetationIndexImageFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..66badd3a8850bba72e67919b25d2c7481f8c5260 --- /dev/null +++ b/Testing/Code/Radiometry/otbTSARVIRAndBAndNIRVegetationIndexImageFilter.cxx @@ -0,0 +1,89 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "itkExceptionObject.h" + +#include "otbRAndBAndNIRVegetationIndexImageFilter.h" +#include "otbImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbVegetationIndex.h" + + +int otbTSARVIRAndBAndNIRVegetationIndexImageFilter(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::Image<PixelType,Dimension> InputRImageType; + typedef otb::Image<PixelType,Dimension> InputBImageType; + typedef otb::Image<PixelType,Dimension> InputNIRImageType; + typedef otb::Image<double,Dimension> OutputImageType; + + typedef otb::ImageFileReader<InputRImageType> RReaderType; + typedef otb::ImageFileReader<InputBImageType> BReaderType; + typedef otb::ImageFileReader<InputNIRImageType> NIRReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + + typedef otb::Functor::TSARVI< InputRImageType::PixelType, + InputBImageType::PixelType, + InputNIRImageType::PixelType, + OutputImageType::PixelType > FunctorType; + + typedef otb::RAndBAndNIRVegetationIndexImageFilter< InputRImageType, + InputBImageType, + InputNIRImageType, + OutputImageType, + FunctorType > RAndBAndNIRVegetationIndexImageFilterType; + + // Instantiating object + RAndBAndNIRVegetationIndexImageFilterType::Pointer filter = RAndBAndNIRVegetationIndexImageFilterType::New(); + RReaderType::Pointer readerR = RReaderType::New(); + BReaderType::Pointer readerB = BReaderType::New(); + NIRReaderType::Pointer readerNIR = NIRReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + + + const char * inputFilenameR = argv[1]; + const char * inputFilenameB = argv[2]; + const char * inputFilenameNIR = argv[3]; + const char * outputFilename = argv[4]; + + double a(::atof(argv[5])); + double b(::atof(argv[6])); + double x(::atof(argv[7])); + double gamma(::atof(argv[8])); + + readerR->SetFileName( inputFilenameR ); + readerB->SetFileName( inputFilenameB ); + readerNIR->SetFileName( inputFilenameNIR ); + writer->SetFileName( outputFilename ); + filter->SetInputR( readerR->GetOutput() ); + filter->SetInputB( readerB->GetOutput() ); + filter->SetInputNIR( readerNIR->GetOutput() ); + + filter->GetFunctor().SetA(a); + filter->GetFunctor().SetB(b); + filter->GetFunctor().SetX(x); + filter->GetFunctor().SetGamma(gamma); + + writer->SetInput( filter->GetOutput() ); + writer->Update(); + + + return EXIT_SUCCESS; +} +