diff --git a/wavelet2.c b/wavelet2.c index 8978960..ec20343 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,834 +1,856 @@ /* ========================================================================= 2-dimensional Daubechies D4 and Haar wavelet transforms on (2^l) x (2^l) sized arrays of 3-tuples, where l > 1. Array boundaries are extended as periodic, which does not require padding coefficients. Compile with -DWAVELET_TEST_2D 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 /* The following defs are const, but strict compilers pretend to be dumber than they are, and refuse to init from a func (in this case sqrt(2) and sqrt(3)), or other consts */ /* Haar wavelet coeffs */ const WAVELET_COEFF h2 = 1 / SQRT2; /* Daubechies D4 wavelet coeffs */ const WAVELET_COEFF h4 [4] = { H4NORM * (1 + SQRT3), H4NORM * (3 + SQRT3), H4NORM * (3 - SQRT3), H4NORM * (1 - SQRT3) }; const WAVELET_COEFF g4 [4] = { H4NORM * (1 - SQRT3), -H4NORM * (3 - SQRT3), H4NORM * (3 + SQRT3), -H4NORM * (1 + SQRT3) }; /* g4 [4] = { h4 [3], -h4 [2], h4 [1], -h4 [0] }; */ WaveletMatrix2 allocWaveletMatrix2 (unsigned len) /* Allocate and return a 2D coefficient array of size (len x len). Returns NULL if allocation failed. */ { unsigned i; WaveletMatrix2 y = NULL; if (len >= 2) { 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; } void freeWaveletMatrix2 (WaveletMatrix2 y, unsigned len) /* Free previously allocated 2D coefficient array y of size (len x len) */ { unsigned i, j; if (y) { for (i = 0; i < len; i++) free(y [i]); free(y); } } -#ifdef WAVELET_DBG - void zeroCoeffs2 (WaveletMatrix2 y, unsigned len) - /* Zero 2D output array to facilitate debugging */ +#ifdef WAVELET_DBG + void clearCoeffs2 (WaveletMatrix2 y, unsigned len) + /* Clear 2D output array to facilitate debugging */ { - unsigned i; + unsigned i, j, k; for (i = 0; i < len; i++) + #if 0 memset(y [i], 0, len * WAVELET_COEFFSIZE); + #else + for (j = 0; j < len; j++) + for (k = 0; k < 3; k++) + y [i] [j] [k] = NAN; + #endif } + static char *coeffStr (const WaveletCoeff3 coeff) + /* Format coefficient in string for dump. Cleared coeffs are represented + * as dots to facilitate debugging. */ + { + static char str [9]; + + if (isnan(coeff [0]) && isnan(coeff [1]) && isnan(coeff [2])) + /* Coeff is blank */ + return " ... "; + + snprintf(str, sizeof(str), "% 5.2f", coeffAvg(coeff)); + return str; + } + + + void dumpCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned y1Len, unsigned y2Len, float thresh ) /* Multipurpose routine to dump 2D arrays y1 and y2 of dims y1Len resp y2Len to stdout. These are presented side-by-side, separated by an arrow to indicate the input (y1) and output (y2) of a transform step. If either y1 or y2 is NULL, it will be omitted from the output, enabling dumping a single array. If thresh > 0, coefficients with absolute magnitude below this value are marked with brackets ('[]') as thresholded. */ { unsigned i, j, k; for (i = 0; i < max(y1Len, y2Len); i++) { if (y1) { for (j = 0; j < y1Len; j++) if (i < y1Len) printf((i | j) && coeffThresh(y1 [i][j], thresh) - ? "[% 5.2f]\t" : " % 5.2f\t", coeffAvg(y1 [i][j]) + ? "[%s]\t" : " %s\t", coeffStr(y1 [i][j]) ); else putchar('\t'); } if (y2) { printf(" -->\t"); for (j = 0; j < y2Len; j++) if (i < y2Len) printf((i | j) && coeffThresh(y2 [i][j], thresh) - ? "[% 5.2f]\t" : " % 5.2f\t", coeffAvg(y2 [i][j]) + ? "[%s]\t" : " %s\t", coeffStr(y2 [i][j]) ); else putchar('\t'); } putchar('\n'); } } float rmseCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned len ) /* Calculate RMSE between 2D matrices y1 and y2 */ { unsigned i, j; float d, rmse = 0; for (i = 0; i < len; i++) for (j = 0; j < len; j++) { #if 0 d = (coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j])) / coeffAvg(y1 [i][j]); #else d = coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j]); #endif rmse += d * d; } return sqrt(rmse / ((float)len * len)); } #endif static int haarStep2 (WaveletMatrix2 y, WaveletMatrix2 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; #ifdef WAVELET_DBG static unsigned axis = 0; #endif if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] = h2 * ( y [i] [j ] [k] + y [i] [j + 1] [k] ); /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = h2 * ( y [i] [j ] [k] - y [i] [j + 1] [k] ); } } } #ifdef WAVELET_DBG printf("%s FWD HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int haarInvStep2 (WaveletMatrix2 y, WaveletMatrix2 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; #ifdef WAVELET_DBG static unsigned axis = 1; #endif if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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++) { for (j = 0; j < len; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h2 * ( y [h ] [i] [k] + /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); yt [i] [j + 1] [k] = h2 * ( y [h ] [i] [k] - /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); } } } #ifdef WAVELET_DBG printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int d4Step2 (WaveletMatrix2 y, WaveletMatrix2 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; #ifdef WAVELET_DBG static unsigned axis = 0; #endif if (len < 4 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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 [h ] [i] [k] = h4 [0] * y [i] [j ] [k] + h4 [1] * y [i] [j + 1] [k] + h4 [2] * y [i] [j + 2] [k] + h4 [3] * y [i] [j + 3] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g4 [0] * y [i] [j ] [k] + g4 [1] * y [i] [j + 1] [k] + g4 [2] * y [i] [j + 2] [k] + g4 [3] * y [i] [j + 3] [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 [h ] [i] [k] = h4 [0] * y [i] [j ] [k] + h4 [1] * y [i] [j + 1] [k] + h4 [2] * y [i] [0 ] [k] + h4 [3] * y [i] [1 ] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g4 [0] * y [i] [j ] [k] + g4 [1] * y [i] [j + 1] [k] + g4 [2] * y [i] [0 ] [k] + g4 [3] * y [i] [1 ] [k]; } } #ifdef WAVELET_DBG printf("%s FWD D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return 0; } static int d4InvStep2 (WaveletMatrix2 y, WaveletMatrix2 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; #ifdef WAVELET_DBG static unsigned axis = 1; #endif if (len < 4 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return -1; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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 [i] [0] [k] = h4 [2] * y [hlen - 1] [i] [k] + /* Last avg */ g4 [2] * y [len - 1] [i] [k] + /* Last diff */ h4 [0] * y [0 ] [i] [k] + /* First avg */ g4 [0] * y [hlen ] [i] [k]; /* First diff */ yt [i] [1] [k] = h4 [3] * y [hlen - 1] [i] [k] + g4 [3] * y [len - 1] [i] [k] + h4 [1] * y [0 ] [i] [k] + g4 [1] * y [hlen ] [i] [k]; } /* Invert until upper boundary */ for (j = 2; j < len; j += 2) { h = (j >> 1) - 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h4 [2] * y [h ] [i] [k] + /* Avg */ g4 [2] * y [hlen + h ] [i] [k] + /* Diff */ h4 [0] * y [h + 1 ] [i] [k] + /* Next avg */ g4 [0] * y [hlen + h + 1] [i] [k]; /* Next diff */ yt [i] [j + 1] [k] = h4 [3] * y [h ] [i] [k] + g4 [3] * y [hlen + h ] [i] [k] + h4 [1] * y [h + 1 ] [i] [k] + g4 [1] * y [hlen + h + 1] [i] [k]; } } } #ifdef WAVELET_DBG printf("%s INV D4 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return 0; } int waveletXform2 (WaveletMatrix2 y, WaveletMatrix2 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; WaveletMatrix2 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 = allocWaveletMatrix2(len))) { 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 (d4Step2(y, yt, li) || d4Step2(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 (haarStep2(y, yt, li) || haarStep2(yt, y, li)) return -1; /* NOTE: All coefficients now in y */ if (ytloc) /* Free yt if allocated on demand */ freeWaveletMatrix2(ytloc, len); return 0; } int waveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 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; WaveletMatrix2 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 = allocWaveletMatrix2(len))) { 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 (haarInvStep2(y, yt, 1) || haarInvStep2(yt, y, 1)) return -1; for (li = 2; li <= l; li++) { /* Invert horizontal & vertical Daubechies D4 transform, swapping input and transposed output arrays */ if (d4InvStep2(y, yt, li) || d4InvStep2(yt, y, li)) return -1; } /* NOTE: Reconstructed signal now in y */ if (ytloc) /* Free yt if allocated on demand */ freeWaveletMatrix2(ytloc, len); return 0; } #ifdef WAVELET_TEST_2D #include #ifdef WAVELET_TEST_mRGBE #include "mrgbe.h" #endif #define WAVELET_TEST_INIT 0 int main (int argc, char *argv []) { int i, j, k, l; unsigned len, numThresh = 0; WaveletMatrix2 y0 = NULL, y = NULL; FILE *dataFile = NULL; WAVELET_COEFF inData, thresh = 0; #ifdef WAVELET_TEST_mRGBE #define HUGE 1e10 mRGBE mrgbeCoeff; mRGBERange mrgbeRange; WaveletMatrix2 ymrgbe = NULL; #endif if (argc < 2) { fprintf(stderr, "%s [threshold] [dataFile]\n", argv [0]); fputs("Missing array resolution l > 1, " "compression threshold >= 0\n", stderr ); return -1; } if (!(l = atoi(argv [1])) || l < 1) { fputs("Invalid array resolution l\n", stderr); return -1; } else len = 1 << l; if (argc > 2 && (thresh = atof(argv [2])) < 0) { fprintf(stderr, "Invalid threshold %.3f\n", thresh); return -1; } /* Allocate arrays for original and reconstruction */ if (!(y0 = allocWaveletMatrix2(len)) || !(y = allocWaveletMatrix2(len)) #ifdef WAVELET_TEST_mRGBE || !(ymrgbe = allocWaveletMatrix2(len)) #endif ) { fprintf(stderr, "Failed allocating %dx%d array\n", len, len); return -1; } if (argc > 3) { /* Load data from file; length must not exceed allocated */ if (!(dataFile = fopen(argv [3], "r"))) { fprintf(stderr, "Failed opening data file %s\n", argv [3]); 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, " %lf", &inData)) { y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = y [i][j][0] = y [i][j][1] = y [i][j][2] = inData; } else { fprintf(stderr, "Error reading from data file %s\n", argv [2] ); fclose(dataFile); return -1; } } } fclose(dataFile); } else { /* Init input */ srand48(111); for (i = 0; i < len; i++) { for (j = 0; j < len; j++) { for (k = 0; k < 3; k++) { y0 [i][j][k] = y [i][j][k] = #if WAVELET_TEST_INIT == 0 /* Random data, channel-independent */ drand48(); #elif WAVELET_TEST_INIT == 1 /* Random data, indentical for all channels */ k ? y [i][j][k - 1] : drand48(); #elif WAVELET_TEST_INIT == 2 /* Monotonically increasing along axis 0 */ i; #elif WAVELET_TEST_INIT == 3 /* Monotonically increasing along axis 1 */ j; #elif WAVELET_TEST_INIT == 4 /* Monotonically increasing along both axes */ i * j; #elif WAVELET_TEST_INIT == 5 /* Monotonically increasing by linear index */ i * len + j; #endif } } } } #ifdef WAVELET_DBG puts("----------------------- FWD XFORM -------------------------\n"); #endif /* Forward xform */ if (waveletXform2(y, NULL, l)) { fputs("Forward xform failed\n", stderr); return -1; } /* Threshold coefficients; we use hard thresholding as it's easier to * implement than soft thresholding, which requires sorting the * coefficients. NOTE: y [0][0] is omitted it contains the coarsest * approximation coefficient, which is essential for the * reconstruction! */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) { if ((i | j) && coeffThresh(y [i][j], thresh)) { y [i][j][0] = y [i][j][1] = y [i][j][2] = 0; numThresh++; #if 0 /* Replace thresholded values with random noise in range [-threshold, threshold] */ y [i][j][0] = y [i][j][1] = y [i][j][2] = thresh * ( 2 * drand48() - 1 ); #endif } } #ifdef WAVELET_DBG /* Dump coefficients */ puts("----------------------- COEFFICIENTS -------------------------\n"); dumpCoeffs2(y, NULL, len, 0, thresh); if (numThresh) printf("\n%d/%d coefficients thresholded = %.1f%% compression", numThresh, len * len, 100. * numThresh / (len * len) ); #endif #ifdef WAVELET_TEST_mRGBE /* Test mRGBE coefficient encoding */ mrgbeRange.min [0] = mrgbeRange.min [1] = mrgbeRange.min [2] = HUGE; mrgbeRange.max [0] = mrgbeRange.max [1] = mrgbeRange.max [2] = 0; /* Find min/max coeff and init mRGBE range */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) for (k = 0; k < 3; k++) { inData = fabs(y [i][j][k]); if (inData < mrgbeRange.min [k]) mrgbeRange.min [k] = inData; if (inData > mrgbeRange.max [k]) mrgbeRange.max [k] = inData; }; mRGBEinit(&mrgbeRange, mrgbeRange.min, mrgbeRange.max); /* Encode coeffs to mRGBE and back, but preserve approximation * coefficient at y [0][0] */ for (i = 0; i < len; i++) for (j = 0; j < len; j++) if (i | j) { mrgbeCoeff = mRGBEencode(y [i][j], &mrgbeRange, 0); mRGBEdecode(mrgbeCoeff, &mrgbeRange, ymrgbe [i][j]); } else for (k = 0; k < 3; k++) ymrgbe [i][j][k] = y [i][j][k]; #ifdef WAVELET_DBG /* Dump mRGBE-decoded coefficients */ puts("\n\n-------------------- mRGBE COEFFICIENTS ----------------------\n"); dumpCoeffs2(y, ymrgbe, len, len, 0); #endif #endif #ifdef WAVELET_DBG puts("\n----------------------- INV XFORM -------------------------\n"); #endif /* Inverse xform (also using mRGBE coeffs if enabled) */ if (waveletInvXform2(y, NULL, l) #ifdef WAVELET_TEST_mRGBE || waveletInvXform2(ymrgbe, NULL, l) #endif ) { fputs("\nInverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("\n--------------------- ORIG vs. INV XFORM ------------------------\n"); dumpCoeffs2(y0, y, len, len, 0); #endif printf("\nAvg RMSE = %.2f\n", rmseCoeffs2(y0, y, len)); #ifdef WAVELET_TEST_mRGBE #ifdef WAVELET_DBG puts("\n------------------ ORIG vs. INV XFORM + mRGBE -------------------\n"); dumpCoeffs2(y0, ymrgbe, len, len, 0); #endif printf("\nAvg RMSE with mRGBE enc = %.2f\n", rmseCoeffs2(y0, ymrgbe, len) ); #endif freeWaveletMatrix2(y0, len); freeWaveletMatrix2(y, len); #ifdef WAVELET_TEST_mRGBE freeWaveletMatrix2(ymrgbe, len); #endif return 0; } #endif /* WAVELET_TEST_2D */ diff --git a/wavelet2.h b/wavelet2.h index dd47a43..40a647f 100644 --- a/wavelet2.h +++ b/wavelet2.h @@ -1,200 +1,200 @@ /* ========================================================================= Definitions for 2D wavelet transforms on arrays of 3-tuples. 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. */ #ifndef WAVELET_COEFF #define WAVELET_COEFF double #endif /* Convenience macros for handling coefficients */ #define min(a, b) ((a) < (b) ? (a) : (b)) #define max(a, b) ((a) > (b) ? (a) : (b)) #define coeffAvg(a) (((a) [0] + (a) [1] + (a) [2]) / 3) #define coeffThresh(c, t) ( \ fabs((c) [0]) < (t) || fabs((c) [1]) < (t) || fabs((c) [2]) < (t) \ ) /* Haar / Daubechies D4 coefficients */ #define SQRT2 1.41421356 #define SQRT3 1.73205081 #define H4NORM (0.25 / SQRT2) /* Haar wavelet coeffs */ extern const WAVELET_COEFF h2; /* Daubechies D4 wavelet coeffs */ extern const WAVELET_COEFF h4 [4]; extern const WAVELET_COEFF g4 [4]; /* 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 WAVELET_COEFF WaveletCoeff3 [3]; #define WAVELET_COEFFSIZE (sizeof(WaveletCoeff3)) typedef WaveletCoeff3 **WaveletMatrix2; WaveletMatrix2 allocWaveletMatrix2 (unsigned len); /* Allocate and return a 2D coefficient array of size (len x len). Returns NULL if allocation failed. */ void freeWaveletMatrix2 (WaveletMatrix2 y, unsigned len); /* Free previously allocated 2D coefficient array y of size (len x len) */ /* ---------- 2D WAVELET TRANSFORM OPTIMISED FOR SIZES 2^l ---------- ------ THESE USE PERIODIC (WRAP-AROUND) BOUNDARY EXTENSION ------- */ int waveletXform2 (WaveletMatrix2 y, WaveletMatrix2 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. */ int waveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 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. */ /* ------------ 2D WAVELET TRANSFORM FOR ARBITRARY SIZES ------------- ---- THESE USE CONSTANT BOUNDARY EXTENSION WITH PADDING COEFFS----- */ unsigned padWaveletXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0 ); /* Perform full 2D multiresolution forward wavelet transform with padded boundary coefficients on array y containing (yLen0 x yLen0) samples of 3-tuples. Note no intra-tuple transform occurs. An additional array yt of identical dimension to y is required as buffer for intermediate results. The wavelet coefficients are returned in array y, containing the coarsest approximation coefficients in y [0][0] followed by horizontal/vertical detail coefficients in order of increasing resolution/frequency. NOTE: Both y and yt must be large enough to accommodate the additional padding coefficients introduced at the array boundaries during the transform. It is the caller's responsibility to preallocate and dimension these arrays accordingly with allocWaveletMatrix2(). The maximum expected number of coefficients (which generally exceeds the number of input samples yLen0) must be interrogated beforehand by calling the function with y or yt set to NULL. The returned value is the maximum number of output coefficients, or 0 if the transform failed. */ unsigned padWaveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0, unsigned yLenOrg ); /* Perform full 2D multiresolution inverse wavelet transform on array y containing (yLen0 x yLen0) boundary padded wavelet coefficients as 3-tuples. The original number of samples (yLenOrg x yLenOrg) is required for the reconstruction. An additional array yt of identical dimension to y is required as buffer for intermediate results. The reconstructed samples are returned in y. This is smaller than (yLen0 x yLen0) since the boundary coefficients from the forward transform are collapsed. The returned value is the reconstructed number of samples per dimension (yLenOrg), or 0 if the inverse transform failed. */ #ifdef WAVELET_DBG /* ----------------- SUPPORT ROUTINES FOR DEBUGGING ----------------- */ - void zeroCoeffs2 (WaveletMatrix2 y, unsigned len); - /* Zero 2D output array to facilitate debugging */ + void clearCoeffs2 (WaveletMatrix2 y, unsigned len); + /* Clear 2D output array to facilitate debugging */ void dumpCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned y1Len, unsigned y2Len, float thresh ); /* Dump 2D arrays y1 and y2 of dims y1Len resp. y2Len side-by-side to stdout (skip y2 if NULL). Coefficients with absolute magnitude less than thresh are marked with brackets ('[]') as thresholded. */ float rmseCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned len ); /* Calculate RMSE between 2D matrices y1 and y2 */ #endif #endif diff --git a/wavelet3.c b/wavelet3.c index c8ded57..eec430a 100644 --- a/wavelet3.c +++ b/wavelet3.c @@ -1,807 +1,811 @@ /* ========================================================================= 2-dimensional Daubechies D4 and Haar wavelet transforms on arbitrary sized arrays of 3-tuples. Array boundaries are extended as constants, which introduces padding coefficients. Compile with -DWAVELET_TEST_2DPAD to build standalone unit tests. The implementation here is mainly distilled from the following sources: -- Filip Wasilewski's pywavelets: https://github.com/PyWavelets/pywt/tree/master/pywt/_extensions/c [D4 with boundary condition handling via padding, arbitrary sizes] -- 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: 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 "wavelet2.h" #include #include #include #include /* Aggregated Dabuechies D4 coeffs for padding at left boundary */ static const WAVELET_COEFF h4Bound = H4NORM * (7 + SQRT3), g4Bound = H4NORM * (1 + SQRT3); /* h4Bound = h4 [0] + h4 [1] + h4 [2]; g4Bound = g4 [0] + g4 [1] + g4 [2]; */ static unsigned padHaarStep2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned len ) /* Single iteration of forward 2D Haar wavelet transform with padding on array y of arbitrary size (len x len) containing 3-tuples. This variant uses constant boundary extension. 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 array dimension len on success, else 0. */ { const unsigned hLen = -(-(signed)len >> 1); int h, i, j, k; #ifdef WAVELET_DBG static unsigned axis = 0; #endif if (len < 2 || !y || !yt) /* Input shorter than wavelet support, or no input/output */ return 0; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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++) { for (j = 0; j <= len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] = h2 * ( y [i] [j ] [k] + y [i] [j + 1] [k] ); /* Detail/diff/highpass */ yt [h + hLen] [i] [k] = h2 * ( y [i] [j ] [k] - y [i] [j + 1] [k] ); } } if (len & 1) { /* Odd length; pad with last value as additional approx term */ for (k = 0; k < 3; k++) /* Smooth/approx/avg/lowpass */ yt [hLen - 1] [i] [k] = y [i] [len - 1] [k]; } } #ifdef WAVELET_DBG printf("%s FWD HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return len; } static unsigned padHaarInvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned len ) /* Single iteration of inverse 2D Haar wavelet transform with padding on array y of arbitrary size (len x len) containing 3-tuples. 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 array dimension len on success, else 0. */ { const unsigned hLen = -(-(signed)len >> 1); int h, i, j, k; #ifdef WAVELET_DBG static unsigned axis = 1; #endif if (len < 2 || !y || !yt) /* Too few coeffs for reconstruction, or no input/output */ return 0; #ifdef WAVELET_DBG - zeroCoeffs2(yt, len); + clearCoeffs2(yt, len); #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++) { for (j = 0; j <= len - 2; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h2 * ( y [h ] [i] [k] + /* Approx */ y [h + hLen] [i] [k] /* Detail */ ); yt [i] [j + 1] [k] = h2 * ( y [h ] [i] [k] - /* Approx */ y [h + hLen] [i] [k] /* Detail */ ); } } - + if (len & 1) { /* Odd length; pad additional value from last approximation term */ for (k = 0; k < 3; k++) yt [i] [len - 1] [k] = y [hLen - 1] [i] [k]; } } #ifdef WAVELET_DBG printf("%s INV HAAR (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0); putchar('\n'); axis ^= 1; #endif return len; } static unsigned padD4Step2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen, unsigned ytOffset ) /* Single iteration of forward 2D Daubechies D4 wavelet transform with boundary padding on array y containing (yLen x yLen) samples of 3-tuples. This variant uses constant boundary extension, which results in padding coefficients. 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. 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. Both the input array y and output array yt must be large enough to also accommodate the additional padding coefficients introduced at the boundaries. The number of output coefficient pairs ytLen returned in yt is obtained as return value, and can be interrogated beforehand without performing a transform by passing either y or yt as NULL. The approximation coefficients are returned in yt [0..ytLen-1] [i], while the details are returned in yt [ytOffset..ytOffset+ytLen-1] [i]. The offset ytOffset must be supplied by the caller to accommodate the detail coefficients (plus padding!) without clobbering the approximation coefficients (again, plus padding!). Consequently, there is a gap between the approximation coefficients in the upper left of the matrix, and the detail coefficients in the lower right. This gap closes with each iteration. */ { unsigned h, i, j, k, l; /* Number of output coefficient pairs ytLen incl. padding at boundaries, offset hyt for detail coeffs */ const unsigned ytLen = (yLen + 3) >> 1, totalLen = ytOffset + ytLen; #ifdef WAVELET_DBG static unsigned axis = 0; #endif if (y && yt && yLen >= 2) { #ifdef WAVELET_DBG - /* Zero output array to facilitate debugging */ - zeroCoeffs2(yt, totalLen); + /* TODO; DO WE NEED THIS??? */ + /* Clear output array to facilitate debugging */ + clearCoeffs2(yt, totalLen); #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 transforms over the alternate axis. This is done by simply * swapping the indices during assignment */ for (i = 0; i < yLen; i++) { /* First iteration (j = -2); pad coeffs at left edge with * constant boundary extension (note multiple factoring of * y [i][0] via aggregated coefficients h4Bound and g4Bound) */ for (k = 0; k < 3; k++) { /* Smooth/approx/avg/lowpass */ yt [0 ] [i] [k] = h4Bound * y [i] [0] [k] + h4 [3] * y [i] [1] [k]; /* Detail/diff/highpass */ yt [ytOffset] [i] [k] = g4Bound * y [i] [0] [k] + g4 [3] * y [i] [1] [k]; } /* Intermediate and last iterations (j >= 0) */ for (j = 0; j < yLen; j += 2) { h = (j >> 1) + 1; for (k = 0; k < 3; k++) { yt [h] [i] [k] = yt [h + ytOffset] [i] [k] = 0; /* Convolve up to boundary of wavelet or input, whichever * comes first */ for (l = j; l < min(j + 4, yLen); l++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] += h4 [l - j] * y [i] [l] [k]; /* Detail/diff/highpass */ yt [h + ytOffset] [i] [k] += g4 [l - j] * y [i] [l] [k]; } /* If necessary, convolve beyond input boundary with * constant extension (note multiple factoring of * y [i] [yLen-1] via aggregated coefficients h4Bound and * g4Bound). This is skipped if the wavelet boundary was * already reached */ for (; l < j + 4; l++) { /* Smooth/approx/avg/lowpass */ yt [h ] [i] [k] += h4 [l - j] * y [i] [yLen - 1] [k]; /* Detail/diff/highpass */ yt [h + ytOffset] [i] [k] += g4 [l - j] * y [i] [yLen - 1] [k]; } } } } + /* TODO: IS THIS NECESSARY??? */ if (axis) { /* Copy padding coeffs from preceding horizontal iteration */ for (i = yLen; i < totalLen; i++) for (j = 0; j < yLen; j++) for (k = 0; k < 3; k++) yt [j] [i] [k] = y [i] [j] [k]; } #ifdef WAVELET_DBG printf("%s FWD D4 (%d x %d) -> (%d x %d)\n", axis ? "VERT" : "HORIZ", yLen, yLen, ytLen << 1, ytLen << 1 ); dumpCoeffs2(y, yt, totalLen, totalLen, 0); putchar('\n'); axis ^= 1; #endif } return ytLen; } static unsigned padD4InvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen, unsigned yOffset ) /* Single iteration of inverse 2D Daubechies D4 wavelet transform with boundary padding on array y containg (yLen x yLen) approximation and detail coefficients in 3-tuples. 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. 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. The number of output coefficients is smaller than the input as padding coefficients are collapsed at the boundaries during the inverse transform. The number of output coefficients per dimension ytLen is obtained as returned value, and can be interrogated beforehand without performing a transform by passing either y or yt as NULL. The inversion expects approximation coefficients at y [0..yLen/2-1] and detail coefficients at y [yOffset..yOffset+yLen/2-1]. The offset yOffset must be supplied by the caller to account for the fact that the number of output coefficients ytLen is smaller than yLen due to the collapsed padding coefficients. Consequently there is a gap between the inverted coefficients in the upper left of the matrix, and the remaining detail coefficients in the lower right. This gap opens with each iteration, until all detail coefficients have been consumed by the inversion. */ { unsigned h, i, j, k, l, ytLen = yLen - 2; /* Max num output coeffs in yt */ const unsigned totalLen = yOffset + (yLen >> 1); #ifdef WAVELET_DBG static unsigned axis = 1; #endif if (y && yt && yLen >= 2) { #ifdef WAVELET_DBG - /* Zero output array to facilitate debugging */ - zeroCoeffs2(yt, totalLen); + /* Clear output array to facilitate debugging */ + /* TODO; DO WE NEED THIS??? */ + clearCoeffs2(yt, totalLen); #endif /* Reconstruct from approximation and detail coefficients */ for (i = 0; i < ytLen; i++) for (j = 0; j < ytLen - 1; j += 2) { h = j >> 1; for (k = 0; k < 3; k++) { yt [i] [j ] [k] = h4 [2] * y [h ] [i] [k] + /* Approx */ g4 [2] * y [h + yOffset ] [i] [k] + /* Detail */ h4 [0] * y [h + 1 ] [i] [k] + /* Approx */ g4 [0] * y [h + 1 + yOffset] [i] [k]; /* Detail */ yt [i] [j + 1] [k] = h4 [3] * y [h ] [i] [k] + /* Approx */ g4 [3] * y [h + yOffset ] [i] [k] + /* Detail */ h4 [1] * y [h + 1 ] [i] [k] + /* Approx */ g4 [1] * y [h + 1 + yOffset] [i] [k]; /* Detail */ } } ytLen = j; - if (axis) { + /* TODO: IS THIS NECESSARY??? */ + if (0 && axis) { /* Copy padding coeffs for subsequent horizontal iteration */ for (i = ytLen; i < totalLen; i++) for (j = 0; j < ytLen; j++) for (k = 0; k < 3; k++) yt [i] [j] [k] = y [j] [i] [k]; } #ifdef WAVELET_DBG printf("%s INV D4 (%d x %d) -> (%d x %d)\n", axis ? "VERT" : "HORIZ", yLen, yLen, ytLen, ytLen ); dumpCoeffs2(y, yt, totalLen, totalLen, 0); putchar('\n'); axis ^= 1; #endif } return ytLen; } unsigned padWaveletXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0 ) /* Perform full 2D multiresolution forward wavelet transform with padded boundary coefficients on array y containing (yLen0 x yLen0) samples of 3-tuples. Note no intra-tuple transform occurs. An additional array yt of identical dimension to y is required as buffer for intermediate results. The wavelet coefficients are returned in array y, containing the coarsest approximation coefficients in y [0][0] followed by horizontal/vertical detail coefficients in order of increasing resolution/frequency. NOTE: Both y and yt must be large enough to accommodate the additional padding coefficients introduced at the array boundaries during the transform. It is the caller's responsibility to preallocate and dimension these arrays accordingly with allocWaveletMatrix2(). The maximum expected number of coefficients (which generally exceeds the number of input samples yLen0) must be interrogated beforehand by calling the function with y or yt set to NULL. The returned value is the maximum number of output coefficients, or 0 if the transform failed. */ { /* Output size (number of coeffs, including boundary padding coeffs), output array offset and pointers */ unsigned yLen, ytLen, ytOffset, totalLen = 0; /* Transform if enough samples, else just return output size. NOTE: the D4 transform will be skipped for input sizes yLen0 < 4. */ if (yLen0 >= 2) { /* Figure out total number of coefficients (= maximum size of output array including boundary padding) by summing the number of detail coeffs returned by each xform iteration. */ for (yLen = yLen0, totalLen = 0; yLen > 3; totalLen += yLen = padD4Step2(NULL, NULL, yLen, 0) ); /* Finally add number of approximation coeffs at coarsest resolution */ totalLen += yLen; if (y && yt) { /* Apply Daubechies D4 transform, accounting for boundary padding */ for (yLen = yLen0, ytOffset = totalLen; yLen > 3; yLen = ytLen ) { /* Get num of coeff pairs returned by the next xform iteration, and calculate the offset of detail coeffs in output array */ ytOffset -= ytLen = padD4Step2(NULL, NULL, yLen, 0); /* Apply horizontal & vertical Daubechies D4 transform, swapping input and transposed output array. The resulting detail coefficients are prepended to those from the previous iteration from position ytOffset; consequently, the returned number of coefficient pairs is the input size for next iteration */ if (!padD4Step2(y, yt, yLen, ytOffset) || !padD4Step2(yt, y, yLen, ytOffset) ) return 0; } /* Apply horizontal & vertical Haar transform in penultimate and final iterations (since yLen is odd) to decompose to a single approximation coefficient y [0] [0] at the coarsest resolution; all other coeffs are details. Obviously no boundary padding is necessary here. */ for (; yLen >= 2; yLen--) { if (!padHaarStep2(y, yt, yLen) || !padHaarStep2(yt, y, yLen)) return 0; } } } /* All coeffs now in y; return output size */ return totalLen; } unsigned padWaveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0, unsigned yLenOrg ) /* Perform full 2D multiresolution inverse wavelet transform on array y containing (yLen0 x yLen0) boundary padded wavelet coefficients as 3-tuples. The original number of samples (yLenOrg x yLenOrg) is required for the reconstruction. An additional array yt of identical dimension to y is required as buffer for intermediate results. The reconstructed samples are returned in y. This is smaller than (yLen0 x yLen0) since the boundary coefficients from the forward transform are collapsed. The returned value is the reconstructed number of samples per dimension (yLenOrg), or 0 if the inverse transform failed. */ { unsigned yLen = 0, nextyLen, ytLen, ytOffset, i, j, d4Iter; if (yLenOrg < 2) /* Too few coefficients for reconstruction */ return 0; /* Figure out num Daubechies iterations based on original num coeffs */ for (d4Iter = 0, yLen = yLenOrg; yLen > 3; d4Iter++, yLen = padD4Step2(NULL, NULL, yLen, 0) ); /* Invert horizontal & vertical Haar transform at two coarsest levels */ for (yLen = 2; yLen <= min(3, yLen0); yLen++) if (!padHaarInvStep2(y, yt, yLen) || !padHaarInvStep2(yt, y, yLen)) return 0; /* Invert horiz & vert Daubechies D4 transform on subsequent levels */ for (i = d4Iter; i; i--) { /* Get num detail coeffs and their offset for current iteration */ for (j = 0, ytOffset = yLen0, yLen = nextyLen = yLenOrg; j < i; j++) { nextyLen = yLen; ytOffset -= yLen = padD4Step2(NULL, NULL, yLen, 0); } if (!(ytLen = padD4InvStep2(y, yt, yLen << 1, ytOffset)) || !(ytLen = padD4InvStep2(yt, y, yLen << 1, ytOffset)) ) return 0; #if 0 /* XXX: Drop excess reconstructed approximation coeffs to match number of detail coeffs in next iteration. */ /* ??? */ yLen = min(ytLen, nextyLen); #endif } return yLenOrg; } #ifdef WAVELET_TEST_2DPAD #include #ifdef WAVELET_TEST_mRGBE #include "mrgbe.h" #endif #define WAVELET_TEST_INIT 1 int main (int argc, char *argv []) { int i, j, k; unsigned y0Len, yLen, numThresh = 0; WaveletMatrix2 y0 = NULL, y = NULL, yt = NULL; FILE *dataFile = NULL; WAVELET_COEFF inData, thresh = 0; #ifdef WAVELET_TEST_mRGBE #define HUGE 1e10 mRGBE mrgbeCoeff; mRGBERange mrgbeRange; WaveletMatrix2 ymrgbe = NULL; #endif if (argc < 2) { fprintf(stderr, "%s [threshold] [dataFile]\n", argv [0]); fputs("Missing array dimension dim >= 2, " "compression threshold >= 0\n", stderr ); return -1; } if (!(y0Len = atoi(argv [1])) || y0Len < 2) { fprintf(stderr, "Invalid array dimension %d\n", y0Len); return -1; } if (argc > 2 && (thresh = atof(argv [2])) < 0) { fprintf(stderr, "Invalid threshold %.3f\n", thresh); return -1; } /* Get size of (padded) output array */ if (!(yLen = padWaveletXform2(NULL, NULL, y0Len))) { fputs("Can't determine number of coefficients " "(input too short?)\n", stderr ); return -1; } /* Allocate & clear arrays for original and reconstruction */ if (!(y0 = allocWaveletMatrix2(y0Len)) || !(y = allocWaveletMatrix2(yLen)) || !(yt = allocWaveletMatrix2(yLen)) #ifdef WAVELET_TEST_mRGBE || !(ymrgbe = allocWaveletMatrix2(yLen)) #endif ) { fprintf(stderr, "Failed allocating %dx%d array\n", yLen, yLen); return -1; } - zeroCoeffs2(y0, y0Len); - zeroCoeffs2(y, yLen); - zeroCoeffs2(yt, yLen); + clearCoeffs2(y0, y0Len); + clearCoeffs2(y, yLen); + clearCoeffs2(yt, yLen); if (argc > 3) { /* Load data from file; length must not exceed allocated */ if (!(dataFile = fopen(argv [3], "r"))) { fprintf(stderr, "Failed opening data file %s\n", argv [3]); return -1; } for (i = 0; i < y0Len; i++) { for (j = 0; j < y0Len; 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, " %lf", &inData)) { y0 [i][j][0] = y0 [i][j][1] = y0 [i][j][2] = y [i][j][0] = y [i][j][1] = y [i][j][2] = inData; } else { fprintf(stderr, "Error reading from data file %s\n", argv [2] ); fclose(dataFile); return -1; } } } fclose(dataFile); } else { /* Init input */ srand48(111); for (i = 0; i < y0Len; i++) { for (j = 0; j < y0Len; j++) { for (k = 0; k < 3; k++) { y0 [i][j][k] = y [i][j][k] = #if WAVELET_TEST_INIT == 0 /* Random data, channel-independent */ drand48(); #elif WAVELET_TEST_INIT == 1 /* Random data, indentical for all channels */ k ? y [i][j][k - 1] : drand48(); #elif WAVELET_TEST_INIT == 2 /* Monotonically increasing along axis 0 */ i; #elif WAVELET_TEST_INIT == 3 /* Monotonically increasing along axis 1 */ j; #elif WAVELET_TEST_INIT == 4 /* Monotonically increasing along both axes */ i * j; #elif WAVELET_TEST_INIT == 5 /* Monotonically increasing by linear index */ i * y0Len + j; #endif } } } } #ifdef WAVELET_DBG puts("----------------------- FWD XFORM -------------------------\n"); #endif /* Forward xform */ if (!padWaveletXform2(y, yt, y0Len)) { fputs("Forward xform failed\n", stderr); return -1; } /* Threshold coefficients; we use hard thresholding as it's easier to * implement than soft thresholding, which requires sorting the * coefficients. NOTE: the coarsest approximation coefficient at * y [0][0] is omitted as it is essential for reconstruction! */ for (i = 0; i < yLen; i++) for (j = 0; j < yLen; j++) { if ((i | j) && coeffThresh(y [i][j], thresh)) { y [i][j][0] = y [i][j][1] = y [i][j][2] = 0; numThresh++; #if 0 /* Replace thresholded values with random noise in range [-threshold, threshold] */ y [i][j][0] = y [i][j][1] = y [i][j][2] = thresh * ( 2 * drand48() - 1 ); #endif } } #ifdef WAVELET_DBG /* Dump coefficients */ puts("----------------------- COEFFICIENTS -------------------------\n"); dumpCoeffs2(y, NULL, yLen, 0, thresh); if (numThresh) printf("\n%d/%d coefficients thresholded = %.1f%% compression", numThresh, yLen * yLen, 100. * numThresh / (yLen * yLen) ); #endif #ifdef WAVELET_TEST_mRGBE /* Test mRGBE coefficient encoding */ mrgbeRange.min [0] = mrgbeRange.min [1] = mrgbeRange.min [2] = HUGE; mrgbeRange.max [0] = mrgbeRange.max [1] = mrgbeRange.max [2] = 0; - zeroCoeffs2(ymrgbe, yLen); + clearCoeffs2(ymrgbe, yLen); /* Find min/max coeff and init mRGBE range */ for (i = 0; i < yLen; i++) for (j = 0; j < yLen; j++) for (k = 0; k < 3; k++) { inData = fabs(y [i][j][k]); if (inData < mrgbeRange.min [k]) mrgbeRange.min [k] = inData; if (inData > mrgbeRange.max [k]) mrgbeRange.max [k] = inData; }; mRGBEinit(&mrgbeRange, mrgbeRange.min, mrgbeRange.max); /* Encode coeffs to mRGBE and back, but preserve approximation * coefficient at y [0][0] */ for (i = 0; i < yLen; i++) for (j = 0; j < yLen; j++) if (i | j) { mrgbeCoeff = mRGBEencode(y [i][j], &mrgbeRange, 0); mRGBEdecode(mrgbeCoeff, &mrgbeRange, ymrgbe [i][j]); } else for (k = 0; k < 3; k++) ymrgbe [i][j][k] = y [i][j][k]; #ifdef WAVELET_DBG /* Dump mRGBE-decoded coefficients */ puts("\n\n-------------------- mRGBE COEFFICIENTS ----------------------\n"); dumpCoeffs2(y, ymrgbe, yLen, yLen, 0); #endif #endif #ifdef WAVELET_DBG puts("\n----------------------- INV XFORM -------------------------\n"); #endif /* Inverse xform (also using mRGBE coeffs if enabled) */ if (!padWaveletInvXform2(y, yt, yLen, y0Len) #ifdef WAVELET_TEST_mRGBE || !padWaveletInvXform2(ymrgbe, yt, yLen, y0Len) #endif ) { fputs("\nInverse xform failed\n", stderr); return -1; } #ifdef WAVELET_DBG puts("\n--------------------- ORIG vs. INV XFORM ------------------------\n"); dumpCoeffs2(y0, y, y0Len, y0Len, 0); #endif printf("\nAvg RMSE = %.2f\n", rmseCoeffs2(y0, y, y0Len)); #ifdef WAVELET_TEST_mRGBE #ifdef WAVELET_DBG puts("\n------------------ ORIG vs. INV XFORM + mRGBE -------------------\n"); dumpCoeffs2(y0, ymrgbe, y0Len, y0Len, 0); #endif printf("\nAvg RMSE with mRGBE enc = %.2f\n", rmseCoeffs2(y0, ymrgbe, y0Len) ); #endif freeWaveletMatrix2(y0, y0Len); freeWaveletMatrix2(y, yLen); freeWaveletMatrix2(yt, yLen); #ifdef WAVELET_TEST_mRGBE freeWaveletMatrix2(ymrgbe, yLen); #endif return 0; } #endif /* WAVELET_TEST_2DPAD */