diff --git a/Code/Hyperspectral/otbVcaImageFilter.h b/Code/Hyperspectral/otbVcaImageFilter.h
index 2622906cf58d711823039947b190273eeafe27b7..ed4f8e9064cb70db696d6eac8b4fd6adbab31926 100644
--- a/Code/Hyperspectral/otbVcaImageFilter.h
+++ b/Code/Hyperspectral/otbVcaImageFilter.h
@@ -27,9 +27,12 @@
#include "otbProjectiveProjectionImageFilter.h"
#include "otbMatrixMultiplyImageFilter.h"
#include "otbVectorImageToMatrixImageFilter.h"
-#include "otbStreamingMinMaxImageFilter2.h"
+#include "otbStreamingMinMaxImageFilter.h"
#include "otbStreamingStatisticsImageFilter.h"
-#include "otbStreamingStatisticsVectorImageFilter2.h"
+#include "otbStreamingStatisticsVectorImageFilter.h"
+#include "otbPCAImageFilter.h"
+#include "otbVectorImageToAmplitudeImageFilter.h"
+#include "otbConcatenateScalarValueImageFilter.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "vnl/algo/vnl_svd.h"
@@ -68,17 +71,21 @@ public:
typedef typename VectorImageType::InternalPixelType InternalPixelType;
typedef InternalPixelType PrecisionType;
- typedef otb::Image<InternalPixelType> ImageType;
+ typedef otb::Image<InternalPixelType> ImageType;
typedef itk::AbsImageFilter<ImageType, ImageType> AbsImageFilterType;
typedef otb::ProjectiveProjectionImageFilter<VectorImageType,VectorImageType,PrecisionType> ProjectiveProjectionImageFilterType;
typedef otb::DotProductImageFilter<VectorImageType,ImageType> DotProductImageFilterType;
typedef otb::MatrixMultiplyImageFilter<VectorImageType,VectorImageType,PrecisionType> MatrixMultiplyImageFilterType;
typedef otb::VectorImageToMatrixImageFilter<VectorImageType> VectorImageToMatrixImageFilterType;
- typedef otb::StreamingMinMaxImageFilter2<ImageType> StreamingMinMaxImageFilterType;
- typedef otb::StreamingStatisticsVectorImageFilter2<VectorImageType,PrecisionType> StreamingStatisticsVectorImageFilterType;
+ typedef otb::StreamingMinMaxImageFilter<ImageType> StreamingMinMaxImageFilterType;
+ typedef otb::StreamingStatisticsVectorImageFilter<VectorImageType,PrecisionType> StreamingStatisticsVectorImageFilterType;
typedef otb::StreamingStatisticsImageFilter<ImageType> StreamingStatisticsImageFilterType;
typedef itk::Statistics::MersenneTwisterRandomVariateGenerator RandomVariateGeneratorType;
+ typedef otb::PCAImageFilter< VectorImageType, VectorImageType, otb::Transform::FORWARD > ForwardPCAImageFilterType;
+ typedef otb::PCAImageFilter< VectorImageType, VectorImageType, otb::Transform::INVERSE > InversePCAImageFilterType;
+ typedef otb::VectorImageToAmplitudeImageFilter< VectorImageType, ImageType > VectorImageToAmplitudeImageFilterType;
+ typedef otb::ConcatenateScalarValueImageFilter< VectorImageType, VectorImageType > ConcatenateScalarValueImageFilterType;
// creation of SmartPointer
itkNewMacro(Self);
diff --git a/Code/Hyperspectral/otbVcaImageFilter.txx b/Code/Hyperspectral/otbVcaImageFilter.txx
index 1e99e1755062d32604f53ca2902f2ba7c88c532c..f639c21474bbebea8076442bb841a8287b47fc90 100644
--- a/Code/Hyperspectral/otbVcaImageFilter.txx
+++ b/Code/Hyperspectral/otbVcaImageFilter.txx
@@ -21,11 +21,14 @@
#ifndef __otbVCAImageFilter_txx
#define __otbVCAImageFilter_txx
-#define otbMsgDevMacroVCA(p) std::cout << p << std::endl;
#include "otbVcaImageFilter.h"
#include "otbStandardWriterWatcher.h"
+
+#define otbMsgDevMacroVCA(p) std::cout << p << std::endl;
+//#define otbMsgDevMacroVCA(p) ;
+
namespace otb {
template<class TImage>
@@ -75,7 +78,8 @@ void VCAImageFilter<TImage>::GenerateData()
VectorImageType* input = const_cast<VectorImageType*>(this->GetInput());
const unsigned int nbBands = this->GetInput()->GetNumberOfComponentsPerPixel();
- otbMsgDevMacroVCA( "Computing image stats" );
+
+ otbMsgDevMacroVCA( "Computing image stats" )
typename StreamingStatisticsVectorImageFilterType::Pointer statsInput =
StreamingStatisticsVectorImageFilterType::New();
@@ -84,7 +88,8 @@ void VCAImageFilter<TImage>::GenerateData()
statsInput->SetEnableMinMax(false);
statsInput->Update();
- otbMsgDevMacroVCA( "Computing SVD of correlation matrix" );
+/*
+// otbMsgDevMacroVCA( "Computing SVD of correlation matrix" );
// Take the correlation matrix
vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
@@ -93,16 +98,131 @@ void VCAImageFilter<TImage>::GenerateData()
vnl_matrix<PrecisionType> U = svd.U();
vnl_matrix<PrecisionType> Ud = U.get_n_columns(0, m_NumberOfEndmembers);
vnl_matrix<PrecisionType> UdT = Ud.transpose();
+*/
+ vnl_matrix<PrecisionType> Ud;
+
+ // TODO : temporary hack. implement SNR estimation formula here
+ double SNR = 40;
+ double SNRth = 30;
+ typename VectorImageType::Pointer Xd;
+ typename VectorImageType::Pointer Y;
+
+ std::vector<itk::ProcessObject::Pointer> refHolder;
+ typename ForwardPCAImageFilterType::Pointer pca = ForwardPCAImageFilterType::New();
+ typedef typename ForwardPCAImageFilterType::NormalizeFilterType NormalizeFilterType;
+ typename NormalizeFilterType::Pointer normalize = NormalizeFilterType::New();
+
+ if (SNR > SNRth)
+ {
+ otbMsgDevMacroVCA( "Using projective projection for dimensionnality reduction" )
+
+ otbMsgDevMacroVCA( "Computing SVD of correlation matrix" )
+ // Take the correlation matrix
+ vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
+
+ // Apply SVD
+ vnl_svd<PrecisionType> svd(R);
+ vnl_matrix<PrecisionType> U = svd.U();
+ Ud = U.get_n_columns(0, m_NumberOfEndmembers);
+ vnl_matrix<PrecisionType> UdT = Ud.transpose();
+
+ otbMsgDevMacroVCA( "Apply dimensionality reduction" )
+ // Xd = Ud.'*M;
+ typename MatrixMultiplyImageFilterType::Pointer mulUd = MatrixMultiplyImageFilterType::New();
+ mulUd->SetInput(this->GetInput());
+ mulUd->SetMatrix(UdT);
+ mulUd->UpdateOutputInformation();
+ refHolder.push_back(mulUd.GetPointer());
+
+ Xd = mulUd->GetOutput();
+
+ // Compute mean(Xd)
+ otbMsgDevMacroVCA( "Compute mean(Xd)" )
+ typename StreamingStatisticsVectorImageFilterType::Pointer
+ statsXd = StreamingStatisticsVectorImageFilterType::New();
+ statsXd->SetEnableSecondOrderStats(false);
+ statsXd->SetInput(Xd);
+ // statsXd->GetStreamer()->SetNumberOfDivisionsStrippedStreaming(10);
+ statsXd->Update();
+ typename VectorImageType::PixelType Xdmean = statsXd->GetMean();
+ otbMsgDevMacroVCA( "mean(Xd) = " << Xdmean)
+
+ // Projective projection
+ // Xd ./ repmat( sum( Xd .* repmat(u,[1 N]) ) ,[d 1]);
+ otbMsgDevMacroVCA( "Compute projective projection" )
+ typename ProjectiveProjectionImageFilterType::Pointer proj = ProjectiveProjectionImageFilterType::New();
+ proj->SetInput(Xd);
+ proj->SetProjectionDirection(Xdmean);
+ refHolder.push_back(proj.GetPointer());
+
+ Xd = proj->GetOutput();
+ Y = Xd;
+
+ }
+ else
+ {
+ otbMsgDevMacroVCA( "Using PCA for dimensionnality reduction" )
+
+ normalize->SetInput(input);
+ normalize->SetMean(statsInput->GetMean());
+ normalize->SetUseMean(true);
+ normalize->SetUseStdDev(false);
+
+ // Take the correlation matrix
+ vnl_matrix<PrecisionType> R = statsInput->GetCovariance().GetVnlMatrix();
+
+ // Apply SVD
+ vnl_svd<PrecisionType> svd(R);
+ vnl_matrix<PrecisionType> U = svd.U();
+ Ud = U.get_n_columns(0, m_NumberOfEndmembers - 1);
+ vnl_matrix<PrecisionType> UdT = Ud.transpose();
+
+ typename MatrixMultiplyImageFilterType::Pointer mulUd = MatrixMultiplyImageFilterType::New();
+ mulUd->SetInput(normalize->GetOutput());
+ mulUd->SetMatrix(UdT);
+ mulUd->UpdateOutputInformation();
+ refHolder.push_back(mulUd.GetPointer());
+
+ Xd = mulUd->GetOutput();
+
+ typename VectorImageToAmplitudeImageFilterType::Pointer normComputer = VectorImageToAmplitudeImageFilterType::New();
+ normComputer->SetInput(Xd);
+
+ typename StreamingMinMaxImageFilterType::Pointer maxNormComputer = StreamingMinMaxImageFilterType::New();
+ maxNormComputer->SetInput(normComputer->GetOutput());
+ maxNormComputer->Update();
+ typename ImageType::PixelType maxNorm = maxNormComputer->GetMaximum();
+
+ typename ConcatenateScalarValueImageFilterType::Pointer concat = ConcatenateScalarValueImageFilterType::New();
+ concat->SetInput(Xd);
+ concat->SetScalarValue(maxNorm);
+ refHolder.push_back(concat.GetPointer());
+
+ Y = concat->GetOutput();
+ Y->UpdateOutputInformation();
+
+ /*
+ pca->SetInput(normalize->GetOutput());
+ pca->SetNumberOfPrincipalComponentsRequired( m_NumberOfEndmembers - 1 );
+ pca->SetCovarianceMatrix(statsInput->GetCovariance());
+ // refHolder.push_back(pca.GetPointer());
+ // pca->GetTransformationMatrix() is not up to date at this point...
+ // Ud = pca->GetTransformationMatrix().GetVnlMatrix().transpose();
+ Xd = pca->GetOutput();
+ */
+
+ }
+/*
otbMsgDevMacroVCA( "Apply dimensionality reduction" );
// Xd = Ud.'*M;
typename MatrixMultiplyImageFilterType::Pointer mulUd = MatrixMultiplyImageFilterType::New();
mulUd->SetInput(this->GetInput());
mulUd->SetMatrix(UdT);
mulUd->UpdateOutputInformation();
-
typename VectorImageType::Pointer Xd = mulUd->GetOutput();
+
// Compute mean(Xd)
otbMsgDevMacroVCA( "Compute mean(Xd)" );
typename StreamingStatisticsVectorImageFilterType::Pointer statsXd = \
@@ -110,7 +230,7 @@ void VCAImageFilter<TImage>::GenerateData()
statsXd->SetEnableMinMax(false);
statsXd->SetEnableSecondOrderStats(false);
statsXd->SetInput(Xd);
-// statsXd->GetStreamer()->SetBufferNumberOfLinesDivisions(10);
+// statsXd->GetStreamer()->SetNumberOfDivisionsStrippedStreaming(10);
statsXd->Update();
typename VectorImageType::PixelType Xdmean = statsXd->GetMean();
@@ -121,7 +241,7 @@ void VCAImageFilter<TImage>::GenerateData()
proj->SetInput(Xd);
proj->SetProjectionDirection(Xdmean);
typename VectorImageType::Pointer Y = proj->GetOutput();
-
+*/
// E : result, will contain the endmembers
vnl_matrix<PrecisionType> E(nbBands, m_NumberOfEndmembers);
@@ -134,21 +254,22 @@ void VCAImageFilter<TImage>::GenerateData()
for (unsigned int i = 0; i < m_NumberOfEndmembers; ++i)
{
- otbMsgDevMacroVCA( "Iteration " << i );
+ otbMsgDevMacroVCA( "----------------------------------------" )
+ otbMsgDevMacroVCA( "Iteration " << i )
// w = rand(q,1)
- otbMsgDevMacroVCA( "Random vector generation " );
+ otbMsgDevMacroVCA( "Random vector generation " )
vnl_vector<PrecisionType> w(m_NumberOfEndmembers);
- for (unsigned int j = 0; j < m_NumberOfEndmembers; ++j)
+ for (unsigned int j = 0; j < w.size(); ++j)
{
w(j) = randomGen->GetVariateWithOpenRange();
}
// f = ((I - A*pinv(A))*w) / (norm(I - A*pinv(A))*w))
- otbMsgDevMacroVCA( "f = ((I - A*pinv(A))*w) /(norm(I - A*pinv(A))*w))" );
- vnl_matrix<PrecisionType> tmpMat(m_NumberOfEndmembers,m_NumberOfEndmembers);
+ otbMsgDevMacroVCA( "f = ((I - A*pinv(A))*w) /(norm(I - A*pinv(A))*w))" )
+ vnl_matrix<PrecisionType> tmpMat(m_NumberOfEndmembers, m_NumberOfEndmembers);
tmpMat.set_identity();
- otbMsgDevMacroVCA( "A" << std::endl << A );
+ otbMsgDevMacroVCA( "A" << std::endl << A )
vnl_svd<PrecisionType> Asvd(A);
tmpMat -= A * Asvd.inverse();
@@ -156,9 +277,11 @@ void VCAImageFilter<TImage>::GenerateData()
vnl_vector<PrecisionType> f = tmpNumerator / tmpNumerator.two_norm();
// v = f.'*Y
+ otbMsgDevMacroVCA( "f = " << f );
otbMsgDevMacroVCA( "v = f.'*Y" );
typename DotProductImageFilterType::Pointer dotfY = DotProductImageFilterType::New();
dotfY->SetInput(Y);
+
typename VectorImageType::PixelType fV(f.data_block(), f.size());
dotfY->SetVector(typename VectorImageType::PixelType(fV));
typename ImageType::Pointer v = dotfY->GetOutput();
@@ -172,17 +295,17 @@ void VCAImageFilter<TImage>::GenerateData()
otbMsgDevMacroVCA( "max(abs(v))" );
typename StreamingMinMaxImageFilterType::Pointer maxAbs = StreamingMinMaxImageFilterType::New();
maxAbs->SetInput(absVFilter->GetOutput());
-// maxAbs->GetStreamer()->SetBufferNumberOfLinesDivisions(10);
+// maxAbs->GetStreamer()->SetNumberOfDivisionsStrippedStreaming(10);
maxAbs->Update();
// k = arg_max( max(abs(v)) )
- otbMsgDevMacroVCA( "k = arg_max( max(abs(v)) )" );
+ otbMsgDevMacroVCA( "k = arg_max( max(abs(v)) )" )
IndexType maxIdx = maxAbs->GetMaximumIndex();
- otbMsgDevMacroVCA( "max : " << maxAbs->GetMaximum() );
- otbMsgDevMacroVCA( "maxIdx : " << maxIdx );
+ otbMsgDevMacroVCA( "max : " << maxAbs->GetMaximum() )
+ otbMsgDevMacroVCA( "maxIdx : " << maxIdx )
// extract Y(:,k)
- otbMsgDevMacroVCA( "Y(:,k)" );
+ otbMsgDevMacroVCA( "Y(:,k)" )
RegionType region;
region.SetIndex( maxIdx );
SizeType size;
@@ -193,33 +316,49 @@ void VCAImageFilter<TImage>::GenerateData()
// store new endmember in A
// A(:,i) = Y(:,k)
- otbMsgDevMacroVCA( "A(:,i) = Y(:,k)" );
+ otbMsgDevMacroVCA( "A(:,i) = Y(:,k)" )
typename VectorImageType::PixelType e = Y->GetPixel(maxIdx);
otbMsgDevMacroVCA( "e = " << e )
A.set_column(i, e.GetDataPointer());
- otbMsgDevMacroVCA( "A" << std::endl << A );
- // reproject in original space
- // u = Ud * Xd(:,k)
- otbMsgDevMacroVCA( "u = Ud * Xd(:,k)" );
- Xd->SetRequestedRegion(region);
- Xd->Update();
- typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
- vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
- vnl_vector<PrecisionType> u = Ud * xdV;
+ otbMsgDevMacroVCA( "A" << std::endl << A )
+
+ // reproject new endmember in original space
+ vnl_vector<PrecisionType> u;
+ if (SNR > SNRth)
+ {
+ // u = Ud * Xd(:,k)
+ otbMsgDevMacroVCA( "u = Ud * Xd(:,k)" )
+ Xd->SetRequestedRegion(region);
+ Xd->Update();
+ typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
+ vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
+ u = Ud * xdV;
+ }
+ else
+ {
+ // u = invPCA( Xd(:,k) )
+ otbMsgDevMacroVCA( "u = Ud * Xd(:,k)" )
+ Xd->SetRequestedRegion(region);
+ Xd->Update();
+ typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
+ vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
+ //Ud = pca->GetTransformationMatrix().GetVnlMatrix().transpose();
+
+ vnl_vector<PrecisionType> mean(statsInput->GetMean().GetDataPointer(), statsInput->GetMean().GetSize());
+ u = Ud * xdV + mean;
+ }
// E(:, i) = u
- otbMsgDevMacroVCA( "E(:, i) = u" );
+ otbMsgDevMacroVCA( "E(:, i) = u" )
otbMsgDevMacroVCA( "u = " << u )
E.set_column(i, u);
}
typename VectorImageType::Pointer output = this->GetOutput();
- otbMsgDevMacroVCA( "output = " << output )
output->SetRegions(output->GetLargestPossibleRegion());
//output->SetNumberOfComponentsPerPixel(input->GetNumberOfComponentsPerPixel());
output->Allocate();
- otbMsgDevMacroVCA( "output = " << output )
itk::ImageRegionIteratorWithIndex<VectorImageType> it(output, output->GetLargestPossibleRegion());
unsigned int i;