1 /***************************************************************************
2                           smooth1dnans.hpp -- nan version of smooth 1d
3                              -------------------
4     begin                : May 30 2017
5     copyright            : (C) 2017 by G. Duvert
6     email                : m_schellens@users.sourceforge.net
7  ***************************************************************************/
8 
9 /***************************************************************************
10  *                                                                         *
11  *   This program is free software; you can redistribute it and/or modify  *
12  *   it under the terms of the GNU General Public License as published by  *
13  *   the Free Software Foundation; either version 2 of the License, or     *
14  *   (at your option) any later version.                                   *
15  *                                                                         *
16  ***************************************************************************/
17 
18 // to be included from datatypes.cpp
19 #ifdef INCLUDE_SMOOTH_1D_NAN
20 //initiate mean of Width first values:
21 DDouble mean = 0;
22 DDouble z;
23 DDouble n = 0;
24 //initial mean. may be Nan
25 for (SizeT i = 0; i < (2 * w + 1); ++i) {
26  DDouble v = data[i];
27  if (isfinite(v)) {
28   n += 1.0;
29   z = 1. / n;
30   mean = (1. - z) * mean + z * v;
31  }
32 }
33 #if defined(USE_EDGE)
34 //start: use mean1, n1
35 DDouble mean1 = mean;
36 DDouble n1 = n;
37 for (SizeT i = 0; i < w; ++i) {
38  if (n1 > 0) res[w - i] = mean1;
39  DDouble v = data[2 * w - i];
40  if (isfinite(v)) {
41   mean1 *= n1;
42   mean1 -= v;
43   n1 -= 1.0;
44   mean1 /= n1;
45  }
46 
47  if (n1 <= 0) mean1 = 0;
48 
49 #if defined (EDGE_WRAP)
50  v = data[dimx - i - 1];
51  if (isfinite(v)) {
52 #elif defined (EDGE_TRUNCATE)
53  v = data[0];
54  if (isfinite(v)) {
55 #elif defined (EDGE_MIRROR)
56  v = data[i];
57  if (isfinite(v)) {
58 #elif defined (EDGE_ZERO)
59  v = 0;
60  {
61 #endif
62   mean1 *= n1;
63   if (n1 < 2 * w + 1) n1 += 1.0;
64   mean1 += v;
65   mean1 /= n1;
66  }
67 }
68 if (n1 > 0) res[0] = mean1;
69 
70 #endif
71 
72 //middle: use mean & n
73 for (SizeT i = w, ivm = 0, ivp = 2 * w + 1; i < dimx - w - 1; ++i, ++ivm, ++ivp) {
74  if (n > 0) res[i] = mean; //otherwise contains already data[i] since res is initialized to data.
75  DDouble v = data[ivm];
76  if (isfinite(v)) {
77   mean *= n;
78   mean -= v;
79   n -= 1.0;
80   mean /= n;
81  }
82 
83  if (n <= 0) mean = 0;
84 
85  v = data[ivp];
86  if (isfinite(v)) {
87   mean *= n;
88   if (n < 2 * w + 1) n += 1.0;
89   mean += v;
90   mean /= n;
91  }
92 }
93 if (n > 0) res[dimx - 1 - w] = mean;
94 
95 #ifdef USE_EDGE
96 //end: use mean2 & n as we go in the same direction...
97 DDouble mean2 = mean;
98 
99 for (SizeT i = dimx - 1 - w; i < dimx - 1; ++i) {
100  if (n > 0) res[i] = mean2;
101  DDouble v = data[i - w];
102  if (isfinite(v)) {
103   mean2 *= n;
104   mean2 -= v;
105   n -= 1.0;
106   mean2 /= n;
107  }
108 
109  if (n <= 0) mean2 = 0;
110 
111 #if defined (EDGE_WRAP)
112  v = data[i - dimx + w + 1];
113  if (isfinite(v)) {
114 #elif defined (EDGE_TRUNCATE)
115  v = data[dimx - 1];
116  if (isfinite(v)) {
117 #elif defined (EDGE_MIRROR)
118  v = data[2 * dimx - i - w - 2];
119  if (isfinite(v)) {
120 #elif defined (EDGE_ZERO)
121  v = 0;
122  {
123 #endif
124   mean2 *= n;
125   if (n < 2 * w + 1) n += 1.0;
126   mean2 += v;
127   mean2 /= n;
128  }
129 }
130 if (n > 0) res[dimx - 1] = mean2;
131 #endif
132 
133 #endif
134