From 6359020bd1256a5847c70e714e0c5a70bd3bea2b Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Fri, 5 Oct 2007 15:26:45 +0000
Subject: [PATCH] =?UTF-8?q?Cr=C3=A9ation=20d'une=20nouvelle=20th=C3=A9mati?=
=?UTF-8?q?que=20:=20Fusion.=20Cr=C3=A9ation=20de=20la=20classe=20de=20bas?=
=?UTF-8?q?e=20FusionBase.=20Cr=C3=A9ation=20de=20la=20classe=20de=20fusio?=
=?UTF-8?q?n=20BayesianFusionImageFilter=20(contribution=20Julien=20Radoux?=
=?UTF-8?q?)=20avec=20les=20classes=20MatrixTransposeMatrix=20et=20Streami?=
=?UTF-8?q?ngStatisticsVectorImageFilter=20dans=20BasicFilter.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../otbMatrixTransposeMatrixImageFilter.h | 214 +++++++
.../otbMatrixTransposeMatrixImageFilter.txx | 374 ++++++++++++
.../otbStreamingStatisticsVectorImageFilter.h | 198 +++++++
...tbStreamingStatisticsVectorImageFilter.txx | 554 ++++++++++++++++++
Code/CMakeLists.txt | 1 +
Code/Fusion/CMakeLists.txt | 18 +
Code/Fusion/foo.cxx | 0
Code/Fusion/otbBayesianFusionFilter.h | 278 +++++++++
Code/Fusion/otbBayesianFusionFilter.txx | 277 +++++++++
Code/Fusion/otbFusionImageBase.h | 119 ++++
Testing/Code/BasicFilters/CMakeLists.txt | 39 ++
.../BasicFilters/otbBasicFiltersTests.cxx | 4 +
.../otbMatrixTransposeMatrixImageFilter.cxx | 86 +++
...otbMatrixTransposeMatrixImageFilterNew.cxx | 59 ++
...tbStreamingStatisticsVectorImageFilter.cxx | 77 +++
...treamingStatisticsVectorImageFilterNew.cxx | 50 ++
Testing/Code/CMakeLists.txt | 1 +
Testing/Code/Fusion/CMakeLists.txt | 62 ++
.../Code/Fusion/otbBayesianFusionFilter.cxx | 83 +++
.../Fusion/otbBayesianFusionFilterNew.cxx | 63 ++
Testing/Code/Fusion/otbFusionImageBaseNew.cxx | 86 +++
Testing/Code/Fusion/otbFusionTests.cxx | 34 ++
otbIncludeDirectories.cmake | 1 +
23 files changed, 2678 insertions(+)
create mode 100644 Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.h
create mode 100644 Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.txx
create mode 100644 Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
create mode 100644 Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
create mode 100755 Code/Fusion/CMakeLists.txt
create mode 100644 Code/Fusion/foo.cxx
create mode 100755 Code/Fusion/otbBayesianFusionFilter.h
create mode 100755 Code/Fusion/otbBayesianFusionFilter.txx
create mode 100755 Code/Fusion/otbFusionImageBase.h
create mode 100644 Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
create mode 100644 Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilterNew.cxx
create mode 100755 Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
create mode 100644 Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilterNew.cxx
create mode 100755 Testing/Code/Fusion/CMakeLists.txt
create mode 100755 Testing/Code/Fusion/otbBayesianFusionFilter.cxx
create mode 100644 Testing/Code/Fusion/otbBayesianFusionFilterNew.cxx
create mode 100644 Testing/Code/Fusion/otbFusionImageBaseNew.cxx
create mode 100755 Testing/Code/Fusion/otbFusionTests.cxx
diff --git a/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.h b/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.h
new file mode 100644
index 0000000000..061ea7896d
--- /dev/null
+++ b/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.h
@@ -0,0 +1,214 @@
+/*=========================================================================
+
+ 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 __otbMatrixTransposeMatrixImageFilter_h
+#define __otbMatrixTransposeMatrixImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkNumericTraits.h"
+#include "itkArray.h"
+#include "itkSimpleDataObjectDecorator.h"
+#include "otbStreamingTraits.h"
+#include "itkImageRegionSplitter.h"
+#include "itkVariableSizeMatrix.h"
+#include "itkVariableLengthVector.h"
+
+
+namespace otb {
+
+ /** \class MatrixTransposeMatrixImageFilter
+ * \brief Compute \f[X^T.Y \f]. Allow a padding of ones.
+ *
+ * \f[X\f] and \f[Y\f] are the input images.
+ * The padding has the effect of adding a component filled with ones to the image
+ *
+ * \sa StreamingTraits
+ * \sa StatisticsImageFilter
+ * \ingroup MathematicalStatisticsImageFilters
+ */
+ template<class TInputImage, class TInputImage2>
+ class ITK_EXPORT MatrixTransposeMatrixImageFilter :
+ public itk::ImageToImageFilter<TInputImage, TInputImage>
+ {
+ public:
+ /** Standard Self typedef */
+ typedef MatrixTransposeMatrixImageFilter Self;
+ typedef itk::ImageToImageFilter<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(MatrixTransposeMatrixImageFilter, ImageToImageFilter);
+
+ /** Image related typedefs. */
+ // First Input
+ typedef TInputImage ImageType;
+ typedef typename TInputImage::Pointer InputImagePointer;
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename TInputImage::SizeType SizeType;
+ typedef typename TInputImage::IndexType IndexType;
+ typedef typename TInputImage::PixelType PixelType;
+ typedef typename TInputImage::InternalPixelType InternalPixelType;
+
+ typedef typename TInputImage2::IndexType IndexType2;
+ typedef typename TInputImage2::PixelType PixelType2;
+ typedef typename TInputImage2::InternalPixelType InternalPixelType2;
+
+ typedef StreamingTraits<TInputImage> StreamingTraitsType;
+ typedef StreamingMode StreamingModeType;
+ typedef itk::ImageRegionSplitter<itkGetStaticConstMacro(InputImageDimension)> SplitterType;
+ typedef typename SplitterType::Pointer SplitterPointer;
+
+ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+
+
+ /** Streaming-related accessors */
+ itkSetMacro(BufferMemorySize, unsigned long);
+ itkGetMacro(BufferMemorySize, unsigned long);
+ itkSetMacro(BufferNumberOfLinesDivisions, unsigned long);
+ itkGetMacro(BufferNumberOfLinesDivisions, unsigned long);
+ itkSetMacro(NumberOfStreamDivisions,unsigned long);
+ itkGetMacro(NumberOfStreamDivisions,unsigned long);
+ itkSetMacro(StreamingMode,StreamingModeType);
+ itkGetMacro(StreamingMode,StreamingModeType);
+
+
+ itkSetMacro(UsePadFirstInput,bool);
+ itkGetMacro(UsePadFirstInput,bool);
+ itkSetMacro(UsePadSecondInput,bool);
+ itkGetMacro(UsePadSecondInput,bool);
+
+
+ itkSetObjectMacro(RegionSplitter, SplitterType);
+ itkGetObjectMacro(RegionSplitter, SplitterType);
+
+ /** Image related typedefs. */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension );
+
+
+
+ /** Type to use for computations. */
+ // First Input
+ typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
+ typedef itk::VariableLengthVector<RealType> RealPixelType;
+
+
+ /** Smart Pointer type to a DataObject. */
+ typedef typename itk::DataObject::Pointer DataObjectPointer;
+
+ /** Type of DataObjects used for scalar outputs */
+ typedef typename itk::Array<long> ArrayLongPixelType;
+ typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
+ typedef typename std::vector<MatrixType> ArrayMatrixType;
+ typedef typename std::vector<RealPixelType> ArrayRealPixelType;
+ typedef typename std::vector<PixelType> ArrayPixelType;
+ typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<MatrixType> MatrixObjectType;
+
+
+ /** Return the computed transpose(Image1)*Image2. */
+ MatrixType GetResult() const { return this->GetResultOutput()->Get(); };
+ MatrixObjectType* GetResultOutput();
+ const MatrixObjectType* GetResultOutput() 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.
+ */
+ void AllocateOutputs();
+ void GenerateOutputInformation();
+
+
+ /** Input wrapper */
+ void SetFirstInput(const TInputImage *input1){ this->SetInput(0, input1 ); };
+ void SetSecondInput(const TInputImage2 *input2){ this->SetInput(1, input2 ); };
+
+ const TInputImage* GetFirstInput()
+ {
+ if( this->GetNumberOfInputs() < 1 )
+ {
+ return 0;
+ }
+ else
+ return( static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(0)) );
+ }
+
+ const TInputImage2* GetSecondInput()
+ {
+ if( this->GetNumberOfInputs() < 2 )
+ {
+ return 0;
+ }
+ else
+ return( static_cast<const TInputImage2 *>(this->itk::ProcessObject::GetInput(1)) );
+ }
+
+ protected:
+ MatrixTransposeMatrixImageFilter();
+ ~MatrixTransposeMatrixImageFilter(){};
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+ /** Initialize some accumulators before the threads run. */
+ virtual void BeforeThreadedGenerateData ();
+ /** Do final mean and variance computation from data accumulated in threads. */
+ virtual void AfterThreadedGenerateData ();
+ /** Multi-thread version GenerateData. */
+ virtual void ThreadedGenerateData (const RegionType& outputRegionForThread,int threadId);
+ // Override since the filter needs all the data for the algorithm
+ virtual void GenerateInputRequestedRegion();
+
+ private:
+ MatrixTransposeMatrixImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ ArrayMatrixType m_ThreadSum;
+ bool m_UsePadFirstInput;
+ bool m_UsePadSecondInput;
+
+ /** Nulber Of Component per Pixel. Change for padding */
+ unsigned int m_NumberOfComponents1;
+ unsigned int m_NumberOfComponents2;
+
+ /** Use to define the method used to calculate number of divisions */
+ unsigned long m_BufferMemorySize;
+ unsigned long m_BufferNumberOfLinesDivisions;
+ unsigned long m_NumberOfStreamDivisions;
+
+ SplitterPointer m_RegionSplitter;
+
+ /** Use to determine method of calculation number of divisions */
+ StreamingModeType m_StreamingMode;
+ }; // end of class
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMatrixTransposeMatrixImageFilter.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.txx b/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.txx
new file mode 100644
index 0000000000..b6b5e8ee0d
--- /dev/null
+++ b/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.txx
@@ -0,0 +1,374 @@
+/*=========================================================================
+
+ 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 _otbMatrixTransposeMatrixImageFilter_txx
+#define _otbMatrixTransposeMatrixImageFilter_txx
+#include "otbMatrixTransposeMatrixImageFilter.h"
+
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkNumericTraits.h"
+#include "itkProgressReporter.h"
+
+
+namespace otb {
+
+template<class TInputImage, class TInputImage2>
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::MatrixTransposeMatrixImageFilter()
+{
+ this->SetNumberOfRequiredInputs(2);
+
+ // 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
+
+
+ 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());
+
+ // false means no pad added
+ m_UsePadFirstInput = false;
+ m_UsePadSecondInput = false;
+ // Streaming initialization
+ m_BufferMemorySize = 0;
+ m_BufferNumberOfLinesDivisions = 0;
+ m_NumberOfStreamDivisions = 0;
+ // default to AUTOMATIC_NUMBER_OF_DIVISIONS
+ m_StreamingMode = SET_AUTOMATIC_NUMBER_OF_STREAM_DIVISIONS;
+ // create default region splitter
+ m_RegionSplitter = itk::ImageRegionSplitter<InputImageDimension>::New();
+
+ // Number of component initialisation
+ m_NumberOfComponents1 = 0;
+ m_NumberOfComponents2 = 0;
+}
+
+
+template<class TInputImage, class TInputImage2>
+itk::DataObject::Pointer
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::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, class TInputImage2>
+typename MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>::MatrixObjectType*
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::GetResultOutput()
+{
+ return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage, class TInputImage2>
+const typename MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>::MatrixObjectType*
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::GetResultOutput() const
+{
+ return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::GenerateInputRequestedRegion()
+{
+ Superclass::GenerateInputRequestedRegion();
+ if ( this->GetInput() )
+ {
+ InputImagePointer image = const_cast< typename Superclass::InputImageType * >( this->GetFirstInput() );
+ InputImagePointer image2 = const_cast< typename Superclass::InputImageType * >( this->GetSecondInput() );
+
+ IndexType index = image->GetLargestPossibleRegion().GetIndex();
+ SizeType size;
+ size.Fill(0);
+ RegionType region;
+ region.SetSize(size);
+ region.SetIndex(index);
+ image->SetRequestedRegion(region);
+
+ IndexType index2 = image2->GetLargestPossibleRegion().GetIndex();
+ SizeType size2;
+ size2.Fill(0);
+ RegionType region2;
+ region2.SetSize(size2);
+ region2.SetIndex(index2);
+ image2->SetRequestedRegion(region2);
+
+ }
+}
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation();
+ if ( this->GetInput() )
+ {
+ this->GetOutput()->CopyInformation(this->GetFirstInput());
+ this->GetOutput()->SetLargestPossibleRegion(this->GetFirstInput()->GetLargestPossibleRegion());
+ }
+}
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::AllocateOutputs()
+{
+
+ // Pass the input through as the output
+ InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
+ this->GraftOutput( image );
+ // Nothing that needs to be allocated for the remaining outputs
+}
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::BeforeThreadedGenerateData()
+{
+ if(this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()==0)
+ {
+ this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
+ }
+
+ if ( this->GetFirstInput()->GetLargestPossibleRegion().GetSize() != this->GetSecondInput()->GetLargestPossibleRegion().GetSize() )
+ {
+ itkExceptionMacro( <<" Can't multiply the transposed matrix of a "
+ << this->GetFirstInput()->GetLargestPossibleRegion().GetSize()
+ << " and a "
+ << this->GetSecondInput()->GetLargestPossibleRegion().GetSize()
+ << " matrix " );
+ }
+
+ m_NumberOfComponents1 = this->GetFirstInput()->GetNumberOfComponentsPerPixel();
+ m_NumberOfComponents2 = this->GetSecondInput()->GetNumberOfComponentsPerPixel();
+ unsigned int numberOfThreads = this->GetNumberOfThreads();
+
+ if ( m_UsePadFirstInput == true )
+ {
+ m_NumberOfComponents1++;
+ }
+ if ( m_UsePadSecondInput == true )
+ {
+ m_NumberOfComponents2++;
+ }
+
+ MatrixType tempMatrix, initMatrix;
+ tempMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
+ tempMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+ m_ThreadSum = ArrayMatrixType(numberOfThreads, tempMatrix);
+
+ initMatrix.SetSize(m_NumberOfComponents2, m_NumberOfComponents2);
+ initMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+ this->GetResultOutput()->Set( initMatrix );
+}
+
+
+
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::AfterThreadedGenerateData()
+{
+ unsigned int numberOfThreads = this->GetNumberOfThreads();
+ MatrixType resultMatrix;
+ resultMatrix.SetSize(m_NumberOfComponents1, m_NumberOfComponents2);
+ resultMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+
+
+ for( unsigned int thread = 0; thread < numberOfThreads; thread++)
+ {
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * resultMatrix += m_ThreadSum[thread];
+ **/
+ for (unsigned int i=0; i<resultMatrix.Rows(); i++)
+ {
+ for (unsigned int j=0; j<resultMatrix.Cols(); j++)
+ {
+ resultMatrix(i, j) += m_ThreadSum[thread](i, j);
+ }
+ }
+ /********END TODO ******/
+ }
+ this->GetResultOutput()->Set( resultMatrix );
+}
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId)
+{
+ /**
+ * Grab the input
+ */
+
+ InputImagePointer input1Ptr = const_cast< TInputImage * >( this->GetFirstInput() );
+ InputImagePointer input2Ptr = const_cast< TInputImage2 * >( this->GetSecondInput() );
+
+ // support progress methods/callbacks
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ // Here we will divide the input region into pieces to process one afther the other.
+ unsigned int numDivisions = static_cast<unsigned int>(StreamingTraitsType
+ ::CalculateNumberOfStreamDivisions(this->GetFirstInput(),
+ outputRegionForThread,
+ m_StreamingMode,
+ m_NumberOfStreamDivisions,
+ m_BufferMemorySize,
+ m_BufferNumberOfLinesDivisions));
+
+
+ otbMsgDebugMacro(<<"ThreadedGeneraeData() - thread "<<threadId<<" - Thread region: "<<outputRegionForThread);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId<<" - Streaming configuration: "<<m_StreamingMode<<" "<<m_NumberOfStreamDivisions<<" "<<m_BufferMemorySize<<" "<<m_BufferNumberOfLinesDivisions);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - nb of divisions from StreamingTraits: "<<numDivisions);
+
+
+ unsigned int numDivisionsFromSplitter = m_RegionSplitter->GetNumberOfSplits(outputRegionForThread, numDivisions);
+
+
+ if (numDivisionsFromSplitter < numDivisions)
+ {
+ numDivisions = numDivisionsFromSplitter;
+ }
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - nb of divisions: "<<numDivisions);
+
+ /*
+ unsigned int numberOfComponents1 = this->GetFirstInput()->GetNumberOfComponentsPerPixel();
+ unsigned int numberOfComponents2 = this->GetSecondInput()->GetNumberOfComponentsPerPixel();
+
+ MatrixType result;
+ result.SetSize(numberOfComponents1, numberOfComponents2);
+ result.Fill(itk::NumericTraits<RealType>::Zero);
+ */
+
+ /**
+ * Loop over the number of pieces, execute the upstream pipeline on each
+ * piece, and copy the results into the output image.
+ */
+ unsigned int piece;
+ for (piece = 0;
+ piece < numDivisions && !this->GetAbortGenerateData();
+ piece++)
+ {
+ RegionType streamRegion = m_RegionSplitter->GetSplit(piece,numDivisions,outputRegionForThread);
+
+
+ input1Ptr->SetRequestedRegion(streamRegion);
+ input2Ptr->SetRequestedRegion(streamRegion);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - streaming region: "<<streamRegion);
+ input1Ptr->PropagateRequestedRegion();
+ input1Ptr->UpdateOutputData();
+ input2Ptr->PropagateRequestedRegion();
+ input2Ptr->UpdateOutputData();
+
+ itk::ImageRegionConstIteratorWithIndex<TInputImage> it1 (input1Ptr, streamRegion);
+ itk::ImageRegionConstIteratorWithIndex<TInputImage> it2 (input2Ptr, streamRegion);
+ it1.GoToBegin();
+ it2.GoToBegin();
+
+ // loop the second image and get one pixel a time
+ while (!it1.IsAtEnd())
+ {
+ //IndexType index1 = it1.GetIndex();
+ PixelType vectorValue1 = it1.Get();
+ //IndexType2 index2 = it2.GetIndex();
+ PixelType2 vectorValue2 = it2.Get();
+
+ // Add a first component to vectorValue2 and vectorValue1 filled with ones.
+ if (m_UsePadFirstInput == true)
+ {
+ PixelType vectortemp1(vectorValue1.Size()+1);
+ vectortemp1[0] = 1;
+ for (unsigned int n=0; n<vectorValue1.Size(); n++)
+ {
+ vectortemp1[n+1] = vectorValue1[n];
+
+ }
+ vectorValue1.SetSize(vectortemp1.Size());
+ vectorValue1 = vectortemp1;
+ }
+
+ if (m_UsePadSecondInput == true)
+ {
+ PixelType2 vectortemp2(vectorValue2.Size()+1);
+ vectortemp2[0] = 1;
+ for (unsigned int m=0; m<vectorValue2.Size(); m++)
+ {
+ vectortemp2[m+1] = vectorValue2[m];
+
+ }
+ vectorValue2.SetSize(vectortemp2.Size());
+ vectorValue2 = vectortemp2;
+ }
+
+ for (unsigned int i=0; i<vectorValue1.Size(); i++)
+ {
+ for (unsigned int j=0; j<vectorValue2.Size(); j++)
+ {
+ m_ThreadSum[threadId](i, j) += vectorValue1[i]*static_cast<InternalPixelType>(vectorValue2[j]);
+ }
+
+ }
+ ++it1;
+ ++it2;
+ progress.CompletedPixel();
+ }
+ }
+ //otbMsgDebugMacro(<<"Leaving ThreadedGenerateData()");
+}
+
+template<class TInputImage, class TInputImage2>
+void
+MatrixTransposeMatrixImageFilter<TInputImage, TInputImage2>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os,indent);
+
+ os << indent << "Result: " <<this->GetResultOutput()->Get()<< std::endl;
+}
+
+
+}// end namespace otb
+#endif
diff --git a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
new file mode 100644
index 0000000000..9d3af2769e
--- /dev/null
+++ b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
@@ -0,0 +1,198 @@
+/*=========================================================================
+
+ 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 __otbStreamingStatisticsVectorImageFilter_h
+#define __otbStreamingStatisticsVectorImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkNumericTraits.h"
+#include "itkArray.h"
+#include "itkSimpleDataObjectDecorator.h"
+#include "otbStreamingTraits.h"
+#include "itkImageRegionSplitter.h"
+#include "itkVariableSizeMatrix.h"
+#include "itkVariableLengthVector.h"
+
+
+namespace otb {
+
+/** \class StreamingStatisticsVectorImageFilter
+ * \brief Compute min. max, covariance of a large image using streaming
+ *
+ * This filter Computes the same statistics as the StatisticsImageFilter, but it
+ * supports large image since it processes reasonnable pieces of the input image one
+ * afther the other. It supports the same streaming mode than the StreamingImageVectorFileWriter.
+ *
+ * Of course streaming at the end of a pipeline is only available when each filter in the pipeline
+ * supports streaming. This filter will also perform multithreading if possible.
+ *
+ * \note The output image has no sense at all and should not be used.
+ *
+ * \sa StreamingTraits
+ * \sa StreamingImageVectorFileWriter
+ * \sa StreamingImageFileWriter
+ * \sa StatisticsImageFilter
+ * \ingroup Streamed
+ * \ingroup Multithreaded
+ * \ingroup MathematicalStatisticsImageFilters
+ */
+template<class TInputImage>
+class ITK_EXPORT StreamingStatisticsVectorImageFilter :
+ public itk::ImageToImageFilter<TInputImage, TInputImage>
+ {
+ public:
+ /** Standard Self typedef */
+ typedef StreamingStatisticsVectorImageFilter Self;
+ typedef itk::ImageToImageFilter<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(StreamingStatisticsVectorImageFilter, ImageToImageFilter);
+
+ /** Image related typedefs. */
+ typedef typename TInputImage::Pointer InputImagePointer;
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename TInputImage::SizeType SizeType;
+ typedef typename TInputImage::IndexType IndexType;
+ typedef typename TInputImage::PixelType PixelType;
+ typedef typename TInputImage::InternalPixelType InternalPixelType;
+ typedef StreamingTraits<TInputImage> StreamingTraitsType;
+ typedef StreamingMode StreamingModeType;
+ typedef itk::ImageRegionSplitter<itkGetStaticConstMacro(InputImageDimension)> SplitterType;
+ typedef typename SplitterType::Pointer SplitterPointer;
+
+ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ /** Streaming-related accessors */
+ itkSetMacro(BufferMemorySize, unsigned long);
+ itkGetMacro(BufferMemorySize, unsigned long);
+ itkSetMacro(BufferNumberOfLinesDivisions, unsigned long);
+ itkGetMacro(BufferNumberOfLinesDivisions, unsigned long);
+ itkSetMacro(NumberOfStreamDivisions,unsigned long);
+ itkGetMacro(NumberOfStreamDivisions,unsigned long);
+ itkSetMacro(StreamingMode,StreamingModeType);
+ itkGetMacro(StreamingMode,StreamingModeType);
+
+ itkSetObjectMacro(RegionSplitter, SplitterType);
+ itkGetObjectMacro(RegionSplitter, SplitterType);
+
+ /** Image related typedefs. */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension );
+
+ /** Type to use for computations. */
+ typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
+ typedef itk::VariableLengthVector<RealType> RealPixelType;
+
+ /** Smart Pointer type to a DataObject. */
+ typedef typename itk::DataObject::Pointer DataObjectPointer;
+
+ /** Type of DataObjects used for scalar outputs */
+ typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
+ typedef typename std::vector<MatrixType> ArrayMatrixType;
+ typedef typename itk::Array<long> ArrayLongPixelType;
+ typedef typename std::vector<RealPixelType> ArrayRealPixelType;
+ typedef typename std::vector<PixelType> ArrayPixelType;
+ typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<MatrixType> MatrixObjectType;
+
+
+ /** Return the computed Minimum. */
+ PixelType GetMinimum() const { return this->GetMinimumOutput()->Get(); };
+ PixelObjectType* GetMinimumOutput();
+ const PixelObjectType* GetMinimumOutput() const;
+
+ /** Return the computed Maximum. */
+ PixelType GetMaximum() const { return this->GetMaximumOutput()->Get(); };
+ PixelObjectType* GetMaximumOutput();
+ const PixelObjectType* GetMaximumOutput() const;
+
+ /** Return the computed Mean. */
+ RealPixelType GetMean() const { return this->GetMeanOutput()->Get(); };
+ RealPixelObjectType* GetMeanOutput();
+ const RealPixelObjectType* GetMeanOutput() const;
+
+ /** Return the compute Sum. */
+ RealPixelType GetSum() const { return this->GetSumOutput()->Get(); };
+ RealPixelObjectType* GetSumOutput();
+ const RealPixelObjectType* GetSumOutput() const;
+
+ /** Return the computed Covariance. */
+ MatrixType GetCovariance() const { return this->GetCovarianceOutput()->Get(); };
+ MatrixObjectType* GetCovarianceOutput();
+ const MatrixObjectType* GetCovarianceOutput() 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.
+ */
+ void AllocateOutputs();
+ void GenerateOutputInformation();
+
+ protected:
+ StreamingStatisticsVectorImageFilter();
+ ~StreamingStatisticsVectorImageFilter(){};
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+ /** Initialize some accumulators before the threads run. */
+ void BeforeThreadedGenerateData ();
+ /** Do final mean and variance computation from data accumulated in threads. */
+ void AfterThreadedGenerateData ();
+ /** Multi-thread version GenerateData. */
+ void ThreadedGenerateData (const RegionType& outputRegionForThread,int threadId);
+ // Override since the filter needs all the data for the algorithm
+ void GenerateInputRequestedRegion();
+
+ private:
+ StreamingStatisticsVectorImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ ArrayRealPixelType m_ThreadSum;
+ ArrayLongPixelType m_Count;
+ ArrayPixelType m_ThreadMin;
+ ArrayPixelType m_ThreadMax;
+ ArrayMatrixType m_XX;
+
+ /** Use to define the method used to calculate number of divisions */
+ unsigned long m_BufferMemorySize;
+ unsigned long m_BufferNumberOfLinesDivisions;
+ unsigned long m_NumberOfStreamDivisions;
+
+ SplitterPointer m_RegionSplitter;
+
+ /** Use to determine method of calculation number of divisions */
+ StreamingModeType m_StreamingMode;
+ }; // end of class
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbStreamingStatisticsVectorImageFilter.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
new file mode 100644
index 0000000000..aef0e7f995
--- /dev/null
+++ b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
@@ -0,0 +1,554 @@
+/*=========================================================================
+
+ 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 _otbStreamingStatisticsVectorImageFilter_txx
+#define _otbStreamingStatisticsVectorImageFilter_txx
+#include "otbStreamingStatisticsVectorImageFilter.h"
+
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkNumericTraits.h"
+#include "itkProgressReporter.h"
+
+
+namespace otb {
+
+template<class TInputImage>
+StreamingStatisticsVectorImageFilter<TInputImage>
+::StreamingStatisticsVectorImageFilter()
+{
+
+ // 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
+
+ for (int i=1; i < 3; ++i)
+ {
+ typename PixelObjectType::Pointer output = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
+ }
+
+ // allocate the data objects for the outputs which are
+ // just decorators around real types
+
+ for (int i=3; i < 5; ++i)
+ {
+ typename RealPixelObjectType::Pointer output = static_cast<RealPixelObjectType*>(this->MakeOutput(i).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
+ }
+
+ // allocate the data objects for the outputs which are
+ // just decorators around matrix types
+ typename MatrixObjectType::Pointer output = static_cast<MatrixObjectType*>(this->MakeOutput(5).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(5, output.GetPointer());
+
+ // Streaming initialization
+ m_BufferMemorySize = 0;
+ m_BufferNumberOfLinesDivisions = 0;
+ m_NumberOfStreamDivisions = 0;
+ // default to AUTOMATIC_NUMBER_OF_DIVISIONS
+ m_StreamingMode = SET_AUTOMATIC_NUMBER_OF_STREAM_DIVISIONS;
+
+ // create default region splitter
+ m_RegionSplitter = itk::ImageRegionSplitter<InputImageDimension>::New();
+
+}
+
+
+template<class TInputImage>
+itk::DataObject::Pointer
+StreamingStatisticsVectorImageFilter<TInputImage>
+::MakeOutput(unsigned int output)
+{
+ switch (output)
+ {
+ case 0:
+ return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
+ break;
+ case 1:
+ case 2:
+ return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
+ break;
+ case 3:
+ case 4:
+ return static_cast<itk::DataObject*>(RealPixelObjectType::New().GetPointer());
+ break;
+ case 5:
+ 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 StreamingStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMinimumOutput()
+{
+ return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage>
+const typename StreamingStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMinimumOutput() const
+{
+ return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+
+template<class TInputImage>
+typename StreamingStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMaximumOutput()
+{
+ return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
+}
+
+template<class TInputImage>
+const typename StreamingStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMaximumOutput() const
+{
+ return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
+}
+
+
+template<class TInputImage>
+typename StreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMeanOutput()
+{
+ return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
+}
+
+template<class TInputImage>
+const typename StreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetMeanOutput() const
+{
+ return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
+}
+
+template<class TInputImage>
+typename StreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetSumOutput()
+{
+ return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
+}
+
+template<class TInputImage>
+const typename StreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetSumOutput() const
+{
+ return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
+}
+
+template<class TInputImage>
+typename StreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetCovarianceOutput()
+{
+ return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
+}
+
+template<class TInputImage>
+const typename StreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GetCovarianceOutput() const
+{
+ return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
+}
+
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GenerateInputRequestedRegion()
+{
+ Superclass::GenerateInputRequestedRegion();
+ if ( this->GetInput() )
+ {
+ InputImagePointer image = const_cast< typename Superclass::InputImageType * >( this->GetInput() );
+ IndexType index = image->GetLargestPossibleRegion().GetIndex();
+ SizeType size;
+ size.Fill(0);
+ RegionType region;
+ region.SetSize(size);
+ region.SetIndex(index);
+ image->SetRequestedRegion(region);
+ }
+}
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation();
+ if ( this->GetInput() )
+ {
+ this->GetOutput()->CopyInformation(this->GetInput());
+ this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
+ }
+}
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::AllocateOutputs()
+{
+
+ // Pass the input through as the output
+ InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
+ this->GraftOutput( image );
+ // Nothing that needs to be allocated for the remaining outputs
+}
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::BeforeThreadedGenerateData()
+{
+ if(this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()==0)
+ {
+ this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
+ }
+
+ unsigned int numberOfThreads = this->GetNumberOfThreads();
+ unsigned int numberOfComponent = this->GetInput()->GetNumberOfComponentsPerPixel();
+ // Variable Initialisation
+ PixelType tempPixel;
+ tempPixel.SetSize(numberOfComponent);
+ tempPixel.Fill(itk::NumericTraits<InternalPixelType>::Zero);
+ this->GetMaximumOutput()->Set( tempPixel );
+
+ tempPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
+ this->GetMinimumOutput()->Set( tempPixel );
+
+ RealPixelType tempRealPixel;
+ tempRealPixel.SetSize(numberOfComponent);
+ tempRealPixel.Fill(itk::NumericTraits<RealType>::max());
+ this->GetMeanOutput()->Set( tempRealPixel );
+
+ tempRealPixel.Fill(itk::NumericTraits<RealType>::Zero);
+ this->GetMeanOutput()->Set( tempRealPixel );
+
+ MatrixType tempMatrix;
+ tempMatrix.SetSize(numberOfComponent, numberOfComponent);
+ tempMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+
+
+ // Initialize the tempories
+ m_Count.SetSize(numberOfThreads);
+ m_Count.Fill( itk::NumericTraits<long>::Zero );
+
+ PixelType tempTemporiesPixel;
+ tempTemporiesPixel.SetSize(numberOfComponent);
+ tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
+ m_ThreadMin = ArrayPixelType(numberOfThreads, tempTemporiesPixel);
+
+ tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::Zero);
+ m_ThreadMax = ArrayPixelType(numberOfThreads, tempTemporiesPixel);
+
+ RealPixelType tempTemporiesRealPixel;
+ tempTemporiesRealPixel.SetSize(numberOfComponent);
+ tempTemporiesRealPixel.Fill(itk::NumericTraits<RealType>::Zero);
+ m_ThreadSum = ArrayRealPixelType(numberOfThreads, tempTemporiesRealPixel);
+
+ m_XX = ArrayMatrixType(numberOfThreads, tempMatrix);
+
+ otbMsgDevMacro(<<"BeforeThreadedGenerateData() - output requested region: "<<this->GetOutput()->GetRequestedRegion());
+ otbMsgDevMacro(<<"Leaving BeforeThreadedGenerateData() - nb threads: "<<numberOfThreads);
+}
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::AfterThreadedGenerateData()
+{
+ //otbMsgDebugMacro(<<"Entering AfterThreadedGenerateData()");
+ int i;
+ long count;
+
+ int numberOfThreads = this->GetNumberOfThreads();
+ unsigned int numberOfComponent = this->GetInput()->GetNumberOfComponentsPerPixel();
+
+ PixelType minimumVector;
+ minimumVector.SetSize( numberOfComponent );
+ minimumVector.Fill(itk::NumericTraits<InternalPixelType>::max());
+
+ PixelType maximumVector;
+ maximumVector.SetSize( numberOfComponent );
+ maximumVector.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
+
+ RealPixelType sumVector;
+ sumVector.SetSize( numberOfComponent );
+ sumVector.Fill(itk::NumericTraits<InternalPixelType>::Zero);
+
+ RealPixelType meanVector;
+ meanVector.SetSize( numberOfComponent );
+ MatrixType crossedMatrix;
+ crossedMatrix.SetSize(numberOfComponent,numberOfComponent);
+ crossedMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+ count = 0;
+
+ // Find the min/max over all threads and accumulate count, sum and
+ // sum of squares
+ for( i = 0; i < numberOfThreads; i++)
+ {
+ count += m_Count[i];
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * crossedMatrix += m_XX[i];
+ **/
+ if( (m_XX[i].Rows() != crossedMatrix.Rows()) || (m_XX[i].Cols() != crossedMatrix.Cols()))
+ {
+ itkExceptionMacro( << "Matrix with size (" << m_XX[i].Rows() << "," <<
+ m_XX[i].Cols() << ") cannot be subtracted from matrix with size (" <<
+ crossedMatrix.Rows() << "," << crossedMatrix.Cols() );
+ }
+
+ for( unsigned int r=0; r<m_XX[i].Rows(); r++)
+ {
+ for( unsigned int c=0; c<m_XX[i].Cols(); c++ )
+ {
+ crossedMatrix(r,c) += m_XX[i](r,c);
+ }
+ }
+ //**** END TODO ****
+
+ for (unsigned int j=0; j<numberOfComponent; j++)
+ {
+ sumVector[j] += m_ThreadSum[i][j];
+ if (m_ThreadMin[i][j] < minimumVector[j])
+ {
+ minimumVector[j] = m_ThreadMin[i][j];
+ }
+ if (m_ThreadMax[i][j] > maximumVector[j])
+ {
+ maximumVector[j] = m_ThreadMax[i][j];
+ }
+ }
+ }
+ for (unsigned int j=0; j<numberOfComponent; j++)
+ {
+ // compute statistics
+ meanVector[j] = sumVector[j] / static_cast<RealType>(count);
+ }
+
+ // Compute Matrix Covariance
+ MatrixType pixelSumMatrix;
+ pixelSumMatrix.SetSize(static_cast<unsigned int>(numberOfComponent), 1);
+ pixelSumMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+ for( unsigned int j = 0; j < numberOfComponent; j++)
+ {
+ pixelSumMatrix(j, 0) = sumVector[j];
+ }
+
+ MatrixType covMatrix, covMatrixTemp, tempTranspose;
+ covMatrix.SetSize(static_cast<unsigned int>(numberOfComponent), static_cast<unsigned int>(numberOfComponent));
+ covMatrixTemp.SetSize(static_cast<unsigned int>(numberOfComponent), static_cast<unsigned int>(numberOfComponent));
+ tempTranspose.SetSize(1, static_cast<unsigned int>(numberOfComponent));
+
+ covMatrix = crossedMatrix/static_cast<RealType>(count);
+ pixelSumMatrix/=static_cast<RealType>(count);
+ tempTranspose = pixelSumMatrix.GetTranspose();
+ covMatrixTemp = pixelSumMatrix*tempTranspose;
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ *covMatrix -= covMatrixTemp;
+ **/
+ if( (covMatrix.Rows() != covMatrixTemp.Rows()) || (covMatrix.Cols() != covMatrixTemp.Cols()))
+ {
+ itkExceptionMacro( << "Matrix with size (" << covMatrix.Rows() << "," <<
+ covMatrix.Cols() << ") cannot be subtracted from matrix with size (" <<
+ covMatrixTemp.Rows() << "," << covMatrixTemp.Cols() );
+ }
+
+ for( unsigned int r=0; r<covMatrix.Rows(); r++)
+ {
+ for( unsigned int c=0; c<covMatrix.Cols(); c++ )
+ {
+ covMatrix(r,c) -= covMatrixTemp(r,c);
+ }
+ }
+ //**** END TODO ****/
+
+ // Set the outputs
+ this->GetMinimumOutput()->Set( minimumVector );
+ this->GetMaximumOutput()->Set( maximumVector );
+ this->GetMeanOutput()->Set( meanVector );
+ this->GetSumOutput()->Set( sumVector);
+ this->GetCovarianceOutput()->Set( covMatrix );
+
+ //otbMsgDebugMacro(<<"Leaving AfterThreadedGenerateData()");
+}
+
+template<class TInputImage>
+void
+StreamingStatisticsVectorImageFilter<TInputImage>
+::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId)
+{
+ //otbMsgDebugMacro(<<"Entering ThreadedGenerateData()");
+ /**
+ * Grab the input
+ */
+
+ InputImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput(0) );
+
+
+ // support progress methods/callbacks
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ // Here we will divide the input region into pieces to process one afther the other.
+ unsigned int numDivisions = static_cast<unsigned int>(StreamingTraitsType
+ ::CalculateNumberOfStreamDivisions(this->GetInput(),
+ outputRegionForThread,
+ m_StreamingMode,
+ m_NumberOfStreamDivisions,
+ m_BufferMemorySize,
+ m_BufferNumberOfLinesDivisions)); // calcul en combien fau diviser idealement l'image
+
+ otbMsgDebugMacro(<<"ThreadedGeneraeData() - thread "<<threadId <<" - Thread region: "<<outputRegionForThread);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId<<" - Streaming configuration: "<<m_StreamingMode<<" "<<m_NumberOfStreamDivisions<<" "<<m_BufferMemorySize<<" "<<m_BufferNumberOfLinesDivisions);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - nb of divisions from StreamingTraits: "<<numDivisions);
+
+ unsigned int numDivisionsFromSplitter = m_RegionSplitter->GetNumberOfSplits(outputRegionForThread, numDivisions);
+ if (numDivisionsFromSplitter < numDivisions)
+ {
+ numDivisions = numDivisionsFromSplitter;
+ }
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - nb of divisions: "<<numDivisions);
+
+
+ MatrixType pixelVector, pixelTransposeVector, pixelSumVector, tempMatrix;
+
+ pixelVector.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(), 1);
+ pixelVector.Fill(itk::NumericTraits<RealType>::Zero);
+ pixelTransposeVector.SetSize(1, this->GetInput()->GetNumberOfComponentsPerPixel());
+ pixelTransposeVector.Fill(itk::NumericTraits<RealType>::Zero);
+ pixelSumVector.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(), this->GetInput()->GetNumberOfComponentsPerPixel());
+ pixelSumVector.Fill(itk::NumericTraits<RealType>::Zero);
+ tempMatrix.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(), this->GetInput()->GetNumberOfComponentsPerPixel());
+
+ /**
+ * Loop over the number of pieces, execute the upstream pipeline on each
+ * piece, and copy the results into the output image.
+ */
+ unsigned int piece;
+ for (piece = 0;
+ piece < numDivisions && !this->GetAbortGenerateData();
+ piece++)
+ {
+ //otbMsgDebugMacro(<<"ThreadedGenerateData() - processing piece: "<<piece<<"/"<<numDivisions);
+ RegionType streamRegion = m_RegionSplitter->GetSplit(piece,numDivisions,outputRegionForThread);
+ //otbMsgDebugMacro(<<"ThreadedGenerateData() - piece region: "<<streamRegion);
+
+ inputPtr->SetRequestedRegion(streamRegion);
+ otbMsgDebugMacro(<<"ThreadedGenerateData() - thread "<<threadId <<" - streaming region: "<<streamRegion);
+ inputPtr->PropagateRequestedRegion();
+ inputPtr->UpdateOutputData();
+
+
+
+ itk::ImageRegionConstIteratorWithIndex<TInputImage> it (inputPtr, streamRegion);
+ it.GoToBegin();
+ // do the work
+ while (!it.IsAtEnd())
+ {
+
+ IndexType index = it.GetIndex();
+ PixelType vectorValue = it.Get();
+ for (unsigned int j=0; j<vectorValue.GetSize(); j++)
+ {
+ InternalPixelType value = vectorValue[j];
+
+ RealType realValue = static_cast<RealType>( value );
+ if (value < m_ThreadMin[threadId][j])
+ {
+ m_ThreadMin[threadId][j] = value;
+ }
+ if (value > m_ThreadMax[threadId][j])
+ {
+ m_ThreadMax[threadId][j] = value;
+ }
+ m_ThreadSum[threadId][j] += realValue;
+ pixelVector(j, 0) = realValue;
+ }
+
+ ++it;
+ progress.CompletedPixel();
+ pixelTransposeVector = pixelVector.GetTranspose();
+
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * m_XX[threadId]+=pixelVector*pixelTransposeVector;
+ **/
+ tempMatrix = pixelVector*pixelTransposeVector;
+ if( (m_XX[threadId].Rows() != tempMatrix.Rows()) || (m_XX[threadId].Cols() != tempMatrix.Cols()))
+ {
+ itkExceptionMacro( << "Matrix with size (" << m_XX[threadId].Rows() << "," <<
+ m_XX[threadId].Cols() << ") cannot be subtracted from matrix with size (" <<
+ tempMatrix.Rows() << "," << tempMatrix.Cols() );
+ }
+
+ for( unsigned int r=0; r<m_XX[threadId].Rows(); r++)
+ {
+ for( unsigned int c=0; c<m_XX[threadId].Cols(); c++ )
+ {
+ m_XX[threadId](r,c) += tempMatrix(r,c);
+ }
+ }
+ //**** END TODO ****
+ m_Count[threadId]++;
+ }
+ }
+
+ //otbMsgDebugMacro(<<"Leaving ThreadedGenerateData()");
+}
+
+template <class TImage>
+void
+StreamingStatisticsVectorImageFilter<TImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os,indent);
+
+ os << indent << "Minimum: " <<this->GetMinimumOutput()->Get()<< std::endl;
+ os << indent << "Maximum: "<< this->GetMaximumOutput()->Get() << std::endl;
+ os << indent << "Sum: " << this->GetSumOutput()->Get() << std::endl;
+ os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
+ os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
+
+}
+
+
+}// end namespace otb
+#endif
diff --git a/Code/CMakeLists.txt b/Code/CMakeLists.txt
index d493a648f8..167aaa0451 100644
--- a/Code/CMakeLists.txt
+++ b/Code/CMakeLists.txt
@@ -10,6 +10,7 @@ DisparityMap
SpatialReasoning
Projections
Radiometry
+Fusion
)
IF(OTB_USE_VISU_GUI)
diff --git a/Code/Fusion/CMakeLists.txt b/Code/Fusion/CMakeLists.txt
new file mode 100755
index 0000000000..004e778918
--- /dev/null
+++ b/Code/Fusion/CMakeLists.txt
@@ -0,0 +1,18 @@
+# Sources of non-templated classes.
+
+FILE(GLOB OTBFusion_SRCS "*.cxx" )
+
+ADD_LIBRARY(OTBFusion ${OTBFusion_SRCS})
+TARGET_LINK_LIBRARIES (OTBFusion ITKCommon)
+
+INSTALL(TARGETS OTBFusion
+RUNTIME DESTINATION ${OTB_INSTALL_BIN_DIR} COMPONENT RuntimeLibraries
+LIBRARY DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
+ARCHIVE DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT Development)
+
+FILE(GLOB __files1 "${CMAKE_CURRENT_SOURCE_DIR}/*.h")
+FILE(GLOB __files2 "${CMAKE_CURRENT_SOURCE_DIR}/*.txx")
+
+INSTALL(FILES ${__files1} ${__files2}
+ DESTINATION ${OTB_INSTALL_INCLUDE_DIR}/Fusion
+ COMPONENT Development)
diff --git a/Code/Fusion/foo.cxx b/Code/Fusion/foo.cxx
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/Code/Fusion/otbBayesianFusionFilter.h b/Code/Fusion/otbBayesianFusionFilter.h
new file mode 100755
index 0000000000..68e76858b9
--- /dev/null
+++ b/Code/Fusion/otbBayesianFusionFilter.h
@@ -0,0 +1,278 @@
+/*=========================================================================
+
+ 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 code is from a Julien Radoux contribution.
+
+ 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 __otbBayesianFusionFilter_h
+#define __otbBayesianFusionFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "otbVectorImage.h"
+#include "itkNumericTraits.h"
+#include "otbStreamingStatisticsVectorImageFilter.h"
+#include "otbFusionImageBase.h"
+#include "otbMatrixTransposeMatrixImageFilter.h"
+#include "otbImageToVectorImageCastFilter.h"
+#include "itkVariableSizeMatrix.h"
+
+
+namespace otb
+{
+ namespace Functor
+ {
+ template <class TInputMultiSpectral,
+ class TInputMultiSpectralInterp,
+ class TInputPanchro,
+ class TOutput>
+ class BayesianFunctor
+ {
+ public:
+ BayesianFunctor() {};
+ ~BayesianFunctor() {};
+ typedef typename TInputMultiSpectral::RealValueType RealType;
+ typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
+
+ void SetLambda(float lambda){ m_Lambda = lambda;};
+ void SetS(float S){ m_S = S;};
+ void SetAlpha(float alpha){ m_Alpha = alpha;};
+ void SetBeta(MatrixType matrix){ m_Beta = matrix;};
+ void SetCovarianceInvMatrix(MatrixType matrix){ m_CovarianceInvMatrix = matrix;};
+ void SetVcondopt(MatrixType matrix){ m_Vcondopt = matrix;};
+ float GetLambda(){ return m_Lambda;};
+ float GetAlpha(){ return m_Alpha;};
+ float GetS(){ return m_S;};
+ MatrixType GetBeta(){ return m_Beta;};
+ MatrixType GetCovarianceInvMatrix(){ return m_CovarianceInvMatrix;};
+ MatrixType GetVcondopt(){ return m_Vcondopt;};
+
+
+
+ inline TOutput operator() (const TInputMultiSpectral & ms, const TInputMultiSpectralInterp & msi, const TInputPanchro & p)
+ {
+ TOutput obs;
+ obs.SetSize(msi.GetSize());
+ MatrixType obsMat, msiVect;
+ obsMat.SetSize(1, obs.GetSize());
+ msiVect.SetSize(1, msi.GetSize());
+ for (unsigned int i=0; i<msi.GetSize();i++)
+ {
+ msiVect(0, i) = msi[i];
+ }
+ obsMat = msiVect*m_CovarianceInvMatrix;
+ obsMat *= 2*(1-m_Lambda);
+ MatrixType PanVect;
+ PanVect = m_Beta.GetTranspose();
+ PanVect *= (p-m_Alpha);
+ PanVect /= m_S;
+ PanVect *= 2*m_Lambda;
+
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * obsMat += PanVect;
+ **/
+ if( (obsMat.Cols() != PanVect.Cols()) || (obsMat.Rows() != PanVect.Rows()) )
+ {
+ itkGenericExceptionMacro( << "Matrix with size (" << obsMat.Rows() << "," <<
+ obsMat.Cols() << ") cannot be subtracted from matrix with size (" <<
+ PanVect.Rows() << "," << PanVect.Cols() << " )" );
+ }
+
+ for( unsigned int r=0; r<obsMat.Rows(); r++)
+ {
+ for( unsigned int c=0; c<obsMat.Cols(); c++ )
+ {
+ obsMat(r,c) += PanVect(r,c);
+ }
+ }
+ //**** END TODO ****/
+ obsMat *= m_Vcondopt;
+ for (unsigned int i=0; i<obs.GetSize();i++)
+ {
+ obs[i] = obsMat(0,i);
+ }
+ return obs;
+ }
+
+ private:
+ float m_Lambda;
+ float m_S;
+ float m_Alpha;
+ MatrixType m_CovarianceInvMatrix;
+ MatrixType m_Beta;
+ MatrixType m_Vcondopt;
+ };
+ }
+
+ /***** TODO ***
+ * Complete the description with J. Radoux text
+ */
+ /***** END TODO ***/
+
+ /** \class StreamingStatisticsVectorImageFilter
+ * \brief Baesian fusion filter. Contribution of Julien Radoux
+ *
+ *
+ * \sa FusionImageBase
+ * \sa MatrixTransposeMatrix
+ * \sa StreamingStatisticsVectorImageFilter
+ * \ingroup Streamed
+ * \ingroup Multithreaded
+ * \ingroup MathematicalStatisticsImageFilters
+ */
+
+
+
+
+template <class TInputMultiSpectralImage,
+ class TInputMultiSpectralInterpImage,
+ class TInputPanchroImage,
+ class TOutputImage>
+class ITK_EXPORT BayesianFusionFilter
+ : public FusionImageBase<TInputMultiSpectralImage,
+ TInputMultiSpectralInterpImage,
+ TInputPanchroImage,
+ TOutputImage,
+ Functor::BayesianFunctor<ITK_TYPENAME TInputMultiSpectralImage::PixelType,
+ ITK_TYPENAME TInputMultiSpectralInterpImage::PixelType,
+ ITK_TYPENAME TInputPanchroImage::PixelType,
+ ITK_TYPENAME TOutputImage::PixelType> >
+{
+public:
+/** Extract input and output images dimensions.*/
+ itkStaticConstMacro( InputImageDimension, unsigned int, TInputMultiSpectralImage::ImageDimension);
+ itkStaticConstMacro( OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
+
+ /** "typedef" to simplify the variables definition and the declaration. */
+ typedef TInputMultiSpectralImage InputMultiSpectralImageType;
+ typedef TInputMultiSpectralInterpImage InputMultiSpectralInterpImageType;
+ typedef TInputPanchroImage InputPanchroImageType;
+ typedef TOutputImage OutputImageType;
+
+
+
+ /** "typedef" for standard classes. */
+ typedef BayesianFusionFilter Self;
+ typedef FusionImageBase< InputMultiSpectralImageType,
+ InputMultiSpectralInterpImageType,
+ InputPanchroImageType,
+ OutputImageType,
+ Functor::BayesianFunctor<ITK_TYPENAME InputMultiSpectralImageType::PixelType,
+ ITK_TYPENAME InputMultiSpectralInterpImageType::PixelType,
+ ITK_TYPENAME InputPanchroImageType::PixelType,
+ ITK_TYPENAME OutputImageType::PixelType> > Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** object factory method. */
+ itkNewMacro(Self);
+
+ /** return class name. */
+ itkTypeMacro(BayesianFusionFilter, FusionImageBase);
+
+ /** Supported images definition. */
+ typedef typename InputMultiSpectralImageType::PixelType InputMultiSpectralPixelType;
+ typedef typename InputMultiSpectralImageType::InternalPixelType InputMultiSpectralInternalPixelType;
+ typedef typename InputMultiSpectralInterpImageType::PixelType InputMultiSpectralInterpPixelType;
+ typedef typename InputMultiSpectralInterpImageType::InternalPixelType InputMultiSpectralInterpInternalPixelType;
+ typedef typename InputPanchroImageType::PixelType InputPanchroPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
+
+ /** Real class typedef definition. */
+ typedef typename itk::NumericTraits<InputPanchroPixelType>::RealType InputPanchroRealType;
+ typedef typename itk::NumericTraits<InputMultiSpectralInternalPixelType>::RealType InputMultiSpectralRealType;
+ typedef typename InputMultiSpectralImageType::RegionType InputMultiSpectralImageRegionType;
+ typedef typename itk::NumericTraits<InputMultiSpectralInterpInternalPixelType>::RealType InputMultiSpectralInterpRealType;
+ typedef typename InputMultiSpectralInterpImageType::RegionType InputMultiSpectralInterpImageRegionType;
+ typedef typename InputPanchroImageType::RegionType InputPanchroImageRegionType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+
+ /** Image size "typedef" definition. */
+ typedef typename InputMultiSpectralImageType::SizeType SizeType;
+
+
+ /** Typedef for statistic computing. */
+ typedef StreamingStatisticsVectorImageFilter<InputMultiSpectralImageType> StreamingStatisticsVectorImageFilterType;
+ typedef typename StreamingStatisticsVectorImageFilterType::MatrixType MatrixType;
+ typedef MatrixTransposeMatrixImageFilter<InputMultiSpectralImageType, InputMultiSpectralImageType> MSTransposeMSType;
+ typedef ImageToVectorImageCastFilter<InputPanchroImageType, InputMultiSpectralImageType> CasterType;
+
+
+ /** Set the ponderation value. */
+ itkSetMacro(Lambda, float);
+ /** Give the ponderation value. */
+ itkGetConstReferenceMacro(Lambda, float);
+
+ /** Set the Beta matrix. */
+ itkSetMacro(Beta, MatrixType);
+ /** Give the Beta matrix. */
+ itkGetConstReferenceMacro(Beta, MatrixType);
+
+ /** Set the Covariance matrix. */
+ itkSetMacro(CovarianceMatrix, MatrixType);
+ /** Give the Covariance matrix. */
+ itkGetConstReferenceMacro(CovarianceMatrix, MatrixType);
+
+ /** Set the Covariance inverse matrix. */
+ itkSetMacro(CovarianceInvMatrix, MatrixType);
+ /** Give the Covariance inverse matrix. */
+ itkGetConstReferenceMacro(CovarianceInvMatrix, MatrixType);
+
+ /** Set the Bayesian Data Fusion matrix. */
+ itkSetMacro(Vcondopt, MatrixType);
+ /** Give the Bayesian Data Fusion matrix. */
+ itkGetConstReferenceMacro(Vcondopt, MatrixType);
+
+ /** Set the S coefficient. */
+ itkSetMacro(S, float);
+ /** Give the S coefficient. */
+ itkGetConstReferenceMacro(S, float);
+
+
+protected:
+ BayesianFusionFilter();
+ virtual ~BayesianFusionFilter();
+ /** Initialize some accumulators before the threads run. */
+ void BeforeThreadedGenerateData ();
+
+
+private:
+ /** Ponderation declaration*/
+ float m_Lambda;
+ float m_S;
+ /** Multispectral covariance matrix */
+ MatrixType m_CovarianceMatrix;
+ /** Multispectral inverse covariance matrix */
+ MatrixType m_CovarianceInvMatrix;
+ /** Regression coefficients matrix */
+ MatrixType m_Beta;
+ /** Optimisation matrix */
+ MatrixType m_Vcondopt;
+
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbBayesianFusionFilter.txx"
+#endif
+
+#endif
diff --git a/Code/Fusion/otbBayesianFusionFilter.txx b/Code/Fusion/otbBayesianFusionFilter.txx
new file mode 100755
index 0000000000..60f991d370
--- /dev/null
+++ b/Code/Fusion/otbBayesianFusionFilter.txx
@@ -0,0 +1,277 @@
+/*=========================================================================
+
+ 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.
+
+=========================================================================*/
+#include "otbBayesianFusionFilter.h"
+
+
+
+namespace otb
+{
+
+ template <class TInputMultiSpectralImage,
+ class TInputMultiSpectralInterpImage,
+ class TInputPanchroImage,
+ class TOutputImage>
+ BayesianFusionFilter<TInputMultiSpectralImage,
+ TInputMultiSpectralInterpImage,
+ TInputPanchroImage,
+ TOutputImage>
+ ::BayesianFusionFilter()
+ {
+ m_Lambda = 0.9999;
+ m_S = 1;
+ }
+
+
+ template <class TInputMultiSpectralImage,
+ class TInputMultiSpectralInterpImage,
+ class TInputPanchroImage,
+ class TOutputImage>
+ BayesianFusionFilter<TInputMultiSpectralImage,
+ TInputMultiSpectralInterpImage,
+ TInputPanchroImage,
+ TOutputImage>
+ ::~BayesianFusionFilter()
+ {
+
+ }
+
+
+ template <class TInputMultiSpectralImage,
+ class TInputMultiSpectralInterpImage,
+ class TInputPanchroImage,
+ class TOutputImage>
+ void
+ BayesianFusionFilter<TInputMultiSpectralImage,
+ TInputMultiSpectralInterpImage,
+ TInputPanchroImage,
+ TOutputImage>
+ ::BeforeThreadedGenerateData ()
+ {
+ // Allocate output
+ typename OutputImageType::Pointer output = this->GetOutput();
+ typename InputMultiSpectralImageType::ConstPointer multiSpec = this->GetMultiSpect();
+ typename InputMultiSpectralInterpImageType::ConstPointer multiSpecInterp = this->GetMultiSpectInterp();
+ typename InputPanchroImageType::ConstPointer panchro = this->GetPanchro();
+
+ /** Variable Initialisaton */
+ m_Beta.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel()+1, 1);
+ m_Beta.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+
+ m_CovarianceMatrix.SetSize( multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel() );
+ m_CovarianceMatrix.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+
+ m_CovarianceInvMatrix.SetSize( multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel() );
+ m_CovarianceInvMatrix.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+
+ m_Vcondopt.SetSize( multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel() );
+ m_Vcondopt.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+
+
+ /** Compute the inverse of the multispectral interpolated image covariance matrix */
+ typename StreamingStatisticsVectorImageFilterType::Pointer covComputefilter = StreamingStatisticsVectorImageFilterType::New();
+ covComputefilter->SetInput(multiSpecInterp);
+ covComputefilter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
+ covComputefilter->SetNumberOfStreamDivisions(200);
+ covComputefilter->Update();
+ MatrixType m_CovarianceMatrix = covComputefilter->GetCovariance();
+ m_CovarianceInvMatrix = m_CovarianceMatrix.GetInverse();
+
+ /** Beta computation : Regression model coefficient */
+ // MatrixTransform only support vectorimage input
+ typename CasterType::Pointer caster = CasterType::New();
+ caster->SetInput(panchro);
+ caster->Update();
+ // Compute the transpose multispectral image multiplied by itself
+ typename MSTransposeMSType::Pointer msTransposeMs = MSTransposeMSType::New();
+ // Compute the transpose multispectral image multiplied by the panchromatic one
+ typename MSTransposeMSType::Pointer msTransposePan = MSTransposeMSType::New();
+
+ // Add a dimension filled with ones to the images
+ msTransposeMs->SetUsePadFirstInput(true);
+ msTransposeMs->SetUsePadSecondInput(true);
+ msTransposePan->SetUsePadFirstInput(true);
+
+ msTransposeMs->SetFirstInput(multiSpec);
+ msTransposeMs->SetSecondInput(multiSpec);
+
+ msTransposePan->SetFirstInput(multiSpec);
+ msTransposePan->SetSecondInput( caster->GetOutput() );
+
+ msTransposeMs->Update();
+ msTransposePan->Update();
+
+ MatrixType temp;
+ temp = msTransposeMs->GetResultOutput()->Get().GetInverse();
+ m_Beta = temp*msTransposePan->GetResultOutput()->Get();
+
+ // S computation : quadratique mean of the regression residue
+ // Compute the transpose panchromatic image multiplied by itself
+ typename MSTransposeMSType::Pointer panTransposePan = MSTransposeMSType::New();
+ panTransposePan->SetFirstInput(caster->GetOutput());
+ panTransposePan->SetSecondInput(caster->GetOutput());
+ panTransposePan->Update();
+
+ MatrixType S, tempS, tempS2;
+ S = panTransposePan->GetResultOutput()->Get();
+ tempS = msTransposePan->GetResultOutput()->Get().GetTranspose();
+ tempS = tempS*m_Beta;
+
+ /** TODO
+ * A modifier en utilisant l'opérateur - de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * S = S-tempS;
+ **/
+ if( (S.Rows() != tempS.Rows()) || (S.Cols() != tempS.Cols()))
+ {
+ itkExceptionMacro( << "Matrix with size (" << S.Rows() << "," <<
+ S.Cols() << ") cannot be subtracted from matrix with size (" <<
+ tempS.Rows() << "," << tempS.Cols() <<" )" );
+ }
+ for( unsigned int r=0; r<S.Rows(); r++)
+ {
+ for( unsigned int c=0; c<S.Cols(); c++ )
+ {
+ S(r,c) -= tempS(r,c);
+ }
+ }
+ //**** END TODO ****/
+
+ tempS = m_Beta.GetTranspose();
+ tempS2 = msTransposePan->GetResultOutput()->Get();
+ tempS = tempS*tempS2;
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * S = S-tempS;
+ **/
+ if( (S.Rows() != tempS.Rows()) || (S.Cols() != tempS.Cols()) )
+ {
+ itkExceptionMacro( << "Matrix with size (" << S.Rows() << "," <<
+ S.Cols() << ") cannot be subtracted from matrix with size (" <<
+ tempS.Rows() << "," << tempS.Cols() << " )" );
+
+ }
+
+ for( unsigned int r=0; r<S.Rows(); r++)
+ {
+ for( unsigned int c=0; c<S.Cols(); c++ )
+ {
+ S(r,c) -= tempS(r,c);
+ }
+ }
+ //**** END TODO ****/
+
+
+ MatrixType xxT, xxTb, xxTbT, xxTbTb;
+ xxT = msTransposeMs->GetResultOutput()->Get().GetTranspose();
+ xxTb = xxT*m_Beta;
+ xxTbT = xxTb.GetTranspose();
+ xxTbTb = xxTbT*m_Beta;
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * S = S-xxTbTb;
+ **/
+ if( (S.Cols() != xxTbTb.Cols()) || (S.Cols() != xxTbTb.Cols()) )
+ {
+ itkExceptionMacro( << "Matrix with size (" << S.Rows() << "," <<
+ S.Cols() << ") cannot be subtracted from matrix with size (" <<
+ xxTbTb.Rows() << "," << xxTbTb.Cols() << " )" );
+ }
+
+ for( unsigned int r=0; r<S.Rows(); r++)
+ {
+ for( unsigned int c=0; c<S.Cols(); c++ )
+ {
+ S(r,c) += xxTbTb(r,c);
+ }
+ }
+ //**** END TODO ****/
+
+ unsigned int size1 = multiSpec->GetLargestPossibleRegion().GetSize()[0]*multiSpec->GetLargestPossibleRegion().GetSize()[1];
+ unsigned int size2 = multiSpec->GetNumberOfComponentsPerPixel()+1;
+ m_S = S(0,0);
+ m_S /= static_cast<float>(size1-size2);
+
+
+ // cutBeta is the N-1 last m_Beta element matrix.
+ // varPan contains transpose(cutBeta)*cutBeta/S
+ MatrixType varPan, cutBeta;
+ varPan.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), 1);
+ varPan.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+ cutBeta.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), 1);
+ cutBeta.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+ // Take the N-1 m_Beta last elements
+ for( unsigned int r=1; r<m_Beta.Rows(); r++ )
+ {
+ cutBeta(r-1,0) = m_Beta(r,0);
+ }
+ varPan = cutBeta;
+
+ MatrixType tempvarPan;
+ tempvarPan = varPan.GetTranspose();
+ varPan *= tempvarPan;
+ varPan /= m_S;
+
+
+ // Compute the optimization matrix : m_Vcondopt
+ // eye is the identical matrix which size is the number of components of the multispectral image
+ MatrixType eye;
+ eye.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel());
+ eye.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
+ for( unsigned int r=1; r<eye.Rows(); r++)
+ {
+ eye(r,r) = vcl_pow(10., -12.);
+ }
+
+ /** TODO
+ * A modifier en utilisant l'opérateur + de la méthode. Pour le moment probleme avec exceptionmacro (pas de GetClassName...)
+ * m_Vcondopt = 2 *m_Lambda*varPan+2*m_CovarianceInvMatrix*(1-m_Lambda)+eye;
+ **/
+ if( (m_Vcondopt.Cols() != varPan.Cols()) || (m_Vcondopt.Cols() != varPan.Cols())
+ || (m_Vcondopt.Cols() != m_CovarianceInvMatrix.Cols()) || (m_Vcondopt.Cols() != m_CovarianceInvMatrix.Cols()))
+ {
+ itkExceptionMacro( << "Matrix with size (" << m_Vcondopt.Rows() << "," <<
+ m_Vcondopt.Cols() << ") cannot be subtracted from matrix with size (" <<
+ varPan.Rows() << "," << varPan.Cols() << " ) or ( " <<
+ m_CovarianceInvMatrix.Rows() << "," << m_CovarianceInvMatrix.Cols()<<")" );
+ }
+ for( unsigned int r=0; r<m_Vcondopt.Rows(); r++)
+ {
+ for( unsigned int c=0; c<m_Vcondopt.Cols(); c++ )
+ {
+ m_Vcondopt(r,c) = 2 *m_Lambda*varPan(r,c)
+ +2*m_CovarianceInvMatrix(r,c)*(1-m_Lambda)
+ +eye(r,c);
+ }
+ }
+ //**** END TODO ****/
+ m_Vcondopt = m_Vcondopt.GetInverse();
+
+
+ // Functor initialisation
+ this->GetFunctor().SetVcondopt(m_Vcondopt);
+ this->GetFunctor().SetBeta(cutBeta);
+ this->GetFunctor().SetAlpha(m_Beta(0,0));
+ this->GetFunctor().SetCovarianceInvMatrix(m_CovarianceInvMatrix);
+ this->GetFunctor().SetLambda(m_Lambda);
+ this->GetFunctor().SetS(m_S);
+ }
+
+
+} // end namespace otb
diff --git a/Code/Fusion/otbFusionImageBase.h b/Code/Fusion/otbFusionImageBase.h
new file mode 100755
index 0000000000..ae9aa10a6f
--- /dev/null
+++ b/Code/Fusion/otbFusionImageBase.h
@@ -0,0 +1,119 @@
+/*=========================================================================
+
+ Program: Insight Segmentation & Registration Toolkit
+ Module: $RCSfile: itkAddImageFilter.h,v $
+ Language: C++
+ Date: $Date: 2006/03/24 16:03:15 $
+ Version: $Revision: 1.11 $
+
+ Copyright (c) Insight Software Consortium. All rights reserved.
+ See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __otbFusionImageBase_h
+#define __otbFusionImageBase_h
+
+#include "itkTernaryFunctorImageFilter.h"
+#include "itkNumericTraits.h"
+
+
+namespace otb
+{
+ /** \class FusionImageBase
+ * Basic class for every Fusion classes.
+ * \sa TernaryFunctorImageFilter
+ */
+ template <class TInputMultiSpectralImage, class TInputMultiSpectralInterpImage, class TInputPanchroImage, class TOutputImage, class TFunctor>
+ class ITK_EXPORT FusionImageBase : public itk::TernaryFunctorImageFilter<TInputMultiSpectralImage, TInputMultiSpectralInterpImage, TInputPanchroImage, TOutputImage, TFunctor>
+ {
+ public:
+ /** Extract input and output images dimensions.*/
+ itkStaticConstMacro( InputImageDimension,unsigned int, TInputMultiSpectralImage::ImageDimension);
+ itkStaticConstMacro( OutputImageDimension,unsigned int, TOutputImage::ImageDimension);
+
+ /** "typedef" to simplify the variables definition and the declaration. */
+ typedef TInputMultiSpectralImage InputMultiSpectralImageType;
+ typedef TInputMultiSpectralInterpImage InputMultiSpectralInterpImageType;
+ typedef TInputPanchroImage InputPanchroImageType;
+ typedef TOutputImage OutputImageType;
+ typedef TFunctor FunctorType;
+
+ /** "typedef" for standard classes. */
+ typedef FusionImageBase Self;
+ typedef itk::TernaryFunctorImageFilter< InputMultiSpectralImageType,
+ InputMultiSpectralInterpImageType,
+ InputPanchroImageType,
+ OutputImageType,
+ FunctorType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** object factory method. */
+ itkNewMacro(Self);
+
+ /** return class name. */
+ itkTypeMacro(FusionImageBase, TernaryFunctorImageFilter);
+
+ /** Supported images definition. */
+ typedef typename InputMultiSpectralImageType::PixelType InputMultiSpectralPixelType;
+ typedef typename InputMultiSpectralInterpImageType::PixelType InputMultiSpectralInterpPixelType;
+ typedef typename InputPanchroImageType::PixelType InputPanchroPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ /** Real class typedef definition. */
+ typedef typename itk::NumericTraits<InputMultiSpectralPixelType>::RealType InputMultiSpectralRealType;
+ typedef typename itk::NumericTraits<InputMultiSpectralInterpPixelType>::RealType InputMultiSpectralInterpRealType;
+ typedef typename itk::NumericTraits<InputPanchroPixelType>::RealType InputPanchroRealType;
+ typedef typename InputMultiSpectralImageType::RegionType InputMultiSpectralImageRegionType;
+ typedef typename InputMultiSpectralInterpImageType::RegionType InputMultiSpectralInterpImageRegionType;
+ typedef typename InputPanchroImageType::RegionType InputPanchroImageRegionType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+ /** Image size "typedef" definition. */
+ typedef typename InputMultiSpectralImageType::SizeType SizeType;
+
+
+ void SetMultiSpect(const InputMultiSpectralImageType *multiSpect){ this->SetInput1( multiSpect ); };
+ void SetMultiSpectInterp(const InputMultiSpectralInterpImageType *multiSpectInterp){ this->SetInput2( multiSpectInterp ); };
+ void SetPanchro(const InputPanchroImageType *panchro){ this->SetInput3( panchro ); };
+
+ const InputMultiSpectralImageType* GetMultiSpect()
+ {
+ if( this->GetNumberOfInputs() < 1 )
+ {
+ return 0;
+ }
+ else
+ return( static_cast<const InputMultiSpectralImageType *>(this->itk::ProcessObject::GetInput(0)) );
+ }
+
+ const InputMultiSpectralInterpImageType* GetMultiSpectInterp()
+ {
+ if( this->GetNumberOfInputs() < 2 )
+ {
+ return 0;
+ }
+ else
+ return( static_cast<const InputMultiSpectralInterpImageType *>(this->itk::ProcessObject::GetInput(1)) );
+ }
+
+ const InputPanchroImageType* GetPanchro()
+ {
+ if( this->GetNumberOfInputs() < 3 )
+ {
+ return 0;
+ }
+ else
+ return( static_cast<const InputPanchroImageType *>(this->itk::ProcessObject::GetInput(2)) );
+ }
+
+
+ };
+
+} // end namespace otb
+
+
+#endif
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 807f6d6ee7..6dd56d127b 100755
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -379,6 +379,41 @@ ADD_TEST(StreamingResampleImageFilterCompareWithITK ${BASICFILTERS_TESTS}
#)
+
+# ------- otb::StreamingStatisticsVectorImageFilter ----------------------------
+
+ADD_TEST(bfTuStreamingStatisticsVectorImageFilterNew ${BASICFILTERS_TESTS}
+ otbStreamingStatisticsVectorImageFilterNew)
+
+ADD_TEST(bfTvStreamingStatisticsVectorImageFilter ${BASICFILTERS_TESTS}
+ --compare-ascii ${TOL}
+ ${BASELINE_FILES}/bfStreamingStatisticsVectorImageFilterResults.txt
+ ${TEMP}/bfTvStreamingStatisticsVectorImageFilterResults.txt
+ otbStreamingStatisticsVectorImageFilter
+ #${IMAGEDATA}/SPOT5_SCENE01/IMAGERY.TIF
+ ${INPUTDATA}/sbuv.png
+ ${TEMP}/bfTvStreamingStatisticsVectorImageFilterResults.txt
+ )
+
+
+# ------- otb::MatrixTransposeMatrixImageFilter ----------------------------
+
+ADD_TEST(bfTuMatrixTransposeMatrixImageFilterNew ${BASICFILTERS_TESTS}
+ otbMatrixTransposeMatrixImageFilterNew)
+
+ADD_TEST(bfTvMatrixTransposeMatrixImageFilter ${BASICFILTERS_TESTS}
+ --compare-ascii ${TOL}
+ ${BASELINE_FILES}/bfTvMatrixTransposeMatrixImageFilterResults.txt
+ ${TEMP}/bfMatrixTransposeMatrixImageFilterResults.txt
+ otbMatrixTransposeMatrixImageFilter
+ ${INPUTDATA}/multiSpect.tif
+ ${INPUTDATA}/multiSpectInterp.tif
+ ${TEMP}/bfTvMatrixTransposeMatrixImageFilterResults.txt
+ )
+
+
+
+
# A enrichir
SET(BasicFilters_SRCS
otbLeeFilter.cxx
@@ -424,6 +459,10 @@ otbVectorImageTo3DScalarImageFilter.cxx
otbStreamingResampleImageFilterNew.cxx
otbStreamingResampleImageFilter.cxx
otbStreamingResampleImageFilterCompareWithITK.cxx
+otbStreamingStatisticsVectorImageFilterNew.cxx
+otbStreamingStatisticsVectorImageFilter.cxx
+otbMatrixTransposeMatrixImageFilterNew.cxx
+otbMatrixTransposeMatrixImageFilter.cxx
)
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
index 4cc88ae3de..d169dfc6ff 100755
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
@@ -70,4 +70,8 @@ REGISTER_TEST(otbVectorImageTo3DScalarImageFilter);
REGISTER_TEST(otbStreamingResampleImageFilterNew);
REGISTER_TEST(otbStreamingResampleImageFilter);
REGISTER_TEST(otbStreamingResampleImageFilterCompareWithITK);
+REGISTER_TEST(otbStreamingStatisticsVectorImageFilterNew);
+REGISTER_TEST(otbStreamingStatisticsVectorImageFilter);
+REGISTER_TEST(otbMatrixTransposeMatrixImageFilterNew);
+REGISTER_TEST(otbMatrixTransposeMatrixImageFilter);
}
diff --git a/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
new file mode 100644
index 0000000000..755a763b5d
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
@@ -0,0 +1,86 @@
+/*=========================================================================
+
+ 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 "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbMatrixTransposeMatrixImageFilter.h"
+
+
+int otbMatrixTransposeMatrixImageFilter( int argc, char * argv[] )
+{
+ try
+ {
+
+ const char * infname1 = argv[1];
+ const char * infname2 = argv[2];
+ const char * outfname = argv[3];
+
+ const unsigned int Dimension = 2;
+ typedef unsigned char InputPixelType;
+ typedef unsigned char OutputPixelType;
+
+ typedef otb::VectorImage<InputPixelType,Dimension> InputImage1Type;
+ typedef otb::VectorImage<InputPixelType,Dimension> InputImage2Type;
+ typedef otb::MatrixTransposeMatrixImageFilter<InputImage1Type, InputImage2Type > MatrixTransposeMatrixImageFilterType;
+ typedef otb::ImageFileReader<InputImage1Type> ReaderType1;
+ typedef otb::ImageFileReader<InputImage2Type> ReaderType2;
+
+ // Instantiation
+ MatrixTransposeMatrixImageFilterType::Pointer filter = MatrixTransposeMatrixImageFilterType::New();
+ ReaderType1::Pointer reader1 = ReaderType1::New();
+ ReaderType2::Pointer reader2 = ReaderType2::New();
+
+ reader1->SetFileName(infname1);
+ reader2->SetFileName(infname2);
+
+ filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
+ filter->SetNumberOfStreamDivisions(200);
+ filter->SetFirstInput(reader1->GetOutput());
+ filter->SetSecondInput(reader2->GetOutput());
+ filter->SetUsePadFirstInput(true);
+ filter->SetUsePadSecondInput(true);
+
+ filter->Update();
+
+ std::ofstream file;
+ file.open(outfname);
+ file<<"transpose : "<<filter->GetResult()<<std::endl;
+ file.close();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+ catch( ... )
+ {
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilterNew.cxx b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilterNew.cxx
new file mode 100644
index 0000000000..0895616f67
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilterNew.cxx
@@ -0,0 +1,59 @@
+/*=========================================================================
+
+ 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 "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbMatrixTransposeMatrixImageFilter.h"
+
+
+int otbMatrixTransposeMatrixImageFilterNew( int argc, char * argv[] )
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef unsigned char InputPixelType;
+ typedef unsigned char OutputPixelType;
+
+ typedef otb::VectorImage<InputPixelType,Dimension> InputImage1Type;
+ typedef otb::VectorImage<InputPixelType,Dimension> InputImage2Type;
+ typedef otb::MatrixTransposeMatrixImageFilter<InputImage1Type, InputImage2Type > MatrixTransposeMatrixImageFilterType;
+
+ // Instantiation
+ MatrixTransposeMatrixImageFilterType::Pointer filter = MatrixTransposeMatrixImageFilterType::New();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+ catch( ... )
+ {
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
new file mode 100755
index 0000000000..c86b02e690
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
@@ -0,0 +1,77 @@
+/*=========================================================================
+
+ 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 "otbStreamingStatisticsVectorImageFilter.h"
+#include "otbImageFileReader.h"
+#include "otbVectorImage.h"
+#include <fstream>
+#include "otbStreamingTraits.h"
+
+int otbStreamingStatisticsVectorImageFilter(int argc, char * argv[])
+{
+ try
+ {
+
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+ typedef otb::VectorImage<PixelType,Dimension> ImageType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef otb::StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType;
+
+ // Instantiating object
+ StreamingStatisticsVectorImageFilterType::Pointer filter = StreamingStatisticsVectorImageFilterType::New();
+
+
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(infname);
+
+ filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
+ filter->SetNumberOfStreamDivisions(200);
+ filter->SetInput(reader->GetOutput());
+ filter->Update();
+
+
+ std::ofstream file;
+ file.open(outfname);
+ file<<"Minimum: "<<filter->GetMinimum()<<std::endl;
+ file<<"Maximum: "<<filter->GetMaximum()<<std::endl;
+ file<<"Sum: "<<filter->GetSum()<<std::endl;
+ file<<"Mean: "<<filter->GetMean()<<std::endl;
+ file<<"Covariance: "<<filter->GetCovariance()<<std::endl;
+ file.close();
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilterNew.cxx b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilterNew.cxx
new file mode 100644
index 0000000000..e8ca952935
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilterNew.cxx
@@ -0,0 +1,50 @@
+/*=========================================================================
+
+ 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 "otbStreamingStatisticsVectorImageFilter.h"
+#include "otbVectorImage.h"
+
+int otbStreamingStatisticsVectorImageFilterNew(int argc, char * argv[])
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef unsigned char PixelType;
+
+ typedef otb::VectorImage<PixelType,Dimension> ImageType;
+ typedef otb::StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType;
+
+ // Instantiating object
+ StreamingStatisticsVectorImageFilterType::Pointer object = StreamingStatisticsVectorImageFilterType::New();
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/CMakeLists.txt b/Testing/Code/CMakeLists.txt
index 803a321f39..2168fb5bad 100644
--- a/Testing/Code/CMakeLists.txt
+++ b/Testing/Code/CMakeLists.txt
@@ -12,6 +12,7 @@ DisparityMap
SpatialReasoning
Projections
Radiometry
+Fusion
)
IF(OTB_USE_VISU_GUI)
diff --git a/Testing/Code/Fusion/CMakeLists.txt b/Testing/Code/Fusion/CMakeLists.txt
new file mode 100755
index 0000000000..b72f6fd187
--- /dev/null
+++ b/Testing/Code/Fusion/CMakeLists.txt
@@ -0,0 +1,62 @@
+
+IF( NOT OTB_DISABLE_CXX_TESTING )
+
+SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images)
+SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files)
+SET(INPUTDATA ${OTB_DATA_ROOT}/Input)
+#Images de teledetection (grosses images )
+SET(IMAGEDATA ${OTB_DATA_ROOT}/LargeInput )
+SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary)
+
+#Tolerance de base
+SET(TOL 0.0)
+
+#Tolerance sur la difference des valeurs numeriques dans le --compare-ascii
+SET(EPS 0.001)
+
+#Tolerance sur diff pixel image
+#EPSILON est different de 0.0 pour les tests multiplateformes avec différentes options de compilation.
+SET(EPSILON 0.00000001)
+
+SET(FUSION_TESTS ${CXX_TEST_PATH}/otbFusionTests)
+
+
+# ------- otb::ImageFusionBase ------------------------------
+ADD_TEST(fuTuFusionImageBaseNew ${FUSION_TESTS}
+ otbFusionImageBaseNew
+)
+# ------- otb::BayesianFusionFilter ------------------------------
+ADD_TEST(fuTuBayesianFusionFilterNew ${FUSION_TESTS}
+ otbBayesianFusionFilterNew
+)
+
+ADD_TEST(fuTvBayesianFusionFilter ${FUSION_TESTS}
+ --compare-image ${TOL} ${BASELINE}/fuTvBayesianFusion.tif
+ ${TEMP}/fuTvBayesianFusion.tif
+ otbBayesianFusionFilter
+ ${INPUTDATA}/multiSpect.tif
+ ${INPUTDATA}/multiSpectInterp.tif
+ ${INPUTDATA}/panchro.tif
+ ${TEMP}/fuTvBayesianFusion.tif
+)
+
+
+# A enrichir
+SET(Fusion_SRCS
+otbFusionImageBaseNew.cxx
+otbBayesianFusionFilterNew.cxx
+otbBayesianFusionFilter.cxx
+)
+
+
+INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
+
+ADD_EXECUTABLE(otbFusionTests otbFusionTests.cxx ${Fusion_SRCS})
+TARGET_LINK_LIBRARIES(otbFusionTests OTBFusion OTBCommon OTBIO gdal ITKIO ITKCommon ITKAlgorithms)
+
+
+ENDIF( NOT OTB_DISABLE_CXX_TESTING )
+
+
+
+
diff --git a/Testing/Code/Fusion/otbBayesianFusionFilter.cxx b/Testing/Code/Fusion/otbBayesianFusionFilter.cxx
new file mode 100755
index 0000000000..625b3ed3a4
--- /dev/null
+++ b/Testing/Code/Fusion/otbBayesianFusionFilter.cxx
@@ -0,0 +1,83 @@
+/*=========================================================================
+
+ 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 "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbVectorImage.h"
+#include "otbImage.h"
+
+#include "otbBayesianFusionFilter.h"
+
+int otbBayesianFusionFilter( int argc, char * argv[] )
+{
+ try
+ {
+ const char * multispect = argv[1];
+ const char * multispectinterp = argv[2];
+ const char * panchro = argv[3];
+ const char * output = argv[4];
+
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+
+ typedef otb::VectorImage<PixelType,Dimension> VectorImageType;
+ typedef otb::Image<PixelType,Dimension> PanchroImageType;
+ typedef otb::ImageFileReader<VectorImageType> VectorReaderType;
+ typedef otb::ImageFileReader<PanchroImageType> ImageReaderType;
+ typedef otb::ImageFileWriter<VectorImageType> VectorImageWriterType;
+ typedef otb::BayesianFusionFilter<VectorImageType, VectorImageType, PanchroImageType, VectorImageType> FilterType;
+
+ VectorReaderType::Pointer multiSpectReader = VectorReaderType::New();
+ VectorReaderType::Pointer multiSpectInterpReader = VectorReaderType::New();
+ ImageReaderType::Pointer panchroReader = ImageReaderType::New();
+ FilterType::Pointer filter = FilterType::New();
+ VectorImageWriterType::Pointer writer = VectorImageWriterType::New();
+
+ multiSpectReader->SetFileName(multispect);
+ multiSpectInterpReader->SetFileName(multispectinterp);
+ panchroReader->SetFileName(panchro);
+
+ filter->SetMultiSpect(multiSpectReader->GetOutput());
+ filter->SetMultiSpectInterp(multiSpectInterpReader->GetOutput());
+ filter->SetPanchro(panchroReader->GetOutput());
+ writer->SetInput( filter->GetOutput() );
+ writer->SetFileName(output);
+ writer->Update();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+ catch( ... )
+ {
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/Fusion/otbBayesianFusionFilterNew.cxx b/Testing/Code/Fusion/otbBayesianFusionFilterNew.cxx
new file mode 100644
index 0000000000..614056dfbd
--- /dev/null
+++ b/Testing/Code/Fusion/otbBayesianFusionFilterNew.cxx
@@ -0,0 +1,63 @@
+/*=========================================================================
+
+ 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 "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbVectorImage.h"
+#include "otbImage.h"
+
+#include "otbBayesianFusionFilter.h"
+
+int otbBayesianFusionFilterNew( int argc, char * argv[] )
+{
+ try
+ {
+
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+
+ typedef otb::VectorImage<PixelType,Dimension> VectorImageType;
+ typedef otb::Image<PixelType,Dimension> PanchroImageType;
+ typedef otb::ImageFileReader<VectorImageType> VectorReaderType;
+ typedef otb::ImageFileReader<PanchroImageType> ImageReaderType;
+ typedef otb::BayesianFusionFilter<VectorImageType, VectorImageType, PanchroImageType, VectorImageType> FilterType;
+
+ FilterType::Pointer filter = FilterType::New();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+ catch( ... )
+ {
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/Fusion/otbFusionImageBaseNew.cxx b/Testing/Code/Fusion/otbFusionImageBaseNew.cxx
new file mode 100644
index 0000000000..7b983ea2b5
--- /dev/null
+++ b/Testing/Code/Fusion/otbFusionImageBaseNew.cxx
@@ -0,0 +1,86 @@
+/*=========================================================================
+
+ 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 "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbFusionImageBase.h"
+
+namespace Functor
+{
+ template <class TInputMultiSpectral, class TInputMultiSpectralInterp, class TInputPanchro, class TOutput>
+ class NewFunctorTest
+ {
+ public:
+ NewFunctorTest() {};
+ ~NewFunctorTest() {};
+
+ inline TOutput operator() (const TInputMultiSpectral & A, const TInputMultiSpectralInterp & B, const TInputPanchro & C)
+ {
+ return(static_cast<TOutput>(A[0]) + static_cast<TOutput>(B[0]) + static_cast<TOutput>(C));
+
+ }
+ };
+}
+
+int otbFusionImageBaseNew( int argc, char * argv[] )
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef unsigned char InputPixelType;
+ typedef unsigned char OutputPixelType;
+
+ typedef otb::Image<InputPixelType,Dimension> InputPanchroImageType;
+ typedef otb::VectorImage<InputPixelType,Dimension> InputMultiSpectralImageType;
+ typedef otb::VectorImage<InputPixelType,Dimension> InputMultiSpectralInterpImageType;
+ typedef otb::Image<OutputPixelType,Dimension> OutputImageType;
+
+ typedef otb::FusionImageBase<InputMultiSpectralImageType,
+ InputMultiSpectralInterpImageType,
+ InputPanchroImageType,
+ OutputImageType,
+ Functor::NewFunctorTest<InputMultiSpectralImageType::PixelType,
+ InputMultiSpectralInterpImageType::PixelType,
+ InputPanchroImageType::PixelType,
+ OutputImageType::PixelType>
+ > FusionImageBaseType;
+
+ // Instantiation
+ FusionImageBaseType::Pointer base = FusionImageBaseType::New();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+ catch( ... )
+ {
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/Fusion/otbFusionTests.cxx b/Testing/Code/Fusion/otbFusionTests.cxx
new file mode 100755
index 0000000000..6fa5dc19c4
--- /dev/null
+++ b/Testing/Code/Fusion/otbFusionTests.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.
+
+=========================================================================*/
+
+// this file defines the otbCommonTest for the test driver
+// and all it expects is that you have a function called RegisterTests
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include <iostream>
+#include "otbTestMain.h"
+
+
+void RegisterTests()
+{
+REGISTER_TEST(otbBayesianFusionFilterNew);
+REGISTER_TEST(otbBayesianFusionFilter);
+REGISTER_TEST(otbFusionImageBaseNew);
+}
diff --git a/otbIncludeDirectories.cmake b/otbIncludeDirectories.cmake
index 4f4f0ee87f..96af8bb6fc 100644
--- a/otbIncludeDirectories.cmake
+++ b/otbIncludeDirectories.cmake
@@ -20,6 +20,7 @@ SET(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
${OTB_SOURCE_DIR}/Code/DisparityMap
${OTB_SOURCE_DIR}/Code/Visu
${OTB_SOURCE_DIR}/Code/Gui
+ ${OTB_SOURCE_DIR}/Code/Fusion
${OTB_SOURCE_DIR}/Code/Projections
${OTB_SOURCE_DIR}/Code/Radiometry
${OTB_SOURCE_DIR}/Utilities/BGL
--
GitLab