diff --git a/wavelet.c b/wavelet.c index 453c24e..c790734 100644 --- a/wavelet.c +++ b/wavelet.c @@ -1,645 +1,499 @@ /* ========================================================================= Lifting scheme implementation of Daubechies D4 and Haar wavelet transforms on arbitrary sized arrays (i.e. not restricted to powers of 2) of 3-tuples. Array boundaries are extended cyclically. 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 #include -#define COEFSIZE (3 * sizeof(double)) +#define COEFSIZE (3 * sizeof(double)) +#define min(a, b) ((a) < (b) ? (a) : (b)) +/* Daubechies D4 wavelet coeffs */ +static const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3), + h [4] = { + norm * (1 + sqrt3), norm * (3 + sqrt3), + norm * (3 - sqrt3), norm * (1 - sqrt3) + }, + g [4] = { h [3], -h [2], h [1], -h [0] }, + /* Summed coeffs for padding at left boundary */ + hBound = h [0] + h [1] + h [2], + gBound = g [0] + g [1] + g [2]; +/* Haar wavelet coeffs */ +static const double c0 = 1 / sqrt(2); -#if 0 /* THIS STUFF IS BULLSHIT AND DOESN'T WORK! */ -/* 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 void d4LiftStep (double x [][3], int len) -/* Single step (at fixed resolution) of a forward Daubechies D4 wavelet - transform on array x 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 uses the lifting scheme to reduce the number of FLOPS - by leveraging redundancies inherent in the transform compared to a - "naive" algebraic implementation using convolution. However, it is not - performed in-place, but rather permutes coefficients on-the-fly to an - auxiliary array y for efficiency, as shuffling data within large arrays - is expected to be more expensive than FLOPS. */ -{ - /* Aux array, update/predict for current/previous/next iteration */ - double y [len][3], updt [3], nextUpdt [3], pred [3], nextPred [3]; - int i, j; - - /* Haar wavelet for 1st iteration to handle boundary effects */ - for (j = 0; j < 3; j++) { - /* Predict (= detail = high pass) */ - y [len >> 1][j] = pred [j] = x [1][j] - x [0][j]; - /* Update (= average = low pass) */ - y [0][j] = updt [j] = x [0][j] + 0.5 * pred [j]; - } - - for (i = 0; i < len - 3; i += 2) { - /* Iterate over triplet */ - for (j = 0; j < 3; j++) { - /* Next update/predict (= average/detail = low/high pass) */ - nextUpdt [j] = x [i + 2][j] + c0 * x [i + 3][j]; - nextPred [j] = x [i + 3][j] - c1 * nextUpdt [j] - c2 * updt [j]; - - if (i > 0) { - /* D4 wavelet for subsequent iterations */ - /* Current 2nd update (= average = low pass); - normalise and store in lower half of y */ - y [i >> 1][j] = c3 * (updt [j] - nextPred [j]); - /* Current predict (= detail = high pass); - normalise and store in upper half of y */ - y [i + len >> 1][j] = c4 * pred [j]; - } - - /* Prepare next iteration */ - updt [j] = nextUpdt [j]; - pred [j] = nextPred [j]; - } - } - - if (i < len) { - /* Haar wavelet for last iteration to handle boundary effects */ - for (j = 0; j < 3; j++) { - /* Predict (= detail = high pass) */ - y [i + len >> 1][j] = pred [j] = x [i + 1][j] - x [i][j]; - /* Update (= average = low pass) */ - y [i >> 1][j] = updt [j] = x [i][j] + 0.5 * pred [j]; - } - } - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); -} - - - -static void d4InvLiftStep (double x [][3], int 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. */ -{ - /* Aux array, update/predict for current/previous/next iteration */ - double y [len][3], updt [3], prevUpdt [3], nextUpdt [3], - pred [3], nextPred [3]; - int i, j; - - /* Haar wavelet for 1st iteration to handle boundary effects */ - for (j = 0; j < 3; j++) { - /* Invert predict/update from upper/lower half of x */ - pred [j] = x [len >> 1][j]; - updt [j] = x [0][j]; - y [0][j] = updt [j] - 0.5 * pred [j]; - y [1][j] = pred [j] + y [0][j]; - } - - for (i = 0; i < len - 3; i += 2) { - /* Iterate over triplet */ - for (j = 0; j < 3; j++) { - /* Invert next predict/update from upper/lower half of x */ - nextPred [j] = c3 * x [(i + len >> 1) + 1][j]; - nextUpdt [j] = c4 * x [(i >> 1) + 1][j]; - - if (i > 0) { - /* D4 wavelet for subsequent iterations */ - /* Invert current 2nd update */ - updt [j] += nextPred [j]; - /* Invert current predict, store linearly in y */ - y [i + 1][j] = pred [j] + c1 * updt [j] + c2 * prevUpdt [j]; - /* Invert current 1st update, store linearly in y */ - y [i][j] = updt [j] - c0 * y [i + 1][j]; - } - - /* Prepare next iteration */ - prevUpdt [j] = updt [j]; - updt [j] = nextUpdt [j]; - pred [j] = nextPred [j]; - } - } +static void haarStep (double (*x) [3], int len, double (*y) [3]) +/* + Single step (at fixed resolution) of a forward Haar wavelet transform on + array x of arbitrary size len containing 3-tuples. - if (i < len) { - /* Haar wavelet for last iteration to handle boundary effects */ - for (j = 0; j < 3; j++) { - /* Invert predict/update from upper/lower half of x */ - pred [j] = x [i + len >> 1][j]; - updt [j] = x [i >> 1][j]; - y [i][j] = updt [j] - 0.5 * pred [j]; - y [i + 1][j] = pred [j] + y [i][j]; - } - } + Note this is a 1D transform, and the triplets per array element are _not_ + decorrelated, but transformed independently. - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); -} -#endif /* BULLSHIT */ - - - -static void haarConvStep (double x [][3], int len) -/* - Single step (at fixed resolution) of a forward Haar wavelet - transform on array x 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. + The wavelet coefficients are returned in the preallocated array y. */ { - /* Aux array, coeffs */ - double y [len][3]; - const double c0 = 1 / sqrt(2); const unsigned hl = -(-len >> 1); int i, hi, j; if (len >= 2) { /* First and intermediate iterations */ for (i = 0; i <= len - 2; i += 2) { hi = i >> 1; for (j = 0; j < 3; j++) { /* Smooth */ y [hi][j] = c0 * (x [i][j] + x [i + 1][j]); /* Detail */ y [hl + hi][j] = c0 * (x [i][j] - x [i + 1][j]); } } if (len & 1) { - /* Odd length; pad with last value as additional smoothing term */ + /* Odd length; pad with last value as additional approx term */ for (j = 0; j < 3; j++) y [hl - 1][j] = x [len - 1][j]; } #ifdef _WAVELET_DBG printf("FWD HAAR XFORM (len = %d)\n", len); for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); + printf("%.3f\t-->\t%.3f\n", x [i][0], y [i][0]); putchar('\n'); #endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); } } -static void haarInvConvStep (double x [][3], int len) +static void haarInvStep (double (*y) [3], int len, double (*x) [3]) /* Single step (at fixed resolution) of an inverse Haar wavelet - transform on coefficient array xIn of arbitrary size len containing + transform on coefficient array y of size len containing 3-tuples. This is essentially reverses the forward transform above; - see additional comments there for implementation details. + see additional comments there for implementation details. + + The reconstructed signal is return in the preallocated array x. */ { - /* Aux array, coeffs */ - double y [len][3]; - const double c0 = 1 / sqrt(2); const unsigned hl = -(-len >> 1); int i, hi, j; if (len >= 2) { for (i = 0; i <= len - 2; i += 2) { hi = i >> 1; for (j = 0; j < 3; j++) { - y [i][j] = c0 * (x [hi][j] + x [hl + hi][j]); - y [i + 1][j] = c0 * (x [hi][j] - x [hl + hi][j]); + x [i][j] = c0 * (y [hi][j] + y [hl + hi][j]); + x [i + 1][j] = c0 * (y [hi][j] - y [hl + hi][j]); } } if (len & 1) { - /* Odd length; pad additional value from last smoothing term */ + /* Odd length; pad additional value from last approx term */ for (j = 0; j < 3; j++) - y [len - 1][j] = x [hl - 1][j]; + x [len - 1][j] = y [hl - 1][j]; } #ifdef _WAVELET_DBG printf("INV HAAR XFORM (len = %d)\n", len); for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); + printf("%.3f\t-->\t%.3f\n", y [i][0], x [i][0]); putchar('\n'); #endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); } } -static void qHaarConvStep (double x [][3], int len) +static unsigned d4Step (double (*x) [3], int xLen, double (*y) [3]) /* - Single step (at fixed resolution) of a forward quadratic Haar wavelet - transform on array x of arbitrary size len containing 3-tuples. + Single step (at fixed resolution) of forward Daubechies D4 wavelet + transform on array x of arbitrary size xLen containing 3-tuples. Note this is a 1D transform, and the triplets per array element - are _not_ decorrelated, but transformed independently. + are _not_ decorrelated, but transformed independently. + + The approximation and detail coefficients are returned in the + preallocated output array y. This array is larger than len to + accommodate padding coefficients to handle boundary effects. This size + is returned by the function, and can be interrogated beforehand by + passing either x or y as NULL pointer. */ -{ - /* Aux array, coeffs */ - double y [len][3]; - const double h0 = 0.5, g0 = sqrt(3) * h0; - const unsigned hl = -(-len >> 1); - int i, hi, j; - - if (len >= 3) { - /* First and intermediate iterations */ - for (i = 0; i <= len - 3; i += 2) { - hi = i >> 1; - for (j = 0; j < 3; j++) { - /* Smooth */ - y [hi][j] = - h0 * (x [i][j] + x [i + 1][j] + x [i + 2][j]); - /* Detail */ - y [hl + hi][j] = - g0 * x [i + 1][j] - h0 * (x [i][j] + x [i + 2][j]); - } - } - - if (i <= len - 2) { - /* Last iteration with boundary extension */ - hi = i >> 1; +{ + unsigned i, j, k, hi; + const unsigned yLen = xLen + 3 >> 1, /* Output length + boundary pad */ + hy = yLen; /* Offset for detail coeffs */ + + if (x && y) { + if (xLen >= 2) { +#ifdef _WAVELET_DBG + /* Zero output array to facilitate debugging */ + memset(y, 0, 2 * yLen * COEFSIZE); +#endif + + /* First iteration (i = -2); pad coeffs at left edge with constant + * boundary extension (note repeated x [0]) */ for (j = 0; j < 3; j++) { - /* Smooth */ - y [hi][j] = - h0 * (x [i][j] + x [i + 1][j] + x [0][j]); + /* Approximation */ + y [0][j] = hBound * x [0][j] + h [3] * x [1][j]; /* Detail */ - y [hl + hi][j] = - g0 * x [i + 1][j] - h0 * (x [i][j] + x [0][j]); + y [hy][j] = gBound * x [0][j] + g [3] * x [1][j]; + } + + /* Intermediate and last iterations (i >= 0) */ + for (i = 0; i < xLen; i += 2) { + hi = (i >> 1) + 1; + + for (j = 0; j < 3; j++) { + y [hi][j] = y [hy + hi][j] = 0; + + /* Convolve up to wavelet or input boundary, whichever comes + * first */ + for (k = i; k < min(i + 4, xLen); k++) { + /* Approximation */ + y [hi][j] += h [k - i] * x [k][j]; + /* Detail */ + y [hy + hi][j] += g [k - i] * x [k][j]; + } + + /* If necessary, convolve beyond input boundary with constant + * extension (note repeated x [xLen - 1]). This is skipped if + * the wavelet boundary was already reached */ + for (; k < i + 4; k++) { + /* Approximation */ + y [hi][j] += h [k - i] * x [xLen - 1][j]; + /* Detail */ + y [hy + hi][j] += g [k - i] * x [xLen - 1][j]; + } + } } - } - - if (len & 1) { - /* Odd length; pad with last value as additional smoothing term */ - for (j = 0; j < 3; j++) - y [hl - 1][j] = x [len - 1][j]; - } #ifdef _WAVELET_DBG - printf("FWD QHAAR XFORM (len = %d)\n", len); - for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); - putchar('\n'); -#endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); - } -} - - - -static void qHaarInvConvStep (double x [][3], int len) -/* - Single step (at fixed resolution) of an inverse quadratic Haar 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. -*/ -{ - /* Aux array, coeffs */ - double y [len][3]; - const double h0 = 0.5, g0 = sqrt(3) * h0; - const unsigned hl = -(-len >> 1); - int i, hi, j; - - if (len >= 3) { - /* First iteration with boundary extension */ - hi = len >> 1; - for (j = 0; j < 3; j++) { - y [0][j] = h0 * ( - x [0][j] - x [hl][j] + - x [hi - 1][j] - x [len - 1][j] + printf( + "D4 XFORM + boundary padding (len = %d -> %d)\n", + xLen, yLen << 1 ); - y [1][j] = h0 * x [0][j] + g0 * x [hl][j]; - } - - /* Intermediate and last iterations */ - for (i = 0; i <= len - 4; i += 2) { - hi = i >> 1; - for (j = 0; j < 3; j++) { - y [i + 2][j] = h0 * ( - x [hi][j] - x [hl + hi][j] + - x [hi + 1][j] - x [hl + hi + 1][j] - ); - y [i + 3][j] = h0 * x [hi][j] + g0 * x [hl + hi][j]; + for (i = 0; i < yLen << 1; i++) { + if (i < xLen) + printf("%.3f\t-->\t", x [i][0]); + else + printf("\t\t"); + printf("%.3f\n", y [i][0]); } - } - - if (len & 1) { - /* Odd length; pad additional value from last smoothing term */ - for (j = 0; j < 3; j++) - y [len - 1][j] = x [hl - 1][j]; - } - -#ifdef _WAVELET_DBG - printf("INV QHAAR XFORM (len = %d)\n", len); - for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); - putchar('\n'); + putchar('\n'); #endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); + } } + + return yLen; } -static void d4ConvStep (double x [][3], int len) +static unsigned d4InvStep (double (*y) [3], int yLen, double (*x) [3]) /* - Single step (at fixed resolution) of a forward Daubechies D4 wavelet - transform on array x 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. + Single step (at fixed resolution) of inverse Daubechies D4 wavelet + transform on coefficient array y of size 2 * yLen (i.e. yLen approx and + detail coeffs each) containing 3-tuples. This essentially reverses the + forward transform above; see additional comments there for implementation + details. + + The reconstructed signal is returned in the preallocated array x. */ { - /* Aux array, coeffs */ - double y [len][3]; - const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3), - h0 = norm * (1 + sqrt3), h1 = norm * (3 + sqrt3), - h2 = norm * (3 - sqrt3), h3 = norm * (1 - sqrt3), - g0 = h3, g1 = -h2, g2 = h1, g3 = -h0; - int i, hi, j; - const unsigned hl = -(-len >> 1); - - if (len >= 4) { - /* First and intermediate iterations */ - for (i = 0; i <= len - 4; i += 2) { - hi = i >> 1; - for (j = 0; j < 3; j++) { - /* Smooth */ - y [hi][j] = - h0 * x [i][j] + h1 * x [i + 1][j] + - h2 * x [i + 2][j] + h3 * x [i + 3][j]; - /* Detail */ - y [hl + hi][j] = - g0 * x [i][j] + g1 * x [i + 1][j] + - g2 * x [i + 2][j] + g3 * x [i + 3][j]; - } - } - - if (i <= len - 2) { - /* Last iteration with wrapped boundary extension */ - hi = i >> 1; - for (j = 0; j < 3; j++) { - /* Smooth */ - y [hi][j] = - h0 * x [i][j] + h1 * x [i + 1][j] + - h2 * x [0][j] + h3 * x [1][j]; - /* Detail */ - y [hl + hi][j] = - g0 * x [i][j] + g1 * x [i + 1][j] + - g2 * x [0][j] + g3 * x [1][j]; + unsigned xLen = 0, i, j, k, hi; + const unsigned hy = yLen; /* Offset for detail coeffs */ + + if (x && y) { + if (yLen >= 2) { + /* Reconstruct from approximation and detail coefficients */ + for (i = 0; i < yLen; i += 2) { + hi = i >> 1; + for (j = 0; j < 3; j++) { + x [i][j] = + h [2] * y [hi][j] + g [2] * y [hy + hi][j] + + h [0] * y [hi + 1][j] + g [0] * y [hy + hi + 1][j]; + x [i + 1][j] = + h [3] * y [hi][j] + g [3] * y [hy + hi][j] + + h [1] * y [hi + 1][j] + g [1] * y [hy + hi + 1][j]; + } } - } - - if (len & 1) { - /* Odd length; pad with last value as additional smoothing term */ - for (j = 0; j < 3; j++) - y [hl - 1][j] = x [len - 1][j]; - } + + /* Last reconstruction for even number of coefficients */ + if (!(yLen & 1)) { + hi = i >> 1; + for (j = 0; j < 3; j++) { + x [i][j] = + h [2] * y [hi][j] + g [2] * y [hy + hi][j] + + h [0] * y [hi + 1][j] + g [0] * y [hy + hi + 1][j]; + } + i++; + } + + xLen = i; #ifdef _WAVELET_DBG - printf("FWD D4 XFORM (len = %d)\n", len); - for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); - putchar('\n'); + printf("D4 INV XFORM (len = %d --> %d)\n", yLen << 1, xLen); + for (i = 0; i < yLen << 1; i++) { + printf("%.3f\t", y [i][0]); + if (i < xLen) + printf("-->\t%.3f\n", x [i][0]); + else putchar('\n'); + } + putchar('\n'); #endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); + } } + + return xLen; } -static void d4InvConvStep (double x [][3], int 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. +unsigned waveletXform (double (*x) [3], unsigned xLen, double (*y) [3]) +/* Perform full multiresolution forward wavelet transform + on array x of length xLen containing 3-tuples. Note no intra-tuple + transform occurs (1-D only). + + The wavelet coefficients are returned in preallocated array y, containing + the coarsest approximation in y [0] followed by details in order of + increasing resolution/frequency. + + Note y must accommodate more than xLen coefficients to accommodate + padding for additional boundary coefficients. The number of coefficients + can be interrogated beforehand by calling the function with either x or y + set to NULL. + + Note also that x is used as a buffer for intermediate results, and is + therefore modified. + + The returned value is the number of output coefficients in y. */ { - /* Aux array, coeffs */ - double y [len][3]; - const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3), - h0 = norm * (1 + sqrt3), h1 = norm * (3 + sqrt3), - h2 = norm * (3 - sqrt3), h3 = norm * (1 - sqrt3), - g0 = h3, g1 = -h2, g2 = h1, g3 = -h0; - int i, hi, j; - const unsigned hl = -(-len >> 1); - - if (len >= 4) { - /* First iteration with wrapped boundary extension */ - hi = len >> 1; - for (j = 0; j < 3; j++) { - y [0][j] = - h0 * x [0][j] + g0 * x [hl][j] + - h2 * x [hi - 1][j] + g2 * x [len - 1][j]; - y [1][j] = - h1 * x [0][j] + g1 * x [hl][j] + - h3 * x [hi - 1][j] + g3 * x [len - 1][j]; - } - - /* Intermediate and last iterations */ - for (i = 0; i <= len - 4; i += 2) { - hi = i >> 1; - for (j = 0; j < 3; j++) { - y [i + 2][j] = - h2 * x [hi][j] + g2 * x [hl + hi][j] + - h0 * x [hi + 1][j] + g0 * x [hl + hi + 1][j]; - y [i + 3][j] = - h3 * x [hi][j] + g3 * x [hl + hi][j] + - h1 * x [hi + 1][j] + g1 * x [hl + hi + 1][j]; - } - } - - if (len & 1) { - /* Odd length; pad additional value from last smoothing term */ - for (j = 0; j < 3; j++) - y [len - 1][j] = x [hl - 1][j]; + /* Output size (number of coeffs, including boundary padding), output + * array offset and pointers */ + unsigned numCoeffs = 0, yOff = 0, yLen = 0, l; + double (*yOffPtr) [3] = NULL; + + /* Figure out total number of coefficients (= minimum size of output + * array including boundary padding). */ + for ( + numCoeffs = l = d4Step(NULL, xLen, NULL); + l > 3; + numCoeffs += l = d4Step(NULL, l, NULL) + ); + numCoeffs += l; + + if (x && y) { + /* Apply Daubechies D4 transform, accounting for boundary padding */ + for (l = xLen, yOff = numCoeffs; l > 3; l = yLen) { + /* Get num coeffs and calculate offset into output array + * (subtracting num coeffs once more to account for approximation + * coeffs, which will be replaced in next iteration) */ + yOff -= yLen = d4Step(NULL, l, NULL); + yOffPtr = y + yOff - yLen; + /* Resulting wavelet coefficients in y are prepended to those from + * previous iteration; consequently, returned num coeffs is input + * size for next iteration */ + d4Step(x, l, yOffPtr); + /* Copy approx coeffs from y into x as input for next iteration */ + memcpy(x, yOffPtr, yLen * COEFSIZE); } -#ifdef _WAVELET_DBG - printf("INV D4 XFORM (len = %d)\n", len); - for (i = 0; i < len; i++) - printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]); - putchar('\n'); +#if 1 + /* Apply Haar transform in penultimate and final (since len is odd) + * iterations to decompose to a single approximation and detail + * coefficient at coarsest resolution at y [0] and y [1], + * respectively. Obviously no boundary padding required here. */ + /* HACK: shifting negated value rounds _up_, thus preserving + * additional iteration if length is odd. */ + for (; l >= 2; l--) { + haarStep(x, l, y); + /* Copy coeffs from y into x as input for next iteration */ + memcpy(x, y, l * COEFSIZE); + } #endif - - /* Return coeffs in x */ - memcpy(x, y, len * COEFSIZE); } + + /* All coeffs now in y; return their number */ + return numCoeffs; } -void waveletXform (double x [][3], unsigned len) -/* Perform full multiresolution forward wavelet transform - on array x of length len containing 3-tuples. Note no intra-tuple - transform occurs (1-D only). - Returns wavelet coefficients in x, containing average in x [0], - followed by details in order or increasing resolution. */ +unsigned waveletInvXform ( + double (*y) [3], unsigned yLen, double (*x) [3], unsigned xLen +) +/* Perform full multiresolution inverse wavelet transform on array y of + * length yLen containing wavelet coefficients as 3-tuples. Note no + * intra-tuple transform occurs (1-D only). The reconstructed signal is + * returned in preallocated array x of length xLen; the latter MUST + * correspond to the original signal length! + * + * Note xLen < yLen since additional boundary coefficients are collapsed. + * X must be either preallocated + * allocated or enlarged (reallocated) if xLen is too small. + * + * The returned value is the reconstructed signal length (=xLen). + */ { - int l, j; - - /* Shifting negated length leverages the fact that shift always rounds - * down; as a result, odd lengths are rounded *up* to account for extra - * low frequency coefficient (cf. Stefan Moebius' "TurboWavelets") */ - for (l = len; l >= 2; l = -(-l >> 1)) { - if (l >= 4) - d4ConvStep(x, l); - else - /* Apply Haar xform in final (and penultimate, if is len odd) - iteration to a obtain single smoothing an detail term */ - haarConvStep(x, l); + unsigned yOff = 0, i, j, l, nIter, numCoeffs; + double (*yOffPtr) [3] = NULL; + + /* Figure out number of iterations */ + for (nIter = 0, l = xLen; l > 3; nIter++, l = d4Step(NULL, l, NULL)); + +#if 1 + /* Invert Haar transform at two coarsest levels */ + for (l = 2; l <= 3; l++) { + haarInvStep(y, l, x); + /* Copy reconstructed signal into y for next iteration */ + memcpy(y, x, l * COEFSIZE); } -} - - - -void waveletInvXform (double x [][3], unsigned len) -/* Perform full multiresolution inverse wavelet transform - on array x of length len containing 3-tuples. Note no intra-tuple - transform occurs (1-D only). Returns reconstructed signal in x. */ -{ - int l, i, j, log2len; - - /* Get log2(len); reconstruct lengths used in forward xform, preserving - * odd lengths by signed shifting (which always rounds down) */ - for (log2len = 0, l = len; l > 1; log2len++, l = -(-l >> 1)); - - for (i = 1; i <= log2len; i++) { - l = -(-(int)len >> (log2len - i)); - if (l >= 4) - d4InvConvStep(x, l); - else - /* Inverse Haar xform in first (and second, if len is odd) - iteration to expand single smoothing and detail term */ - haarInvConvStep(x, l); +#endif + + /* Invert Daubechies D4 transform */ + for (i = nIter; i; i--) { + for ( + j = 0, yOff = yLen, l = xLen; + j < i; + j++, yOff -= l = d4Step(NULL, l, NULL) + ); + +#ifdef _WAVELET_DBG + /* Zero output array to facilitate debugging */ + memset(x, 0, xLen * COEFSIZE); +#endif + + yOffPtr = y + yOff; + numCoeffs = d4InvStep(yOffPtr - l, l, x); + + /* Copy approximation coeffs from x into y for next iteration */ + memcpy(yOffPtr + l - numCoeffs, x, numCoeffs * COEFSIZE); } + + return xLen; } #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]; + int i, j, xLen, yLen; + double (*xOrg) [3], (*x) [3], (*y) [3] = NULL; char errstr [ERRSTRLEN]; if (argc < 2) { fputs("Missing array length\n", stderr); return -1; } - if (!(len = atoi(argv [1]))) { + if (!(xLen = atoi(argv [1]))) { fputs("Invalid array length\n", stderr); return -1; } - /* Allocate arrays for original and reconstruction */ - if (!(x0 = calloc(len, COEFSIZE)) || !(x1 = calloc(len, COEFSIZE))) { + /* Allocate arrays for original and reconstruction (num coeffs is + * obtained from waveletXform() before allocating output array) */ + if ( + !(xOrg = calloc(xLen, COEFSIZE)) || + !(x = calloc(xLen, COEFSIZE)) || + !(y = calloc(waveletXform(NULL, xLen, NULL), 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++) { + for (i = 0; i < xLen; i++) { #if 1 - x0 [i][0] = x0 [i][1] = x0 [i][2] = drand48(); + xOrg [i][0] = xOrg [i][1] = xOrg [i][2] = drand48(); #else for (j = 0; j < 3; j++) - x0 [i][j] = drand48(); + xOrg [i][j] = drand48(); #endif } - qsort(x0, len, COEFSIZE, tupleCmp); - memcpy(x1, x0, (size_t)len * COEFSIZE); + qsort(xOrg, xLen, COEFSIZE, tupleCmp); + memcpy(x, xOrg, (size_t)xLen * COEFSIZE); /* Transform and reconstruct data */ - waveletXform(x1, len); + yLen = waveletXform(x, xLen, y); #ifdef _WAVELET_DBG puts("-----------------------------------------------------------\n"); #endif - waveletInvXform(x1, len); + waveletInvXform(y, yLen, x, xLen); #ifdef _WAVELET_DBG puts("-----------------------------------------------------------\n"); #endif #ifdef _WAVELET_DBG puts("ORIG\t vs.\tINV XFORM"); - for (i = 0; i < len; i++) + for (i = 0; i < xLen; 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] + xOrg [i][0], xOrg [i][1], xOrg [i][2], + x [i][0], x [i][1], x [i][2] ); #else - printf("%.6f\t%.6f\n", x0 [i][0], x1 [i][0]); + printf("%.3f\t%.3f\n", xOrg [i][0], x [i][0]); #endif #endif - free(x0); - free(x1); + free(xOrg); + free(x); + free(y); return 0; } #endif