diff --git a/wavelet2.c b/wavelet2.c index 22a3584..1aafa9e 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,605 +1,617 @@ /* ========================================================================= 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 COEFSIZE (sizeof(Coeff3)) #define min(a, b) ((a) < (b) ? (a) : (b)) #define coeffAvg(a) (((a).coeff [0] + (a).coeff [1] + (a).coeff [2]) / 3) /* Daubechies D4 wavelet coeffs */ static const Coeff h4Norm = 0.25 / sqrt(2), 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] }, /* Summed coeffs for boundary extension */ h4Bound = h4 [2] + h4 [3], g4Bound = g4 [2] + g4 [3], h4InvBound [2] = { h4 [2] + h4 [0], h4 [3] + h4 [1] }, g4InvBound [2] = { g4 [2] + g4 [0], g4 [3] + g4 [1] }; /* Haar wavelet coeffs */ static const Coeff h2 = 1 / sqrt(2); #ifdef WAVELET_DBG static void zeroCoeffs (CoeffArray2 y, unsigned l) /* Zero output array to facilitate debugging */ { const unsigned len = 1 << l; unsigned i; for (i = 0; i < len; i++) memset(y [i], 0, len * COEFSIZE); } static void dumpCoeffs (const CoeffArray2 y, const CoeffArray2 yt, unsigned l ) { 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])); printf("-->\t"); for (j = 0; j < len; j++) printf("%4.3f\t", coeffAvg(yt [i][j])); putchar('\n'); } } #endif static int haarStep (CoeffArray2 y, CoeffArray2 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 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. */ { - const unsigned len = 1 << l, hlen = 1 << (l - 1); - unsigned h, i, j, k; - Coeff3 *s [2], *a, *d; + static unsigned axis = 0; + const unsigned len = 1 << l, hlen = 1 << (l - 1); + unsigned h, i, j, k; + Coeff3 *s [2], *a, *d; if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif for (i = 0; i < len; i++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* yt transposed on the fly by swapping indices in assignment */ /* Smooth/approx/avg/lowpass */ yt [h ][i].coeff [k] = h2 * (y [i][j].coeff [k] + y [i][j + 1].coeff [k]); /* Detail/diff/highpass */ yt [hlen + h][i].coeff [k] = h2 * (y [i][j].coeff [k] - y [i][j + 1].coeff [k]); } } } #ifdef WAVELET_DBG - printf("HAAR XFORM (%d x %d)\n", len, len); + 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 (CoeffArray2 y, CoeffArray2 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. */ { - const unsigned len = 1 << l, hlen = 1 << (l - 1); - unsigned h, i, j, k; - Coeff *s [2], *a, *d; + static unsigned axis = 0; + const unsigned len = 1 << l, hlen = 1 << (l - 1); + unsigned h, i, j, k; + Coeff *s [2], *a, *d; if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif for (i = 0; i < len; i++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* yt transposed on the fly by swapping indices in assignment */ yt [j ][i].coeff [k] = h2 * (y [i][h].coeff [k] + y [i][hlen + h].coeff [k]); yt [j + 1][i].coeff [k] = h2 * (y [i][h].coeff [k] - y [i][hlen + h].coeff [k]); } } } - + #ifdef WAVELET_DBG - printf("INV HAAR XFORM (%d x %d)\n", len, len); + printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); -#endif + putchar('\n'); + axis ^= 1; +#endif return 0; } static int d4Step (CoeffArray2 y, CoeffArray2 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 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. */ { - const unsigned len = 1 << l, hlen = 1 << (l - 1); - unsigned h, i, j, k; - Coeff *s [4], *a, *d; + static unsigned axis = 0; + const unsigned len = 1 << l, hlen = 1 << (l - 1); + unsigned h, i, j, k; + Coeff *s [4], *a, *d; if (len < 4 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif for (i = 0; i < len; i++) { /* Transform up to boundary */ for (j = 0; j < len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* yt transposed on the fly by swapping indices in assignment */ /* Smooth/approx/avg/lowpass */ 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 [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 boundary with constant extension at y [len-1]. Note j is still set to last index from previous loop */ h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ][i].coeff [k] = h4 [0] * y [i][j ].coeff [k] + h4 [1] * y [i][j + 1 ].coeff [k] + h4Bound * y [i][len - 1].coeff [k]; /* Detail/diff/highpass */ yt [hlen + h][i].coeff [k] = g4 [0] * y [i][j ].coeff [k] + g4 [1] * y [i][j + 1 ].coeff [k] + g4Bound * y [i][len - 1].coeff [k]; } } - + #ifdef WAVELET_DBG - printf("D4 XFORM (%d x %d)\n", len, len); + 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 (CoeffArray2 y, CoeffArray2 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). 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. */ { - const unsigned len = 1 << l, hlen = 1 << (l - 1); - unsigned h, i, j, k; - Coeff *s [2], *a, *d; + static unsigned axis = 0; + const unsigned len = 1 << l, hlen = 1 << (l - 1); + unsigned h, i, j, k; + Coeff *s [2], *a, *d; if (len < 4 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG zeroCoeffs(yt, l); #endif for (i = 0; i < len; i++) { /* Invert transform up to boundary */ for (j = 0; j < len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { yt [j ][i].coeff [k] = h4 [2] * y [i][h ].coeff [k] + g4 [2] * y [i][hlen + h ].coeff [k] + h4 [0] * y [i][h + 1 ].coeff [k] + g4 [0] * y [i][hlen + h + 1].coeff [k]; yt [j + 1][i].coeff [k] = h4 [3] * y [i][h ].coeff [k] + g4 [3] * y [i][hlen + h ].coeff [k] + h4 [1] * y [i][h + 1 ].coeff [k] + g4 [1] * y [i][hlen + h + 1].coeff [k]; } } /* Invert transform at boundary with constant extension at y [len/2-1] and y [len-1]. Note j is still set to last index from previous loop */ h = j >> 1; for (k = 0; k < 3; k++) { yt [j ][i].coeff [k] = h4InvBound [0] * y [i][h ].coeff [k] + g4InvBound [0] * y [i][hlen + h].coeff [k]; yt [j + 1][i].coeff [k] = h4InvBound [1] * y [i][h ].coeff [k] + g4InvBound [1] * y [i][hlen + h].coeff [k]; } } - + #ifdef WAVELET_DBG - printf("INV D4 XFORM (%d x %d)\n", len, len); + printf("%s INV D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs(y, yt, l); -#endif + putchar('\n'); + axis ^= 1; +#endif return 0; } CoeffArray2 waveletAllocCoeffArray2 (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 len = 1 << l; unsigned i; CoeffArray2 y = NULL; if (l >= 1) { if (!(y = calloc(len, sizeof(Coeff3*)))) return NULL; for (i = 0; i < len; i++) if (!(y [i] = calloc(len, COEFSIZE))) return NULL; } return y; } void waveletFreeCoeffArray2 (CoeffArray2 y, unsigned l) /* Free previously allocated 2D coefficient array y of size (2^l) x (2^l) */ { unsigned i, j; const unsigned len = 1 << l; if (y) { for (i = 0; i < len; i++) free(y [i]); free(y); } } int waveletXform2 (CoeffArray2 y, CoeffArray2 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][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. 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; CoeffArray2 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 = waveletAllocCoeffArray2(l))) { fprintf(stderr, "ERROR - Failed allocating %dx%d buffer array" " in WaveletXform2()", len, len); return -1; } for (li = l; li > 1; li--) { /* Apply horizontal & vertical Daubechies D4 transform, swapping input and transposed output array */ if (d4Step(y, yt, li) || d4Step(yt, y, li)) return -1; } /* Apply horizontal & vertical Haar transform at coarsest resolution (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 */ waveletFreeCoeffArray2(ytloc, l); return 0; } int waveletInvXform2 (CoeffArray2 y, CoeffArray2 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. 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; CoeffArray2 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 = waveletAllocCoeffArray2(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++) { /* Invert horizontal & vertical Daubechies D4 transform, swapping input and transposed output arrays */ if (d4InvStep(y, yt, li) || d4InvStep(yt, y, li)) return -1; } /* NOTE: Reconstructed signal now in y */ if (ytloc) /* Free yt if allocated on demand */ waveletFreeCoeffArray2(ytloc, l); return 0; } #ifdef WAVELET_TEST #include #define ERRSTRLEN 1024 int main (int argc, char *argv []) { int i, j, k, l, *xIdx; unsigned len; CoeffArray2 y0 = NULL, y = NULL; char errstr [ERRSTRLEN]; if (argc < 2) { 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; } /* Allocate arrays for original and reconstruction */ len = 1 << l; if (!(y0 = waveletAllocCoeffArray2(len)) || !(y = waveletAllocCoeffArray2(len)) ) { fprintf(stderr, "Failed allocating %dx%d array\n", len, len); return -1; } /* Init with random data */ srand48(111); for (i = 0; i < len; i++) for (j = 0; j < len; j++) #if 1 for (k = 0; k < 3; k++) y0 [i][j].coeff [k] = y [i][j].coeff [k] = drand48(); #else 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"); for (i = 0; i < len; i++) { for (j = 0; j < len; j++) printf("%4.3f\t", coeffAvg(y [i][j])); putchar('\n'); } puts("\n-----------------------------------------------------------\n"); #endif /* Inverse xform */ if (waveletInvXform2(y, NULL, l)) { fputs("Inverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("-----------------------------------------------------------\n"); #endif puts("ORIG\tvs.\tINV XFORM"); for (i = 0; i < len; i++) { for (j = 0; j < len; j++) printf("%4.3f\t", coeffAvg(y0 [i][j])); printf("\t\t"); for (j = 0; j < len; j++) printf("%4.3f\t", coeffAvg(y [i][j])); putchar('\n'); } waveletFreeCoeffArray2(y0, l); waveletFreeCoeffArray2(y, l); return 0; } #endif