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/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;
+}