diff --git a/wavelet.c b/wavelet.c index 0cc4c8e..3d6b1cf 100644 --- a/wavelet.c +++ b/wavelet.c @@ -1,292 +1,294 @@ /* ========================================================================= Lifting scheme implementation of a Daubechies D4 wavelet transform on arbitrary sized arrays (i.e. not restricted to powers of 2) of 3-tuples. Array boundaries are extended as constants. The implementation here is mainly distilled from the following sources: -- Ingrid Daubechies and Wim Sweldens, "Factoring Wavelet Tranforms into Lifting Steps", section 7.5: https://services.math.duke.edu/~ingrid/publications/J_Four_Anal_Appl_4_p247.pdf [Lifting scheme on D4] -- Ian Kaplan's webpages on wavelets and signal processing:: http://bearcave.com/misl/misl_tech/wavelets/daubechies/index.html http://bearcave.com/misl/misl_tech/wavelets/lifting/basiclift.html [Lifting scheme on D4] -- Stefan Moebius' TurboWavelets.Net code: https://github.com/codeprof/TurboWavelets.Net/tree/master/TurboWavelets [Boundary condition handling with arbitrary sizes, backreferencing] -- "Numerical Recipes in C", chapter 13.10: http://cacs.usc.edu/education/cs653/NRC-c13.10-Wavelets.pdf [Linear algebra (matrix mult.) implementation of D4] Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id$ */ #include "wavelet.h" #include #include +#include -#define ELEMSIZE (3 * sizeof(double)) +#define COEFSIZE (3 * sizeof(double)) /* Daubechies lifting coefficients; note c3 * c4 = 1 */ -static const double c0 = sqrt(3), - c1 = c0 / 4, - c2 = c1 - 0.5, - norm = 1/sqrt(2), - c3 = norm * (c0 - 1), - c4 = norm * (c0 + 1); +static const double c0 = sqrt(3), + c1 = c0 / 4, + c2 = c1 - 0.5, + norm = 1 / sqrt(2), + c3 = norm * (c0 - 1), + c4 = norm * (c0 + 1); static void d4WaveletXformStep ( double xIn [][3], double xOut [][3], unsigned len ) /* Single step (at fixed resolution) of a forward Daubechies D4 wavelet transform on array xIn of arbitrary size len containing 3-tuples. Note this is a 1D transform, and the triplets per array element are _not_ decorrelated, but transformed independently. This implementation use the lifting scheme to reduce the number of FLOPS by leveraging redundancies inherent in the transform compared to a "naive" algebraic (matrix) implementation. However, it is not performed in-place, but rather permutes coefficients on-the-fly to xOut for efficiency, as shuffling data within large arrays is expected to be more expensive than FLOPS. */ { /* Update/predict terms for current/previous/next iteration */ double updt [3], prevUpdt [3], nextUpdt [3], pred [3], nextPred [3]; unsigned i, j; double *coef (double x [][3], int idx) /* Array access with boundary extension (constant in this case) */ { - return idx < len ? x [idx] : idx >= 0 ? x [len - 1] : x [0]; + return idx >= len ? x [len - 1] : idx < 0 ? x [0] : x [idx]; } /* Init update/predict for previous and current iteration; do so for each 3-tuple item */ for (j = 0; j < 3; j++) { /* Update for iteration -1; assume constant signal beyond boundary as substitute for: xIn [-2][j] + c0 * xIn [-1][j] */ prevUpdt [j] = (1 + c0) * xIn [0] [j]; /* First update (= average, low pass) */ updt [j] = xIn [0][j] + c0 * xIn [1][j]; /* First predict (= detail, high/band pass) */ pred [j] = xIn [1][j] - c1 * updt [j] - c2 * prevUpdt [j]; } for (i = 0; i < len; i += 2) { /* Iterate over triplet */ for (j = 0; j < 3; j++) { /* Next update/predict (= average/detail = low/high pass) */ nextUpdt [j] = coef(xIn, i + 2)[j] + c0 * coef(xIn, i + 3)[j]; nextPred [j] = coef(xIn, i + 3)[j] - c1 * nextUpdt [j] - c2 * updt [j]; /* Current 2nd update (= average = low pass); normalise and store in lower half of xOut */ xOut [i >> 1][j] = c3 * (updt [j] - nextPred [j]); /* Current predict (= detail = high pass); normalise and store in upper half of xOut */ xOut [i + len >> 1][j] = c4 * pred [j]; /* Prepare next iteration */ updt [j] = nextUpdt [j]; pred [j] = nextPred [j]; } } } static void d4WaveletInvXformStep ( double xIn [][3], double xOut [][3], unsigned len ) /* Single step (at fixed resolution) of an inverse Daubechies D4 wavelet transform on coefficient array xIn of arbitrary size len containing 3-tuples. This is essentially reverses the forward transform above; see additional comments there for implementation details. */ { /* Update/predict terms for current/previous/next iteration */ double updt [3], prevUpdt [3], nextUpdt [3], pred [3], nextPred [3]; unsigned i, j; double *coef (double x [][3], int idx) /* Array access with boundary extension (constant in this case) */ { - return idx < len ? x [idx] : idx >= 0 ? x [len - 1] : x [0]; + return idx >= len ? x [len - 1] : idx < 0 ? x [0] : x [idx]; } /* Init update/predict for previous and current iteration; do so for each triplet item */ for (j = 0; j < 3; j++) { /* Invert first normalisation (by swapping reciprocal coeffs) and predict/update from avg/diff coeffs in lower/upper half of xIn */ pred [j] = c3 * xIn [len >> 1][j]; - // updt [j] = c4 * xIn [0][j]; + updt [j] = c4 * xIn [0][j]; /* Invert update for iteration -1; assume constant signal beyond boundary as substitute for: c4 * xIn [-1][j] + pred [j] */ - // prevUpdt [j] = updt [j] + pred [j]; - prevUpdt [j] = c4 * xIn [0][j] + pred [j]; + prevUpdt [j] = updt [j] + pred [j]; + // prevUpdt [j] = c4 * xIn [0][j] + pred [j]; } for (i = 0; i < len; i += 2) { /* Iterate over triplet */ for (j = 0; j < 3; j++) { /* Invert next predict/update from upper/lower half of xIn */ nextPred [j] = c3 * coef(xIn, (i + len >> 1) + 1) [j]; - // nextUpdt [j] = c4 * coef(xIn, (i >> 1) + 1) [j]; + nextUpdt [j] = c4 * coef(xIn, (i >> 1) + 1) [j]; /* Invert current 2nd update */ - updt [j] = c4 * xIn [i >> 1][j] + nextPred [j]; - // updt [j] += nextPred [j]; + // updt [j] = c4 * xIn [i >> 1][j] + nextPred [j]; + updt [j] += nextPred [j]; /* Invert current predict, store linearly in xOut */ coef(xOut, i + 1) [j] = pred [j] + c1 * updt [j] + c2 * prevUpdt [j]; /* Invert current 1st update, store linearly in xOut */ xOut [i][j] = updt [j] - c0 * coef(xOut, i + 1) [j]; /* Prepare next iteration */ prevUpdt [j] = updt [j]; - // updt [j] = nextUpdt [j]; + updt [j] = nextUpdt [j]; pred [j] = nextPred [j]; } } } void d4WaveletXform (double xIn [][3], unsigned len) /* Perform full multiresolution forward Daubechies D4 wavelet transform on array xIn of length len containing 3-tuples. Note no intra-tuple transform occurs (1-D only). Returns wavelet coefficients in xIn, containing average in xIn [0], followed by details in order or increasing resolution. This implementation uses an auxiliary array xOut to permute and "bounce" coefficients between xIn and xOut alternatively for each iteration. */ { /* Auxiliary array, array pointers; the latter alternate per iteration */ double xOut [len][3], (*x0) [3] = xIn, (*x1) [3] = xOut, (*xTmp) [3]; - unsigned l; + unsigned l, h; - for (l = len; l >= 2; l >>= 1) { - d4WaveletXformStep(x0, x1, l); + for (l = len, h = len >> 1; l >= 2; l = h, h >>= 1) { + d4WaveletXformStep(x0, x1, l); + /* Carry over detail coeffs in upper half for next iteration */ + memcpy(x0 [h], x1 [h], h * COEFSIZE); /* Swap array pointers */ xTmp = x0; x0 = x1; x1 = xTmp; } - if (x1 != xIn) + if (x0 != xIn) /* Coeffs not in xIn, so copy them there (is memcpy safe?) */ - memcpy(xIn, xOut, len * ELEMSIZE); + memcpy(xIn, xOut, len * COEFSIZE); } void d4WaveletInvXform (double xIn [][3], unsigned len) /* Perform full multiresolution inverse Daubechies D4 wavelet transform on array xIn of length len containing 3-tuples. Note no intra-tuple transform occurs (1-D only). Returns reconstructed signal in xIn. This implementation uses an auxiliary array xOut to permute and "bounce" coefficients between xIn and xOut alternatively for each iteration. */ { /* Auxiliary array, array pointers; the latter alternate per iteration */ double xOut [len][3], (*x0) [3] = xIn, (*x1) [3] = xOut, (*xTmp) [3]; unsigned l; - for (l = 2; l <= len; l <<= 1) { + for (l = 2; l <= len; l <<= 1) { d4WaveletInvXformStep(x0, x1, l); /* Swap array pointers */ xTmp = x0; x0 = x1; x1 = xTmp; } - if (x1 != xIn) - /* Data not in xIn, so copy it there (is memcpy safe?) */ - memcpy(xIn, xOut, ELEMSIZE); + if (x0 != xIn) + /* Recovered signal not in xIn, so copy it there (is memcpy safe?) */ + memcpy(xIn, xOut, len * COEFSIZE); } #ifdef _WAVELET_TEST #include #include #define ERRSTRLEN 1024 - static int tupleCmp (const void *t0, const void *t1) /* Comparison callback for sorting 3-tuples; reduce by summation */ { int i; double sum = 0; for (i = 0; i < 3; i++) sum += ((double *)t0) [i] - ((double *)t1) [i]; return sum > 0 ? 1 : sum < 0 ? -1 : 0; } int main (int argc, char *argv []) { unsigned i, j, len; double (*x0) [3], (*x1) [3]; char errstr [ERRSTRLEN]; if (argc < 2) { fputs("Missing array length\n", stderr); return -1; } if (!(len = atoi(argv [1]))) { fputs("Invalid array length\n", stderr); return -1; } /* Allocate arrays for original and reconstruction */ - if (!(x0 = calloc(len, ELEMSIZE)) || !(x1 = calloc(len, ELEMSIZE))) { + if (!(x0 = calloc(len, COEFSIZE)) || !(x1 = calloc(len, COEFSIZE))) { fputs("Failed array allocation\n", stderr); return -1; } /* Init with random data, sort by tuple sum */ srand48(111); for (i = 0; i < len; i++) { #if 1 x0 [i][0] = x0 [i][1] = x0 [i][2] = drand48(); #else for (j = 0; j < 3; j++) x0 [i][j] = drand48(); #endif } - qsort(x0, len, ELEMSIZE, tupleCmp); - memcpy(x1, x0, (size_t)len * ELEMSIZE); + qsort(x0, len, COEFSIZE, tupleCmp); + memcpy(x1, x0, (size_t)len * COEFSIZE); /* Transform and reconstruct data */ d4WaveletXform(x1, len); - //d4WaveletInvXform(x1, len); + d4WaveletInvXform(x1, len); for (i = 0; i < len; i++) #if 0 printf( "%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n", x0 [i][0], x0 [i][1], x0 [i][2], x1 [i][0], x1 [i][1], x1 [i][2] ); #else printf("%.6f\t%.6f\n", x0 [i][0], x1 [i][0]); #endif free(x0); free(x1); return 0; } #endif diff --git a/wavelet.h b/wavelet.h index 29d17d3..73dbee1 100644 --- a/wavelet.h +++ b/wavelet.h @@ -1,39 +1,33 @@ /* ========================================================================= Lifting scheme implementation of a Daubechies D4 wavelet transform on - arbitrary sized arrays (i.e. not restricted to powers of 2) of arbitrary - record types containing double 3-tuples. Records are addressed via a - supplied callback. Array boundaries are extended as constants. + arbitrary sized arrays (i.e. not restricted to powers of 2) of double + 3-tuples. Array boundaries are extended as constants. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id$ */ #ifndef _PMAP_WAVELET_H #define _PMAP_WAVELET_H - /* Callback to return coefficients in idx'th record in array x as double - 3-tuples */ - typedef double [3] (*)(const void *x, unsigned idx) WaveletCoeffCallback; - - - void d4WaveletXform (void *xIn, unsigned len, WaveletCoeffCallback coeff); + void d4WaveletXform (double xIn [][3], unsigned len); /* Perform full multiresolution forward Daubechies D4 wavelet transform - on array xIn of length len. Coefficients in records are accessed via - coeff as 3-tuples. Note no intra-tuple transform occurs (1-D only). + on array xIn of length len containing 3-tuples of type double. + Note no intra-tuple transform occurs (1-D only). Returns wavelet coefficients in xIn, containing average in xIn [0], followed by details in order or increasing resolution. */ - void d4WaveletInvXform (void *xIn, unsigned len, WaveletCoeffCallback coeff); + void d4WaveletInvXform (double xIn [][3], unsigned len); /* Perform full multiresolution inverse Daubechies D4 wavelet transform - on array xIn of length len. Coefficients in records are accessed via - coeff as 3-tuples. Note no intra-tuple transform occurs (1-D only). + on array xIn of length len containing 3-tuples of type double. + Note no intra-tuple transform occurs (1-D only). Returns reconstructed data in xIn. */ #endif