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

Generated at Sat Oct 11 2014 16:15:47 for Orfeo Toolbox with doxygen 1.8.3.1