OTB  5.11.0
Orfeo Toolbox
otbKullbackLeiblerDistanceImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  * Copyright (C) 2007-2012 Institut Mines Telecom / Telecom Bretagne
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbKullbackLeiblerDistanceImageFilter_txx
23 #define otbKullbackLeiblerDistanceImageFilter_txx
24 
26 #include <vector>
27 
29 
30 #include "otbMacro.h"
31 
32 namespace otb
33 {
34 
35 /* *******************************************************************
36  * CumulantsForEdgeworth
37  *********************************************************************
38  */
39 template <class TInput>
41 ::CumulantsForEdgeworth (const TInput& input)
42 {
43  MakeSumAndMoments(input);
44  MakeCumulants();
45 }
46 
47 template <class TInput>
50 {
51  MakeSumAndMoments(input);
52  MakeCumulants();
53 }
54 
55 /* ========================== Divergence de KL ======================= */
56 
57 template <class TInput>
58 template <class TInput2>
59 double
62 {
63  double cum1 = GetMean();
64  double cum2 = GetVariance();
65  double cum3 = GetSkewness();
66 
67  double sqrt_cum2 = sqrt(cum2);
68  double cum2_3 = cum2 * cum2 * cum2;
69  double cum3_2 = cum3 * cum3;
70 
71  double tilde_cum1 = cumulants.GetMean();
72  double tilde_cum2 = cumulants.GetVariance();
73  double tilde_cum3 = cumulants.GetSkewness();
74  double tilde_cum4 = cumulants.GetKurtosis();
75 
76  double tilde_cum2_2 = cum2 * cum2;
77  double tilde_cum2_3 = cum2 * tilde_cum2_2;
78  double tilde_cum2_6 = tilde_cum2_3 * tilde_cum2_3;
79  double tilde_cum3_2 = tilde_cum3 * tilde_cum3;
80 
81  double beta = sqrt_cum2 / tilde_cum2;
82  double alpha = (cum1 - tilde_cum1) / tilde_cum2;
83 
84  double alpha_2 = alpha * alpha;
85  double alpha_4 = alpha_2 * alpha_2;
86  double alpha_6 = alpha_2 * alpha_4;
87 
88  double beta_2 = beta * beta;
89  double beta_4 = beta_2 * beta_2;
90  double beta_6 = beta_2 * beta_4;
91 
92  double c2 = alpha_2 + beta_2;
93  double c3 = alpha * (alpha_2 + 3.0 * beta_2);
94  double c4 = alpha_4 + 6.0 * alpha_2 * beta_2 + 3.0 * beta_4;
95  double c6 = alpha_6 + 15.0 * alpha_4 * beta_2 + 45.0 * alpha_2 * beta_4 + 15.0 * beta_6;
96 
97  double a1 = c3 - 3.0 * alpha / tilde_cum2;
98  double a2 = c4 - 6.0 * c2 / tilde_cum2 + 3.0 / tilde_cum2_2; // Watch for tilde_cum2_2 mistake in the article!
99  double a3 = c6 - 15.0 * c4 / tilde_cum2 + 45.0 * c2 / tilde_cum2_2 - 15.0 / tilde_cum2_3;
100 
101  double tmp = cum1 - tilde_cum1 + sqrt_cum2;
102  double resu = cum3_2 / (12.0 * cum2_3)
103  + (log(tilde_cum2 / cum2)
104  - 1.0 + tmp * tmp / tilde_cum2) / 2.0
105  - (tilde_cum3 * a1 / 6.0 + tilde_cum4 * a2 / 24.0 + tilde_cum3_2 * a3 / 72.0)
106  - tilde_cum3_2 * (c6 - 6.0 * c4 / cum2 + 9.0 * c2 / tilde_cum2_2) / 72.0
107  - 10.0 * cum3 * tilde_cum3 * (cum1 - tilde_cum1) * (cum2 - tilde_cum2) / tilde_cum2_6;
108 
109  return resu < 0.0 ? 0.0 : resu;
110 }
111 
112 /* ========== Moment estimation from initial neighborhood ============ */
113 
114 template <class TInput>
115 void
117 ::MakeSumAndMoments(const TInput& input)
118 {
119 
120  fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
121  double pixel, pixel_2;
122 
123  for (unsigned long i = 0; i < input.Size(); ++i)
124  {
125  pixel = static_cast<double> (input.GetPixel(i));
126  pixel_2 = pixel * pixel;
127 
128  fSum0 += 1.0;
129  fSum1 += pixel;
130  fSum2 += pixel_2;
131  fSum3 += pixel_2 * pixel;
132  fSum4 += pixel_2 * pixel_2;
133  }
134 
135  fMu1 = fSum1 / fSum0;
136  fMu2 = fSum2 / fSum0 - fMu1 * fMu1;
137 
138  if (fMu2 <= 0.0)
139  {
140  //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
141  fMu3 = 0.0;
142  fMu4 = 4.0;
143  fDataAvailable = false;
144  return;
145  }
146 
147  double sigma = sqrt(fMu2);
148 
149  itk::VariableLengthVector<double> tab(input.Size());
150  double * x = const_cast<double *> (tab.GetDataPointer());
151  for (unsigned long i = 0; i < input.Size(); ++i)
152  *x++ = (static_cast<double> (input.GetPixel(i)) - fMu1) / sigma;
153 
154  fMu3 = fMu4 = 0.0;
155  x = const_cast<double *> (tab.GetDataPointer());
156  for (unsigned long i = 0; i < input.Size(); ++i)
157  {
158  pixel = *x++;
159  pixel_2 = pixel * pixel;
160 
161  fMu3 += pixel * pixel_2;
162  fMu4 += pixel_2 * pixel_2;
163  }
164 
165  fMu3 /= fSum0;
166  fMu4 /= fSum0;
167 
168  fDataAvailable = true;
169 }
170 
171 /* ================= Moment estimation from raw data ================= */
172 
173 template <class TInput>
174 void
177 {
178 
179  fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
180  double pixel, pixel_2;
181 
183  typedef itk::ImageRegionConstIterator<LocalImageType> ImageRegionConstIteratorType;
184 
185  ImageRegionConstIteratorType inputIter(input, input->GetRequestedRegion());
186  inputIter.GoToBegin();
187  unsigned int inputSize = 0;
188 
189  while (!inputIter.IsAtEnd())
190  {
191  pixel = static_cast<double> (inputIter.Get());
192  pixel_2 = pixel * pixel;
193 
194  fSum0 += 1.0;
195  fSum1 += pixel;
196  fSum2 += pixel_2;
197  fSum3 += pixel_2 * pixel;
198  fSum4 += pixel_2 * pixel_2;
199 
200  ++inputIter;
201  ++inputSize;
202  }
203 
204  fMu1 = fSum1 / fSum0;
205  fMu2 = fSum2 / fSum0 - fMu1 * fMu1;
206 
207  if (fMu2 <= 0.0)
208  {
209  //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
210  fMu3 = 0.0;
211  fMu4 = 4.0;
212 
213  fDataAvailable = false;
214 
215  return;
216  }
217 
218  double sigma = sqrt(fMu2);
219 
220  std::vector<double> tab(inputSize);
221  std::vector<double>::iterator iterTab = tab.begin();
222 
223  inputIter.GoToBegin();
224 
225  while (!inputIter.IsAtEnd())
226  {
227  *iterTab = (static_cast<double> (inputIter.Get()) - fMu1) / sigma;
228 
229  ++inputIter;
230  ++iterTab;
231  }
232 
233  fMu3 = fMu4 = 0.0;
234  for (unsigned int i = 0; i < inputSize; ++i)
235  {
236  pixel = tab[i];
237  pixel_2 = pixel * pixel;
238 
239  fMu3 += pixel * pixel_2;
240  fMu4 += pixel_2 * pixel_2;
241  }
242 
243  fMu3 /= fSum0;
244  fMu4 /= fSum0;
245 
246  otbGenericMsgDebugMacro(<< "Moments: " << fMu1 << ",\t" << fMu2 << ",\t" << fMu3 << ",\t" << fMu4);
247 
248  fDataAvailable = true;
249  //return 0;
250 }
251 
252 /* ================= moments -> cumulants transformation ============= */
253 
254 template <class TInput>
255 void
258 {
259  if (fDataAvailable)
260  {
261  fMean = fMu1;
262  fVariance = fMu2;
263  fSkewness = fMu3;
264  fKurtosis = fMu4 - 3.0;
265  }
266  //return 0;
267 }
268 
269 } // end of namespace otb
270 
271 #endif
#define otbGenericMsgDebugMacro(x)
Definition: otbMacro.h:70
virtual const RegionType & GetRequestedRegion() const
Helper class for KullbackLeiblerDistanceImageFilter. Please refer to KullbackLeiblerDistanceImageFilt...
double Divergence(const CumulantsForEdgeworth< TInput2 > &cumulants)