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