OTB  6.1.0
Orfeo Toolbox
otbLeeImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbLeeImageFilter_txx
22 #define otbLeeImageFilter_txx
23 
24 #include "otbLeeImageFilter.h"
25 
26 #include "itkDataObject.h"
29 #include "itkImageRegionIterator.h"
31 #include "itkOffset.h"
32 #include "itkProgressReporter.h"
33 
34 namespace otb
35 {
36 
40 template <class TInputImage, class TOutputImage>
42 {
43  m_Radius.Fill(1);
44  SetNbLooks(1.0);
45 }
47 
48 template <class TInputImage, class TOutputImage>
49 void LeeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
50  {
51  // call the superclass' implementation of this method
52  Superclass::GenerateInputRequestedRegion();
53 
54  // get pointers to the input and output
55  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
56  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
57 
58  if (!inputPtr || !outputPtr)
59  {
60  return;
61  }
62 
63  // get a copy of the input requested region (should equal the output
64  // requested region)
65  typename TInputImage::RegionType inputRequestedRegion;
66  inputRequestedRegion = inputPtr->GetRequestedRegion();
67 
68  // pad the input requested region by the operator radius
69  inputRequestedRegion.PadByRadius(m_Radius);
70 
71  // crop the input requested region at the input's largest possible region
72  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
73  {
74  inputPtr->SetRequestedRegion(inputRequestedRegion);
75  return;
76  }
77  else
78  {
79  // Couldn't crop the region (requested region is outside the largest
80  // possible region). Throw an exception.
81 
82  // store what we tried to request (prior to trying to crop)
83  inputPtr->SetRequestedRegion(inputRequestedRegion);
84 
85  // build an exception
86  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
87  std::ostringstream msg;
88  msg << static_cast<const char *>(this->GetNameOfClass())
89  << "::GenerateInputRequestedRegion()";
90  e.SetLocation(msg.str().c_str());
91  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
92  e.SetDataObject(inputPtr);
93  throw e;
94  }
95  }
96 
97 template<class TInputImage, class TOutputImage>
99  const OutputImageRegionType& outputRegionForThread,
100  itk::ThreadIdType threadId
101  )
102 {
103  unsigned int i;
105 
108 
109  typename OutputImageType::Pointer output = this->GetOutput();
110  typename InputImageType::ConstPointer input = this->GetInput();
111 
112  // Find the data-set boundary "faces"
115 
117  faceList = bC(input, outputRegionForThread, m_Radius);
118 
119  // support progress methods/callbacks
120  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
121 
122  // InputRealType pixel;
123  InputRealType sum;
124  InputRealType sum2;
125 
126  double Ci2, Cu2, w, E_I, I, Var_I, dPixel;
127 
128  //Compute the ratio using the number of looks
129  Cu2 = 1.0/m_NbLooks;
130 
131  // Process each of the boundary faces. These are N-d regions which border
132  // the edge of the buffer.
133  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
134  {
135  bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
136  const unsigned int neighborhoodSize = bit.Size();
138  bit.OverrideBoundaryCondition(&nbc);
139 
140  bit.GoToBegin();
141  it.GoToBegin();
142 
143 
144  while (!bit.IsAtEnd())
145  {
148 
149  //Parcours du voisinage
150  for (i = 0; i < neighborhoodSize; ++i)
151  {
152  dPixel = static_cast<double>(bit.GetPixel(i));
153  sum += dPixel;
154  }
155  E_I = sum / static_cast<double>(neighborhoodSize);
156 
157  for (i = 0; i < neighborhoodSize; ++i)
158  {
159  dPixel = static_cast<double>(bit.GetPixel(i));
160  sum2 += (dPixel-E_I) * (dPixel-E_I);
161  }
162  Var_I = sum2 / static_cast<double>(neighborhoodSize -1);
163 
164  I = static_cast<double>(bit.GetCenterPixel());
165 
166  Ci2 = Var_I / (E_I * E_I);
167 
168  const double epsilon = 0.0000000001;
169  if (vcl_abs(E_I) < epsilon)
170  {
172  }
173  else if (vcl_abs(Var_I) < epsilon)
174  {
175  dPixel = E_I;
176  }
177  else if (Ci2 < Cu2)
178  {
179  dPixel = E_I;
180  }
181  else
182  {
183  w = 1 - Cu2 / Ci2;
184  dPixel = I*w + E_I*(1-w);
185  }
186 
187  // set the weighted value
188  it.Set(static_cast<OutputPixelType>(dPixel));
189 
190  ++bit;
191  ++it;
192 
193  progress.CompletedPixel();
194  }
195  }
196 }
197 
201 template <class TInputImage, class TOutput>
202 void
204 {
205  Superclass::PrintSelf(os, indent);
206  os << indent << "Radius: " << m_Radius << std::endl;
207 }
209 
210 } // end namespace otb
211 
212 #endif
virtual bool IsAtEnd() const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
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
PixelType GetCenterPixel() const
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
unsigned int ThreadIdType
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
void GenerateInputRequestedRegion() ITK_OVERRIDE
OutputImageType::RegionType OutputImageRegionType