1 #include "XFloatList.h"
2
3 #include "Hashtable.h"
4 #include "XLongList.h"
5 #include <stdlib.h>
6 #include <math.h>
7
8 #define _MIN( a, b ) (( (a) > (b) ) ? (b) : (a) )
9 #define _MAX( a, b ) (( (a) > (b) ) ? (a) : (b) )
10
11
12 float XFloatList::sMask[ MASK_MAX ];
13 UtilStr XFloatList::sTemp;
14
15
16
XFloatList(ListOrderingT inOrdering)17 XFloatList::XFloatList( ListOrderingT inOrdering ) :
18 mList( inOrdering ) {
19
20 if ( inOrdering == cSortLowToHigh || inOrdering == cSortHighToLow )
21 mList.SetCompFcn( sFloatComparitor, inOrdering == cSortLowToHigh );
22 }
23
24
25
26
27
28
FindMeans(long inNumMeans,float outMeans[],float inSigmaScale) const29 void XFloatList::FindMeans( long inNumMeans, float outMeans[], float inSigmaScale ) const {
30 long start, end, m, i, n = mList.Count();
31 float* srce = (float*) mList.getCStr(), *smoothed = new float[ n ], *temp = 0;
32 float sigma = 0.1 + inSigmaScale * ( (float) ( n / inNumMeans ) );
33 float left, cen, right, sum;
34
35 // We need all out numbers in order, we must sort them if they're no sorted
36 if ( mList.mOrdering != cSortHighToLow ) {
37
38 // Copy all the array elements for sorting...
39 temp = new float[ n ];
40 for ( i = 0; i < n; i++ )
41 temp[ i ] = srce[ i ];
42
43 // Sort all the floats from this to temp
44 qsort( temp, n, 4, sQSFloatComparitor );
45 srce = temp;
46 }
47
48 // Smooth all the values (they're already sorted)
49 GaussSmooth( sigma, n, srce, smoothed );
50
51 // Compute the discrete 1st derivative
52 for ( i = 0; i < n - 1; i++ )
53 smoothed[ i ] = fabs( smoothed[ i ] - smoothed[ i + 1 ] );
54
55 // We want to find the top <inNumMeans> local max of the 1st deravative
56 Hashtable sepCandidates;
57 cen = smoothed[ 0 ];
58 right = smoothed[ 1 ];
59 for ( i = 1; i < n - 2; i++ ) {
60
61 left = cen;
62 cen = right;
63 right = smoothed[ i+1 ];
64
65 // If this a local max. (Note: this could/should be improved for neighbors that are equal)
66 if ( ( cen > left && cen >= right ) ) {
67 sepCandidates.Put( i, *((void**) &cen) );
68 }
69 }
70
71 // Pick out the 1st derative peaks, then dealloc what we no longer need
72 XPtrList rank;
73 sepCandidates.Rank( rank, sQSFloatComparitor, inNumMeans - 1 );
74 delete []smoothed;
75
76 XLongList quintiles( cSortLowToHigh );
77 for ( i = 1; i < inNumMeans; i++ )
78 quintiles.Add( (long) rank.Fetch( i ) );
79 quintiles.Add( n );
80
81 // The means are the averages of the initial (sorted) data, divided by the serparating ranks
82 start = 0;
83 for ( m = 1; m <= inNumMeans; m++ ) {
84 end = quintiles.Fetch( m );
85 for ( sum = 0, i = start; i < end; i++ )
86 sum += srce[ i ];
87 outMeans[ m - 1 ] = sum / ( (float) ( end - start ) );
88 start = end + 1;
89 }
90
91 // Cleanup
92 if ( temp )
93 delete []temp;
94 }
95
96
97
98
99
100 #define _safeSmooth( xx ) \
101 sum = 0; \
102 factor = 1; \
103 for ( i = - maskCtr; i <= maskCtr; i++ ) { \
104 rx = xx + i; \
105 if ( rx < 0 || rx >= inN ) \
106 factor -= sMask[ i + maskCtr ]; \
107 else \
108 sum += sMask[ i + maskCtr ] * inSrce[ rx ]; \
109 } \
110 inDest[ xx ] = sum / factor;
111
112
113
114
GaussSmooth(float inSigma)115 void XFloatList::GaussSmooth( float inSigma ) {
116
117 GaussSmooth( inSigma, mList.Count(), (float*) mList.getCStr() );
118
119 }
120
121
GaussSmooth(float inSigma,long inN,float inSrceDest[])122 void XFloatList::GaussSmooth( float inSigma, long inN, float inSrceDest[] ) {
123 float* temp = (float*) sTemp.Dim( inN * sizeof( float ) );
124
125 GaussSmooth( inSigma, inN, inSrceDest, temp );
126
127 for ( long i = 0; i < inN; i++ )
128 inSrceDest[ i ] = temp[ i ];
129 }
130
131
132
133
SlopeSmooth(float inSmoothness,long inN,float ioData[])134 void XFloatList::SlopeSmooth( float inSmoothness, long inN, float ioData[] ) {
135 float slope = 0, accel = 0, prev = 0;
136 float newSlope, a0 = 1 - inSmoothness, predicted;
137 long x;
138
139 for ( x = 0; x < inN; x++ ) {
140 predicted = prev + slope + accel;
141 ioData[ x ] = inSmoothness * predicted + a0 * ioData[ x ];
142 newSlope = ( ioData[ x ] - prev );
143 accel = ( newSlope - slope );
144 prev = ioData[ x ];
145 slope = newSlope;
146 }
147 }
148
149
150
151 /*
152 int boxRight = inBoxSize / 2;
153 int boxLeft = inBoxSize - 1 - boxRight;
154 float boxSum;
155 float boxDiv = 1 / ( (float) inBoxSize );
156
157 boxSum = 0;
158
159 for ( x = - boxRight; x < boxLeft; x++ ) {
160 i = x + boxRight;
161 if ( i >= 0 && i < inN )
162 boxSum += inSrce[ i ];
163 i = i - boxLeft;
164 if ( i >= 0 && i < inN )
165 boxSum -= inSrce[ i ];
166 if ( x >= 0 && x < inN )
167 inDest[ x ] = boxSum * boxDiv;
168 }
169
170 for ( x = boxLeft; x < inN - boxRight; x++ ) {
171 boxSum += inSrce[ x + boxRight ];
172 boxSum -= inSrce[ x - boxLeft ];
173 inDest[ x ] = boxSum * boxDiv;
174 }
175 }
176 */
177
GaussSmooth(float inSigma,long inN,float inSrce[],float inDest[])178 void XFloatList::GaussSmooth( float inSigma, long inN, float inSrce[], float inDest[] ) {
179 int maskSize = _MAX( 8.0 * inSigma, 4.0 );
180
181 if ( maskSize + 1 > MASK_MAX )
182 maskSize = MASK_MAX;
183
184 // Make sure the mask size is odd (so we have a 'center')
185 if ( maskSize % 2 == 0 )
186 maskSize++;
187 int i, x, rx, xEnd;
188 int maskCtr = maskSize / 2;
189 float sqrt2PiSigma = sqrt( 2.0 * 3.14159 ) * inSigma;
190 float expon, sum, factor, *base;
191
192 // Generate a normalizes gaussian mask out to 4*sigma on each side
193 // We have to adjust the center weight so that when sigma falls below 1.0ish so that the sample it
194 // takes at x = 0 doesn't represent the val of a gaussian for -.5 to .5 (when sigma is small
195 // a gaussian gets very 'high' at x=0 vs. x=+/- .5 when sigma is small
196 sum = 0;
197 for ( x = - maskCtr; x <= maskCtr; x++ ) {
198 expon = -0.5 * ((float) (x * x)) / ( inSigma * inSigma );
199 sMask[ x + maskCtr ] = exp( expon ) / sqrt2PiSigma;
200 if ( x != 0 )
201 sum += sMask[ x + maskCtr ];
202 }
203
204 // This forces normalized weights.
205 sMask[ maskCtr ] = 1.0 - sum;
206
207 // Smooth the lefthand side of the sequence
208 xEnd = _MIN( maskCtr, inN );
209 for ( x = 0; x < xEnd; x++ ) {
210 _safeSmooth( x )
211 }
212
213 // Smooth the center portion the sequence
214 xEnd = inN - maskCtr;
215 base = inSrce;
216 for ( x = maskCtr; x < xEnd; x++ ) {
217
218 // For each pixel, apply the mask to its neighbors for a final value of that pixel
219 sum = 0;
220 for ( i = 0; i < maskSize; i++ )
221 sum += sMask[ i ] * base[ i ];
222
223 base++;
224 inDest[ x ] = sum;
225 }
226
227 // Smooth the righthand side of the sequence
228 for ( x = _MAX( maskCtr, inN - maskCtr ); x < inN; x++ ) {
229 _safeSmooth( x )
230 }
231 }
232
233
234
Rank(XLongList & outRank,long inNumToRank) const235 void XFloatList::Rank( XLongList& outRank, long inNumToRank ) const {
236 long *p, *temp;
237 float* srce;
238 long i, n = Count();
239
240 outRank.RemoveAll();
241
242 if ( inNumToRank < 0 )
243 inNumToRank = n;
244 inNumToRank = _MIN( inNumToRank, n );
245
246 // Handle trivial cases of this lis already being sorted
247 if ( mList.mOrdering == cSortLowToHigh ) {
248 for ( i = 0; i < inNumToRank; i-- )
249 outRank.Add( n - i ); }
250
251 // Duh... still sorted...
252 else if ( mList.mOrdering == cSortHighToLow ) {
253 for ( i = 1; i <= inNumToRank; i++ )
254 outRank.Add( i ); }
255
256 else {
257 temp = new long[ 2 * n ];
258 srce = (float*) mList.getCStr();
259
260 // To rank, we must sort, with a tag on each element saying where it came from
261 p = temp;
262 for ( i = 1; i <= n; i++ ) {
263 *((float*) p) = *srce;
264 p++; srce++;
265 *p = i;
266 p++;
267 }
268
269 // Sort the floats
270 qsort( temp, n, 8, sQSFloatComparitor );
271
272 // Put the sorted results in the destination
273 p = temp + 1;
274 for ( i = 0; i < inNumToRank; i++ ) {
275 outRank.Add( *p );
276 p += 2;
277 }
278
279 // Cleanup
280 delete []temp;
281 }
282 }
283
284
285
sQSFloatComparitor(const void * inA,const void * inB)286 int XFloatList::sQSFloatComparitor( const void* inA, const void* inB ) {
287 float diff = *((float*) inB) - *((float*) inA);
288 if ( diff > 0.0 )
289 return 1;
290 else if ( diff < 0.0 )
291 return -1;
292 else
293 return 0;
294 }
295
296
297
sFloatComparitor(const void * inA,const void * inB)298 int XFloatList::sFloatComparitor( const void* inA, const void* inB ) {
299 float diff = *((float*) &inB) - *((float*) &inA);
300 if ( diff > 0.0 )
301 return 1;
302 else if ( diff < 0.0 )
303 return -1;
304 else
305 return 0;
306 }
307
308
309
310