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