Orfeo Toolbox  3.16
otbPolarimetricSynthesisFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbPolarimetricSynthesisFilter_txx
19 #define __otbPolarimetricSynthesisFilter_txx
20 
21 #include "otbPolarimetricSynthesisFilter.h"
22 #include "otbMath.h"
23 
24 namespace otb
25 {
26 
30 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
31  class TFunction>
32 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
33 ::PolarimetricSynthesisFilter()
34 {
35  this->SetNumberOfRequiredInputs(0);
36  this->SetNumberOfInputs(4);
37  SetMode(0);
38  SetGain(1);
39  m_PresentInputImages[0] = false;
40  m_PresentInputImages[1] = false;
41  m_PresentInputImages[2] = false;
42  m_PresentInputImages[3] = false;
43  m_ArchitectureType = PolarimetricData::New();
44 }
45 
46 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
47  class TFunction>
48 void
49 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
50 ::GenerateOutputInformation()
51 {
52  // if HH is set, use HH to generate output information
53  if (m_PresentInputImages[0])
54  {
55  this->GetOutput()->CopyInformation(this->GetInput(0));
56  }
57  // else, use VH
58  else if (m_PresentInputImages[2])
59  {
60  this->GetOutput()->CopyInformation(this->GetInput(2));
61  }
62  else
63  {
64  itkExceptionMacro(<< "Bad input polarization images: neither HH image nor VH image is set!");
65  }
66 }
67 
71 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
72  class TFunction>
73 void
74 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
75 ::SetInputHH(const TInputImageHH * image)
76 {
77  // Process object is not const-correct so the const casting is required.
78  this->SetInput1(image);
79  m_PresentInputImages[0] = true;
80 }
81 
85 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
86  class TFunction>
87 void
88 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
89 ::SetInputHV(const TInputImageHV * image)
90 {
91  this->SetInput2(image);
92  m_PresentInputImages[1] = true;
93 }
94 
98 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
99  class TFunction>
100 void
101 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
102 ::SetInputVH(const TInputImageVH * image)
103 {
104  this->SetInput3(image);
105  m_PresentInputImages[2] = true;
106 }
107 
111 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
112  class TFunction>
113 void
114 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
115 ::SetInputVV(const TInputImageVV * image)
116 {
117  this->SetInput4(image);
118  m_PresentInputImages[3] = true;
119 }
120 
124 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
125  class TFunction>
126 void
127 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
128 ::PrintSelf(std::ostream& os, itk::Indent indent) const
129 {
130  Superclass::PrintSelf(os, indent);
131 
132  os << "PsiI: " << m_PsiI << std::endl;
133  os << "KhiI: " << m_KhiI << std::endl;
134  os << "PsiR: " << m_PsiR << std::endl;
135  os << "KhiR: " << m_KhiR << std::endl;
136 
137  os << "Ei0 im: " << m_Ei[0].imag() << std::endl;
138  os << "Ei0 re: " << m_Ei[0].real() << std::endl;
139  os << "Ei1 im: " << m_Ei[1].imag() << std::endl;
140  os << "Ei1 re: " << m_Ei[1].real() << std::endl;
141 
142  os << "Er0 im: " << m_Er[0].imag() << std::endl;
143  os << "Er0 re: " << m_Er[0].real() << std::endl;
144  os << "Er1 im: " << m_Er[1].imag() << std::endl;
145  os << "Er1 re: " << m_Er[1].real() << std::endl;
146 }
147 
151 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
152  class TFunction>
153 void
154 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
155 ::ForceCoPolar()
156 {
157  SetPsiR(m_PsiI);
158  SetKhiR(m_KhiI);
159  SetMode(1);
160 }
161 
165 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
166  class TFunction>
167 void
168 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
169 ::ForceCrossPolar()
170 {
171  SetPsiR(m_PsiI + 90);
172  SetKhiR(-m_KhiI);
173  SetMode(2);
174 }
175 
179 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
180  class TFunction>
181 void
182 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
183 ::VerifyAndForceInputs()
184 {
185 
186  switch (m_ArchitectureType->GetArchitectureType())
187  {
188  case HH_HV:
189  // Forcing KhiI=0 PsiI=0
190  this->SetKhiI(0);
191  this->SetPsiI(0);
192  break;
193 
194  case VH_VV:
195  // Forcing KhiI=0 PsiI=90
196  this->SetKhiI(0);
197  this->SetPsiI(90);
198  break;
199 
200  case HH_VV:
201  itkExceptionMacro("Only the HH and VV channels are available : Polarimetric synthesis is not supported !");
202  return;
203 
204  case UNKNOWN_ARCHITECTURE:
205  itkExceptionMacro("Unknown architecture : Polarimetric synthesis is impossible !");
206  return;
207 
208  default:
209  break;
210 
211  }
212 
213  if (GetMode() == 1) ForceCoPolar();
214  else if (GetMode() == 2) ForceCrossPolar();
215 
216 }
217 
221 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
222  class TFunction>
223 void
224 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
225 ::ComputeElectromagneticFields()
226 {
227  ComplexArrayType AEi, AEr;
228 
230  double DTOR = CONST_PI_180;
231  double real, imag;
232 
233  real = vcl_cos(DTOR * m_PsiI) * vcl_cos(DTOR * m_KhiI);
234  imag = -vcl_sin(DTOR * m_PsiI) * vcl_sin(DTOR * m_KhiI);
235  ComplexType Ei0(real, imag);
236 
237  real = vcl_sin(DTOR * m_PsiI) * vcl_cos(DTOR * m_KhiI);
238  imag = vcl_cos(DTOR * m_PsiI) * vcl_sin(DTOR * m_KhiI);
239  ComplexType Ei1(real, imag);
240 
241  real = vcl_cos(DTOR * m_PsiR) * vcl_cos(DTOR * m_KhiR);
242  imag = -vcl_sin(DTOR * m_PsiR) * vcl_sin(DTOR * m_KhiR);
243  ComplexType Er0(real, imag);
244 
245  real = vcl_sin(DTOR * m_PsiR) * vcl_cos(DTOR * m_KhiR);
246  imag = vcl_cos(DTOR * m_PsiR) * vcl_sin(DTOR * m_KhiR);
247  ComplexType Er1(real, imag);
248 
249  AEi[0] = Ei0;
250  AEi[1] = Ei1;
251  AEr[0] = Er0;
252  AEr[1] = Er1;
253 
254  this->SetEi(AEi);
255  this->SetEr(AEr);
256 
257 }
258 
262 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
263  class TFunction>
264 void
265 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
266 ::BeforeThreadedGenerateData()
267 {
268 
269  // First Part. Determine the kind of architecture
270  m_ArchitectureType->DetermineArchitecture(m_PresentInputImages);
271 
272  // Second Part. Verify and force the inputs
273  VerifyAndForceInputs();
274 
275  // Third Part. Estimation of the incident field Ei and the reflected field Er
276  ComputeElectromagneticFields();
277 
278 }
279 
283 template <class TInputImageHH, class TInputImageHV, class TInputImageVH, class TInputImageVV, class TOutputImage,
284  class TFunction>
285 void
286 PolarimetricSynthesisFilter<TInputImageHH, TInputImageHV, TInputImageVH, TInputImageVV, TOutputImage, TFunction>
287 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
288  int threadId)
289 {
290 
291  OutputImagePointer outputPtr = this->GetOutput(0);
292  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
293  outputIt.GoToBegin();
294 
295  switch (m_ArchitectureType->GetArchitectureType())
296  {
297  // With 4 channels :
298  case HH_HV_VH_VV:
299  {
300  // We use dynamic_cast since inputs are stored as DataObjects. The
301  // ImageToImageFilter::GetInput(int) always returns a pointer to a
302  // TInputImage1 so it cannot be used for the second or third input.
303  HHInputImagePointer inputPtrHH
304  = dynamic_cast<const TInputImageHH*>((itk::ProcessObject::GetInput(0)));
305  HVInputImagePointer inputPtrHV
306  = dynamic_cast<const TInputImageHV*>((itk::ProcessObject::GetInput(1)));
307  VHInputImagePointer inputPtrVH
308  = dynamic_cast<const TInputImageVH*>((itk::ProcessObject::GetInput(2)));
309  VVInputImagePointer inputPtrVV
310  = dynamic_cast<const TInputImageVV*>((itk::ProcessObject::GetInput(3)));
311 
312  itk::ImageRegionConstIterator<TInputImageHH> inputItHH(inputPtrHH, outputRegionForThread);
313  itk::ImageRegionConstIterator<TInputImageHV> inputItHV(inputPtrHV, outputRegionForThread);
314  itk::ImageRegionConstIterator<TInputImageVH> inputItVH(inputPtrVH, outputRegionForThread);
315  itk::ImageRegionConstIterator<TInputImageVV> inputItVV(inputPtrVV, outputRegionForThread);
316 
317  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
318 
319  inputItHH.GoToBegin();
320  inputItHV.GoToBegin();
321  inputItVH.GoToBegin();
322  inputItVV.GoToBegin();
323 
324  while (!inputItHH.IsAtEnd())
325  {
326  outputIt.Set(m_Gain *
327  Superclass::GetFunctor() (inputItHH.Get(), inputItHV.Get(), inputItVH.Get(), inputItVV.Get()));
328  ++inputItHH;
329  ++inputItHV;
330  ++inputItVH;
331  ++inputItVV;
332  ++outputIt;
333  progress.CompletedPixel(); // potential exception thrown here
334  }
335  break;
336  }
337  // With 3 channels : HH HV VV
338  case HH_HV_VV:
339 
340  {
341  HHInputImagePointer inputPtrHH
342  = dynamic_cast<const TInputImageHH*>((itk::ProcessObject::GetInput(0)));
343  HVInputImagePointer inputPtrHV
344  = dynamic_cast<const TInputImageHV*>((itk::ProcessObject::GetInput(1)));
345  VVInputImagePointer inputPtrVV
346  = dynamic_cast<const TInputImageVV*>((itk::ProcessObject::GetInput(3)));
347 
348  itk::ImageRegionConstIterator<TInputImageHH> inputItHH(inputPtrHH, outputRegionForThread);
349  itk::ImageRegionConstIterator<TInputImageHV> inputItHV(inputPtrHV, outputRegionForThread);
350  itk::ImageRegionConstIterator<TInputImageVV> inputItVV(inputPtrVV, outputRegionForThread);
351  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
352 
353  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
354 
355  inputItHH.GoToBegin();
356  inputItHV.GoToBegin();
357  inputItVV.GoToBegin();
358 
359  while (!inputItHH.IsAtEnd())
360  {
361  outputIt.Set(m_Gain *
362  Superclass::GetFunctor() (inputItHH.Get(), inputItHV.Get(), inputItHV.Get(), inputItVV.Get()));
363  ++inputItHH;
364  ++inputItHV;
365  ++inputItVV;
366  ++outputIt;
367  progress.CompletedPixel(); // potential exception thrown here
368  }
369  break;
370  }
371  // With 3 channels : HH VH VV
372  case HH_VH_VV:
373  {
374 
375  HHInputImagePointer inputPtrHH
376  = dynamic_cast<const TInputImageHH*>((itk::ProcessObject::GetInput(0)));
377  VHInputImagePointer inputPtrVH
378  = dynamic_cast<const TInputImageVH*>((itk::ProcessObject::GetInput(2)));
379  VVInputImagePointer inputPtrVV
380  = dynamic_cast<const TInputImageVV*>((itk::ProcessObject::GetInput(3)));
381  OutputImagePointer outputPtr = this->GetOutput(0);
382 
383  itk::ImageRegionConstIterator<TInputImageHH> inputItHH(inputPtrHH, outputRegionForThread);
384  itk::ImageRegionConstIterator<TInputImageVH> inputItVH(inputPtrVH, outputRegionForThread);
385  itk::ImageRegionConstIterator<TInputImageVV> inputItVV(inputPtrVV, outputRegionForThread);
386 
387  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
388 
389  inputItHH.GoToBegin();
390  inputItVH.GoToBegin();
391  inputItVV.GoToBegin();
392 
393  while (!inputItHH.IsAtEnd())
394  {
395  outputIt.Set(m_Gain *
396  Superclass::GetFunctor() (inputItHH.Get(), inputItVH.Get(), inputItVH.Get(), inputItVV.Get()));
397  ++inputItHH;
398  ++inputItVH;
399  ++inputItVV;
400  ++outputIt;
401  progress.CompletedPixel(); // potential exception thrown here
402  }
403  break;
404  }
405  // With 2 channels : HH HV
406  case HH_HV:
407  {
408 
409  HHInputImagePointer inputPtrHH
410  = dynamic_cast<const TInputImageHH*>((itk::ProcessObject::GetInput(0)));
411  HVInputImagePointer inputPtrHV
412  = dynamic_cast<const TInputImageHV*>((itk::ProcessObject::GetInput(1)));
413 
414  itk::ImageRegionConstIterator<TInputImageHH> inputItHH(inputPtrHH, outputRegionForThread);
415  itk::ImageRegionConstIterator<TInputImageHV> inputItHV(inputPtrHV, outputRegionForThread);
416  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
417 
418  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
419 
420  inputItHH.GoToBegin();
421  inputItHV.GoToBegin();
422 
423  while (!inputItHH.IsAtEnd())
424  {
425  outputIt.Set(m_Gain * Superclass::GetFunctor() (inputItHH.Get(), inputItHV.Get(), 0, 0));
426  ++inputItHH;
427  ++inputItHV;
428  ++outputIt;
429  progress.CompletedPixel(); // potential exception thrown here
430  }
431  break;
432  }
433  // With 2 channels : VH VV
434  case VH_VV:
435  {
436  VHInputImagePointer inputPtrVH
437  = dynamic_cast<const TInputImageVH*>((itk::ProcessObject::GetInput(2)));
438  VVInputImagePointer inputPtrVV
439  = dynamic_cast<const TInputImageVV*>((itk::ProcessObject::GetInput(3)));
440 
441  itk::ImageRegionConstIterator<TInputImageVH> inputItVH(inputPtrVH, outputRegionForThread);
442  itk::ImageRegionConstIterator<TInputImageVV> inputItVV(inputPtrVV, outputRegionForThread);
443  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
444 
445  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
446 
447  inputItVH.GoToBegin();
448  inputItVV.GoToBegin();
449 
450  while (!inputItVH.IsAtEnd())
451  {
452  outputIt.Set(m_Gain * Superclass::GetFunctor() (0, 0, inputItVH.Get(), inputItVV.Get()));
453  ++inputItVH;
454  ++inputItVV;
455  ++outputIt;
456  progress.CompletedPixel(); // potential exception thrown here
457  }
458  break;
459  }
460  default:
461  itkExceptionMacro("Unknown architecture : Polarimetric synthesis is impossible !");
462  return;
463 
464  }
465 } // end namespace otb
466 
467 }
468 
469 #endif
Generated at Sun Jun 16 2013 00:44:53 for Orfeo Toolbox with doxygen 1.8.3.1