OTB  5.0.0
Orfeo Toolbox
otbLeeImageFilter.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 __otbLeeImageFilter_txx
19 #define __otbLeeImageFilter_txx
20 
21 #include "otbLeeImageFilter.h"
22 
23 #include "itkDataObject.h"
26 #include "itkImageRegionIterator.h"
28 #include "itkOffset.h"
29 #include "itkProgressReporter.h"
30 
31 namespace otb
32 {
33 
37 template <class TInputImage, class TOutputImage>
39 {
40  m_Radius.Fill(1);
41  SetNbLooks(1.0);
42 }
44 
45 template <class TInputImage, class TOutputImage>
46 void LeeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
47  {
48  // call the superclass' implementation of this method
49  Superclass::GenerateInputRequestedRegion();
50 
51  // get pointers to the input and output
52  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
53  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
54 
55  if (!inputPtr || !outputPtr)
56  {
57  return;
58  }
59 
60  // get a copy of the input requested region (should equal the output
61  // requested region)
62  typename TInputImage::RegionType inputRequestedRegion;
63  inputRequestedRegion = inputPtr->GetRequestedRegion();
64 
65  // pad the input requested region by the operator radius
66  inputRequestedRegion.PadByRadius(m_Radius);
67 
68  // crop the input requested region at the input's largest possible region
69  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
70  {
71  inputPtr->SetRequestedRegion(inputRequestedRegion);
72  return;
73  }
74  else
75  {
76  // Couldn't crop the region (requested region is outside the largest
77  // possible region). Throw an exception.
78 
79  // store what we tried to request (prior to trying to crop)
80  inputPtr->SetRequestedRegion(inputRequestedRegion);
81 
82  // build an exception
83  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
84  std::ostringstream msg;
85  msg << static_cast<const char *>(this->GetNameOfClass())
86  << "::GenerateInputRequestedRegion()";
87  e.SetLocation(msg.str().c_str());
88  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
89  e.SetDataObject(inputPtr);
90  throw e;
91  }
92  }
93 
94 template<class TInputImage, class TOutputImage>
96  const OutputImageRegionType& outputRegionForThread,
97  itk::ThreadIdType threadId
98  )
99 {
100  unsigned int i;
102 
105 
106  // Allocate output
107  typename OutputImageType::Pointer output = this->GetOutput();
108  typename InputImageType::ConstPointer input = this->GetInput();
109 
110  // Find the data-set boundary "faces"
113 
115  faceList = bC(input, outputRegionForThread, m_Radius);
116 
117  // support progress methods/callbacks
118  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
119 
120  // InputRealType pixel;
121  InputRealType sum;
122  InputRealType sum2;
123 
124  double Cr2, Cv2, E_I, I, Var_I, dPixel;
125 
126  //Compute the ratio using the number of looks
127  Cv2 = 1. / (vcl_sqrt(m_NbLooks));
128  Cv2 *= Cv2;
129 
130  // Process each of the boundary faces. These are N-d regions which border
131  // the edge of the buffer.
132  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
133  {
134  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
135  const unsigned int neighborhoodSize = bit.Size();
137  bit.OverrideBoundaryCondition(&nbc);
138  bit.GoToBegin();
139 
140  while (!bit.IsAtEnd())
141  {
144  //Parcours du voisinage
145  for (i = 0; i < neighborhoodSize; ++i)
146  {
147  dPixel = static_cast<double>(bit.GetPixel(i));
148  sum += dPixel;
149  sum2 += dPixel * dPixel;
150  }
151  E_I = sum / static_cast<double>(neighborhoodSize);
152  Var_I = sum2 / static_cast<double>(neighborhoodSize) - E_I * E_I;
153  I = static_cast<double>(bit.GetCenterPixel());
154 
155  const double epsilon = 0.0000000001;
156  if (vcl_abs(E_I) < epsilon)
157  {
159  }
160  else
161  {
162  Cr2 = Var_I / (E_I * E_I);
163  dPixel = E_I + ((I - E_I) * (Cr2)) / (Cr2 + Cv2);
164 
165  }
166  // get the mean value
167  it.Set(static_cast<OutputPixelType>(dPixel));
168 
169  ++bit;
170  ++it;
171  progress.CompletedPixel();
172  }
173  }
174 }
175 
179 template <class TInputImage, class TOutput>
180 void
182 {
183  Superclass::PrintSelf(os, indent);
184  os << indent << "Radius: " << m_Radius << std::endl;
185 }
187 
188 } // end namespace otb
189 
190 #endif
virtual bool IsAtEnd() const
virtual void GenerateInputRequestedRegion()
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
virtual void OverrideBoundaryCondition(const ImageBoundaryConditionPointerType i)
itk::NumericTraits< InputPixelType >::RealType InputRealType
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
PixelType GetCenterPixel() const
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
OutputImageType::RegionType OutputImageRegionType
void PrintSelf(std::ostream &os, itk::Indent indent) const
unsigned int ThreadIdType