OTB  9.0.0
Orfeo Toolbox
otbMultiChannelsPolarimetricSynthesisFilter.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 otbMultiChannelsPolarimetricSynthesisFilter_hxx
22 #define otbMultiChannelsPolarimetricSynthesisFilter_hxx
23 
24 #include <complex>
25 
27 #include "itkImageRegionIterator.h"
28 #include "itkProgressReporter.h"
29 #include "otbMath.h"
30 
31 namespace otb
32 {
33 
37 template <class TInputImage, class TOutputImage, class TFunction>
39  : m_PsiI(0.0), m_KhiI(0.0), m_PsiR(0.0), m_KhiR(0.0), m_Gain(1.0), m_Mode(0), m_EmissionH(false), m_EmissionV(false)
40 {
41  this->SetNumberOfRequiredInputs(1);
42  this->InPlaceOff();
44 }
46 
50 template <class TInputImage, class TOutputImage, class TFunction>
52 {
53  // do not call the superclass' implementation of this method since
54  // this filter allows the input the output to be of different dimensions
55 
56  // get pointers to the input and output
57  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
58  typename Superclass::InputImageConstPointer inputPtr = this->GetInput();
59 
60  if (!outputPtr || !inputPtr)
61  {
62  return;
63  }
64 
65  // Set the output image largest possible region. Use a RegionCopier
66  // so that the input and output images can be different dimensions.
67  OutputImageRegionType outputLargestPossibleRegion;
68  this->CallCopyInputRegionToOutputRegion(outputLargestPossibleRegion, inputPtr->GetLargestPossibleRegion());
69  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
70 
71  // Set the output spacing and origin
72  const itk::ImageBase<Superclass::InputImageDimension>* phyData;
73 
74  phyData = dynamic_cast<const itk::ImageBase<Superclass::InputImageDimension>*>(this->GetInput());
75 
76  if (phyData)
77  {
78  // Copy what we can from the image from spacing and origin of the input
79  // This logic needs to be augmented with logic that select which
80  // dimensions to copy
81  unsigned int i, j;
82  const typename InputImageType::SpacingType& inputSpacing = inputPtr->GetSignedSpacing();
83  const typename InputImageType::PointType& inputOrigin = inputPtr->GetOrigin();
84  const typename InputImageType::DirectionType& inputDirection = inputPtr->GetDirection();
85 
86  typename OutputImageType::SpacingType outputSpacing;
87  typename OutputImageType::PointType outputOrigin;
88  typename OutputImageType::DirectionType outputDirection;
89 
90  // copy the input to the output and fill the rest of the
91  // output with zeros.
92  for (i = 0; i < Superclass::InputImageDimension; ++i)
93  {
94  outputSpacing[i] = inputSpacing[i];
95  outputOrigin[i] = inputOrigin[i];
96  for (j = 0; j < Superclass::OutputImageDimension; ++j)
97  {
98  if (j < Superclass::InputImageDimension)
99  {
100  outputDirection[j][i] = inputDirection[j][i];
101  }
102  else
103  {
104  outputDirection[j][i] = 0.0;
105  }
106  }
107  }
108  for (; i < Superclass::OutputImageDimension; ++i)
109  {
110  outputSpacing[i] = 1.0;
111  outputOrigin[i] = 0.0;
112  for (j = 0; j < Superclass::OutputImageDimension; ++j)
113  {
114  if (j == i)
115  {
116  outputDirection[j][i] = 1.0;
117  }
118  else
119  {
120  outputDirection[j][i] = 0.0;
121  }
122  }
123  }
124 
125  // set the spacing and origin
126  outputPtr->SetSignedSpacing(outputSpacing);
127  outputPtr->SetOrigin(outputOrigin);
128  outputPtr->SetDirection(outputDirection);
129  }
130  else
131  {
132  // pointer could not be cast back down
133  itkExceptionMacro(<< "otb::MultiChannelsPolarimetricSynthesisFilter::GenerateOutputInformation "
134  << "cannot cast input to " << typeid(itk::ImageBase<Superclass::InputImageDimension>*).name());
135  }
136 }
137 
141 template <class TInputImage, class TOutputImage, class TFunction>
143  itk::ThreadIdType threadId)
144 {
145 
146  InputImagePointer inputPtr = this->GetInput();
147  OutputImagePointer outputPtr = this->GetOutput(0);
148 
149  // Define the portion of the input to walk for this thread, using
150  // the CallCopyOutputRegionToInputRegion method allows for the input
151  // and output images to be different dimensions
152  InputImageRegionType inputRegionForThread;
153  this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
154 
155  // Define the iterators
156  itk::ImageRegionConstIterator<TInputImage> inputIt(inputPtr, inputRegionForThread);
157  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
158 
159  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
160 
161  inputIt.GoToBegin();
162  outputIt.GoToBegin();
163 
164  ArchitectureType val = m_ArchitectureType->GetArchitectureType();
165 
166  // Computation with 4 channels
167  switch (val)
168  {
169  case HH_HV_VH_VV:
170  while (!inputIt.IsAtEnd())
171  {
172  outputIt.Set(m_Gain * GetFunctor()(inputIt.Get()[0], inputIt.Get()[1], inputIt.Get()[2], inputIt.Get()[3]));
173  ++inputIt;
174  ++outputIt;
175  progress.CompletedPixel(); // potential exception thrown here
176  }
177  break;
178 
179  // With 3 channels : HH HV VV or HH VH VV
180  case HH_HV_VV:
181  while (!inputIt.IsAtEnd())
182  {
183  outputIt.Set(m_Gain * GetFunctor()(inputIt.Get()[0], inputIt.Get()[1], inputIt.Get()[1], inputIt.Get()[2]));
184  ++inputIt;
185  ++outputIt;
186  progress.CompletedPixel(); // potential exception thrown here
187  }
188  break;
189 
190  // Only HH and HV are present
191  case HH_HV:
192  while (!inputIt.IsAtEnd())
193  {
194  outputIt.Set(m_Gain * GetFunctor()(inputIt.Get()[0], inputIt.Get()[1], 0, 0));
195  ++inputIt;
196  ++outputIt;
197  progress.CompletedPixel(); // potential exception thrown here
198  }
199  break;
200 
201  // Only VH and VV are present
202  case VH_VV:
203  while (!inputIt.IsAtEnd())
204  {
205  outputIt.Set(m_Gain * GetFunctor()(0, 0, inputIt.Get()[2], inputIt.Get()[3]));
206  ++inputIt;
207  ++outputIt;
208  progress.CompletedPixel(); // potential exception thrown here
209  }
210  break;
211 
212  default:
213  itkExceptionMacro("Unknown architecture : Polarimetric synthesis is impossible !");
214  return;
215  }
216 }
217 
221 template <class TInputImage, class TOutputImage, class TFunction>
223 {
224  ComplexArrayType AEi, AEr;
225 
227  double DTOR = CONST_PI_180;
228  double real, imag;
229 
230  real = std::cos(DTOR * m_PsiI) * std::cos(DTOR * m_KhiI);
231  imag = -std::sin(DTOR * m_PsiI) * std::sin(DTOR * m_KhiI);
232  ComplexType Ei0(real, imag);
233 
234  real = std::sin(DTOR * m_PsiI) * std::cos(DTOR * m_KhiI);
235  imag = std::cos(DTOR * m_PsiI) * std::sin(DTOR * m_KhiI);
236  ComplexType Ei1(real, imag);
237 
238  real = std::cos(DTOR * m_PsiR) * std::cos(DTOR * m_KhiR);
239  imag = -std::sin(DTOR * m_PsiR) * std::sin(DTOR * m_KhiR);
240  ComplexType Er0(real, imag);
241 
242  real = std::sin(DTOR * m_PsiR) * std::cos(DTOR * m_KhiR);
243  imag = std::cos(DTOR * m_PsiR) * std::sin(DTOR * m_KhiR);
244  ComplexType Er1(real, imag);
245 
246  AEi[0] = Ei0;
247  AEi[1] = Ei1;
248  AEr[0] = Er0;
249  AEr[1] = Er1;
250 
251  this->SetEi(AEi);
252  this->SetEr(AEr);
253 }
254 
258 template <class TInputImage, class TOutputImage, class TFunction>
260 {
261 
262  ArchitectureType val = m_ArchitectureType->GetArchitectureType();
263 
264  switch (val)
265  {
266 
267  case HH_HV_VH_VV:
268  break;
269  case HH_HV_VV:
270  break;
271  case HH_VH_VV:
272  break;
273  // Only HH and HV are present
274  case HH_HV:
275 
276  // Forcing KhiI=0 PsiI=0
277  this->SetKhiI(0);
278  this->SetPsiI(0);
279  break;
280 
281  // Only VH and VV are present
282  case VH_VV:
283 
284  // Forcing KhiI=0 PsiI=90
285  this->SetKhiI(0);
286  this->SetPsiI(90);
287  break;
288 
289  default:
290  itkExceptionMacro("Unknown architecture : Polarimetric synthesis is impossible !!");
291  return;
292  }
293 
294  if (GetMode() == 1)
295  ForceCoPolar();
296  else if (GetMode() == 2)
297  ForceCrossPolar();
298 }
299 
303 template <class TInputImage, class TOutputImage, class TFunction>
305 {
306 
307  int NumberOfImages = this->GetInput()->GetNumberOfComponentsPerPixel();
308 
309  // First Part. Determine the kind of architecture of the input picture
310  m_ArchitectureType->DetermineArchitecture(NumberOfImages, GetEmissionH(), GetEmissionV());
311 
312  // Second Part. Verify and force the inputs
313  VerifyAndForceInputs();
314 
315  // Third Part. Estimation of the incident field Ei and the reflected field Er
316  ComputeElectromagneticFields();
317 }
318 
322 template <class TInputImage, class TOutputImage, class TFunction>
324 {
325  SetPsiR(m_PsiI);
326  SetKhiR(m_KhiI);
327 }
329 
333 template <class TInputImage, class TOutputImage, class TFunction>
335 {
336  SetPsiR(m_PsiI + 90);
337  SetKhiR(-m_KhiI);
338  SetMode(2);
339 }
341 
345 template <class TInputImage, class TOutputImage, class TFunction>
347 {
348  this->Superclass::PrintSelf(os, indent);
349  os << indent << "PsiI: " << m_PsiI << std::endl;
350  os << indent << "KhiI: " << m_KhiI << std::endl;
351  os << indent << "PsiR: " << m_PsiR << std::endl;
352  os << indent << "KhiR: " << m_KhiR << std::endl;
354 
355  os << indent << "Ei0 im: " << m_Ei[0].imag() << std::endl;
356  os << indent << "Ei0 re: " << m_Ei[0].real() << std::endl;
357  os << indent << "Ei1 im: " << m_Ei[1].imag() << std::endl;
358  os << indent << "Ei1 re: " << m_Ei[1].real() << std::endl;
359 
360  os << indent << "Er0 im: " << m_Er[0].imag() << std::endl;
361  os << indent << "Er0 re: " << m_Er[0].real() << std::endl;
362  os << indent << "Er1 im: " << m_Er[1].imag() << std::endl;
363  os << indent << "Er1 re: " << m_Er[1].real() << std::endl;
364 }
365 
366 } // end namespace otb
367 
368 #endif
otb::MultiChannelsPolarimetricSynthesisFilter::ForceCoPolar
void ForceCoPolar()
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:323
otb::MultiChannelsPolarimetricSynthesisFilter::OutputImageRegionType
OutputImageType::RegionType OutputImageRegionType
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:74
otb::MultiChannelsPolarimetricSynthesisFilter::GenerateOutputInformation
void GenerateOutputInformation() override
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:51
otb::MultiChannelsPolarimetricSynthesisFilter::MultiChannelsPolarimetricSynthesisFilter
MultiChannelsPolarimetricSynthesisFilter()
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:38
otb::HH_HV_VH_VV
@ HH_HV_VH_VV
Definition: otbPolarimetricData.h:33
otbMath.h
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::MultiChannelsPolarimetricSynthesisFilter::ComplexType
std::complex< double > ComplexType
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:76
otb::MultiChannelsPolarimetricSynthesisFilter::ComputeElectromagneticFields
void ComputeElectromagneticFields()
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:222
otb::MultiChannelsPolarimetricSynthesisFilter::InputImageRegionType
InputImageType::RegionType InputImageRegionType
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:70
otb::MultiChannelsPolarimetricSynthesisFilter::ThreadedGenerateData
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:142
otb::MultiChannelsPolarimetricSynthesisFilter::InputImagePointer
InputImageType::ConstPointer InputImagePointer
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:69
otb::MultiChannelsPolarimetricSynthesisFilter::ForceCrossPolar
void ForceCrossPolar()
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:334
otb::CONST_PI_180
constexpr double CONST_PI_180
Definition: otbMath.h:56
otbMultiChannelsPolarimetricSynthesisFilter.h
otb::ArchitectureType
ArchitectureType
This enumeration describes the different architectures we can find in polarimetry....
Definition: otbPolarimetricData.h:33
otb::MultiChannelsPolarimetricSynthesisFilter::OutputImagePointer
OutputImageType::Pointer OutputImagePointer
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:73
otb::PolarimetricData::New
static Pointer New()
otb::HH_HV_VV
@ HH_HV_VV
Definition: otbPolarimetricData.h:33
otb::HH_VH_VV
@ HH_VH_VV
Definition: otbPolarimetricData.h:33
otb::MultiChannelsPolarimetricSynthesisFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:346
otb::MultiChannelsPolarimetricSynthesisFilter::BeforeThreadedGenerateData
void BeforeThreadedGenerateData() override
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:304
otb::MultiChannelsPolarimetricSynthesisFilter::m_ArchitectureType
PolarimetricData::Pointer m_ArchitectureType
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:235
otb::VH_VV
@ VH_VV
Definition: otbPolarimetricData.h:33
otb::HH_HV
@ HH_HV
Definition: otbPolarimetricData.h:33
otb::MultiChannelsPolarimetricSynthesisFilter::ComplexArrayType
itk::FixedArray< ComplexType, 2 > ComplexArrayType
Definition: otbMultiChannelsPolarimetricSynthesisFilter.h:77
otb::MultiChannelsPolarimetricSynthesisFilter::VerifyAndForceInputs
void VerifyAndForceInputs()
Definition: otbMultiChannelsPolarimetricSynthesisFilter.hxx:259