OTB  5.11.0
Orfeo Toolbox
otbSparseUnmixingImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbSparseUnmixingImageFilter_txx
22 #define otbSparseUnmixingImageFilter_txx
23 
25 
26 #include "otbMath.h"
27 #include "itkProcessObject.h"
28 
29 namespace otb {
30 
31 template < class TInputImage, class TOutputImage,
32  unsigned int VNbInputImage, class TPrecision,
33  Wavelet::Wavelet TMotherWaveletOperator >
36 {
37  this->SetNumberOfRequiredInputs( NumberOfInputImages );
38  this->SetNumberOfRequiredOutputs(1);
39  this->SetNthOutput(0, OutputImageListType::New());
40 
41  m_NumberOfComponentsRequired = 0;
42  m_NumberOfHistogramBins = 100;
43  m_AngleList = AngleListType::New();
44 
45  m_WvltFilterList = WvltFilterListType::New();
46  m_WvltFilterList->Resize( NumberOfInputImages );
47  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
48  {
49  m_WvltFilterList->SetNthElement(i, WvltFilterType::New());
50  m_WvltFilterList->GetNthElement(i)->SetNumberOfDecompositions(2);
51  m_WvltFilterList->GetNthElement(i)->SetSubsampleImageFactor(1);
52  }
53 
54  m_AngleListFilter = AngleListFilterType::New();
55  m_AngleListFilter->SetThresholdValue( 10. );
56 
57  m_Histogram = HistogramType::New();
58  m_Transformer = TransformFilterType::New();
59 }
60 
61 template < class TInputImage, class TOutputImage,
62  unsigned int VNbInputImage, class TPrecision,
63  Wavelet::Wavelet TMotherWaveletOperator >
64 void
66 ::SetInput ( unsigned int i, const InputImageType * img )
67 {
69  const_cast< InputImageType * >( img ) );
70 }
71 
72 template < class TInputImage, class TOutputImage,
73  unsigned int VNbInputImage, class TPrecision,
74  Wavelet::Wavelet TMotherWaveletOperator >
75 const TInputImage *
77 ::GetInput ( unsigned int i ) const
78 {
79  if ( i >= this->GetNumberOfInputs() )
80  {
81  return ITK_NULLPTR;
82  }
83 
84  return static_cast<const InputImageType * >
85  (this->itk::ProcessObject::GetInput(i) );
86 }
87 
88 template < class TInputImage, class TOutputImage,
89  unsigned int VNbInputImage, class TPrecision,
90  Wavelet::Wavelet TMotherWaveletOperator >
91 void
94 {
96  progress->SetMiniPipelineFilter(this);
97 
98  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
99  {
100  m_WvltFilterList->GetNthElement(i)->SetInput( this->GetInput(i) );
101  progress->RegisterInternalFilter( m_WvltFilterList->GetNthElement(i),
102  0.5f/static_cast<float>(NumberOfInputImages) );
103  m_WvltFilterList->GetNthElement(i)->Update();
104 
105  m_AngleListFilter->SetInput( i, m_WvltFilterList->GetNthElement(i)->GetOutput() );
106  }
107 
108  progress->RegisterInternalFilter(m_AngleListFilter, 0.25f);
109  m_AngleListFilter->Update();
110 
111  GenerateNumberOfComponentsRequired();
112  otbMsgDebugMacro( << m_NumberOfComponentsRequired << " sources found\n" );
113 
114  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
115  m_Transformer->SetInput( i, this->GetInput(i) );
116  m_Transformer->SetAngleList( m_AngleList );
117  progress->RegisterInternalFilter(m_Transformer, 0.25f);
118  m_Transformer->Update();
119 
120  this->GetOutput()->Resize( m_Transformer->GetOutput()->Size() );
121  for ( unsigned int i = 0; i < this->GetOutput()->Size(); ++i )
122  this->GetOutput()->SetNthElement( i, m_Transformer->GetOutput()->GetNthElement(i) );
123 }
124 
125 template < class TInputImage, class TOutputImage,
126  unsigned int VNbInputImage, class TPrecision,
127  Wavelet::Wavelet TMotherWaveletOperator >
128 void
131 {
137  HistogramSizeType size;
138  size.Fill( m_NumberOfHistogramBins );
139  MeasurementVectorType theMin (0.);
140  theMin[NumberOfInputImages-2] = -otb::CONST_PI;
142  m_Histogram->Initialize( size, theMin, theMax );
143 
144  typename InternalSampleListType::Iterator angleIter
145  = m_AngleListFilter->GetOutputSampleList()->Begin();
146 
147  if ( m_AngleListFilter->GetOutputSampleList()->Begin()
148  == m_AngleListFilter->GetOutputSampleList()->End() )
149  {
150  throw itk::ExceptionObject( __FILE__, __LINE__,
151  "The value of threshold is too high so that there is no angle value to consider for unmixing",
152  ITK_LOCATION);
153  }
154 
155  while ( angleIter != m_AngleListFilter->GetOutputSampleList()->End() )
156  {
157  if ( !m_Histogram->IncreaseFrequencyOfMeasurement( angleIter.GetMeasurementVector(), 1 ) )
158  std::cerr << "Data out of bounds\n";
159 
160  ++angleIter;
161  }
162 
163 #if 0
164  for ( unsigned int index = 0; index < m_Histogram->GetSize()[0]; ++index )
165  std::cerr << index << "\t" << m_Histogram->GetMeasurementVector( index )[0]
166  << "\t" << m_Histogram->GetFrequency(index) << "\n";
167 #endif
168 
169  AngleListPointerType angles = AngleListType::New();
170 
171  for ( unsigned int id = 0; id < m_Histogram->Size(); ++id )
172  {
173  HistogramIndexType curIdx = m_Histogram->GetIndex(id);
174 
175  // Is this curIdx a local max ?
176  bool isLocalMax = true;
177  for ( unsigned int k = 0; k < m_Histogram->GetSize().Size(); ++k )
178  {
179  HistogramIndexType prevIdx = curIdx;
180  if ( prevIdx[k] == 0 )
181  prevIdx[k] = m_Histogram->GetSize(k)-1;
182  else
183  prevIdx[k] = curIdx[k]-1;
184 
185 
186  HistogramIndexType nextIdx = curIdx;
187  if ( static_cast<PrecisionType>( nextIdx[k] ) == m_Histogram->GetSize(k)-1 )
188  nextIdx[k] = 0;
189  else
190  nextIdx[k] = curIdx[k]+1;
191 
192  double freq_prev = m_Histogram->GetFrequency( prevIdx );
193  double freq_cur = m_Histogram->GetFrequency( curIdx );
194  double freq_next = m_Histogram->GetFrequency( nextIdx );
195 
196  if ( !( freq_prev < freq_cur && freq_cur > freq_next ) )
197  {
198  isLocalMax = false;
199  break;
200  }
201  }
202 
203  if ( isLocalMax )
204  {
205  angles->PushBack( m_Histogram->GetMeasurementVector( curIdx ) );
206  }
207  }
208 
209  m_NumberOfComponentsRequired = angles->Size();
210  m_AngleList = angles;
211 
212 #if 1
213  std::cout << m_NumberOfComponentsRequired << " sources found\n";
214  for ( unsigned int i = 0; i < angles->Size(); i++ )
215  {
216  std::cerr << "Source " << i << ":";
217  for ( unsigned int j = 0; j < m_Histogram->GetSize().Size(); ++j )
218  std::cout << "\t" << angles->GetMeasurementVector(i)[j];
219  std::cout << "\n";
220  }
221 #endif
222 }
223 
224 } // end of namespace otb
225 
226 #endif
const double CONST_PI
Definition: otbMath.h:48
HistogramType::MeasurementVectorType MeasurementVectorType
static Pointer New()
void Fill(TValue const &v)
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:58
SizeValueType GetSize(void) const
DataObject * GetInput(const DataObjectIdentifierType &key)
void SetInput(unsigned int i, const InputImageType *)
virtual void SetNthInput(DataObjectPointerArraySizeType num, DataObject *input)
InputImageType * GetInput(void)