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