OTB  6.7.0
Orfeo Toolbox
otbSparseWvltToAngleMapperListFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 otbSparseWvltToAngleMapperListFilter_hxx
22 #define otbSparseWvltToAngleMapperListFilter_hxx
24 
25 #include <vnl/vnl_math.h>
26 
27 #include "itkProgressReporter.h"
28 
29 namespace otb {
30 
31 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
34 {
35  this->SetNumberOfRequiredInputs(NumberOfInputImages);
36 
37  // Generate the output sample list
38  typename OutputSampleListObjectType::Pointer outputPtr =
39  static_cast< OutputSampleListObjectType * >(this->MakeOutput(0).GetPointer());
40  this->ProcessObject::SetNthOutput(0, outputPtr.GetPointer());
41 
42  m_ThresholdValue = static_cast<ValueType>( 10. );
43 }
44 
45 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
46 void
48 ::SetInput ( unsigned int i, const InputImageListType * img )
49 {
51  const_cast< InputImageListType * >( img ) );
52 }
53 
54 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
55 const TInputImageList *
57 ::GetInput ( unsigned int i ) const
58 {
59  if ( i >= this->GetNumberOfInputs() )
60  {
61  return nullptr;
62  }
63 
64  return static_cast<const InputImageListType * >
65  (this->itk::ProcessObject::GetInput(i) );
66 }
67 
68 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
70 ::DataObjectPointer
73 {
74  typename OutputSampleListObjectType::Pointer outputPtr = OutputSampleListObjectType::New();
75  OutputSampleListPointer outputSampleList = OutputSampleListType::New();
76  outputPtr->Set(outputSampleList);
77  return static_cast<DataObjectPointer>(outputPtr);
78 }
79 
80 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
82 ::OutputSampleListType *
85 {
86  typename OutputSampleListObjectType::Pointer dataObjectPointer
87  = static_cast<OutputSampleListObjectType * > (this->ProcessObject::GetOutput(0) );
88  return const_cast<OutputSampleListType *>(dataObjectPointer->Get());
89 }
90 
91 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
93 ::OutputSampleListObjectType *
96 {
97  return static_cast<OutputSampleListObjectType * >
98  (this->ProcessObject::GetOutput(0) );
99 }
100 
101 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
102 void
104 ::PrintSelf(std::ostream& os, itk::Indent indent) const
105 {
106  Superclass::PrintSelf(os, indent);
107  os << indent << "Threshold : " << m_ThresholdValue << "\n";
108 }
109 
110 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
111 void
114 {
115  InputImageListConstIteratorVectorType it ( NumberOfInputImages );
116  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
117  it[i] = this->GetInput(i)->Begin();
118 
119  OutputSampleListType * outputList = this->GetOutputSampleList();
120 
121  // Set-up progress reporting
122  itk::ProgressReporter progress(this, 0, this->GetInput(0)->Size());
123 
124  bool iteratorsNotAtEnd = true;
125  while ( iteratorsNotAtEnd )
126  {
127  ImageConstIteratorVectorType imgIt ( NumberOfInputImages );
128  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
129  {
130  imgIt[i] = ImageConstIteratorType( it[i].Get(), it[i].Get()->GetLargestPossibleRegion() );
131  imgIt[i].GoToBegin();
132  }
133 
134  bool localIteratorsNotAtEnd = true;
135  while ( localIteratorsNotAtEnd )
136  {
137  if ( IsToGenerate( imgIt ) )
138  {
139  outputList->PushBack( GenerateData( imgIt ) );
140  }
141 
142  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
143  {
144  ++(imgIt[i]);
145  if ( imgIt[i].IsAtEnd() )
146  localIteratorsNotAtEnd = false;
147  }
148  }
149 
150  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
151  {
152  ++(it[i]);
153  if ( it[i] != this->GetInput(i)->End() )
154  iteratorsNotAtEnd = false;
155  }
156 
157  progress.CompletedPixel();
158  }
159 }
160 
161 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
162 bool
165 {
166  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
167  {
168  if ( it[i].Get() < m_ThresholdValue )
169  return false;
170  }
171 
172  return true;
173 }
174 
175 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
177 ::OutputMeasurementVectorType
180 {
181  return FromEuclideanToSphericalSpace( it );
182 }
183 
184 template < class TInputImageList, class TOutputSampleList, unsigned int VNbInputImages >
186 ::OutputMeasurementVectorType
189 ( const ImageConstIteratorVectorType & it ) const
190 {
191  // First, get the modulus of the vector
192  double modulus = 0;
193  for ( unsigned int i = 0; i < NumberOfInputImages; ++i )
194  {
195  modulus += std::pow( static_cast<double>( it[i].Get() ), 2. );
196  }
197  // The modulus cannot be nul since it is over the threshold...
198  modulus = std::sqrt( modulus );
199 
200  // FIXME test if NaN nor infinite
201 
203  angle.Fill( static_cast< OutputValueType >( 0. ) );
204 
205  if ( NumberOfInputImages == 2 )
206  {
207  if ( it[1].Get() < 0 )
208  {
209  angle[0] = std::acos( it[0].Get() / modulus );
210  }
211  else
212  {
213  angle[0] = - std::acos( it[0].Get() / modulus );
214  }
215  }
216  else
217  {
218  for ( unsigned int k = 0; k < angle.Size()-1; ++k )
219  {
220  double phase = std::acos( it[k].Get() / modulus );
221  angle[k] = phase;
222  modulus *= std::sin( phase );
223 
224  // FIXME test also if not finite
225  if ( modulus < 1e-5 )
226  {
227  while ( ++k < angle.Size() )
228  angle[k] = 0.;
229  return angle;
230  }
231  }
232 
233  /*
234  * With this sign modification, we can put the same
235  * images for all the components and recover the good direction
236  */
237  double sign = NumberOfInputImages == 3 ? -1. : 1.;
238  if ( it[ NumberOfInputImages-1 ].Get() < 0 )
239  {
240  angle[ angle.Size()-1 ] = sign * std::acos( it[ NumberOfInputImages-2 ].Get() / modulus );
241  }
242  else
243  {
244  angle[ angle.Size()-1 ] = - sign * std::acos( it[ NumberOfInputImages-2 ].Get() / modulus );
245  }
246  }
247 
248  return angle;
249 }
250 
251 } // end of namespace otb
252 
253 #endif
254 
virtual bool IsToGenerate(const ImageConstIteratorVectorType &) const
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
const InputImageListType * GetInput(unsigned int i) const
OutputSampleListType::MeasurementVectorType OutputMeasurementVectorType
ObjectType * GetPointer() const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
DataObject * GetInput(const DataObjectIdentifierType &key)
std::vector< InputImageListConstIteratorType > InputImageListConstIteratorVectorType
std::vector< ImageConstIteratorType > ImageConstIteratorVectorType
virtual OutputMeasurementVectorType FromEuclideanToSphericalSpace(const ImageConstIteratorVectorType &) const
void SetInput(unsigned int i, const InputImageListType *)
virtual void SetNthInput(DataObjectPointerArraySizeType num, DataObject *input)
This class select N-uple join-wvlt coeff for sparse unmixing.