OTB  6.7.0
Orfeo Toolbox
otbKullbackLeiblerProfileImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 otbKullbackLeiblerProfileImageFilter_hxx
23 #define otbKullbackLeiblerProfileImageFilter_hxx
24 
25 #include <vector>
26 
28 #include "otbMath.h"
29 
30 namespace otb
31 {
32 
33 /* *******************************************************************
34 * Classe CumulantsForEdgeworthProfile
35 * ********************************************************************
36 */
37 template <class TInput>
40  (const TInput& input, std::vector<itk::Array2D<int> >& mask)
41 {
42  m_debug = MakeSumAndMoments(input, mask);
43  MakeCumulants();
44 }
45 
46 /* ===================== Kullback-Leibler Profile ==================== */
47 
48 template <class TInput>
49 template <class TInput2>
53 {
54  itk::VariableLengthVector<double> resu(fCum.size());
55 
56  Iterator iter1 = fCum.begin();
57  Iterator iter2 = cumulants.fCum.begin();
58 
59  for (unsigned int level = 0; level < resu.GetSize(); level++)
60  resu[level] = KL_profile((*iter1++), (*iter2++));
61 
62  return resu;
63 }
64 
65 /* =========== Kullback-Leibler divergence at a given scale ========== */
66 
67 template <class TInput>
68 double
71  (const CumulantType& cumulants1, const CumulantType& cumulants2)
72 {
73  double cum1 = cumulants1[0];
74  double cum2 = cumulants1[1];
75  double cum3 = cumulants1[2];
76  // double cum4 = cumulants1[3]; // unused
77 
78  double sqrt_cum2 = sqrt(cum2);
79  double cum2_3 = cum2 * cum2 * cum2;
80  double cum3_2 = cum3 * cum3;
81 
82  double tilde_cum1 = cumulants2[0];
83  double tilde_cum2 = cumulants2[1];
84  double tilde_cum3 = cumulants2[2];
85  double tilde_cum4 = cumulants2[3];
86 
87  double tilde_cum2_2 = cum2 * cum2;
88  double tilde_cum2_3 = cum2 * tilde_cum2_2;
89  double tilde_cum2_6 = tilde_cum2_3 * tilde_cum2_3;
90  double tilde_cum3_2 = tilde_cum3 * tilde_cum3;
91 
92  double beta = sqrt_cum2 / tilde_cum2;
93  double alpha = (cum1 - tilde_cum1) / tilde_cum2;
94 
95  double alpha_2 = alpha * alpha;
96  double alpha_4 = alpha_2 * alpha_2;
97  double alpha_6 = alpha_2 * alpha_4;
98 
99  double beta_2 = beta * beta;
100  double beta_4 = beta_2 * beta_2;
101  double beta_6 = beta_2 * beta_4;
102 
103  double c2 = alpha_2 + beta_2;
104  double c3 = alpha * (alpha_2 + 3.0 * beta_2);
105  double c4 = alpha_4 + 6.0 * alpha_2 * beta_2 + 3.0 * beta_4;
106  double c6 = alpha_6 + 15.0 * alpha_4 * beta_2 + 45.0 * alpha_2 * beta_4 + 15.0 * beta_6;
107 
108  double a1 = c3 - 3.0 * alpha / tilde_cum2;
109  double a2 = c4 - 6.0 * c2 / tilde_cum2 + 3.0 / tilde_cum2_2; // Watch for tilde_cum2_2 mistake in the article!
110  double a3 = c6 - 15.0 * c4 / tilde_cum2 + 45.0 * c2 / tilde_cum2_2 - 15.0 / tilde_cum2_3;
111 
112  double tmp = cum1 - tilde_cum1 + sqrt_cum2;
113  double resu = cum3_2 / (12.0 * cum2_3)
114  + (log(tilde_cum2 / cum2)
115  - 1.0 + tmp * tmp / tilde_cum2) / 2.0
116  - (tilde_cum3 * a1 / 6.0 + tilde_cum4 * a2 / 24.0 + tilde_cum3_2 * a3 / 72.0)
117  - tilde_cum3_2 * (c6 - 6.0 * c4 / cum2 + 9.0 * c2 / tilde_cum2_2) / 72.0
118  - 10.0 * cum3 * tilde_cum3 * (cum1 - tilde_cum1) * (cum2 - tilde_cum2) / tilde_cum2_6;
119 
120  if (vnl_math_isnan(resu) || resu > 1e3) resu = 1e3;
121 
122  return resu < 0.0 ? 0.0 : resu;
123 }
124 
125 /* ====== Moments estimation from nested neighborhoods ==== */
126 
127 template <class TInput>
128 int
131  (const TInput& input, std::vector<itk::Array2D<int> >& mask)
132 {
133  fMu.resize(mask.size());
134  std::vector<itk::Array2D<int> >::iterator iter = mask.begin();
135 
136  if (InitSumAndMoments (input, (*iter++))) return 1;
137 
138  for (unsigned int level = 1; level < mask.size(); level++)
139  if (ReInitSumAndMoments(input, (*iter++), level)) return 1;
140 
141  return 0;
142 }
143 
144 /* ===================== Moments estimation ====================== */
145 /* =============== from the small window size =========== */
146 
147 template <class TInput>
148 int
151  (const TInput& input, itk::Array2D<int>& mask)
152 {
153  fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
154  fMu[0].Fill(0.0);
155 
156  unsigned int i, j;
157  unsigned long k = 0;
158  double pixel, pixel_2;
159 
160  // for ( unsigned long i = 0; i < input.Size(); ++i )
161  for (i = 0; i < mask.rows(); ++i)
162  {
163  for (j = 0; j < mask.cols(); ++j)
164  {
165  // std::cerr << "(" << i << "," << j << ") k=" << k << " ";
166  if (mask.get(i, j) == 1)
167  {
168  pixel = static_cast<double> (input.GetPixel(k));
169  pixel_2 = pixel * pixel;
170 
171  fSum0 += 1.0;
172  fSum1 += pixel;
173  fSum2 += pixel_2;
174  fSum3 += pixel_2 * pixel;
175  fSum4 += pixel_2 * pixel_2;
176  // std::cerr << "*\n";
177  }
178  ++k;
179  }
180  }
181  if (fSum0 == 0.0)
182  {
183  fDataAvailable = false;
184  return 1;
185  }
186 
187  double mu1, mu2;
188 
189  mu1 = fSum1 / fSum0;
190  mu2 = fSum2 / fSum0 - mu1 * mu1;
191 
192  if (mu2 == 0.0)
193  {
194  fDataAvailable = false;
195  return 1;
196  }
197 
198  double sigma = sqrt(mu2);
199 
200  itk::VariableLengthVector<double> tab(input.Size());
201  double * x = const_cast<double *> (tab.GetDataPointer());
202  for (unsigned long cp = 0; cp < input.Size(); ++cp)
203  *x++ = (static_cast<double> (input.GetPixel(cp)) - mu1) / sigma;
204 
205  double mu3 = 0.0;
206  double mu4 = 0.0;
207  x = const_cast<double *> (tab.GetDataPointer());
208 
209  // for ( unsigned long i = 0; i < input.Size(); ++i )
210  for (i = 0; i < mask.rows(); ++i)
211  {
212  for (j = 0; j < mask.cols(); ++j)
213  {
214  if (mask.get(i, j) == 1)
215  {
216  pixel = *x++;
217  pixel_2 = pixel * pixel;
218 
219  mu3 += pixel * pixel_2;
220  mu4 += pixel_2 * pixel_2;
221  }
222  else ++x;
223  }
224  }
225 
226  mu3 /= fSum0;
227  mu4 /= fSum0;
228 
229  if (vnl_math_isnan(mu3) || vnl_math_isnan(mu4))
230  {
231  fDataAvailable = false;
232  return 1;
233  }
234 
235  fMu[0][0] = mu1;
236  fMu[0][1] = mu2;
237  fMu[0][2] = mu3;
238  fMu[0][3] = mu4;
239 
240  fDataAvailable = true;
241 
242  return 0;
243 }
244 
245 /* ================ Window size growth ============ */
246 
247 template <class TInput>
248 int
251  (const TInput& input, itk::Array2D<int>& mask, int level)
252 {
253  fMu[level].Fill(0.0);
254  // mise a jour du comptage...
255  double sum0 = 0.0,
256  sum1 = 0.0,
257  sum2 = 0.0,
258  sum3 = 0.0,
259  sum4 = 0.0;
260 
261  double pixel, pixel_2;
262 
263  unsigned int i, j;
264  unsigned long k = 0L;
265 
266  // for ( unsigned long i = 0; i < input.Size(); ++i )
267  for (i = 0; i < mask.rows(); ++i)
268  {
269  for (j = 0; j < mask.cols(); ++j)
270  {
271  if (mask.get(i, j) == 1)
272  {
273  pixel = static_cast<double> (input.GetPixel(k));
274  pixel_2 = pixel * pixel;
275 
276  sum0 += 1.0;
277  sum1 += pixel;
278  sum2 += pixel_2;
279  sum3 += pixel * pixel_2;
280  sum4 += pixel_2 * pixel_2;
281  }
282  ++k;
283  }
284  }
285 
286  fSum0 += sum0;
287  fSum1 += sum1;
288  fSum2 += sum2;
289  fSum3 += sum3;
290  fSum4 += sum4;
291 
292  double mu = fSum1 / fSum0;
293  double mu_2 = mu * mu;
294  double mu_3 = mu_2 * mu;
295  double mu_4 = mu_2 * mu_2;
296 
297  fMu[level][0] = mu;
298  fMu[level][1] = fSum2 / fSum0 - mu_2;
299 
300  double sigma = sqrt(fSum2);
301  double sigma_2 = fSum2;
302  double sigma_3 = sigma * sigma_2;
303  double sigma_4 = sigma_2 * sigma_2;
304 
305  fMu[level][2] = (fSum3 - 3.0 * mu * fSum2 + 3.0 * mu_2 * fSum1 - fSum0 * mu_3)
306  / (sigma_3 * fSum0);
307  fMu[level][3] = (fSum4 - 4.0 * mu * fSum3 + 6.0 * mu_2 * fSum2
308  - 4.0 * mu_3 * fSum1 + fSum0 * mu_4)
309  / (sigma_4 * fSum0);
310 
311  return 0;
312 }
313 
314 /* =========== transformation moment -> cumulants ==================== */
315 
316 template <class TInput>
317 int
320 {
321  if (!IsDataAvailable()) return 1;
322 
323  fCum.resize(fMu.size());
324  fCum = fMu;
325 
326  for (unsigned int i = 0; i < fCum.size(); ++i)
327  fCum[i][3] -= 3.0;
328 
329  return 0;
330 }
331 
332 /* *******************************************************************
333 *
334 * Functor
335 *
336 * ********************************************************************
337 */
338 
339 namespace Functor
340 {
341 
342 template<class TInput1, class TInput2, class TOutput>
345 {
346  m_RadiusMin = 1;
347  m_RadiusMax = 3;
348 }
349 
350 /* =========== Gives the radius min and max of neighborhood ========== */
351 
352 template<class TInput1, class TInput2, class TOutput>
353 void
355 ::SetRadius(const unsigned char& min, const unsigned char& max)
356 {
357  m_RadiusMin = min < max ? min : max;
358  m_RadiusMax = max > min ? max : min;
359  MakeMultiscaleProfile();
360 }
361 
362 template<class TInput1, class TInput2, class TOutput>
363 unsigned char
366 {
367  return m_RadiusMin;
368 }
369 
370 template<class TInput1, class TInput2, class TOutput>
371 unsigned char
374 {
375  return m_RadiusMax;
376 }
377 
378 /* ====== Make the set of masks to play the increase in window size == */
379 
380 template<class TInput1, class TInput2, class TOutput>
381 void
384 {
385  m_mask.resize(m_RadiusMax - m_RadiusMin + 1);
386  int lenMax = 2 * m_RadiusMax + 1;
387  int i, j, middle = m_RadiusMax;
388 
389  // let's begin by the smaller neighborhood
390  std::vector<itk::Array2D<int> >::iterator outer_iter = m_mask.begin();
391  (*outer_iter).SetSize(lenMax, lenMax);
392  (*outer_iter).fill(0);
393  for (i = middle - m_RadiusMin; i <= middle + m_RadiusMin; ++i)
394  for (j = middle - m_RadiusMin; j <= middle + m_RadiusMin; ++j)
395  (*outer_iter).put(i, j, 1);
396 
397  // std::cerr << "outerIter = " << (*outer_iter) << std::endl;
398 
399  // let's continue with increasing neighborhoods
400  ++outer_iter;
401  for (int radius = m_RadiusMin + 1; radius <= m_RadiusMax; ++radius)
402  {
403  (*outer_iter).SetSize(lenMax, lenMax);
404  (*outer_iter).fill(0);
405 
406  for (i = middle - radius; i <= middle + radius; ++i)
407  {
408  (*outer_iter).put(i, middle - radius, 1);
409  (*outer_iter).put(i, middle + radius, 1);
410  (*outer_iter).put(middle - radius, i, 1);
411  (*outer_iter).put(middle + radius, i, 1);
412  }
413 
414  // std::cerr << "outerIter = " << (*outer_iter) << std::endl;
415  ++outer_iter;
416  }
417 }
418 
419 /* ========================== Functor ================================ */
420 
421 template<class TInput1, class TInput2, class TOutput>
422 TOutput
424 ::operator()
425  (const TInput1 &it1, const TInput2 &it2)
426  {
427  CumulantsForEdgeworthProfile<TInput1> cum1(it1, m_mask);
428 
429  if (cum1.m_debug)
430  {
431  itk::VariableLengthVector<double> resu(m_RadiusMax - m_RadiusMin + 1);
432  resu.Fill(1e3);
433  return static_cast<TOutput> (resu);
434  }
435 
436  CumulantsForEdgeworthProfile<TInput2> cum2(it2, m_mask);
437 
438  if (cum2.m_debug)
439  {
440  itk::VariableLengthVector<double> resu(m_RadiusMax - m_RadiusMin + 1);
441  resu.Fill(1e3);
442  return static_cast<TOutput> (resu);
443  }
444 
445  return static_cast<TOutput> (cum1.KL_profile(cum2) + cum2.KL_profile(cum1));
446  }
447 
448 } // Functor
449 
450 } // namespace otb
451 
452 #endif
int InitSumAndMoments(const TInput &input, itk::Array2D< int > &mask)
int MakeSumAndMoments(const TInput &input, std::vector< itk::Array2D< int > > &mask)
void SetRadius(const unsigned char &min, const unsigned char &max)
void Fill(TValue const &v) noexcept
itk::VariableLengthVector< double > KL_profile(CumulantsForEdgeworthProfile< TInput2 > &cumulants)
Functor for KullbackLeiblerProfileImageFilter. Please refer to KullbackLeiblerProfileImageFilter.
Helper class for KullbackLeiblerProfileImageFilter. Please refer to KullbackLeibleProfileImageFilter...
int ReInitSumAndMoments(const TInput &input, itk::Array2D< int > &mask, int level)