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