Orfeo Toolbox  4.2
otbSparseUnmixingImageFilter.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 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbSparseUnmixingImageFilter_txx
19 #define __otbSparseUnmixingImageFilter_txx
20 
22 
23 #include "otbMath.h"
24 #include "itkProcessObject.h"
25 
26 namespace otb {
27 
28 template < class TInputImage, class TOutputImage,
29  unsigned int VNbInputImage, class TPrecision,
30  Wavelet::Wavelet TMotherWaveletOperator >
33 {
34  this->SetNumberOfRequiredInputs( NumberOfInputImages );
35  this->SetNumberOfRequiredOutputs(1);
36  this->SetNthOutput(0, OutputImageListType::New());
37 
38  m_NumberOfComponentsRequired = 0;
39  m_NumberOfHistogramBins = 100;
40  m_AngleList = AngleListType::New();
41 
42  m_WvltFilterList = WvltFilterListType::New();
43  m_WvltFilterList->Resize( NumberOfInputImages );
44  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
45  {
46  m_WvltFilterList->SetNthElement(i, WvltFilterType::New());
47  m_WvltFilterList->GetNthElement(i)->SetNumberOfDecompositions(2);
48  m_WvltFilterList->GetNthElement(i)->SetSubsampleImageFactor(1);
49  }
50 
51  m_AngleListFilter = AngleListFilterType::New();
52  m_AngleListFilter->SetThresholdValue( 10. );
53 
54  m_Histogram = HistogramType::New();
55  m_Transformer = TransformFilterType::New();
56 }
57 
58 template < class TInputImage, class TOutputImage,
59  unsigned int VNbInputImage, class TPrecision,
60  Wavelet::Wavelet TMotherWaveletOperator >
61 void
63 ::SetInput ( unsigned int i, const InputImageType * img )
64 {
66  const_cast< InputImageType * >( img ) );
67 }
68 
69 template < class TInputImage, class TOutputImage,
70  unsigned int VNbInputImage, class TPrecision,
71  Wavelet::Wavelet TMotherWaveletOperator >
72 const TInputImage *
74 ::GetInput ( unsigned int i ) const
75 {
76  if ( i >= this->GetNumberOfInputs() )
77  {
78  return 0;
79  }
80 
81  return static_cast<const InputImageType * >
82  (this->itk::ProcessObject::GetInput(i) );
83 }
84 
85 template < class TInputImage, class TOutputImage,
86  unsigned int VNbInputImage, class TPrecision,
87  Wavelet::Wavelet TMotherWaveletOperator >
88 void
91 {
93  progress->SetMiniPipelineFilter(this);
94 
95  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
96  {
97  m_WvltFilterList->GetNthElement(i)->SetInput( this->GetInput(i) );
98  progress->RegisterInternalFilter( m_WvltFilterList->GetNthElement(i),
99  0.5f/static_cast<float>(NumberOfInputImages) );
100  m_WvltFilterList->GetNthElement(i)->Update();
101 
102  m_AngleListFilter->SetInput( i, m_WvltFilterList->GetNthElement(i)->GetOutput() );
103  }
104 
105  progress->RegisterInternalFilter(m_AngleListFilter, 0.25f);
106  m_AngleListFilter->Update();
107 
108  GenerateNumberOfComponentsRequired();
109  otbMsgDebugMacro( << m_NumberOfComponentsRequired << " sources found\n" );
110 
111  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
112  m_Transformer->SetInput( i, this->GetInput(i) );
113  m_Transformer->SetAngleList( m_AngleList );
114  progress->RegisterInternalFilter(m_Transformer, 0.25f);
115  m_Transformer->Update();
116 
117  this->GetOutput()->Resize( m_Transformer->GetOutput()->Size() );
118  for ( unsigned int i = 0; i < this->GetOutput()->Size(); ++i )
119  this->GetOutput()->SetNthElement( i, m_Transformer->GetOutput()->GetNthElement(i) );
120 }
121 
122 template < class TInputImage, class TOutputImage,
123  unsigned int VNbInputImage, class TPrecision,
124  Wavelet::Wavelet TMotherWaveletOperator >
125 void
128 {
134  HistogramSizeType size;
135  size.Fill( m_NumberOfHistogramBins );
136  MeasurementVectorType theMin (0.);
137  theMin[NumberOfInputImages-2] = -otb::CONST_PI;
139  m_Histogram->Initialize( size, theMin, theMax );
140 
141  typename InternalSampleListType::Iterator angleIter
142  = m_AngleListFilter->GetOutputSampleList()->Begin();
143 
144  if ( m_AngleListFilter->GetOutputSampleList()->Begin()
145  == m_AngleListFilter->GetOutputSampleList()->End() )
146  {
147  throw itk::ExceptionObject( __FILE__, __LINE__,
148  "The value of threshold is too high so that there is no angle value to consider for unmixing",
149  ITK_LOCATION);
150  }
151 
152  while ( angleIter != m_AngleListFilter->GetOutputSampleList()->End() )
153  {
154  if ( !m_Histogram->IncreaseFrequencyOfMeasurement( angleIter.GetMeasurementVector(), 1 ) )
155  std::cerr << "Data out of bounds\n";
156 
157  ++angleIter;
158  }
159 
160 #if 0
161  for ( unsigned int index = 0; index < m_Histogram->GetSize()[0]; ++index )
162  std::cerr << index << "\t" << m_Histogram->GetMeasurementVector( index )[0]
163  << "\t" << m_Histogram->GetFrequency(index) << "\n";
164 #endif
165 
166  AngleListPointerType angles = AngleListType::New();
167 
168  for ( unsigned int id = 0; id < m_Histogram->Size(); ++id )
169  {
170  HistogramIndexType curIdx = m_Histogram->GetIndex(id);
171 
172  // Is this curIdx a local max ?
173  bool isLocalMax = true;
174  for ( unsigned int k = 0; k < m_Histogram->GetSize().Size(); ++k )
175  {
176  HistogramIndexType prevIdx = curIdx;
177  if ( prevIdx[k] == 0 )
178  prevIdx[k] = m_Histogram->GetSize(k)-1;
179  else
180  prevIdx[k] = curIdx[k]-1;
181 
182 
183  HistogramIndexType nextIdx = curIdx;
184  if ( static_cast<PrecisionType>( nextIdx[k] ) == m_Histogram->GetSize(k)-1 )
185  nextIdx[k] = 0;
186  else
187  nextIdx[k] = curIdx[k]+1;
188 
189  double freq_prev = m_Histogram->GetFrequency( prevIdx );
190  double freq_cur = m_Histogram->GetFrequency( curIdx );
191  double freq_next = m_Histogram->GetFrequency( nextIdx );
192 
193  if ( !( freq_prev < freq_cur && freq_cur > freq_next ) )
194  {
195  isLocalMax = false;
196  break;
197  }
198  }
199 
200  if ( isLocalMax )
201  {
202  angles->PushBack( m_Histogram->GetMeasurementVector( curIdx ) );
203  }
204  }
205 
206  m_NumberOfComponentsRequired = angles->Size();
207  m_AngleList = angles;
208 
209 #if 1
210  std::cout << m_NumberOfComponentsRequired << " sources found\n";
211  for ( unsigned int i = 0; i < angles->Size(); i++ )
212  {
213  std::cerr << "Source " << i << ":";
214  for ( unsigned int j = 0; j < m_Histogram->GetSize().Size(); ++j )
215  std::cout << "\t" << angles->GetMeasurementVector(i)[j];
216  std::cout << "\n";
217  }
218 #endif
219 }
220 
221 } // end of namespace otb
222 
223 #endif

Generated at Sat Aug 16 2014 16:20:23 for Orfeo Toolbox with doxygen 1.8.3.1