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