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