diff --git a/wavelet2.c b/wavelet2.c index 2fd46e8..152260f 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,670 +1,668 @@ /* ========================================================================= 2-dimensional wavelet transform on (2^l) x (2^l) sized arrays of 3-tuples, where l > 1. Compile with -DWAVELET_TEST to build standalone unit tests. 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 "wavelet2.h" #include #include #include #include #define min(a, b) ((a) < (b) ? (a) : (b)) -#define coeffAvg(a) (0.333333 * ((a).coeff [0] + (a).coeff [1] + (a).coeff [2])) +#define coeffAvg(a) (((a).coeff [0] + (a).coeff [1] + (a).coeff [2]) / 3) + /* Haar wavelet coeffs */ static const WAVELET_COEFF h2 = 1 / sqrt(2); /* Daubechies D4 wavelet coeffs */ static const WAVELET_COEFF h4Norm = 0.25 * h2, sqrt3 = sqrt(3), h4 [4] = { h4Norm * (1 + sqrt3), h4Norm * (3 + sqrt3), h4Norm * (3 - sqrt3), h4Norm * (1 - sqrt3) }, g4 [4] = { h4 [3], -h4 [2], h4 [1], -h4 [0] }; -static void zeroCoeffs (WaveletMatrix y, unsigned l) -/* Zero output array to facilitate debugging */ +WaveletMatrix allocWaveletMatrix (unsigned l) +/* + Allocate and return a 2D coefficient array of size (2^l) x (2^l), + where l >= 1. Returns NULL if allocation failed. +*/ { - const unsigned lenSq = 1 << (l << 1); - memset(y, 0, lenSq * COEFFSIZE); + const unsigned len = 1 << l; + unsigned i; + WaveletMatrix y = NULL; + + if (l >= 1) { + if (!(y = calloc(len, sizeof(WaveletCoeff3*)))) + return NULL; + + for (i = 0; i < len; i++) + if (!(y [i] = calloc(len, WAVELET_COEFFSIZE))) + return NULL; + } + + return y; } - -static void dumpCoeffs (const WaveletMatrix y, const WaveletMatrix yt, - unsigned l -) -/* Dump arrays y and yt side-by-side to stdout (skip yt if NULL) */ + +void freeWaveletMatrix (WaveletMatrix y, unsigned l) +/* + Free previously allocated 2D coefficient array y of size (2^l) x (2^l) +*/ { - const unsigned len = 1 << l; unsigned i, j; + const unsigned len = 1 << l; + + if (y) { + for (i = 0; i < len; i++) + free(y [i]); + + free(y); + } +} - for (i = 0; i < len; i++) { - for (j = 0; j < len; j++) - printf("%4.3f\t", coeffAvg(y [WAVELET_IDX(i, j, l)])); + + +#ifdef WAVELET_DBG + static void zeroCoeffs (WaveletMatrix y, unsigned l) + /* Zero output array to facilitate debugging */ + { + const unsigned len = 1 << l; + unsigned i; - if (yt) { - printf("--->\t"); - for (j = 0; j < len; j++) - printf("%4.3f\t", coeffAvg(yt [WAVELET_IDX(i, j, l)])); - } + for (i = 0; i < len; i++) + memset(y [i], 0, len * WAVELET_COEFFSIZE); + } + + + + static void dumpCoeffs (const WaveletMatrix y, const WaveletMatrix yt, + unsigned l + ) + /* Dump arrays y and yt side-by-side to stdout (skip yt if NULL) */ + { + const unsigned len = 1 << l; + unsigned i, j, k; + + for (i = 0; i < len; i++) { + for (j = 0; j < len; j++) + printf("%4.3f\t", coeffAvg(y [i][j])); - putchar('\n'); + if (yt) { + printf("-->\t"); + for (j = 0; j < len; j++) + printf("%4.3f\t", coeffAvg(yt [i][j])); + } + + putchar('\n'); + } } -} +#endif static int haarStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of forward 2D Haar wavelet transform on array y of size - (2^l)x(2^l) containing 3-tuples, where l >= 1. The transform is + (2^l) x (2^l) containing 3-tuples, where l >= 1. The transform is performed over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). Note the triplets per array element are _not_ decorrelated, but transformed independently. The wavelet coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the generated coefficients vertically, and prepares the array for another transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been decorrelated. Returns 0 on success, else -1. */ { static unsigned axis = 0; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: yt is transposed on the fly such that the next function call * transforms over the alternate axis. This is done by simply swapping * the indices during assignment */ - WaveletCoeff3 *y0 = y, *y1 = y + 1; - for (i = 0; i < len; i++) { -#if 0 for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ - yt [WAVELET_IDX(h, i, l)].coeff [k] = h2 * ( - y [WAVELET_IDX(i, j, l)].coeff [k] + - y [WAVELET_IDX(i, j + 1, l)].coeff [k] + yt [h ][i].coeff [k] = h2 * ( + y [i][j ].coeff [k] + + y [i][j + 1].coeff [k] ); /* Detail/diff/highpass */ - yt [WAVELET_IDX(hlen + h, i, l)].coeff [k] = h2 * ( - y [WAVELET_IDX(i, j, l)].coeff [k] - - y [WAVELET_IDX(i, j + 1, l)].coeff [k] - ); - } - } -#else - /* Optimised access through pointers to minimise indexing */ - WaveletCoeff3 *yt0 = yt + WAVELET_IDX(0, i, l), - *yt1 = yt + WAVELET_IDX(hlen, i, l); - - for (j = 0; j < len; - j += 2, yt0 += len, yt1 += len, y0 += 2, y1 += 2 - ) { - h = j >> 1; - for (k = 0; k < 3; k++) { - /* Smooth/approx/avg/lowpass */ - yt0 -> coeff [k] = h2 * (y0 -> coeff [k] + y1 -> coeff [k]); - - /* Detail/diff/highpass */ - yt1 -> coeff [k] = h2 * (y0 -> coeff [k] - y1 -> coeff [k]); + yt [hlen + h][i].coeff [k] = h2 * ( + y [i][j ].coeff [k] - + y [i][j + 1].coeff [k] + ); } } -#endif } #ifdef WAVELET_DBG printf("%s FWD HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); putchar('\n'); axis ^= 1; #endif return 0; } static int haarInvStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of inverse 2D Haar wavelet transform on coefficient array y of size (2^l) x (2^l) containing 3-tuples, where l >= 1. This reverses the forward transform above. The transform is inverted over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). The inverted coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the inverted coefficients vertically, and prepares the array for another inverse transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent inverse transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been inversely transformed. Returns 0 on success, else -1. */ { static unsigned axis = 1; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: i, j are swapped relative to the forward transform, as axis * order is now reversed. */ /* NOTE: yt is transposed on the fly such that the next function call * inverts over the alternate axis. This is done by simply swapping * the indices during assignment */ - WaveletCoeff3 *yt0 = yt, *yt1 = yt + 1; - for (i = 0; i < len; i++) { -#if 0 for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { - yt [WAVELET_IDX(i, j, l)].coeff [k] = h2 * ( - y [WAVELET_IDX(h, i, l)].coeff [k] + /* Avg */ - y [WAVELET_IDX(hlen + h, i, l)].coeff [k] /* Diff */ + yt [i][j ].coeff [k] = h2 * ( + y [h ][i].coeff [k] + /* Avg */ + y [hlen + h][i].coeff [k] /* Diff */ ); - yt [WAVELET_IDX(i, j + 1, l)].coeff [k] = h2 * ( - y [WAVELET_IDX(h, i, l)].coeff [k] - /* Avg */ - y [WAVELET_IDX(hlen + h, i, l)].coeff [k] /* Diff */ + yt [i][j + 1].coeff [k] = h2 * ( + y [h ][i].coeff [k] - /* Avg */ + y [hlen + h][i].coeff [k] /* Diff */ ); } } -#else - /* Optimised access through pointers to minimise indexing */ - WaveletCoeff3 *y0 = y + WAVELET_IDX(0, i, l), - *y1 = y + WAVELET_IDX(hlen, i, l); - - for (j = 0; j < len; - j += 2, yt0 += 2, yt1 += 2, y0 += len, y1 += len - ) { - h = j >> 1; - for (k = 0; k < 3; k++) { - yt0 -> coeff [k] = h2 * ( - y0 -> coeff [k] + /* Avg */ - y1 -> coeff [k] /* Diff */ - ); - - yt1 -> coeff [k] = h2 * ( - y0 -> coeff [k] - /* Avg */ - y1 -> coeff [k] /* Diff */ - ); - } - } -#endif } #ifdef WAVELET_DBG printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); putchar('\n'); axis ^= 1; #endif return 0; } static int d4Step (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Single step of forward 2D Daubechies D4 wavelet transform on array y of - size (2^l)x(2^l) containing 3-tuples, where l >= 2. The transform is + size (2^l) x (2^l) containing 3-tuples, where l >= 2. The transform is performed over y's 2nd axis (i.e. horizontally, assuming row-major addressing as per C convention). Note the triplets per array element are _not_ decorrelated, but transformed independently. The wavelet coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the generated coefficients vertically, and prepares the array for another transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been decorrelated. Returns 0 on success, else -1. */ { static unsigned axis = 0; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 4 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: yt is transposed on the fly such that the next function call * transforms over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { /* Transform until upper boundary */ for (j = 0; j < len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ - yt [WAVELET_IDX(h, i, l)].coeff [k] = - h4 [0] * y [WAVELET_IDX(i, j, l)].coeff [k] + - h4 [1] * y [WAVELET_IDX(i, j + 1, l)].coeff [k] + - h4 [2] * y [WAVELET_IDX(i, j + 2, l)].coeff [k] + - h4 [3] * y [WAVELET_IDX(i, j + 3, l)].coeff [k]; + yt [h ][i].coeff [k] = + h4 [0] * y [i][j ].coeff [k] + + h4 [1] * y [i][j + 1].coeff [k] + + h4 [2] * y [i][j + 2].coeff [k] + + h4 [3] * y [i][j + 3].coeff [k]; /* Detail/diff/highpass */ - yt [WAVELET_IDX(hlen + h, i, l)].coeff [k] = - g4 [0] * y [WAVELET_IDX(i, j, l)].coeff [k] + - g4 [1] * y [WAVELET_IDX(i, j + 1, l)].coeff [k] + - g4 [2] * y [WAVELET_IDX(i, j + 2, l)].coeff [k] + - g4 [3] * y [WAVELET_IDX(i, j + 3, l)].coeff [k]; + yt [hlen + h][i].coeff [k] = + g4 [0] * y [i][j ].coeff [k] + + g4 [1] * y [i][j + 1].coeff [k] + + g4 [2] * y [i][j + 2].coeff [k] + + g4 [3] * y [i][j + 3].coeff [k]; } } /* Transform at upper boundary with wraparound. Note j is set to last index from previous loop */ h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ - yt [WAVELET_IDX(h, i, l)].coeff [k] = - h4 [0] * y [WAVELET_IDX(i, j, l)].coeff [k] + - h4 [1] * y [WAVELET_IDX(i, j + 1, l)].coeff [k] + - h4 [2] * y [WAVELET_IDX(i, 0, l)].coeff [k] + - h4 [3] * y [WAVELET_IDX(i, 1, l)].coeff [k]; + yt [h ][i].coeff [k] = + h4 [0] * y [i][j ].coeff [k] + + h4 [1] * y [i][j + 1].coeff [k] + + h4 [2] * y [i][0 ].coeff [k] + + h4 [3] * y [i][1 ].coeff [k]; /* Detail/diff/highpass */ - yt [WAVELET_IDX(hlen + h, i, l)].coeff [k] = - g4 [0] * y [WAVELET_IDX(i, j, l)].coeff [k] + - g4 [1] * y [WAVELET_IDX(i, j + 1, l)].coeff [k] + - g4 [2] * y [WAVELET_IDX(i, 0, l)].coeff [k] + - g4 [3] * y [WAVELET_IDX(i, 1, l)].coeff [k]; + yt [hlen + h][i].coeff [k] = + g4 [0] * y [i][j ].coeff [k] + + g4 [1] * y [i][j + 1].coeff [k] + + g4 [2] * y [i][0 ].coeff [k] + + g4 [3] * y [i][1 ].coeff [k]; } } #ifdef WAVELET_DBG printf("%s FWD D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); putchar('\n'); axis ^= 1; #endif return 0; } static int d4InvStep (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* - Single step of inverse 2D Daubechies D4 wavelet transform on coefficient array y - of size (2^l) x (2^l) containing 3-tuples, where l >= 2. This reverses - the forward transform above. The transform is inverted over y's 2nd axis - (i.e. horizontally, assuming row-major addressing as per C convention). + Single step of inverse 2D Daubechies D4 wavelet transform on coefficient + array y of size (2^l) x (2^l) containing 3-tuples, where l >= 2. This + reverses the forward transform above. The transform is inverted over y's + 2nd axis (i.e. horizontally, assuming row-major addressing as per C + convention). The inverted coefficients are returned in the *TRANSPOSED* output array yt, of identical dimensions to y. The transpose arranges the inverted coefficients vertically, and prepares the array for another inverse transform over its 2nd axis in a subsequent call (this time as input y). The result of a subsequent inverse transform restores yt to y's original orientation, in which case both horizontal and vertical axes have been inversely transformed. Returns 0 on success, else -1. */ { static unsigned axis = 1; const unsigned len = 1 << l, hlen = 1 << (l - 1); unsigned h, i, j, k; if (len < 4 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif /* NOTE: i, j are swapped relative to the forward transform, as axis * order is now reversed. */ /* NOTE: yt is transposed on the fly such that the next function call * inverts over the alternate axis. This is done by simply swapping * the indices during assignment */ for (i = 0; i < len; i++) { /* Invert at lower boundary with wraparound */ for (k = 0; k < 3; k++) { - yt [WAVELET_IDX(i, 0, l)].coeff [k] = - h4 [2] * y [WAVELET_IDX(hlen - 1, i, l)].coeff [k] + /* Last avg */ - g4 [2] * y [WAVELET_IDX(len - 1, i, l)].coeff [k] + /* Last diff */ - h4 [0] * y [WAVELET_IDX(0, i, l)].coeff [k] + /* First avg */ - g4 [0] * y [WAVELET_IDX(hlen, i, l)].coeff [k]; /* First diff */ + yt [i][0].coeff [k] = + h4 [2] * y [hlen - 1][i].coeff [k] + /* Last avg */ + g4 [2] * y [len - 1][i].coeff [k] + /* Last diff */ + h4 [0] * y [0 ][i].coeff [k] + /* First avg */ + g4 [0] * y [hlen ][i].coeff [k]; /* First diff */ - yt [WAVELET_IDX(i, 1, l)].coeff [k] = - h4 [3] * y [WAVELET_IDX(hlen - 1, i, l)].coeff [k] + - g4 [3] * y [WAVELET_IDX(len - 1, i, l)].coeff [k] + - h4 [1] * y [WAVELET_IDX(0, i, l)].coeff [k] + - g4 [1] * y [WAVELET_IDX(hlen, i, l)].coeff [k]; + yt [i][1].coeff [k] = + h4 [3] * y [hlen - 1][i].coeff [k] + + g4 [3] * y [len - 1][i].coeff [k] + + h4 [1] * y [0 ][i].coeff [k] + + g4 [1] * y [hlen ][i].coeff [k]; } /* Invert until upper boundary */ for (j = 2; j < len; j += 2) { h = (j >> 1) - 1; for (k = 0; k < 3; k++) { - yt [WAVELET_IDX(i, j, l)].coeff [k] = - h4 [2] * y [WAVELET_IDX(h, i, l)].coeff [k] + /* Avg */ - g4 [2] * y [WAVELET_IDX(hlen + h, i, l)].coeff [k] + /* Diff */ - h4 [0] * y [WAVELET_IDX(h + 1, i, l)].coeff [k] + /* Next avg */ - g4 [0] * y [WAVELET_IDX(hlen + h + 1, i, l)].coeff [k]; /* Next diff */ + yt [i][j ].coeff [k] = + h4 [2] * y [h ][i].coeff [k] + /* Avg */ + g4 [2] * y [hlen + h ][i].coeff [k] + /* Diff */ + h4 [0] * y [h + 1 ][i].coeff [k] + /* Next avg */ + g4 [0] * y [hlen + h + 1][i].coeff [k]; /* Next diff */ - yt [WAVELET_IDX(i, j + 1, l)].coeff [k] = - h4 [3] * y [WAVELET_IDX(h, i, l)].coeff [k] + - g4 [3] * y [WAVELET_IDX(hlen + h, i, l)].coeff [k] + - h4 [1] * y [WAVELET_IDX(h + 1, i, l)].coeff [k] + - g4 [1] * y [WAVELET_IDX(hlen + h + 1, i, l)].coeff [k]; + yt [i][j + 1].coeff [k] = + h4 [3] * y [h ][i].coeff [k] + + g4 [3] * y [hlen + h ][i].coeff [k] + + h4 [1] * y [h + 1 ][i].coeff [k] + + g4 [1] * y [hlen + h + 1][i].coeff [k]; } } } #ifdef WAVELET_DBG printf("%s INV D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); putchar('\n'); axis ^= 1; #endif return 0; } int waveletXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Perform full 2D multiresolution forward wavelet transform on array y of size (2^l) x (2^l) containing original signal as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. The wavelet coefficients are returned in array y, containing the coarsest - approximation in y [0] followed by horizontal/vertical details in + approximation in y [0][0] followed by horizontal/vertical details in order of increasing resolution/frequency. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum - expected size. - WARNING: The dimensions of yt are not checked; this is the caller's - responsibility. + expected size. The dimensions of yt are not checked; this is the + caller's responsibility. Returns 0 on success, else -1. */ { const unsigned len = 1 << l; unsigned li; WaveletMatrix ytloc = NULL; /* Skip transform if input too short or missing */ if (l < 1 || !y) return -1; if (!yt) /* No buffer supplied; allocate one on demand */ - if (!(yt = ytloc = WAVELET_ALLOC(l))) { + if (!(yt = ytloc = allocWaveletMatrix(l))) { fprintf(stderr, "ERROR - Failed allocating %dx%d buffer array" " in WaveletXform2()", len, len); return -1; } - + for (li = l; li > 1; li--) { #if 0 /* Apply horizontal & vertical Daubechies D4 transform, swapping input and transposed output array */ if (d4Step(y, yt, li) || d4Step(yt, y, li)) #else if (haarStep(y, yt, li) || haarStep(yt, y, li)) #endif return -1; } /* Apply horizontal & vertical Haar transform at coarsest resolution - (li==1) to obtain single approximation coefficient at y [0]; all + (li==1) to obtain single approximation coefficient at y [0][0]; all other coeffs are details. */ if (haarStep(y, yt, li) || haarStep(yt, y, li)) return -1; /* NOTE: All coefficients now in y */ if (ytloc) /* Free yt if allocated on demand */ - free(ytloc); + freeWaveletMatrix(ytloc, l); return 0; } int waveletInvXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l) /* Perform full 2D multiresolution inverse wavelet transform on array y of size (2^l) x (2^l) containing wavelet coefficients as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum - expected size. - WARNING: The dimensions of yt are not checked; this is the caller's - responsibility. + expected size. The dimensions of yt are not checked; this is the + caller's responsibility. The reconstructed signal is returned in array y. Returns 0 on success, else -1. */ { const unsigned len = 1 << l; unsigned li; WaveletMatrix ytloc = NULL; /* Skip inverse transform if input too short or missing */ if (l < 1 || !y) return -1; if (!yt) /* No buffer supplied; allocate one on demand */ - if (!(yt = ytloc = WAVELET_ALLOC(l))) { + if (!(yt = ytloc = allocWaveletMatrix(l))) { fprintf(stderr, "ERROR - Failed allocating %dx%d buffer array" " in WaveletInvXform2()", len, len); return -1; } - /* Invert horizontal & vertical Haar transform at coarsest level (li==1), swapping input and transposed output array */ if (haarInvStep(y, yt, 1) || haarInvStep(yt, y, 1)) return -1; for (li = 2; li <= l; li++) { #if 0 /* Invert horizontal & vertical Daubechies D4 transform, swapping input and transposed output arrays */ if (d4InvStep(y, yt, li) || d4InvStep(yt, y, li)) #else if (haarInvStep(y, yt, li) || haarInvStep(yt, y, li)) #endif return -1; } /* NOTE: Reconstructed signal now in y */ if (ytloc) /* Free yt if allocated on demand */ - free(ytloc); + freeWaveletMatrix(ytloc, l); return 0; } #ifdef WAVELET_TEST #include #define ERRSTRLEN 1024 int main (int argc, char *argv []) { int i, j, k, l; - unsigned len, idx; + unsigned len; WaveletMatrix y0 = NULL, y = NULL; char errstr [ERRSTRLEN]; FILE *dataFile = NULL; float inData; if (argc < 2) { fprintf(stderr, "%s [dataFile]\n", argv [0]); fputs("Missing array resolution l > 1\n", stderr); return -1; } if (!(l = atoi(argv [1])) || l < 1) { fputs("Invalid array resolution l\n", stderr); return -1; } else len = 1 << l; /* Allocate arrays for original and reconstruction */ - if (!(y0 = WAVELET_ALLOC(l)) || !(y = WAVELET_ALLOC(l))) { + if (!(y0 = allocWaveletMatrix(l)) || !(y = allocWaveletMatrix(l))) { fprintf(stderr, "Failed allocating %dx%d array\n", len, len); return -1; } if (argc > 2) { /* Load data from file; length must not exceed allocated */ if (!(dataFile = fopen(argv [2], "r"))) { fprintf(stderr, "Failed opening data file %s\n", argv [2]); return -1; } for (i = 0; i < len; i++) { for (j = 0; j < len; j++) { if (feof(dataFile)) { fprintf(stderr, "Premature end of file reading data from %s\n", argv [2] ); fclose(dataFile); return -1; } - + /* Read next float, skipping any leading whitespace */ if (fscanf(dataFile, " %f", &inData)) { - idx = WAVELET_IDX(i, j, l); - - y0 [idx].coeff [0] = y0 [idx].coeff [1] = y0 [idx].coeff [2] = - y [idx].coeff [0] = y [idx].coeff [1] = y [idx].coeff [2] = + y0 [i][j].coeff [0] = y0 [i][j].coeff [1] = + y0 [i][j].coeff [2] = y [i][j].coeff [0] = + y [i][j].coeff [1] = y [i][j].coeff [2] = (WAVELET_COEFF)inData; } else { fprintf(stderr, "Error reading from data file %s\n", argv [2] ); fclose(dataFile); return -1; } } } fclose(dataFile); } else { /* Init with random data */ srand48(111); for (i = 0; i < len; i++) - for (j = 0; j < len; j++) { - idx = WAVELET_IDX(i, j, l); + for (j = 0; j < len; j++) #if 1 for (k = 0; k < 3; k++) - y0 [idx].coeff [k] = y [idx].coeff [k] = drand48(); + y0 [i][j].coeff [k] = y [i][j].coeff [k] = drand48(); #else - y0 [idx].coeff [0] = y0 [idx].coeff [1] = y0 [idx].coeff [2] = - y [idx].coeff [0] = y [idx].coeff [1] = y [idx].coeff [2] = drand48(); + y0 [i][j].coeff [0] = y0 [i][j].coeff [1] = + y0 [i][j].coeff [2] = y [i][j].coeff [0] = + y [i][j].coeff [1] = y [i][j].coeff [2] = drand48(); #endif - } } /* Forward xform */ if (waveletXform2(y, NULL, l)) { fputs("Forward xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG /* Dump coefficients */ puts("-----------------------------------------------------------\n"); dumpCoeffs(y, NULL, l); puts("\n-----------------------------------------------------------\n"); #endif /* Inverse xform */ if (waveletInvXform2(y, NULL, l)) { fputs("Inverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("-----------------------------------------------------------\n"); puts("ORIG vs. INV XFORM"); dumpCoeffs(y0, y, l); #endif - free(y0); - free(y); + freeWaveletMatrix(y0, l); + freeWaveletMatrix(y, l); return 0; } #endif diff --git a/wavelet2.h b/wavelet2.h index 5494e2e..b413bed 100644 --- a/wavelet2.h +++ b/wavelet2.h @@ -1,89 +1,94 @@ /* ========================================================================= 2-dimensional wavelet transform on (2^l) x (2^l) sized arrays of 3-tuples, where l >= 1. 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 _WAVELET2_H #define _WAVELET2_H - /* Wavelet coefficient type; defaults to double if not already defined - * (may be overridden in compiler command line, for example). - * NOTE: single precision float may improve cache performance during - * 2D transform with large arrays. */ + + /* Wavelet coefficient type; defaults to double if not already defined. + (may be overridden in compiler command line, for example). + NOTE: single precision float may improve cache performance during + 2D transform with large arrays. */ #ifndef WAVELET_COEFF - #define WAVELET_COEFF double + #define WAVELET_COEFF double #endif - /* Wavelet coefficient matrices defs. These are stored as flattened 2D - * arrays to ensure they are allocated in a contiguous block. This - * simplifies interfacing to calling routines, and improves cache - * performance through localised access. */ + + /* Wavelet matrix defs; these are stored as arrays of pointers (a.k.a. + * Iliffe vectors), as this proved to be more performant than a flattened + * representation. Encapsulating the coefficient's triplets in a struct + * prevents the compiler from introducing another layer of indirection. + * */ typedef struct { - WAVELET_COEFF coeff [3]; + WAVELET_COEFF coeff [3]; } WaveletCoeff3; - typedef WaveletCoeff3 *WaveletMatrix; - #define COEFFSIZE (sizeof(WaveletCoeff3)) - - /* Allocate wavelet coefficient matrix of size 2^l per dimension. - NOTE: Matrices are stored as contiguous flattened 2D arrays for cache - optimisation */ - #define WAVELET_ALLOC(l) calloc(1L << (l << 1), COEFFSIZE) + #define WAVELET_COEFFSIZE (sizeof(WaveletCoeff3)) + typedef WaveletCoeff3 **WaveletMatrix; + + + + WaveletMatrix allocWaveletMatrix (unsigned l); + /* + Allocate and return a 2D coefficient array of size (2^l) x (2^l), + where l >= 1. Returns NULL if allocation failed. + */ - - /* Map 2D matrix indices to linear 1D index for size 2^l per dimension */ - #define WAVELET_IDX(y, x, l) (((y) << l) + (x)) - /* Map linear 1D index to 2D matrix indices for size 2^l per dimension */ - #define WAVELET_Y(i, l) ((i) >> l) - #define WAVELET_X(i, l) ((i) & (1 << l) - 1) + void freeWaveletMatrix (WaveletMatrix y, unsigned l); + /* + Free previously allocated 2D coefficient array y of size (2^l) x (2^l) + */ + int waveletXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l); - /* + /* Perform full 2D multiresolution forward wavelet transform on array y of size (2^l) x (2^l) containing original signal as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. The wavelet coefficients are returned in array y, containing the - coarsest approximation in y [0] followed by horizontal/vertical + coarsest approximation in y [0][0] followed by horizontal/vertical details in order of increasing resolution/frequency. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. Ideally, yt should be preallocated to the maximum - expected size. - WARNING: The dimensions of a preallocated yt are not checked; this is - the caller's responsibility. + expected size. The dimensions of a preallocated yt are not checked; + this is the caller's responsibility. Returns 0 on success, else -1. */ + int waveletInvXform2 (WaveletMatrix y, WaveletMatrix yt, unsigned l); /* Perform full 2D multiresolution inverse wavelet transform on array y of size (2^l) x (2^l) containing wavelet coefficients as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. Ideally, yt should be preallocated to the maximum expected size. The dimensions of a preallocated yt are not checked; this is the caller's responsibility. The reconstructed signal is returned in array y. Returns 0 on success, else -1. */ #endif