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