Orfeo Toolbox  4.0
otbVcaImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12  Some parts of this code are derived from ITK. See ITKCopyright.txt
13  for details.
14 
15 
16  This software is distributed WITHOUT ANY WARRANTY; without even
17  the implied warranty of MERCHANT2ABILITY or FITNESS FOR A PARTICULAR
18  PURPOSE. See the above copyright notices for more information.
19 
20 =========================================================================*/
21 #ifndef __otbVCAImageFilter_txx
22 #define __otbVCAImageFilter_txx
23 
24 #include "otbVcaImageFilter.h"
26 
27 namespace otb {
28 
29 template<class TImage>
31  : m_NumberOfEndmembers(0)
32 {
33 }
34 
35 template<class TImage>
37 {
38 }
39 
40 template<class TImage>
42 {
43  // call the superclass' implementation of this method
44  Superclass::GenerateOutputInformation();
45 
46  // get pointers to the input and output
47  typename VectorImageType::ConstPointer inputPtr = this->GetInput();
48  typename VectorImageType::Pointer outputPtr = this->GetOutput();
49 
50  if ( !inputPtr || !outputPtr )
51  {
52  return;
53  }
54 
55  const unsigned int numberOfBands = inputPtr->GetNumberOfComponentsPerPixel();
56  RegionType outputRegion;
57  IndexType outputOrigin;
58  SizeType outputSize;
59 
60  outputOrigin.Fill(0);
61  outputSize[0] = m_NumberOfEndmembers;
62  outputSize[1] = 1;
63  outputRegion.SetIndex(outputOrigin);
64  outputRegion.SetSize(outputSize);
65 
66  outputPtr->SetLargestPossibleRegion( outputRegion );
67  outputPtr->SetNumberOfComponentsPerPixel( numberOfBands );
68 }
69 
70 template<class TImage>
72 {
73  typedef typename ForwardPCAImageFilterType::NormalizeFilterType NormalizeFilterType;
74 
75  VectorImageType* input = const_cast<VectorImageType*>(this->GetInput());
76  const unsigned int nbBands = this->GetInput()->GetNumberOfComponentsPerPixel();
77 
78 
79  otbMsgDevMacro( "Computing image stats" )
81  StreamingStatisticsVectorImageFilterType::New();
82 
83  statsInput->SetInput(input);
84  statsInput->SetEnableMinMax(false);
85  statsInput->Update();
86 
87  double SNR, SNRth;
88  vnl_matrix<PrecisionType> Ud;
89 
90  // SNR computation
91  {
92  vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
93  vnl_svd<PrecisionType> svd(R);
94  vnl_matrix<PrecisionType> U = svd.U();
95  vnl_matrix<PrecisionType> Ud = U.get_n_columns(0, m_NumberOfEndmembers);
96  vnl_matrix<PrecisionType> Udt = Ud.transpose();
97 
98  // To remove the mean
99  typename NormalizeFilterType::Pointer normalize = NormalizeFilterType::New();
100  normalize->SetInput(input);
101  normalize->SetMean(statsInput->GetMean());
102  normalize->SetUseMean(true);
103  normalize->SetUseStdDev(false);
104 
105  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
106  mulUd->MatrixByVectorOn();
107  mulUd->SetInput(normalize->GetOutput());
108  mulUd->SetMatrix(Udt);
109  mulUd->UpdateOutputInformation();
110 
111  typename StreamingStatisticsVectorImageFilterType::Pointer transformedStats =
112  StreamingStatisticsVectorImageFilterType::New();
113 
114  transformedStats->SetInput(mulUd->GetOutput());
115  transformedStats->SetEnableMinMax(false);
116  transformedStats->Update();
117 
118  double P_R = nbBands * statsInput->GetComponentCorrelation();
119  double P_Rp = m_NumberOfEndmembers * transformedStats->GetComponentCorrelation()
120  + statsInput->GetMean().GetSquaredNorm();
121  //const double qL = static_cast<double>(m_NumberOfEndmembers) / nbBands;
122  SNR = vcl_abs(10*vcl_log10( (P_Rp - (m_NumberOfEndmembers/nbBands)*P_R) / (P_R - P_Rp) ));
123  }
124 
125  SNRth = 15.0 + 10.0 * vcl_log( static_cast<double>(m_NumberOfEndmembers) ) + 8.0;
126 
127  typename VectorImageType::Pointer Xd;
128  typename VectorImageType::Pointer Y;
129 
130  std::vector<itk::ProcessObject::Pointer> refHolder;
131  typename ForwardPCAImageFilterType::Pointer pca = ForwardPCAImageFilterType::New();
132  typename NormalizeFilterType::Pointer normalize = NormalizeFilterType::New();
133 
134  if (SNR > SNRth)
135  {
136  otbMsgDevMacro( "Using projective projection for dimensionnality reduction" )
137 
138  otbMsgDevMacro( "Computing SVD of correlation matrix" )
139  // Take the correlation matrix
140  vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
141 
142  // Apply SVD
143  vnl_svd<PrecisionType> svd(R);
144  vnl_matrix<PrecisionType> U = svd.U();
145  Ud = U.get_n_columns(0, m_NumberOfEndmembers);
146  vnl_matrix<PrecisionType> UdT = Ud.transpose();
147 
148  otbMsgDevMacro( "Apply dimensionality reduction" )
149  // Xd = Ud.'*M;
150  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
151  mulUd->MatrixByVectorOn();
152  mulUd->SetInput(this->GetInput());
153  mulUd->SetMatrix(UdT);
154  mulUd->UpdateOutputInformation();
155  refHolder.push_back(mulUd.GetPointer());
156 
157  Xd = mulUd->GetOutput();
158 
159  // Compute mean(Xd)
160  otbMsgDevMacro( "Compute mean(Xd)" )
162  statsXd = StreamingStatisticsVectorImageFilterType::New();
163  statsXd->SetEnableSecondOrderStats(false);
164  statsXd->SetInput(Xd);
165  // statsXd->GetStreamer()->SetNumberOfDivisionsStrippedStreaming(10);
166  statsXd->Update();
167  typename VectorImageType::PixelType Xdmean = statsXd->GetMean();
168  otbMsgDevMacro( "mean(Xd) = " << Xdmean)
169 
170  // Projective projection
171  // Xd ./ repmat( sum( Xd .* repmat(u, [1 N]) ) , [d 1]);
172  otbMsgDevMacro( "Compute projective projection" )
173  typename ProjectiveProjectionImageFilterType::Pointer proj = ProjectiveProjectionImageFilterType::New();
174  proj->SetInput(Xd);
175  proj->SetProjectionDirection(Xdmean);
176  refHolder.push_back(proj.GetPointer());
177 
178  Xd = proj->GetOutput();
179  Y = Xd;
180 
181  }
182  else
183  {
184  otbMsgDevMacro( "Using PCA for dimensionnality reduction" )
185 
186  normalize->SetInput(input);
187  normalize->SetMean(statsInput->GetMean());
188  normalize->SetUseMean(true);
189  normalize->SetUseStdDev(false);
190 
191  // Take the correlation matrix
192  vnl_matrix<PrecisionType> R = statsInput->GetCovariance().GetVnlMatrix();
193 
194  // Apply SVD
195  vnl_svd<PrecisionType> svd(R);
196  vnl_matrix<PrecisionType> U = svd.U();
197  Ud = U.get_n_columns(0, m_NumberOfEndmembers - 1);
198  vnl_matrix<PrecisionType> UdT = Ud.transpose();
199 
200  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
201  mulUd->MatrixByVectorOn();
202  mulUd->SetInput(normalize->GetOutput());
203  mulUd->SetMatrix(UdT);
204  mulUd->UpdateOutputInformation();
205  refHolder.push_back(mulUd.GetPointer());
206 
207  Xd = mulUd->GetOutput();
208 
209  typename VectorImageToAmplitudeImageFilterType::Pointer normComputer = VectorImageToAmplitudeImageFilterType::New();
210  normComputer->SetInput(Xd);
211 
212  typename StreamingMinMaxImageFilterType::Pointer maxNormComputer = StreamingMinMaxImageFilterType::New();
213  maxNormComputer->SetInput(normComputer->GetOutput());
214  maxNormComputer->Update();
215  typename ImageType::PixelType maxNorm = maxNormComputer->GetMaximum();
216  otbMsgDevMacro( "maxNorm : " << maxNorm)
217 
218  typename ConcatenateScalarValueImageFilterType::Pointer concat = ConcatenateScalarValueImageFilterType::New();
219  concat->SetInput(Xd);
220  concat->SetScalarValue(maxNorm);
221  refHolder.push_back(concat.GetPointer());
222 
223  Y = concat->GetOutput();
224  Y->UpdateOutputInformation();
225  }
226 
227  // E : result, will contain the endmembers
228  vnl_matrix<PrecisionType> E(nbBands, m_NumberOfEndmembers);
229 
230  // A = zeros(q, q)
231  // A(q, 1) = 1
232  vnl_matrix<PrecisionType> A(m_NumberOfEndmembers, m_NumberOfEndmembers);
233  A.fill(0);
234  A(m_NumberOfEndmembers - 1, 0) = 1;
235  typename RandomVariateGeneratorType::Pointer randomGen = RandomVariateGeneratorType::New();
236 
237  for (unsigned int i = 0; i < m_NumberOfEndmembers; ++i)
238  {
239  otbMsgDevMacro( "----------------------------------------" )
240  otbMsgDevMacro( "Iteration " << i )
241 
242  // w = rand(q, 1)
243  otbMsgDevMacro( "Random vector generation " )
244  vnl_vector<PrecisionType> w(m_NumberOfEndmembers);
245  for (unsigned int j = 0; j < w.size(); ++j)
246  {
247  w(j) = randomGen->GetVariateWithOpenRange();
248  }
249 
250  // f = ((I - A*pinv(A))*w) / (norm(I - A*pinv(A))*w))
251  otbMsgDevMacro( "f = ((I - A*pinv(A))*w) /(norm(I - A*pinv(A))*w))" )
252  vnl_matrix<PrecisionType> tmpMat(m_NumberOfEndmembers, m_NumberOfEndmembers);
253  tmpMat.set_identity();
254  otbMsgDevMacro( "A" << std::endl << A )
255  vnl_svd<PrecisionType> Asvd(A);
256  tmpMat -= A * Asvd.inverse();
257 
258  vnl_vector<PrecisionType> tmpNumerator = tmpMat * w;
259  vnl_vector<PrecisionType> f = tmpNumerator / tmpNumerator.two_norm();
260 
261  // v = f.'*Y
262  otbMsgDevMacro( "f = " << f );
263  otbMsgDevMacro( "v = f.'*Y" );
264  typename DotProductImageFilterType::Pointer dotfY = DotProductImageFilterType::New();
265  dotfY->SetInput(Y);
266 
267  typename VectorImageType::PixelType fV(f.data_block(), f.size());
268  dotfY->SetVector(typename VectorImageType::PixelType(fV));
269  typename ImageType::Pointer v = dotfY->GetOutput();
270 
271  // abs(v)
272  otbMsgDevMacro( "abs(v)" );
273  typename AbsImageFilterType::Pointer absVFilter = AbsImageFilterType::New();
274  absVFilter->SetInput(v);
275 
276  // max(abs(v))
277  otbMsgDevMacro( "max(abs(v))" );
278  typename StreamingMinMaxImageFilterType::Pointer maxAbs = StreamingMinMaxImageFilterType::New();
279  maxAbs->SetInput(absVFilter->GetOutput());
280  maxAbs->Update();
281 
282  // k = arg_max( max(abs(v)) )
283  otbMsgDevMacro( "k = arg_max( max(abs(v)) )" )
284  IndexType maxIdx = maxAbs->GetMaximumIndex();
285  otbMsgDevMacro( "max : " << maxAbs->GetMaximum() )
286  otbMsgDevMacro( "maxIdx : " << maxIdx )
287 
288  // extract Y(:, k)
289  otbMsgDevMacro( "Y(:, k)" )
290  RegionType region;
291  region.SetIndex( maxIdx );
292  SizeType size;
293  size.Fill(1);
294  region.SetSize(size);
295  Y->SetRequestedRegion(region);
296  Y->Update();
297 
298  // store new endmember in A
299  // A(:, i) = Y(:, k)
300  otbMsgDevMacro( "A(:, i) = Y(:, k)" )
301  typename VectorImageType::PixelType e = Y->GetPixel(maxIdx);
302  otbMsgDevMacro( "e = " << e )
303  A.set_column(i, e.GetDataPointer());
304 
305  otbMsgDevMacro( "A" << std::endl << A )
306 
307  // reproject new endmember in original space
308  vnl_vector<PrecisionType> u;
309  if (SNR > SNRth)
310  {
311  // u = Ud * Xd(:, k)
312  otbMsgDevMacro( "u = Ud * Xd(:, k)" )
313  Xd->SetRequestedRegion(region);
314  Xd->Update();
315  typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
316  vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
317  u = Ud * xdV;
318  }
319  else
320  {
321  // u = invPCA( Xd(:, k) )
322  otbMsgDevMacro( "u = Ud * Xd(:, k)" )
323  Xd->SetRequestedRegion(region);
324  Xd->Update();
325  typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
326  vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
327  //Ud = pca->GetTransformationMatrix().GetVnlMatrix().transpose();
328 
329  vnl_vector<PrecisionType> mean(statsInput->GetMean().GetDataPointer(), statsInput->GetMean().GetSize());
330  u = Ud * xdV + mean;
331  }
332 
333  // E(:, i) = u
334  otbMsgDevMacro( "E(:, i) = u" )
335  otbMsgDevMacro( "u = " << u )
336  E.set_column(i, u);
337  }
338 
339  typename VectorImageType::Pointer output = this->GetOutput();
340  output->SetRegions(output->GetLargestPossibleRegion());
341  //output->SetNumberOfComponentsPerPixel(input->GetNumberOfComponentsPerPixel());
342  output->Allocate();
343 
344  itk::ImageRegionIteratorWithIndex<VectorImageType> it(output, output->GetLargestPossibleRegion());
345  unsigned int i;
346  for (it.GoToBegin(), i = 0; !it.IsAtEnd(); ++it, ++i)
347  {
348  vnl_vector<PrecisionType> e = E.get_column(i);
349  typename VectorImageType::PixelType pixel(input->GetNumberOfComponentsPerPixel());
350  for (unsigned int j = 0; j < e.size(); ++j)
351  {
352  pixel[j] = E(j, i);
353  }
354  it.Set(pixel);
355  }
356 }
357 
361 template<class TImage>
362 void VCAImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
363 {
364  Superclass::PrintSelf( os, indent );
365 }
366 
367 } //end namespace otb
368 
369 #endif

Generated at Sat Mar 8 2014 16:24:25 for Orfeo Toolbox with doxygen 1.8.3.1