OTB  5.7.0
Orfeo Toolbox
otbKullbackLeiblerDistanceImageFilter.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  Copyright (c) Institut Mines-Telecom. All rights reserved.
13  See IMTCopyright.txt for details.
14 
15  This software is distributed WITHOUT ANY WARRANTY; without even
16  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
17  PURPOSE. See the above copyright notices for more information.
18 
19 =========================================================================*/
20 #ifndef otbKullbackLeiblerDistanceImageFilter_txx
21 #define otbKullbackLeiblerDistanceImageFilter_txx
22 
24 #include <vector>
25 
27 
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 /* *******************************************************************
34  * CumulantsForEdgeworth
35  *********************************************************************
36  */
37 template <class TInput>
39 ::CumulantsForEdgeworth (const TInput& input)
40 {
41  MakeSumAndMoments(input);
42  MakeCumulants();
43 }
44 
45 template <class TInput>
48 {
49  MakeSumAndMoments(input);
50  MakeCumulants();
51 }
52 
53 /* ========================== Divergence de KL ======================= */
54 
55 template <class TInput>
56 template <class TInput2>
57 double
60 {
61  double cum1 = GetMean();
62  double cum2 = GetVariance();
63  double cum3 = GetSkewness();
64 
65  double sqrt_cum2 = sqrt(cum2);
66  double cum2_3 = cum2 * cum2 * cum2;
67  double cum3_2 = cum3 * cum3;
68 
69  double tilde_cum1 = cumulants.GetMean();
70  double tilde_cum2 = cumulants.GetVariance();
71  double tilde_cum3 = cumulants.GetSkewness();
72  double tilde_cum4 = cumulants.GetKurtosis();
73 
74  double tilde_cum2_2 = cum2 * cum2;
75  double tilde_cum2_3 = cum2 * tilde_cum2_2;
76  double tilde_cum2_6 = tilde_cum2_3 * tilde_cum2_3;
77  double tilde_cum3_2 = tilde_cum3 * tilde_cum3;
78 
79  double beta = sqrt_cum2 / tilde_cum2;
80  double alpha = (cum1 - tilde_cum1) / tilde_cum2;
81 
82  double alpha_2 = alpha * alpha;
83  double alpha_4 = alpha_2 * alpha_2;
84  double alpha_6 = alpha_2 * alpha_4;
85 
86  double beta_2 = beta * beta;
87  double beta_4 = beta_2 * beta_2;
88  double beta_6 = beta_2 * beta_4;
89 
90  double c2 = alpha_2 + beta_2;
91  double c3 = alpha * (alpha_2 + 3.0 * beta_2);
92  double c4 = alpha_4 + 6.0 * alpha_2 * beta_2 + 3.0 * beta_4;
93  double c6 = alpha_6 + 15.0 * alpha_4 * beta_2 + 45.0 * alpha_2 * beta_4 + 15.0 * beta_6;
94 
95  double a1 = c3 - 3.0 * alpha / tilde_cum2;
96  double a2 = c4 - 6.0 * c2 / tilde_cum2 + 3.0 / tilde_cum2_2; // Watch for tilde_cum2_2 mistake in the article!
97  double a3 = c6 - 15.0 * c4 / tilde_cum2 + 45.0 * c2 / tilde_cum2_2 - 15.0 / tilde_cum2_3;
98 
99  double tmp = cum1 - tilde_cum1 + sqrt_cum2;
100  double resu = cum3_2 / (12.0 * cum2_3)
101  + (log(tilde_cum2 / cum2)
102  - 1.0 + tmp * tmp / tilde_cum2) / 2.0
103  - (tilde_cum3 * a1 / 6.0 + tilde_cum4 * a2 / 24.0 + tilde_cum3_2 * a3 / 72.0)
104  - tilde_cum3_2 * (c6 - 6.0 * c4 / cum2 + 9.0 * c2 / tilde_cum2_2) / 72.0
105  - 10.0 * cum3 * tilde_cum3 * (cum1 - tilde_cum1) * (cum2 - tilde_cum2) / tilde_cum2_6;
106 
107  return resu < 0.0 ? 0.0 : resu;
108 }
109 
110 /* ========== Moment estimation from initial neighborhood ============ */
111 
112 template <class TInput>
113 void
115 ::MakeSumAndMoments(const TInput& input)
116 {
117 
118  fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
119  double pixel, pixel_2;
120 
121  for (unsigned long i = 0; i < input.Size(); ++i)
122  {
123  pixel = static_cast<double> (input.GetPixel(i));
124  pixel_2 = pixel * pixel;
125 
126  fSum0 += 1.0;
127  fSum1 += pixel;
128  fSum2 += pixel_2;
129  fSum3 += pixel_2 * pixel;
130  fSum4 += pixel_2 * pixel_2;
131  }
132 
133  fMu1 = fSum1 / fSum0;
134  fMu2 = fSum2 / fSum0 - fMu1 * fMu1;
135 
136  if (fMu2 <= 0.0)
137  {
138  //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
139  fMu3 = 0.0;
140  fMu4 = 4.0;
141  fDataAvailable = false;
142  return;
143  }
144 
145  double sigma = sqrt(fMu2);
146 
147  itk::VariableLengthVector<double> tab(input.Size());
148  double * x = const_cast<double *> (tab.GetDataPointer());
149  for (unsigned long i = 0; i < input.Size(); ++i)
150  *x++ = (static_cast<double> (input.GetPixel(i)) - fMu1) / sigma;
151 
152  fMu3 = fMu4 = 0.0;
153  x = const_cast<double *> (tab.GetDataPointer());
154  for (unsigned long i = 0; i < input.Size(); ++i)
155  {
156  pixel = *x++;
157  pixel_2 = pixel * pixel;
158 
159  fMu3 += pixel * pixel_2;
160  fMu4 += pixel_2 * pixel_2;
161  }
162 
163  fMu3 /= fSum0;
164  fMu4 /= fSum0;
165 
166  fDataAvailable = true;
167 }
168 
169 /* ================= Moment estimation from raw data ================= */
170 
171 template <class TInput>
172 void
175 {
176 
177  fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
178  double pixel, pixel_2;
179 
181  typedef itk::ImageRegionConstIterator<LocalImageType> ImageRegionConstIteratorType;
182 
183  ImageRegionConstIteratorType inputIter(input, input->GetRequestedRegion());
184  inputIter.GoToBegin();
185  unsigned int inputSize = 0;
186 
187  while (!inputIter.IsAtEnd())
188  {
189  pixel = static_cast<double> (inputIter.Get());
190  pixel_2 = pixel * pixel;
191 
192  fSum0 += 1.0;
193  fSum1 += pixel;
194  fSum2 += pixel_2;
195  fSum3 += pixel_2 * pixel;
196  fSum4 += pixel_2 * pixel_2;
197 
198  ++inputIter;
199  ++inputSize;
200  }
201 
202  fMu1 = fSum1 / fSum0;
203  fMu2 = fSum2 / fSum0 - fMu1 * fMu1;
204 
205  if (fMu2 <= 0.0)
206  {
207  //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
208  fMu3 = 0.0;
209  fMu4 = 4.0;
210 
211  fDataAvailable = false;
212 
213  return;
214  }
215 
216  double sigma = sqrt(fMu2);
217 
218  std::vector<double> tab(inputSize);
219  std::vector<double>::iterator iterTab = tab.begin();
220 
221  inputIter.GoToBegin();
222 
223  while (!inputIter.IsAtEnd())
224  {
225  *iterTab = (static_cast<double> (inputIter.Get()) - fMu1) / sigma;
226 
227  ++inputIter;
228  ++iterTab;
229  }
230 
231  fMu3 = fMu4 = 0.0;
232  for (unsigned int i = 0; i < inputSize; ++i)
233  {
234  pixel = tab[i];
235  pixel_2 = pixel * pixel;
236 
237  fMu3 += pixel * pixel_2;
238  fMu4 += pixel_2 * pixel_2;
239  }
240 
241  fMu3 /= fSum0;
242  fMu4 /= fSum0;
243 
244  otbGenericMsgDebugMacro(<< "Moments: " << fMu1 << ",\t" << fMu2 << ",\t" << fMu3 << ",\t" << fMu4);
245 
246  fDataAvailable = true;
247  //return 0;
248 }
249 
250 /* ================= moments -> cumulants transformation ============= */
251 
252 template <class TInput>
253 void
256 {
257  if (fDataAvailable)
258  {
259  fMean = fMu1;
260  fVariance = fMu2;
261  fSkewness = fMu3;
262  fKurtosis = fMu4 - 3.0;
263  }
264  //return 0;
265 }
266 
267 } // end of namespace otb
268 
269 #endif
#define otbGenericMsgDebugMacro(x)
Definition: otbMacro.h:67
virtual const RegionType & GetRequestedRegion() const
Helper class for KullbackLeiblerDistanceImageFilter. Please refer to KullbackLeiblerDistanceImageFilt...
double Divergence(const CumulantsForEdgeworth< TInput2 > &cumulants)