ADD: add a GPU implementation of StreamingStatisticsVectorImageFilter based on CuBLAS

CMakeLists.txt

@@ -11,6 +11,8 @@ ENDIF(NOT EXECUTABLE_OUTPUT_PATH)
 
 INCLUDE_REGULAR_EXPRESSION("^.*$")
 
+OPTION(BUILD_SHARED_LIBS "Build OTB with shared libraries." ON)
+
 # Find OTB and load its settings.
 FIND_PACKAGE(OTB)
   IF(OTB_FOUND)
@@ -19,8 +21,16 @@ FIND_PACKAGE(OTB)
     MESSAGE(FATAL_ERROR "OTB not found.  Please set OTB_DIR")
 ENDIF(OTB_FOUND)
 
+OPTION(OTB_USE_CUBLAS "Use Cuda CUBLAS library." ON)
+IF(OTB_USE_CUBLAS)
+  FIND_PACKAGE(CUDA)
+  include_directories(${CUDA_INCLUDE_DIRS})
+ENDIF(OTB_USE_CUBLAS)
+
 INCLUDE_DIRECTORIES(
 ${Hyperspectral_SOURCE_DIR}/Code/BasicFilters
+${Hyperspectral_SOURCE_DIR}/Code/BasicFilters/GPUCuda
+${Hyperspectral_SOURCE_DIR}/Code/CudaMath
 ${Hyperspectral_SOURCE_DIR}/Code/Hyperspectral
 ${Hyperspectral_SOURCE_DIR}/Code/Vahine
 ${Hyperspectral_SOURCE_DIR}/Code/RX
@@ -31,6 +41,3 @@ ADD_SUBDIRECTORY(Application)
 
 ENABLE_TESTING()
 ADD_SUBDIRECTORY(Testing)
-
-
-

Code/BasicFilters/CMakeLists.txt

+ADD_SUBDIRECTORY(GPUCuda)

Code/BasicFilters/GPUCuda/CMakeLists.txt

+# Sources of non-templated classes.
+FILE(GLOB OTBBasicFiltersGPUCuda_SRCS "*.cxx" )
+
+ADD_LIBRARY(OTBBasicFiltersGPUCuda ${OTBBasicFiltersGPUCuda_SRCS})
+TARGET_LINK_LIBRARIES (OTBBasicFiltersGPUCuda OTBCommon OTBBasicFilters ITKBasicFilters )
+
+IF(OTB_LIBRARY_PROPERTIES)
+  SET_TARGET_PROPERTIES(OTBBasicFiltersGPUCuda PROPERTIES ${OTB_LIBRARY_PROPERTIES})
+ENDIF(OTB_LIBRARY_PROPERTIES)
+
+IF(OTB_USE_CUBLAS)
+  CUDA_ADD_CUBLAS_TO_TARGET( OTBBasicFiltersGPUCuda )
+ENDIF(OTB_USE_CUBLAS)

Code/BasicFilters/GPUCuda/otbCublasStreamingStatisticsVectorImageFilter.h

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+  Some parts of this code are derived from ITK. See ITKCopyright.txt
+  for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANT2ABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbCublasStreamingStatisticsVectorImageFilter_h
+#define __otbCublasStreamingStatisticsVectorImageFilter_h
+
+#include "otbPersistentImageFilter.h"
+#include "otbPersistentFilterStreamingDecorator.h"
+#include "itkNumericTraits.h"
+#include "itkArray.h"
+#include "itkSimpleDataObjectDecorator.h"
+#include "itkImageRegionSplitter.h"
+#include "itkVariableSizeMatrix.h"
+#include "itkVariableLengthVector.h"
+
+#include "cublas.h"
+
+namespace otb
+{
+
+/** \class PersistentCublasStreamingStatisticsVectorImageFilter
+ * \brief Compute mean, covariance & correlation of a large image using streaming
+ *
+ *  This filter persists its temporary data. It means that if you Update it n times on n different
+ * requested regions, the output statistics will be the statitics of the whole set of n regions.
+ *
+ * To reset the temporary data, one should call the Reset() function.
+ *
+ * To get the statistics once the regions have been processed via the pipeline, use the Synthetize() method.
+ *
+ * \sa PersistentImageFilter
+ * \ingroup Streamed
+ * \ingroup Multithreaded
+ * \ingroup MathematicalStatisticsImageFilters
+ *
+ */
+template<class TInputImage>
+class ITK_EXPORT PersistentCublasStreamingStatisticsVectorImageFilter :
+  public PersistentImageFilter<TInputImage, TInputImage>
+{
+public:
+  /** Standard Self typedef */
+  typedef PersistentCublasStreamingStatisticsVectorImageFilter           Self;
+  typedef PersistentImageFilter<TInputImage, TInputImage> Superclass;
+  typedef itk::SmartPointer<Self>                         Pointer;
+  typedef itk::SmartPointer<const Self>                   ConstPointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Runtime information support. */
+  itkTypeMacro(PersistentCublasStreamingStatisticsVectorImageFilter, PersistentImageFilter);
+
+  /** Image related typedefs. */
+  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;
+
+  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+
+  /** Image related typedefs. */
+  itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
+
+  /** Type to use for computations. */
+  //typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
+  typedef float 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 itk::VariableSizeMatrix<RealType>    MatrixType;
+  typedef std::vector<MatrixType>              ArrayMatrixType;
+  typedef itk::Array<long>                     ArrayLongPixelType;
+  typedef std::vector<RealPixelType>           ArrayRealPixelType;
+  typedef std::vector<PixelType>               ArrayPixelType;
+  typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
+  typedef itk::SimpleDataObjectDecorator<PixelType>     PixelObjectType;
+  typedef itk::SimpleDataObjectDecorator<MatrixType>    MatrixObjectType;
+
+  /** Return the computed Mean. */
+  RealPixelType GetMean() const
+  {
+    return this->GetMeanOutput()->Get();
+  }
+  RealPixelObjectType* GetMeanOutput();
+  const RealPixelObjectType* GetMeanOutput() const;
+
+  /** Return the computed Covariance. */
+  MatrixType GetCovariance() const
+  {
+    return this->GetCovarianceOutput()->Get();
+  }
+  MatrixObjectType* GetCovarianceOutput();
+  const MatrixObjectType* GetCovarianceOutput() const;
+
+  /** Return the computed Covariance. */
+  MatrixType GetCorrelation() const
+  {
+    return this->GetCorrelation()->Get();
+  }
+  MatrixObjectType* GetCorrelationOutput();
+  const MatrixObjectType* GetCorrelationOutput() const;
+
+  /** Make a DataObject of the correct type to be used as the specified
+   * output.
+   */
+  virtual DataObjectPointer MakeOutput(unsigned int idx);
+
+  /** Pass the input through unmodified. Do this by Grafting in the
+   *  AllocateOutputs method.
+   */
+  virtual void AllocateOutputs();
+  virtual void GenerateOutputInformation();
+  virtual void Synthetize(void);
+  virtual void Reset(void);
+
+protected:
+  PersistentCublasStreamingStatisticsVectorImageFilter();
+  virtual ~PersistentCublasStreamingStatisticsVectorImageFilter();
+
+  virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+  /** Multi-thread version GenerateData. */
+  //void  ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId);
+  void  GenerateData();
+
+private:
+  PersistentCublasStreamingStatisticsVectorImageFilter(const Self &); //purposely not implemented
+  void operator =(const Self&); //purposely not implemented
+
+  RealPixelType m_FirstOrderAccumulator;
+  MatrixType    m_SecondOrderAccumulator;
+
+  float*   m_GPUImage;
+  SizeType m_GPUImageSize;
+  float*   m_GPUSecondOrderAccumulator;
+
+}; // end of class PersistentCublasStreamingStatisticsVectorImageFilter
+
+
+/** \class CublasStreamingStatisticsVectorImageFilter
+ * \brief This class streams the whole input image through the PersistentCublasStreamingStatisticsVectorImageFilter.
+ *
+ * This way, it allows to compute the first and second order global statistics of this image. It calls the
+ * Reset() method of the PersistentCublasStreamingStatisticsVectorImageFilter before streaming the image and the
+ * Synthetize() method of the PersistentCublasStreamingStatisticsVectorImageFilter after having streamed the image
+ * to compute the statistics. The accessor on the results are wrapping the accessors of the
+ * internal PersistentCublasStreamingStatisticsVectorImageFilter.
+ *
+ * \sa PersistentCublasStreamingStatisticsVectorImageFilter
+ * \sa PersistentImageFilter
+ * \sa PersistentFilterStreamingDecorator
+ * \sa StreamingImageVirtualWriter
+ * \ingroup Streamed
+ * \ingroup Multithreaded
+ * \ingroup MathematicalStatisticsImageFilters
+ */
+
+template<class TInputImage>
+class ITK_EXPORT CublasStreamingStatisticsVectorImageFilter :
+  public PersistentFilterStreamingDecorator<PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage> >
+{
+public:
+  /** Standard Self typedef */
+  typedef CublasStreamingStatisticsVectorImageFilter Self;
+  typedef PersistentFilterStreamingDecorator
+  <PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage> > Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  /** Type macro */
+  itkNewMacro(Self);
+
+  /** Creation through object factory macro */
+  itkTypeMacro(CublasStreamingStatisticsVectorImageFilter, PersistentFilterStreamingDecorator);
+
+  typedef TInputImage                                 InputImageType;
+  typedef typename Superclass::FilterType             StatFilterType;
+  typedef typename StatFilterType::PixelType          PixelType;
+  typedef typename StatFilterType::RealType           RealType;
+  typedef typename StatFilterType::RealPixelType      RealPixelType;
+  typedef typename StatFilterType::MatrixType         MatrixType;
+  typedef typename StatFilterType::ArrayMatrixType    ArrayMatrixType;
+  typedef typename StatFilterType::ArrayLongPixelType ArrayLongPixelType;
+  typedef typename StatFilterType::ArrayRealPixelType ArrayRealPixelType;
+  typedef typename StatFilterType::ArrayPixelType     ArrayPixelType;
+
+  /** Type of DataObjects used for scalar outputs */
+  typedef typename StatFilterType::RealPixelObjectType RealPixelObjectType;
+  typedef typename StatFilterType::PixelObjectType     PixelObjectType;
+  typedef typename StatFilterType::MatrixObjectType    MatrixObjectType;
+
+  void SetInput(InputImageType * input)
+  {
+    this->GetFilter()->SetInput(input);
+  }
+  const InputImageType * GetInput()
+  {
+    return this->GetFilter()->GetInput();
+  }
+
+  /** Return the computed Mean. */
+  RealPixelType GetMean() const
+  {
+    return this->GetFilter()->GetMeanOutput()->Get();
+  }
+  RealPixelObjectType* GetMeanOutput()
+  {
+    return this->GetFilter()->GetMeanOutput();
+  }
+  const RealPixelObjectType* GetMeanOutput() const
+  {
+    return this->GetFilter()->GetMeanOutput();
+  }
+
+  /** Return the computed Covariance. */
+  MatrixType GetCovariance() const
+  {
+    return this->GetFilter()->GetCovarianceOutput()->Get();
+  }
+  MatrixObjectType* GetCovarianceOutput()
+  {
+    return this->GetFilter()->GetCovarianceOutput();
+  }
+  const MatrixObjectType* GetCovarianceOutput() const
+  {
+    return this->GetFilter()->GetCovarianceOutput();
+  }
+
+  /** Return the computed Covariance. */
+  MatrixType GetCorrelation() const
+  {
+    return this->GetFilter()->GetCorrelationOutput()->Get();
+  }
+  MatrixObjectType* GetCorrelationOutput()
+  {
+    return this->GetFilter()->GetCorrelationOutput();
+  }
+  const MatrixObjectType* GetCorrelationOutput() const
+  {
+    return this->GetFilter()->GetCorrelationOutput();
+  }
+
+protected:
+  /** Constructor */
+  CublasStreamingStatisticsVectorImageFilter() {};
+  /** Destructor */
+  virtual ~CublasStreamingStatisticsVectorImageFilter() {}
+
+private:
+  CublasStreamingStatisticsVectorImageFilter(const Self &); //purposely not implemented
+  void operator =(const Self&); //purposely not implemented
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbCublasStreamingStatisticsVectorImageFilter.txx"
+#endif
+
+#endif

Code/BasicFilters/GPUCuda/otbCublasStreamingStatisticsVectorImageFilter.txx

+/*=========================================================================
+
+ 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 __otbCublasStreamingStatisticsVectorImageFilter_txx
+#define __otbCublasStreamingStatisticsVectorImageFilter_txx
+
+#include "otbCublasStreamingStatisticsVectorImageFilter.h"
+
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkNumericTraits.h"
+#include "itkProgressReporter.h"
+#include "otbMacro.h"
+
+namespace otb
+{
+
+template<class TInputImage>
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::PersistentCublasStreamingStatisticsVectorImageFilter()
+{
+  // first output is a copy of the image, DataObject created by
+  // superclass
+  //
+  // allocate the data objects for the outputs which are
+  // just decorators around vector/matric types
+
+  this->itk::ProcessObject::SetNthOutput(1, this->MakeOutput(1).GetPointer());
+  this->itk::ProcessObject::SetNthOutput(2, this->MakeOutput(2).GetPointer());
+  this->itk::ProcessObject::SetNthOutput(3, this->MakeOutput(3).GetPointer());
+
+  cublasStatus status = cublasInit();
+  if (status != CUBLAS_STATUS_SUCCESS)
+    {
+    printf("cublasInit failed");
+    }
+
+  m_GPUImage = 0;
+  m_GPUImageSize.Fill(0);
+  m_GPUSecondOrderAccumulator = 0;
+}
+
+template<class TInputImage>
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::~PersistentCublasStreamingStatisticsVectorImageFilter()
+{
+  cublasStatus status = cublasShutdown();
+  if (status != CUBLAS_STATUS_SUCCESS)
+    {
+    printf("cublasShutdown failed");
+    }
+}
+
+template<class TInputImage>
+itk::DataObject::Pointer PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::MakeOutput(
+                                                                                                       unsigned int output)
+{
+  switch (output)
+    {
+    case 0:
+      return static_cast<itk::DataObject*> (TInputImage::New().GetPointer());
+      break;
+    case 1:
+      return static_cast<itk::DataObject*> (RealPixelObjectType::New().GetPointer());
+      break;
+    case 2:
+    case 3:
+      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 PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetMeanOutput()
+{
+  return static_cast<RealPixelObjectType*> (this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage>
+const typename PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetMeanOutput() const
+{
+  return static_cast<const RealPixelObjectType*> (this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage>
+typename PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetCorrelationOutput()
+{
+  return static_cast<MatrixObjectType*> (this->itk::ProcessObject::GetOutput(2));
+}
+
+template<class TInputImage>
+const typename PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetCorrelationOutput() const
+{
+  return static_cast<const MatrixObjectType*> (this->itk::ProcessObject::GetOutput(2));
+}
+
+template<class TInputImage>
+typename PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetCovarianceOutput()
+{
+  return static_cast<MatrixObjectType*> (this->itk::ProcessObject::GetOutput(3));
+}
+
+template<class TInputImage>
+const typename PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
+PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GetCovarianceOutput() const
+{
+  return static_cast<const MatrixObjectType*> (this->itk::ProcessObject::GetOutput(3));
+}
+
+template<class TInputImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GenerateOutputInformation()
+{
+  Superclass::GenerateOutputInformation();
+  if (this->GetInput())
+    {
+    this->GetOutput()->CopyInformation(this->GetInput());
+    this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
+
+    if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
+      {
+      this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
+      }
+    }
+}
+
+template<class TInputImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::AllocateOutputs()
+{
+  // This is commented to prevent the streaming of the whole image for the first stream strip
+  // It shall not cause any problem because the output image of this filter is not intended to be used.
+  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
+  //this->GraftOutput( image );
+  // Nothing that needs to be allocated for the remaining outputs
+}
+
+template<class TInputImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::Reset()
+{
+  TInputImage * inputPtr = const_cast<TInputImage *> (this->GetInput());
+  inputPtr->UpdateOutputInformation();
+
+  unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
+
+  RealPixelType zeroRealPixel;
+  zeroRealPixel.SetSize(numberOfComponent);
+  zeroRealPixel.Fill(itk::NumericTraits<RealType>::Zero);
+  this->GetMeanOutput()->Set(zeroRealPixel);
+
+  MatrixType zeroMatrix;
+  zeroMatrix.SetSize(numberOfComponent, numberOfComponent);
+  zeroMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+  this->GetCovarianceOutput()->Set(zeroMatrix);
+  this->GetCorrelationOutput()->Set(zeroMatrix);
+
+  m_FirstOrderAccumulator = zeroRealPixel;
+  m_SecondOrderAccumulator = zeroMatrix;
+
+  if (m_GPUSecondOrderAccumulator == 0)
+    {
+    otbMsgDevMacro( << "Allocating " << numberOfComponent << " * " << numberOfComponent << " matrix on GPU (" <<
+        numberOfComponent * numberOfComponent * sizeof(float) / 1024 / 1024 << " Mo)" );
+
+    cublasStatus status;
+    status = cublasAlloc(numberOfComponent * numberOfComponent, sizeof(float), (void**) &m_GPUSecondOrderAccumulator);
+
+    // TODO : check status, throw exception on error & clean up GPU memory
+    if (status != CUBLAS_STATUS_SUCCESS)
+      {
+      otbMsgDevMacro( "cublasAlloc m_GPUSecondOrderAccumulator failed" );
+      }
+
+    status = cublasSetMatrix(numberOfComponent, numberOfComponent, sizeof(float),
+                             m_SecondOrderAccumulator.GetVnlMatrix().data_block(), numberOfComponent,
+                             m_GPUSecondOrderAccumulator, numberOfComponent);
+
+    // TODO : check status, throw exception on error & clean up GPU memory
+    if (status != CUBLAS_STATUS_SUCCESS)
+      {
+      otbMsgDevMacro( "cublasSetMatrix m_GPUSecondOrderAccumulator failed" );
+      }
+
+    }
+
+  cublasFree(m_GPUImage);
+  cublasFree(m_GPUSecondOrderAccumulator);
+}
+
+template<class TInputImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::Synthetize()
+{
+  TInputImage * inputPtr = const_cast<TInputImage *> (this->GetInput());
+  const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
+  const unsigned int nbPixels = inputPtr->GetBufferedRegion().GetNumberOfPixels();
+
+  MatrixType cpuSecondOrderAccumulator(numberOfComponent, numberOfComponent);
+  cublasStatus status;
+  status = cublasGetMatrix(numberOfComponent, numberOfComponent, sizeof(float), m_GPUSecondOrderAccumulator,
+                           numberOfComponent, cpuSecondOrderAccumulator.GetVnlMatrix().data_block(), numberOfComponent);
+  // TODO : check status, throw exception on error & clean up GPU memory
+  if (status != CUBLAS_STATUS_SUCCESS)
+    {
+    otbMsgDevMacro( "cublasGetMatrix m_GPUSecondOrderAccumulator failed");
+    }
+
+  MatrixType cor = cpuSecondOrderAccumulator / nbPixels;
+  this->GetCorrelationOutput()->Set(cor);
+
+}
+
+template<class TInputImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TInputImage>::GenerateData()
+{
+  // Support progress methods/callbacks
+  //  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+  // Grab the input
+  InputImagePointer inputPtr = const_cast<TInputImage *> (this->GetInput());
+  const unsigned int nbPixels = inputPtr->GetBufferedRegion().GetNumberOfPixels();
+  otbMsgDevMacro( "buffered region : " << inputPtr->GetBufferedRegion() << " (" << nbPixels << " pixels)" );
+  const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
+
+  cublasStatus status;
+
+  if (inputPtr->GetBufferedRegion().GetSize() != m_GPUImageSize)
+    {
+    otbMsgDevMacro( "Allocating " << nbPixels << " pixels on GPU (" << numberOfComponent * nbPixels * sizeof(float) / 1024
+        / 1024 << " Mo)" );
+
+
+    status = cublasAlloc(numberOfComponent * nbPixels, sizeof(float), (void**) &m_GPUImage);
+
+    // TODO : check status, throw exception on error & clean up GPU memory
+    if (status != CUBLAS_STATUS_SUCCESS)
+      {
+      otbMsgDevMacro( "cublasAlloc m_GPUImage failed" );
+      }
+
+    m_GPUImageSize = inputPtr->GetBufferedRegion().GetSize();
+
+    status = cublasSetMatrix(numberOfComponent, inputPtr->GetBufferedRegion().GetNumberOfPixels(), sizeof(float),
+                             inputPtr->GetBufferPointer(), numberOfComponent, m_GPUImage, numberOfComponent);
+    if (status != CUBLAS_STATUS_SUCCESS)
+      {
+      otbMsgDevMacro( "cublasSetMatrix m_GPUImage failed" );
+      }
+
+    }
+
+  cublasSsyrk('u', 'n', numberOfComponent, inputPtr->GetBufferedRegion().GetNumberOfPixels(), 1.0f, m_GPUImage,
+              numberOfComponent, 1, m_GPUSecondOrderAccumulator, numberOfComponent);
+
+  status = cublasGetError();
+  // TODO : check status, throw exception on error & clean up GPU memory
+  if (status != CUBLAS_STATUS_SUCCESS)
+    {
+    otbMsgDevMacro( "cublasSsyrk failed" );
+    }
+
+}
+
+template<class TImage>
+void PersistentCublasStreamingStatisticsVectorImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+
+  //os << indent << "Mean: "        << this->GetMeanOutput()->Get()        << std::endl;
+  //os << indent << "Covariance: "  << this->GetCovarianceOutput()->Get()  << std::endl;
+  os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
+}
+
+} // end namespace otb
+#endif

Code/CMakeLists.txt

 ADD_SUBDIRECTORY(BasicFilters)
+ADD_SUBDIRECTORY(Hyperspectral)
 ADD_SUBDIRECTORY(Vahine)
 ADD_SUBDIRECTORY(RX)

Testing/CMakeLists.txt

@@ -16,9 +16,7 @@ SET( otbHyperTests1_SRC
 )
 
 ADD_EXECUTABLE(otbHyperTests1 ${otbHyperTests1_SRC} otbHyperTests1.cxx)
-TARGET_LINK_LIBRARIES(otbHyperTests1 OTBCommon OTBIO OTBBasicFilters OTBTesting )
-
-
+TARGET_LINK_LIBRARIES(otbHyperTests1 OTBCommon OTBIO OTBBasicFilters OTBTesting)
 
 ADD_TEST(vahineVCA1
          ${TESTEXE_DIR}/otbHyperTests1
@@ -30,7 +28,7 @@ ADD_TEST(vahineVCA1
 ADD_TEST(vahineVCA2
          ${TESTEXE_DIR}/otbHyperTests1
          vahineVCA
-         ${DATA}/SYNTHETIC/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
          ${TEMP}/synthetic_unmix_5endmembers.tif
          5)
 
@@ -49,18 +47,27 @@ ADD_TEST(bfTvStreamingStatisticsVectorImageFilter
          ${DATA}/AVIRIS/extract128/extract128_goodbands.hdr
          ${TEMP}/extract128_stats_otb.txt)
 
-ADD_TEST(bfStreamingStatisticsVectorImageFilterNew2
+ADD_TEST(bfTuStreamingStatisticsVectorImageFilterNew2
          ${TESTEXE_DIR}/otbHyperTests1
          otbStreamingStatisticsVectorImageFilter2NewTest)
 
+#ADD_TEST(bfTvStreamingStatisticsVectorImageFilter_AVIRIS512
+#         ${TESTEXE_DIR}/otbHyperTests1
+#         otbStreamingStatisticsVectorImageFilterTest
+#         ${DATA}/AVIRIS/extract512/f970619t01p02_r02_sc01_goodbands.hdr
+#         ${TEMP}/aviris512_stats_otb.txt)
+       
 ADD_TEST(bfTvStreamingStatisticsVectorImageFilter2
          ${TESTEXE_DIR}/otbHyperTests1
          otbStreamingStatisticsVectorImageFilter2Test
          ${DATA}/AVIRIS/extract128/extract128_goodbands.hdr
          ${TEMP}/extract128_stats_rt.txt)
 
-
-
+ADD_TEST(bfTvStreamingStatisticsVectorImageFilter2_AVIRIS512
+         ${TESTEXE_DIR}/otbHyperTests1
+         otbStreamingStatisticsVectorImageFilter2Test
+         ${DATA}/AVIRIS/extract512/f970619t01p02_r02_sc01_goodbands.hdr
+         ${TEMP}/aviris512_stats_rt.txt)
 
 ADD_TEST(bfTuMatrixMultiplyImageFilterNew
          ${TESTEXE_DIR}/otbHyperTests1
@@ -73,8 +80,8 @@ ADD_TEST(bfTvMatrixMultiplyImageFilter
 ADD_TEST(bfTvMatrixMultiplyImageFilterRealImageTest
          ${TESTEXE_DIR}/otbHyperTests1
          otbMatrixMultiplyImageFilterTest
-         ${DATA}/SYNTHETIC/hsi_cube.tif
-         ${DATA}/SYNTHETIC/endmembers.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/endmembers.tif
          ${TEMP}/bfMatrixMultiplyImageFilter.tif)
 
 
@@ -86,8 +93,8 @@ ADD_TEST(bfTuUnConstrainedLeastSquareImageFilterNew
 ADD_TEST(bfTvUnConstrainedLeastSquareImageFilter
          ${TESTEXE_DIR}/otbHyperTests1
          otbUnConstrainedLeastSquareImageFilterTest
-         ${DATA}/SYNTHETIC/hsi_cube.tif
-         ${DATA}/SYNTHETIC/endmembers.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/endmembers.tif
          ${TEMP}/bfUnConstrainedLeastSquareImageFilter.tif)
 
 
@@ -125,39 +132,62 @@ ADD_TEST(hyTvEigenvalueLikelihoodMaximization_AVIRIS2
 ADD_TEST(hyTvEigenvalueLikelihoodMaximization2
          ${TESTEXE_DIR}/otbHyperTests1
          otbEigenvalueLikelihoodMaximizationTest
-         ${DATA}/SYNTHETIC/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
          ${TEMP}/synthetic_elm.txt)
 
 ADD_TEST(testProjectiveProjection
          ${TESTEXE_DIR}/otbHyperTests1
          otbProjectiveProjectionTestHighSNR
-         ${DATA}/SYNTHETIC/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
          8
          ${TEMP}/hsi_cube_projectiveprojected_8.tif )
 
 
-
-
-ADD_TEST(hyVCA_Synthetic
+ADD_TEST(hyTvVCA_Synthetic
          ${TESTEXE_DIR}/otbHyperTests1
          otbVCATestHighSNR
-         ${DATA}/SYNTHETIC/hsi_cube.tif
+         ${DATA}/SYNTHETIC100/hsi_cube.tif
          ${TEMP}/synthetic_vca_5.hdr
          5 )
 
 
-ADD_TEST(hyVCA_AVIRIS
+ADD_TEST(hyTvVCA_AVIRIS
          ${TESTEXE_DIR}/otbHyperTests1
          otbVCATestHighSNR
          ${DATA}/AVIRIS/extract512/f970619t01p02_r02_sc01_goodbands.hdr
          ${TEMP}/aviris_vca_5.hdr
          13 )
 
+###################
+## Cublas
+###################
+
+SET( otbHyperCublasTests1_SRC
+  otbCublasStreamingStatisticsVectorImageFilter.cxx
+)
+
+ADD_EXECUTABLE(otbHyperCublasTests1 ${otbHyperCublasTests1_SRC} otbHyperCublasTests1.cxx)
+TARGET_LINK_LIBRARIES(otbHyperCublasTests1 OTBCommon OTBIO OTBBasicFilters OTBTesting OTBBasicFiltersGPUCuda)
 
+ADD_TEST(bfTuCublasStreamingStatisticsVectorImageFilterNew
+         ${TESTEXE_DIR}/otbHyperCublasTests1
+         otbCublasStreamingStatisticsVectorImageFilterNewTest)
+
+ADD_TEST(bfTvCublasStreamingStatisticsVectorImageFilter
+         ${TESTEXE_DIR}/otbHyperCublasTests1
+         otbCublasStreamingStatisticsVectorImageFilterTest
+         ${DATA}/AVIRIS/extract128/extract128_goodbands.hdr
+         ${TEMP}/extract128_stats_cublas.txt)
+
+ADD_TEST(bfTvCublasStreamingStatisticsVectorImageFilter_AVIRIS512
+         ${TESTEXE_DIR}/otbHyperCublasTests1
+         otbCublasStreamingStatisticsVectorImageFilterTest
+         ${DATA}/AVIRIS/extract512/f970619t01p02_r02_sc01_goodbands.hdr
+         ${TEMP}/aviris512_stats_cublas.txt)
 
 
 ###################
-##Anomaly detection
+## Anomaly detection
 ###################
 
 SET( otbAnomalyTests1_SRC

Testing/otbCublasStreamingStatisticsVectorImageFilter.cxx

+/*=========================================================================
+
+  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 "otbCublasStreamingStatisticsVectorImageFilter.h"
+#include "otbImageFileReader.h"
+#include "otbVectorImage.h"
+#include <fstream>
+#include "otbStreamingTraits.h"
+
+const unsigned int Dimension = 2;
+typedef float PixelType;
+
+typedef otb::VectorImage<PixelType, Dimension> ImageType;
+typedef otb::ImageFileReader<ImageType> ReaderType;
+typedef otb::CublasStreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType;
+
+int otbCublasStreamingStatisticsVectorImageFilterNewTest(int argc, char * argv[])
+{
+  StreamingStatisticsVectorImageFilterType::Pointer filter = StreamingStatisticsVectorImageFilterType::New();
+  std::cout << filter << std::endl;
+  return EXIT_SUCCESS;
+}
+
+int otbCublasStreamingStatisticsVectorImageFilterTest(int argc, char * argv[])
+{
+  const char * infname = argv[1];
+  const char * outfname = argv[2];
+
+  StreamingStatisticsVectorImageFilterType::Pointer filter = StreamingStatisticsVectorImageFilterType::New();
+
+  ReaderType::Pointer reader = ReaderType::New();
+  reader->SetFileName(infname);
+
+  //filter->GetStreamer()->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
+  filter->GetStreamer()->SetBufferMemorySize(100 * 1024 * 1024);
+  filter->SetInput(reader->GetOutput());
+  filter->Update();
+
+  std::ofstream file;
+  file.open(outfname);
+  file.precision(5);
+  //file << "Mean: " << filter->GetMean() << std::endl;
+  //file << "Covariance: " << filter->GetCovariance() << std::endl;
+  file << "Correlation: " << filter->GetCorrelation() << std::endl;
+  file.close();
+
+  return EXIT_SUCCESS;
+}

Testing/otbHyperCublasTests1.cxx

+/*=========================================================================
+
+  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
+#include <iostream>
+#include "otbTestMain.h"
+
+void RegisterTests()
+{
+  REGISTER_TEST(otbCublasStreamingStatisticsVectorImageFilterNewTest);
+  REGISTER_TEST(otbCublasStreamingStatisticsVectorImageFilterTest);
+//  REGISTER_TEST(otbCudaCorrelationNoStreamingTest);
+}