Homec4science
Diffusion RADIANCE Photon Map (fixSpecular) Browse

RadiancePmap/wavelet.c
eda1af4994b7fixSpecular

Last Change
Raw File

wavelet.c
View Options

/*
=========================================================================
1-dimensional Daubechies D2 and D1 (Haar) wavelet transforms on arbitrary
sized arrays of 3-tuples. Array boundaries are extended as constants,
which introduces padding coefficients.
Compile with -DWAVELET_TEST_1D 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
[D2 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 D2]
-- 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 D2]
-- 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 D2]
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
=========================================================================
$Id$
*/
#include "wavelet.h"
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#define COEFSIZE (3 * sizeof(Coeff))
#define min(a, b) ((a) < (b) ? (a) : (b))
/* Constants for wavelet coeffs */
#define SQRT2 1.41421356
#define SQRT3 1.73205081
#define H2NORM (0.25 / SQRT2)
/* 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 */
/* Daubechies D1 (Haar) wavelet coeffs */
static const Coeff c0 = (1 / SQRT2);
/* Daubechies D2 wavelet coeffs */
static const Coeff h [4] = {
H2NORM * (1 + SQRT3), H2NORM * (3 + SQRT3),
H2NORM * (3 - SQRT3), H2NORM * (1 - SQRT3)
},
g [4] = {
H2NORM * (1 - SQRT3), -H2NORM * (3 - SQRT3),
H2NORM * (3 + SQRT3), -H2NORM * (1 + SQRT3)
}, /* = { h [3], -h [2], h [1], -h [0] } */
/* Summed coeffs for padding at left boundary */
hBound = H2NORM * (7 + SQRT3),
/* = h [0] + h [1] + h [2] */
gBound = H2NORM * (1 + SQRT3);
/* = g [0] + g [1] + g [2] */
static int tupleCmp (const void *t0, const void *t1)
/* Comparison callback for sorting 3-tuples;
reduce by summation, assuming unsigned. */
{
int i;
Coeff sum0 = 0, sum1 = 0, diff;
for (i = 0; i < 3; i++) {
sum0 += (*(Coeff **)t0) [i];
sum1 += (*(Coeff **)t1) [i];
}
diff = sum0 - sum1;
return diff > 0 ? 1 : diff < 0 ? -1 : 0;
}
void waveletPresort (Coeff (*x) [3], int xIdx [], int xLen)
/* Presort x to obtain a symmetric ("gaussian") distribution. This minimises
discontinuities in the signal, thus minimising the resulting detail
coefficients over the range of the signal. This essentially optimises
the compression.
Returns resorted signal in x and the permuted indices in xIdx to
subsequently reconstruct the original order via waveletPostsort().*/
{
int i, j;
Coeff *xp [xLen], *yp [xLen], y [xLen][3];
/* Initialise array of pointers to sort x via indirection, thus permuting
the pointers rather than the signal itself. */
for (i = 0; i < xLen; i++)
xp [i] = (Coeff *)&x [i];
/* Sort pointers in xp by the values they refer to
(see deref in tupleCmp() above) */
qsort(xp, xLen, sizeof(Coeff *), tupleCmp);
/* Redistribute xp symmetrically from both boundaries towards
centre. Chuck into yp.*/
for (i = 0; i < xLen - 1; i += 2) {
j = i >> 1;
yp [j] = xp [i];
yp [xLen - j - 1] = xp [i + 1];
}
if (xLen & 1)
/* Last partial iteration if odd size */
yp [(i >> 1) + 1] = xp [i];
/* Transfer values referenced by yp into a temp array (to prevent
clobbering x and SEGFAULTing), and transfer their
corresponding indices into xIdx. */
for (i = 0; i < xLen; i++) {
memcpy(y [i], yp [i], COEFSIZE);
/* XXX: Assume tuples in x are allocated consecutively! */
xIdx [i] = (yp [i] - x [0]) / 3;
}
/* Transfer temp to x */
memcpy(x, y, xLen * COEFSIZE);
}
static void d1Step (Coeff (*x) [3], int len, Coeff (*y) [3])
/*
Single step (at fixed resolution) of a forward Daubechies D1 (Haar)
wavelet transform on array x of arbitrary size len containing 3-tuples.
Note this is a 1D transform, and the triplets per array element are _not_
decorrelated, but transformed independently.
The wavelet coefficients are returned in the preallocated array y.
*/
{
const unsigned hl = -(-len >> 1);
int i, hi, j;
if (len >= 2) {
/* First and intermediate iterations */
for (i = 0; i <= len - 2; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
/* Smooth */
y [hi][j] = c0 * (x [i][j] + x [i + 1][j]);
/* Detail */
y [hl + hi][j] = c0 * (x [i][j] - x [i + 1][j]);
}
}
if (len & 1) {
/* Odd length; pad with last value as additional approx term */
for (j = 0; j < 3; j++)
y [hl - 1][j] = x [len - 1][j];
}
#ifdef WAVELET_DBG
printf("FWD HAAR XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.3f\t-->\t%.3f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
}
}
static void d1InvStep (Coeff (*y) [3], int len, Coeff (*x) [3])
/*
Single step (at fixed resolution) of an inverse Daubechies D1 (Haar)
wavelet transform on coefficient array y of size len containing 3-tuples.
This is essentially reverses the forward transform above; see additional
comments there for implementation details.
The reconstructed signal is return in the preallocated array x.
*/
{
const unsigned hl = -(-len >> 1);
int i, hi, j;
if (len >= 2) {
for (i = 0; i <= len - 2; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
x [i][j] = c0 * (y [hi][j] + y [hl + hi][j]);
x [i + 1][j] = c0 * (y [hi][j] - y [hl + hi][j]);
}
}
if (len & 1) {
/* Odd length; pad additional value from last approx term */
for (j = 0; j < 3; j++)
x [len - 1][j] = y [hl - 1][j];
}
#ifdef WAVELET_DBG
printf("INV HAAR XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.3f\t-->\t%.3f\n", y [i][0], x [i][0]);
putchar('\n');
#endif
}
}
static unsigned d2Step (Coeff (*x) [3], int xLen, Coeff (*y) [3])
/*
Single step (at fixed resolution) of forward Daubechies D2 wavelet
transform on array x of arbitrary size xLen containing 3-tuples.
Note this is a 1D transform, and the triplets per array element
are _not_ decorrelated, but transformed independently.
The approximation and detail coefficients are returned in the
preallocated output array y. This array is larger than xLen as the
transform introduces padding coefficients to handle boundary effects.
This size is obtained as return value, and can be interrogated beforehand
without performing a transform by passing either x or y as NULL pointer.
*/
{
unsigned i, j, k, hi;
const unsigned yLen = xLen + 3 >> 1, /* Output length + boundary pad */
hy = yLen; /* Offset for detail coeffs */
if (x && y) {
if (xLen >= 2) {
#ifdef WAVELET_DBG
/* Zero output array to facilitate debugging */
memset(y, 0, 2 * yLen * COEFSIZE);
#endif
/* First iteration (i = -2); pad coeffs at left edge with constant
* boundary extension (note repeated x [0]) */
for (j = 0; j < 3; j++) {
/* Approximation */
y [0][j] = hBound * x [0][j] + h [3] * x [1][j];
/* Detail */
y [hy][j] = gBound * x [0][j] + g [3] * x [1][j];
}
/* Intermediate and last iterations (i >= 0) */
for (i = 0; i < xLen; i += 2) {
hi = (i >> 1) + 1;
for (j = 0; j < 3; j++) {
y [hi][j] = y [hy + hi][j] = 0;
/* Convolve up to wavelet or input boundary, whichever comes
* first */
for (k = i; k < min(i + 4, xLen); k++) {
/* Approximation */
y [hi][j] += h [k - i] * x [k][j];
/* Detail */
y [hy + hi][j] += g [k - i] * x [k][j];
}
/* If necessary, convolve beyond input boundary with constant
* extension (note repeated x [xLen - 1]). This is skipped if
* the wavelet boundary was already reached */
for (; k < i + 4; k++) {
/* Approximation */
y [hi][j] += h [k - i] * x [xLen - 1][j];
/* Detail */
y [hy + hi][j] += g [k - i] * x [xLen - 1][j];
}
}
}
#ifdef WAVELET_DBG
printf(
"D2 XFORM + boundary padding (len = %d -> %d)\n",
xLen, yLen << 1
);
for (i = 0; i < yLen << 1; i++) {
if (i < xLen)
printf("%.3f\t-->\t", x [i][0]);
else
printf("\t\t");
printf("%.3f\n", y [i][0]);
}
putchar('\n');
#endif
}
}
return yLen;
}
static unsigned d2InvStep (Coeff (*y) [3], int yLen, Coeff (*x) [3])
/*
Single step (at fixed resolution) of inverse Daubechies D2 wavelet
transform on coefficient array y of size yLen (i.e. yLen/2 approx and
detail coeffs each) containing 3-tuples. This essentially reverses the
forward transform above; see additional comments there for implementation
details.
The reconstructed signal is returned in the preallocated array x. This
array is smaller than yLen as boundary padding coefficients are collapsed
in the inverse transform. This size is obtained as returned value, and
can be interrogated beforehand without performing a transform by passing
either x or y as NULL pointer.
*/
{
unsigned i, j, k, hi,
xLen = yLen - 2; /* Max size of output array x */
const unsigned hy = yLen >> 1; /* Offset for detail coeffs */
if (x && y) {
if (yLen >= 2) {
/* Reconstruct from approximation and detail coefficients */
for (i = 0; i < xLen - 1; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
x [i][j] =
h [2] * y [hi][j] + g [2] * y [hy + hi][j] +
h [0] * y [hi + 1][j] + g [0] * y [hy + hi + 1][j];
x [i + 1][j] =
h [3] * y [hi][j] + g [3] * y [hy + hi][j] +
h [1] * y [hi + 1][j] + g [1] * y [hy + hi + 1][j];
}
}
xLen = i;
#ifdef WAVELET_DBG
printf("D2 INV XFORM (len = %d --> %d)\n", yLen, xLen);
for (i = 0; i < yLen; i++) {
printf("%.3f\t", y [i][0]);
if (i < xLen)
printf("-->\t%.3f\n", x [i][0]);
else putchar('\n');
}
putchar('\n');
#endif
}
}
return xLen;
}
unsigned waveletXform (Coeff (*x) [3], unsigned xLen, Coeff (*y) [3])
/* Perform full multiresolution forward wavelet transform
on array x of length xLen containing 3-tuples. Note no intra-tuple
transform occurs (1-D only).
The wavelet coefficients are returned in preallocated array y, containing
the coarsest approximation in y [0] followed by details in order of
increasing resolution/frequency.
Note y must accommodate more than xLen coefficients to accommodate
padding for additional boundary coefficients. The number of coefficients
can be interrogated beforehand by calling the function with either x or y
set to NULL.
Note also that x is used as a buffer for intermediate results, and is
therefore modified.
The returned value is the number of output coefficients in y, or 0 if the
transform failed, e.g. if the input is too short.
*/
{
/* Output size (number of coeffs, including boundary padding), output
* array offset and pointers */
unsigned numCoeffs = 0, yOff, yLen, l;
Coeff (*yOffPtr) [3] = NULL;
/* Skip transform if input too short */
if (xLen > 3) {
/* Figure out total number of coefficients (= minimum size of output
* array including boundary padding). */
for (
numCoeffs = 0, l = xLen;
l > 3;
numCoeffs += l = d2Step(NULL, l, NULL)
);
numCoeffs += l;
if (x && y) {
/* Apply Daubechies D2 transform, accounting for boundary padding */
for (l = xLen, yOff = numCoeffs; l > 3; l = yLen) {
/* Get num coeffs and calculate offset into output array
* (subtracting num coeffs once more to account for approx
* coeffs, which will be replaced in next iteration) */
yOff -= yLen = d2Step(NULL, l, NULL);
yOffPtr = y + yOff - yLen;
/* Resulting wavelet coefficients in y are prepended to those
* from previous iteration; consequently, returned num coeffs is
* input size for next iteration */
d2Step(x, l, yOffPtr);
/* Copy approx coeffs from y into x as input for next iteration */
memcpy(x, yOffPtr, yLen * COEFSIZE);
}
/* Apply Daubechies D1 (Haar) transform in penultimate and final
* (since len is odd) iterations to decompose to a single
* approximation and detail coefficient at coarsest resolution at
* y [0] and y [1], respectively. Obviously no boundary padding
* required here. */
for (; l >= 2; l--) {
d1Step(x, l, y);
/* Copy coeffs from y into x as input for next iteration */
if (l > 2)
memcpy(x, y, l * COEFSIZE);
}
}
}
/* All coeffs now in y; return their number */
return numCoeffs;
}
unsigned waveletInvXform (
Coeff (*y) [3], unsigned yLen, Coeff (*x) [3], unsigned xLen0
)
/* Perform full multiresolution inverse wavelet transform on array y of
* length yLen containing wavelet coefficients as 3-tuples. Note no
* intra-tuple transform occurs (1-D only). The original signal length xLen0
* is required for the reconstruction.
*
* The reconstructed signal is returned in the preallocated array x. This
* is shorter than yLen since the additional boundary coefficients in y are
* collapsed. Note y is used as buffer for intermediate results, and is
* therefore modified.
*
* The returned value is the reconstructed signal length, which should match
* xLen, or 0 if the transform failed.
*/
{
unsigned xLen = 0, yOff, i, j, l, nextl, nIter;
/* The temporary output buffer xTmp is passed to d2InvStep() for
* intermediate xform steps. The additional space is required if xLen0
* is odd, as an additional coefficient is reconstructed (which is
* subsequently dropped). If x were to be passed directly to
* d2InvStep() under these circumstances, it would overflow. */
Coeff (*yOffPtr) [3] = NULL, xTmp [xLen0 + 1][3];
/* Figure out number of iterations based on original signal length */
for (nIter = 0, l = xLen0; l > 3; nIter++, l = d2Step(NULL, l, NULL));
/* Invert Daubechies D1 (Haar) transform at two coarsest levels */
for (l = 2; l <= 3; l++) {
d1InvStep(y, l, x);
/* Copy reconstructed signal into y for next iteration */
memcpy(y, x, l * COEFSIZE);
}
/* Invert Daubechies D2 transform */
for (i = nIter; i; i--) {
/* Figure out length of current detail coeffs based on number of
* iterations */
for (
j = 0, yOff = yLen, l = nextl = xLen0;
j < i;
j++, nextl = l, yOff -= l = d2Step(NULL, l, NULL)
);
#ifdef WAVELET_DBG
/* Zero output buffer to facilitate debugging */
memset(xTmp, 0, sizeof(xTmp));
#endif
yOffPtr = y + yOff;
/* Drop excess reconstructed approximation coeffs to match number of
* detail coeffs in next iteration */
xLen = min(d2InvStep(yOffPtr - l, l << 1, xTmp), nextl);
/* Copy reconstructed coeffs/signal from output buffer xTmp into y for
* next iteration, or into x after last iteration */
memcpy(i > 1 ? yOffPtr + l - xLen : x, xTmp, xLen * COEFSIZE);
}
return xLen;
}
#ifdef WAVELET_TEST_1D
#include <stdio.h>
#define ERRSTRLEN 1024
int main (int argc, char *argv [])
{
int i, j, xLen0, yLen, xLen, *xIdx;
Coeff (*xOrg) [3], (*x) [3], (*y) [3] = NULL;
char errstr [ERRSTRLEN];
if (argc < 2) {
fputs("Missing array length\n", stderr);
return -1;
}
if (!(xLen0 = atoi(argv [1]))) {
fputs("Invalid array length\n", stderr);
return -1;
}
/* Allocate arrays for original and reconstruction (num coeffs is
* obtained from waveletXform() before allocating output array) */
if (
!(xOrg = calloc(xLen0, COEFSIZE)) ||
!(xIdx = calloc(xLen0, sizeof(int))) ||
!(x = calloc(xLen0, COEFSIZE)) ||
!(y = calloc(waveletXform(NULL, xLen0, NULL), COEFSIZE))
) {
fputs("Failed array allocation\n", stderr);
return -1;
}
/* Init with random data */
srand48(111);
for (i = 0; i < xLen0; i++) {
#if 1
for (j = 0; j < 3; j++)
xOrg [i][j] = drand48();
#else
xOrg [i][0] = xOrg [i][1] = xOrg [i][2] = drand48();
#endif
}
/* Optimise signal distribution */
waveletPresort(xOrg, xIdx, xLen0);
#ifdef WAVELET_DBG
/* Dump presorted signal */
puts("-----------------------------------------------------------");
puts("Presorted Original & Indices");
for (i = 0; i < xLen0; i++)
printf(
"%.3f\t%d\n",
xOrg [i][0] + xOrg [i][1] + xOrg [i][2],
xIdx [i]
);
puts("-----------------------------------------------------------\n");
#endif
/* Transform and reconstruct */
memcpy(x, xOrg, xLen0 * COEFSIZE);
if (!(yLen = waveletXform(x, xLen0, y))) {
puts("Forward xform failed\n");
return -1;
}
#ifdef WAVELET_DBG
/* Dump coefficients */
puts("-----------------------------------------------------------\n");
printf("%d coefficients:\n\n", yLen);
for (i = 0; i < yLen; i++)
printf("%.3f\t%.3f\t%.3f\n", y [i][0], y [i][1], y [i][2]);
puts("\n-----------------------------------------------------------\n");
#endif
if (!(xLen = waveletInvXform(y, yLen, x, xLen0))) {
puts("Inverse xform failed\n");
return -1;
}
if (xLen != xLen0) {
printf(
"Original length (%d) != reconstructed length (%d)\n",
xLen0, xLen
);
return -1;
}
#ifdef WAVELET_DBG
puts("-----------------------------------------------------------\n");
#endif
puts("ORIG\t vs.\tINV XFORM");
for (i = 0; i < xLen; i++)
#if 1
printf(
"%.3f\t%.3f\t%.3f\t\t%.3f\t%.3f\t%.3f\n",
xOrg [i][0], xOrg [i][1], xOrg [i][2],
x [i][0], x [i][1], x [i][2]
);
#else
printf("%.3f\t\t%.3f\n", xOrg [i][0], x [i][0]);
#endif
free(xOrg);
free(xIdx);
free(x);
free(y);
return 0;
}
#endif /* WAVELET_TEST_1D */
c4science · Help