diff --git a/wavelet2.c b/wavelet2.c index f869245..cf71df3 100644 --- a/wavelet2.c +++ b/wavelet2.c @@ -1,912 +1,916 @@ /* ========================================================================= 2-dimensional Daubechies D2 and D1 (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 /* 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 */ /* 2-tap Haar wavelet/scaling func coeffs */ const WAVELET_COEFF h1 = 1 / SQRT2; /* 4-tap wavelet/scaling func coeffs; default is Daubechies D2 */ #ifdef WAVELET_BIOR /* Biorthogonal 3.1 */ const WAVELET_COEFF h2 [4] = { -0.3535533905932738, 1.0606601717798214, 1.0606601717798214, -0.3535533905932738 }, g2 [4] = { -0.1767766952966369, 0.5303300858899107, -0.5303300858899107, 0.1767766952966369 }, /* Inverse coeffs */ h2i [4] = { 0.1767766952966369, 0.5303300858899107, 0.5303300858899107, 0.1767766952966369 }, /* = { -g2 [0], g2 [1], -g2 [2], g2 [3] } */ g2i [4] = { -0.3535533905932738, -1.0606601717798214, 1.0606601717798214, 0.3535533905932738 }; /* = { h2 [0], -h2 [1], h2 [2], -h2 [3] } */ #else /* Daubechies D2 */ const WAVELET_COEFF h2 [4] = { H2NORM * (1 + SQRT3), H2NORM * (3 + SQRT3), H2NORM * (3 - SQRT3), H2NORM * (1 - SQRT3) }, g2 [4] = { H2NORM * (1 - SQRT3), -H2NORM * (3 - SQRT3), H2NORM * (3 + SQRT3), -H2NORM * (1 + SQRT3) }, /* = { h2 [3], -h2 [2], h2 [1], -h2 [0] } */ /* Inverse coeffs */ h2i [4] = { H2NORM * (1 - SQRT3), H2NORM * (3 - SQRT3), H2NORM * (3 + SQRT3), H2NORM * (1 + SQRT3) }, /* = { h2 [3], h2 [2], h2 [1], h2 [0] } */ g2i [4] = { -H2NORM * (1 + SQRT3), H2NORM * (3 + SQRT3), -H2NORM * (3 - SQRT3), H2NORM * (1 - SQRT3) }; /* = { g2 [3], g2 [2], g2 [1], g2 [0] } */ #endif 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); } } void zeroCoeffs2 (WaveletMatrix2 y, unsigned len) /* Set 2D array coefficients to zero */ { unsigned i; if (y) for (i = 0; i < len; i++) memset(y [i], 0, len * WAVELET_COEFFSIZE); } #if 1 void clearCoeffs2 (WaveletMatrix2 y, unsigned len) /* Clear 2D array to facilitate debugging */ { unsigned i, j, k; if (y) for (i = 0; i < len; i++) for (j = 0; j < len; j++) for (k = 0; k < 3; k++) y [i] [j] [k] = NAN; } - static char *coeffStr (const WaveletCoeff3 coeff) - /* Format coefficient in string for dump. Cleared coeffs are represented - * as dots to facilitate debugging. */ + static char *coeffStr (const WaveletCoeff3 coeff, int compact) + /* Format coefficient in string for dump. If the compact flag is set + (e.g. for dumping to the console), coefficients are printed with + low precision.to reduce clutter, and cleared coefficients are + represented as dots for clarity. + If the compact flag is unset, coefficients are printed at full + precision, and cleared coefficients are printed as-is (zeros). */ { static char str [9]; - if (coeffIsEmpty(coeff)) + if (compact && coeffIsEmpty(coeff)) /* Coeff is blank */ return " .. "; - snprintf(str, sizeof(str), "% 5.2f", coeffAvg(coeff)); + snprintf(str, sizeof(str), compact ? "% 5.2f" : "%f", coeffAvg(coeff)); return str; } void dumpCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned y1Len, unsigned y2Len, float thresh, FILE *stream ) /* Multipurpose routine to dump 2D arrays y1 and y2 of dims y1Len resp y2Len. 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. If stream is not NULL, the coefficients will be dumped to the - specified file, else to stdout. + specified file, else to stdout. */ { unsigned i, j, k; FILE *out = stream ? stream : stdout; for (i = 0; i < max(y1Len, y2Len); i++) { if (y1) { for (j = 0; j < y1Len; j++) if (i < y1Len) fprintf(out, "%s\t", coeffThresh(y1, i, j, thresh) - ? "[ .. ]" : coeffStr(y1 [i][j]) + ? "[ .. ]" : coeffStr(y1 [i][j], out == stdout) ); else fputc('\t', out); } if (y2) { fprintf(out, " -->\t"); for (j = 0; j < y2Len; j++) if (i < y2Len) fprintf(out, "%s\t", coeffThresh(y2, i, j, thresh) - ? "[ .. ]" : coeffStr(y2 [i][j]) + ? "[ .. ]" : coeffStr(y2 [i][j], out == stdout) ); else fputc('\t', out); } fputc('\n', out); } } 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 d1Step2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l) /* Single step of forward 2D Daubechies D1 (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 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] = h1 * ( y [i] [j ] [k] + y [i] [j + 1] [k] ); /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = h1 * ( y [i] [j ] [k] - y [i] [j + 1] [k] ); } } } #ifdef WAVELET_DBG printf("%s FWD D1 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0, NULL); putchar('\n'); axis ^= 1; #endif return 0; } static int d1InvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l) /* Single step of inverse 2D Daubechies D1 (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 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] = h1 * ( y [h ] [i] [k] + /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); yt [i] [j + 1] [k] = h1 * ( y [h ] [i] [k] - /* Avg */ y [hlen + h] [i] [k] /* Diff */ ); } } } #ifdef WAVELET_DBG printf("%s INV D1 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0, NULL); putchar('\n'); axis ^= 1; #endif return 0; } static int d4Step2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l) /* Single step of forward 2D Daubechies D2 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 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] = h2 [0] * y [i] [j ] [k] + h2 [1] * y [i] [j + 1] [k] + h2 [2] * y [i] [j + 2] [k] + h2 [3] * y [i] [j + 3] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g2 [0] * y [i] [j ] [k] + g2 [1] * y [i] [j + 1] [k] + g2 [2] * y [i] [j + 2] [k] + g2 [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] = h2 [0] * y [i] [j ] [k] + h2 [1] * y [i] [j + 1] [k] + h2 [2] * y [i] [0 ] [k] + h2 [3] * y [i] [1 ] [k]; /* Detail/diff/highpass */ yt [hlen + h] [i] [k] = g2 [0] * y [i] [j ] [k] + g2 [1] * y [i] [j + 1] [k] + g2 [2] * y [i] [0 ] [k] + g2 [3] * y [i] [1 ] [k]; } } #ifdef WAVELET_DBG printf("%s FWD D2 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0, NULL); putchar('\n'); axis ^= 1; #endif return 0; } static int d4InvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l) /* Single step of inverse 2D Daubechies D2 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 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] = h2i [1] * y [hlen - 1] [i] [k] + /* Last avg */ g2i [1] * y [len - 1] [i] [k] + /* Last diff */ h2i [3] * y [0 ] [i] [k] + /* First avg */ g2i [3] * y [hlen ] [i] [k]; /* First diff */ yt [i] [1] [k] = h2i [0] * y [hlen - 1] [i] [k] + g2i [0] * y [len - 1] [i] [k] + h2i [2] * y [0 ] [i] [k] + g2i [2] * 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] = h2i [1] * y [h ] [i] [k] + /* Avg */ g2i [1] * y [hlen + h ] [i] [k] + /* Diff */ h2i [3] * y [h + 1 ] [i] [k] + /* Next avg */ g2i [3] * y [hlen + h + 1] [i] [k]; /* Next diff */ yt [i] [j + 1] [k] = h2i [0] * y [h ] [i] [k] + g2i [0] * y [hlen + h ] [i] [k] + h2i [2] * y [h + 1 ] [i] [k] + g2i [2] * y [hlen + h + 1] [i] [k]; } } } #ifdef WAVELET_DBG printf("%s INV D2 (%d x %d)\n", axis ? "VERT" : "HORIZ", len, len); dumpCoeffs2(y, yt, len, len, 0, NULL); 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 D2 transform, swapping input and transposed output array */ if (d4Step2(y, yt, li) || d4Step2(yt, y, li)) return -1; } #ifdef WAVELET_FINAL_HAAR /* Apply horizontal & vertical Daubechies D1 (Haar) transform at coarsest resolution (li==1) to obtain single approximation coefficient at y [0][0]; all other coeffs are details. */ if (d1Step2(y, yt, li) || d1Step2(yt, y, li)) return -1; #endif /* 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; } #ifdef WAVELET_FINAL_HAAR /* Invert horizontal & vertical Haar transform at coarsest level (li==1), swapping input and transposed output array */ if (d1InvStep2(y, yt, 1) || d1InvStep2(yt, y, 1)) return -1; #endif for (li = 2; li <= l; li++) { /* Invert horizontal & vertical Daubechies D2 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 #ifdef WAVELET_TEST_mRGBE #include "mrgbe.h" #endif #ifndef WAVELET_TEST_INIT /* Default test data initialisation */ #define WAVELET_TEST_INIT 0 #endif #if (WAVELET_TEST_INIT < 0 || WAVELET_TEST_INIT > 5) #error Invalid wavelet unit test data initialiation option #endif 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 TINY 1e-6 #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 (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, NULL); 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, NULL); #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, NULL); #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, NULL); #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 */