OTB  9.0.0
Orfeo Toolbox
otbSparseUnmixingImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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_hxx
22 #define otbSparseUnmixingImageFilter_hxx
23 
25 
26 #include "otbMath.h"
27 #include "itkProcessObject.h"
28 
29 namespace otb
30 {
31 
32 template <class TInputImage, class TOutputImage, unsigned int VNbInputImage, class TPrecision, Wavelet::Wavelet TMotherWaveletOperator>
34 {
35  this->SetNumberOfRequiredInputs(NumberOfInputImages);
36  this->SetNumberOfRequiredOutputs(1);
37  this->SetNthOutput(0, OutputImageListType::New());
38 
39  m_NumberOfComponentsRequired = 0;
40  m_NumberOfHistogramBins = 100;
41  m_AngleList = AngleListType::New();
42 
43  m_WvltFilterList = WvltFilterListType::New();
44  m_WvltFilterList->Resize(NumberOfInputImages);
45  for (unsigned int i = 0; i < NumberOfInputImages; ++i)
46  {
47  m_WvltFilterList->SetNthElement(i, WvltFilterType::New());
48  m_WvltFilterList->GetNthElement(i)->SetNumberOfDecompositions(2);
49  m_WvltFilterList->GetNthElement(i)->SetSubsampleImageFactor(1);
50  }
51 
52  m_AngleListFilter = AngleListFilterType::New();
53  m_AngleListFilter->SetThresholdValue(10.);
54 
55  m_Histogram = HistogramType::New();
56  m_Transformer = TransformFilterType::New();
57 }
58 
59 template <class TInputImage, class TOutputImage, unsigned int VNbInputImage, class TPrecision, Wavelet::Wavelet TMotherWaveletOperator>
61  const InputImageType* img)
62 {
63  this->itk::ProcessObject::SetNthInput(i, const_cast<InputImageType*>(img));
64 }
65 
66 template <class TInputImage, class TOutputImage, unsigned int VNbInputImage, class TPrecision, Wavelet::Wavelet TMotherWaveletOperator>
68 {
69  if (i >= this->GetNumberOfInputs())
70  {
71  return nullptr;
72  }
73 
74  return static_cast<const InputImageType*>(this->itk::ProcessObject::GetInput(i));
75 }
76 
77 template <class TInputImage, class TOutputImage, unsigned int VNbInputImage, class TPrecision, Wavelet::Wavelet TMotherWaveletOperator>
79 {
80  itk::ProgressAccumulator::Pointer progress = itk::ProgressAccumulator::New();
81  progress->SetMiniPipelineFilter(this);
82 
83  for (unsigned int i = 0; i < NumberOfInputImages; ++i)
84  {
85  m_WvltFilterList->GetNthElement(i)->SetInput(this->GetInput(i));
86  progress->RegisterInternalFilter(m_WvltFilterList->GetNthElement(i), 0.5f / static_cast<float>(NumberOfInputImages));
87  m_WvltFilterList->GetNthElement(i)->Update();
88 
89  m_AngleListFilter->SetInput(i, m_WvltFilterList->GetNthElement(i)->GetOutput());
90  }
91 
92  progress->RegisterInternalFilter(m_AngleListFilter, 0.25f);
93  m_AngleListFilter->Update();
94 
95  GenerateNumberOfComponentsRequired();
96  otbMsgDebugMacro(<< m_NumberOfComponentsRequired << " sources found\n");
97 
98  for (unsigned int i = 0; i < NumberOfInputImages; ++i)
99  m_Transformer->SetInput(i, this->GetInput(i));
100  m_Transformer->SetAngleList(m_AngleList);
101  progress->RegisterInternalFilter(m_Transformer, 0.25f);
102  m_Transformer->Update();
103 
104  this->GetOutput()->Resize(m_Transformer->GetOutput()->Size());
105  for (unsigned int i = 0; i < this->GetOutput()->Size(); ++i)
106  this->GetOutput()->SetNthElement(i, m_Transformer->GetOutput()->GetNthElement(i));
107 }
108 
109 template <class TInputImage, class TOutputImage, unsigned int VNbInputImage, class TPrecision, Wavelet::Wavelet TMotherWaveletOperator>
111 {
117  HistogramSizeType size;
118  size.Fill(m_NumberOfHistogramBins);
119  MeasurementVectorType theMin(0.);
120  theMin[NumberOfInputImages - 2] = -otb::CONST_PI;
122  m_Histogram->Initialize(size, theMin, theMax);
123 
124  typename InternalSampleListType::Iterator angleIter = m_AngleListFilter->GetOutputSampleList()->Begin();
125 
126  if (m_AngleListFilter->GetOutputSampleList()->Begin() == m_AngleListFilter->GetOutputSampleList()->End())
127  {
128  throw itk::ExceptionObject(__FILE__, __LINE__, "The value of threshold is too high so that there is no angle value to consider for unmixing", ITK_LOCATION);
129  }
130 
131  while (angleIter != m_AngleListFilter->GetOutputSampleList()->End())
132  {
133  if (!m_Histogram->IncreaseFrequencyOfMeasurement(angleIter.GetMeasurementVector(), 1))
134  std::cerr << "Data out of bounds\n";
135 
136  ++angleIter;
137  }
138 
139 #if 0
140  for ( unsigned int index = 0; index < m_Histogram->GetSize()[0]; ++index )
141  std::cerr << index << "\t" << m_Histogram->GetMeasurementVector( index )[0]
142  << "\t" << m_Histogram->GetFrequency(index) << "\n";
143 #endif
144 
145  AngleListPointerType angles = AngleListType::New();
146 
147  for (unsigned int id = 0; id < m_Histogram->Size(); ++id)
148  {
149  HistogramIndexType curIdx = m_Histogram->GetIndex(id);
150 
151  // Is this curIdx a local max ?
152  bool isLocalMax = true;
153  for (unsigned int k = 0; k < m_Histogram->GetSize().Size(); ++k)
154  {
155  HistogramIndexType prevIdx = curIdx;
156  if (prevIdx[k] == 0)
157  prevIdx[k] = m_Histogram->GetSize(k) - 1;
158  else
159  prevIdx[k] = curIdx[k] - 1;
160 
161 
162  HistogramIndexType nextIdx = curIdx;
163  if (static_cast<PrecisionType>(nextIdx[k]) == m_Histogram->GetSize(k) - 1)
164  nextIdx[k] = 0;
165  else
166  nextIdx[k] = curIdx[k] + 1;
167 
168  double freq_prev = m_Histogram->GetFrequency(prevIdx);
169  double freq_cur = m_Histogram->GetFrequency(curIdx);
170  double freq_next = m_Histogram->GetFrequency(nextIdx);
171 
172  if (!(freq_prev < freq_cur && freq_cur > freq_next))
173  {
174  isLocalMax = false;
175  break;
176  }
177  }
178 
179  if (isLocalMax)
180  {
181  angles->PushBack(m_Histogram->GetMeasurementVector(curIdx));
182  }
183  }
184 
185  m_NumberOfComponentsRequired = angles->Size();
186  m_AngleList = angles;
187 
188 #if 1
189  std::cout << m_NumberOfComponentsRequired << " sources found\n";
190  for (unsigned int i = 0; i < angles->Size(); i++)
191  {
192  std::cerr << "Source " << i << ":";
193  for (unsigned int j = 0; j < m_Histogram->GetSize().Size(); ++j)
194  std::cout << "\t" << angles->GetMeasurementVector(i)[j];
195  std::cout << "\n";
196  }
197 #endif
198 }
199 
200 } // end of namespace otb
201 
202 #endif
otbSparseUnmixingImageFilter.h
otb::CONST_PI
constexpr double CONST_PI
Definition: otbMath.h:49
otb::SparseUnmixingImageFilter::PrecisionType
TPrecision PrecisionType
Definition: otbSparseUnmixingImageFilter.h:80
otb::SparseUnmixingImageFilter::GenerateNumberOfComponentsRequired
virtual void GenerateNumberOfComponentsRequired()
Definition: otbSparseUnmixingImageFilter.hxx:110
otbMath.h
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::ImageToImageListFilter::GetInput
InputImageType * GetInput(void)
otb::SparseUnmixingImageFilter::HistogramSizeType
HistogramType::SizeType HistogramSizeType
Definition: otbSparseUnmixingImageFilter.h:104
otb::SparseUnmixingImageFilter::SetInput
void SetInput(unsigned int i, const InputImageType *)
Definition: otbSparseUnmixingImageFilter.hxx:60
otb::SparseUnmixingImageFilter::AngleListPointerType
AngleListType::Pointer AngleListPointerType
Definition: otbSparseUnmixingImageFilter.h:97
otb::SparseUnmixingImageFilter::SparseUnmixingImageFilter
SparseUnmixingImageFilter()
Definition: otbSparseUnmixingImageFilter.hxx:33
otb::SparseUnmixingImageFilter::MeasurementVectorType
HistogramType::MeasurementVectorType MeasurementVectorType
Definition: otbSparseUnmixingImageFilter.h:106
otbMsgDebugMacro
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:62
otb::SparseUnmixingImageFilter::HistogramIndexType
HistogramType::IndexType HistogramIndexType
Definition: otbSparseUnmixingImageFilter.h:105
otb::SparseUnmixingImageFilter::GenerateData
void GenerateData() override
Definition: otbSparseUnmixingImageFilter.hxx:78
otb::ImageToImageListFilter::InputImageType
TInputImage InputImageType
Definition: otbImageToImageListFilter.h:52