OTB  5.4.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  typename OutputImageType::Pointer output = this->GetOutput();
107  typename InputImageType::ConstPointer input = this->GetInput();
108 
109  // Find the data-set boundary "faces"
112 
114  faceList = bC(input, outputRegionForThread, m_Radius);
115 
116  // support progress methods/callbacks
117  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
118 
119  // InputRealType pixel;
120  InputRealType sum;
121  InputRealType sum2;
122 
123  double Ci2, Cu2, w, E_I, I, Var_I, dPixel;
124 
125  //Compute the ratio using the number of looks
126  Cu2 = 1.0/m_NbLooks;
127 
128  // Process each of the boundary faces. These are N-d regions which border
129  // the edge of the buffer.
130  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
131  {
132  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
133  const unsigned int neighborhoodSize = bit.Size();
135  bit.OverrideBoundaryCondition(&nbc);
136 
137  bit.GoToBegin();
138  it.GoToBegin();
139 
140 
141  while (!bit.IsAtEnd())
142  {
145 
146  //Parcours du voisinage
147  for (i = 0; i < neighborhoodSize; ++i)
148  {
149  dPixel = static_cast<double>(bit.GetPixel(i));
150  sum += dPixel;
151  }
152  E_I = sum / static_cast<double>(neighborhoodSize);
153 
154  for (i = 0; i < neighborhoodSize; ++i)
155  {
156  dPixel = static_cast<double>(bit.GetPixel(i));
157  sum2 += (dPixel-E_I) * (dPixel-E_I);
158  }
159  Var_I = sum2 / static_cast<double>(neighborhoodSize -1);
160 
161  I = static_cast<double>(bit.GetCenterPixel());
162 
163  Ci2 = Var_I / (E_I * E_I);
164 
165  const double epsilon = 0.0000000001;
166  if (vcl_abs(E_I) < epsilon)
167  {
169  }
170  else if (vcl_abs(Var_I) < epsilon)
171  {
172  dPixel = E_I;
173  }
174  else if (Ci2 < Cu2)
175  {
176  dPixel = E_I;
177  }
178  else
179  {
180  w = 1 - Cu2 / Ci2;
181  dPixel = I*w + E_I*(1-w);
182  }
183 
184  // set the weighted value
185  it.Set(static_cast<OutputPixelType>(dPixel));
186 
187  ++bit;
188  ++it;
189 
190  progress.CompletedPixel();
191  }
192  }
193 }
194 
198 template <class TInputImage, class TOutput>
199 void
201 {
202  Superclass::PrintSelf(os, indent);
203  os << indent << "Radius: " << m_Radius << std::endl;
204 }
206 
207 } // end namespace otb
208 
209 #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
unsigned int ThreadIdType
OutputImageType::RegionType OutputImageRegionType
void PrintSelf(std::ostream &os, itk::Indent indent) const