Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/wavelet.c
73332c190648pmapmerge

Last Change
Raw File

wavelet.c
View Options

/*
=========================================================================
Lifting scheme implementation of Daubechies D4 and Haar wavelet
transforms on arbitrary sized arrays (i.e. not restricted to powers of 2)
of 3-tuples. Array boundaries are extended cyclically.
The implementation here is mainly distilled from the following sources:
-- 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 code:
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 "wavelet.h"
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#define COEFSIZE (3 * sizeof(double))
#define min(a, b) ((a) < (b) ? (a) : (b))
/* Daubechies D4 wavelet coeffs */
static const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3),
h [4] = {
norm * (1 + sqrt3), norm * (3 + sqrt3),
norm * (3 - sqrt3), norm * (1 - sqrt3)
},
g [4] = { h [3], -h [2], h [1], -h [0] },
/* Summed coeffs for padding at left boundary */
hBound = h [0] + h [1] + h [2],
gBound = g [0] + g [1] + g [2];
/* Haar wavelet coeffs */
static const double c0 = 1 / sqrt(2);
static void haarStep (double (*x) [3], int len, double (*y) [3])
/*
Single step (at fixed resolution) of a forward 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 haarInvStep (double (*y) [3], int len, double (*x) [3])
/*
Single step (at fixed resolution) of an inverse 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 d4Step (double (*x) [3], int xLen, double (*y) [3])
/*
Single step (at fixed resolution) of forward Daubechies D4 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(
"D4 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 d4InvStep (double (*y) [3], int yLen, double (*x) [3])
/*
Single step (at fixed resolution) of inverse Daubechies D4 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("D4 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 (double (*x) [3], unsigned xLen, double (*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;
double (*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 = d4Step(NULL, l, NULL)
);
numCoeffs += l;
if (x && y) {
/* Apply Daubechies D4 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 = d4Step(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 */
d4Step(x, l, yOffPtr);
/* Copy approx coeffs from y into x as input for next iteration */
memcpy(x, yOffPtr, yLen * COEFSIZE);
}
/* Apply 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--) {
haarStep(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 (
double (*y) [3], unsigned yLen, double (*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 d4InvStep() 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
* d4InvStep() under these circumstances, it would overflow. */
double (*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 = d4Step(NULL, l, NULL));
/* Invert Haar transform at two coarsest levels */
for (l = 2; l <= 3; l++) {
haarInvStep(y, l, x);
/* Copy reconstructed signal into y for next iteration */
memcpy(y, x, l * COEFSIZE);
}
/* Invert Daubechies D4 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 = d4Step(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(d4InvStep(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
#include <stdio.h>
#include <stdlib.h>
#define ERRSTRLEN 1024
static int tupleCmp (const void *t0, const void *t1)
/* Comparison callback for sorting 3-tuples; reduce by summation */
{
int i;
double sum0 = 0, sum1 = 0, diff;
for (i = 0; i < 3; i++) {
sum0 += ((double *)t0) [i];
sum1 += ((double *)t1) [i];
}
diff = sum0 - sum1;
return diff > 0 ? 1 : diff < 0 ? -1 : 0;
}
int main (int argc, char *argv [])
{
int i, j, xLen0, yLen, xLen;
double (*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)) ||
!(x = calloc(xLen0, COEFSIZE)) ||
!(y = calloc(waveletXform(NULL, xLen0, NULL), COEFSIZE))
) {
fputs("Failed array allocation\n", stderr);
return -1;
}
/* Init with random data, sort by tuple sum */
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
}
qsort(xOrg, xLen0, COEFSIZE, tupleCmp);
memcpy(x, xOrg, (size_t)xLen0 * COEFSIZE);
/* Transform and reconstruct data */
if (!(yLen = waveletXform(x, xLen0, y))) {
puts("Forward xform failed\n");
return -1;
}
#ifdef _WAVELET_DBG
puts("-----------------------------------------------------------\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(x);
free(y);
return 0;
}
#endif
c4science · Help