- 01
- 02
- 03
- 04
- 05
- 06
- 07
- 08
- 09
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
template < int Order, typename T >
struct PrefixSum {
static inline void update ( T* a ) throw () { *a += *(a-1); PrefixSum < Order-1, T > :: update( a+1 ); }
};
template < typename T >
struct PrefixSum < 1,T > {
static inline void update ( T* a ) throw () { *a += *(a-1); }
};
template < int P, int N, int Condition = 0 >
struct Bpn {
enum { value = N * ( Bpn < P-1, N-1, (N > P) > :: value - Bpn < P-1, N, (N > P-1) > :: value ) };
};
template < int P, int N > struct Bpn < P, N, !0 > { enum { value = 0 }; };
template < int P > struct Bpn < P, 0, 0 > { enum { value = P ? 0 : 1 }; };
template < int P > struct Bpn < P, 1, 0 > { enum { value = P & 1 ? 1 : -1 }; };
template < typename Ta, typename Tb, bool C > struct IfThenElse;
template < typename Ta, typename Tb > struct IfThenElse < Ta, Tb, true > { typedef Ta TRes; };
template < typename Ta, typename Tb > struct IfThenElse < Ta, Tb, false > { typedef Tb TRes; };
template < int K, int I = 1 >
struct MomentSeries {
typedef typename IfThenElse < MomentSeries<K,I+1>, MomentSeries<K,K>, (I<K) >::TRes SubT;
static inline double accumulate ( double const* psum ) throw () {
return Bpn < K, I > :: value * psum [ I + 1 ] + SubT :: accumulate ( psum );
}
};
template < int K >
struct MomentSeries < K, K > {
static inline double accumulate ( double const* psum ) throw () {
return Bpn < K, K > :: value * psum [ K + 1 ];
}
};
template < int Order >
struct MomentLoop {
static inline void assign ( double *moments, size_t momentStride, double const* psum ) throw() {
*(moments - momentStride) = MomentSeries < Order-1 > :: accumulate ( psum );
MomentLoopAssign < Order-1 > :: assign ( moments - momentStride, momentStride, psum );
}
};
template <>
struct MomentLoop < 1 > {
static inline void assign ( double *moments, size_t momentStride, double const* psum ) throw() {
moments [ 0 ] = MomentSeries < 1, 1 > :: accumulate ( psum );
*(moments - momentStride) = psum [ 1 ];
}
};
/**
* Function computes a series of geometric moments by prefix summation method:
* Zhou F., Kornerup P. Computing Moments by Prefix Sums. // VLSI Signal Proc. - 2000. - 25. - P. 5 - 17.
* @param data is first data elemet address.
* @param ndataItems is number of data items.
* @param dataStride is the number of elements between two neighbor data items,
* in case of simple array dataStride is equal to 1.
* @param moments is address of 0-order moment.
* @param momentStride is number of elements between two neigbor moment items,
* in case of consequtive moments placement, momentStride is equal to 1.
*/
template < int Order, class OutPolicy >
inline void psmoment ( double const* data,
size_t const ndataItems,
size_t const dataStride,
double* moments,
size_t const momentStride ) throw() {
double psum [ Order+1 ] = { 0 };
// Initialize prefix sum
for ( size_t i = ndataItems, j = ndataItems * dataStride; i; ) {
--i; j -= dataStride; psum [ 0 ] = data [ j ];
PrefixSum < Order, double > :: update ( psum+1 );
}
// convert psum to moment values
OutPolicy :: assign ( moments + Order * momentStride, momentStride, psum );
}
Follow us!