Homec4science
Diffusion RADIANCE Photon Map (master) Browse

RadiancePmap/wavelet.c
3cb2be88ebbcmaster

Last Change
Raw File

wavelet.c
View Options

/*
=========================================================================
Lifting scheme implementation of a Daubechies D4 wavelet transform on
arbitrary sized arrays (i.e. not restricted to powers of 2) of 3-tuples.
Array boundaries are extended as constants.
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 <stdio.h>
#define COEFSIZE (3 * sizeof(double))
/* Daubechies lifting coefficients; note c3 * c4 = 1 */
static const double c0 = sqrt(3),
c1 = c0 / 4,
c2 = c1 - 0.5,
norm = 1 / sqrt(2),
c3 = norm * (c0 - 1),
c4 = norm * (c0 + 1);
static void d4WaveletXformStep (double x [][3], int len)
/* Single step (at fixed resolution) of a forward Daubechies D4 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.
This implementation uses the lifting scheme to reduce the number of FLOPS
by leveraging redundancies inherent in the transform compared to a
"naive" algebraic implementation using convolution. However, it is not
performed in-place, but rather permutes coefficients on-the-fly to an
auxiliary array y for efficiency, as shuffling data within large arrays
is expected to be more expensive than FLOPS. */
{
/* Aux array, update/predict for current/previous/next iteration */
double y [len][3], updt [3], nextUpdt [3], pred [3], nextPred [3];
int i, j;
/* Init update/predict for previous and current iteration;
do so for each 3-tuple item */
for (j = 0; j < 3; j++) {
/* Update for iteration -1; assume constant signal beyond boundary
as substitute for: x [-2][j] + c0 * x [-1][j] */
/* First update (= average, low pass) */
updt [j] = x [0][j] + c0 * x [1][j];
/* First predict (= detail, high/band pass) */
pred [j] = x [1][j] - (c1 + c2) * updt [j];
}
for (i = 0; i < len - 3; i += 2) {
/* Iterate over triplet */
for (j = 0; j < 3; j++) {
/* Next update/predict (= average/detail = low/high pass) */
nextUpdt [j] = x [i + 2][j] + c0 * x [i + 3][j];
nextPred [j] = x [i + 3][j] - c1 * nextUpdt [j] - c2 * updt [j];
/* Current 2nd update (= average = low pass);
normalise and store in lower half of y */
y [i >> 1][j] = c3 * (updt [j] - nextPred [j]);
/* Current predict (= detail = high pass);
normalise and store in upper half of y */
y [i + len >> 1][j] = c4 * pred [j];
/* Prepare next iteration */
updt [j] = nextUpdt [j];
pred [j] = nextPred [j];
}
}
/* Last iteration with boundary handling */
for (j = 0; j < 3; j++) {
y [i >> 1][j] = c3 * (updt [j] - 0);
y [i + len >> 1][j] = c4 * pred [j];
}
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
static void d4WaveletInvXformStep (double x [][3], int len)
/* Single step (at fixed resolution) of an inverse Daubechies D4 wavelet
transform on coefficient array xIn of arbitrary size len containing
3-tuples. This is essentially reverses the forward transform above;
see additional comments there for implementation details. */
{
/* Aux array, update/predict for current/previous/next iteration */
double y [len][3], updt [3], prevUpdt [3], nextUpdt [3],
pred [3], nextPred [3];
int i, j;
/* Init update/predict for previous and current iteration;
do so for each triplet item */
for (j = 0; j < 3; j++) {
/* Invert first normalisation (by swapping reciprocal coeffs) and
predict/update from avg/diff coeffs in lower/upper half of x */
pred [j] = c3 * x [len >> 1][j];
updt [j] = c4 * x [0][j];
/* Invert update for iteration -1; assume constant signal beyond
boundary as substitute for: c4 * x [-1][j] + pred [j] */
prevUpdt [j] = updt [j] + pred [j];
}
for (i = 0; i < len - 3; i += 2) {
/* Iterate over triplet */
for (j = 0; j < 3; j++) {
/* Invert next predict/update from upper/lower half of x */
nextPred [j] = c3 * x [(i + len >> 1) + 1][j];
nextUpdt [j] = c4 * x [(i >> 1) + 1][j];
/* Invert current 2nd update */
updt [j] += nextPred [j];
/* Invert current predict, store linearly in y */
y [i + 1][j] = pred [j] + c1 * updt [j] + c2 * prevUpdt [j];
/* Invert current 1st update, store linearly in y */
y [i][j] = updt [j] - c0 * y [i + 1][j];
/* Prepare next iteration */
prevUpdt [j] = updt [j];
updt [j] = nextUpdt [j];
pred [j] = nextPred [j];
}
}
/* Last iteration with boundary handling */
for (j = 0; j < 3; j++) {
y [i + 1][j] = pred [j] + c1 * updt [j] + c2 * prevUpdt [j];
y [i][j] = updt [j] - c0 * y [i + 1][j];
}
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
void d4WaveletXform (double x [][3], unsigned len)
/* Perform full multiresolution forward Daubechies D4 wavelet transform
on array x of length len containing 3-tuples. Note no intra-tuple
transform occurs (1-D only).
Returns wavelet coefficients in x, containing average in x [0],
followed by details in order or increasing resolution. */
{
unsigned l;
for (l = len; l >= 2; l >>= 1)
d4WaveletXformStep(x, l);
}
void d4WaveletInvXform (double x [][3], unsigned len)
/* Perform full multiresolution inverse Daubechies D4 wavelet transform
on array x of length len containing 3-tuples. Note no intra-tuple
transform occurs (1-D only). Returns reconstructed signal in x. */
{
unsigned l;
for (l = 2; l <= len; l <<= 1)
d4WaveletInvXformStep(x, l);
}
#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 sum = 0;
for (i = 0; i < 3; i++)
sum += ((double *)t0) [i] - ((double *)t1) [i];
return sum > 0 ? 1 : sum < 0 ? -1 : 0;
}
int main (int argc, char *argv [])
{
unsigned i, j, len;
double (*x0) [3], (*x1) [3];
char errstr [ERRSTRLEN];
if (argc < 2) {
fputs("Missing array length\n", stderr);
return -1;
}
if (!(len = atoi(argv [1]))) {
fputs("Invalid array length\n", stderr);
return -1;
}
/* Allocate arrays for original and reconstruction */
if (!(x0 = calloc(len, COEFSIZE)) || !(x1 = calloc(len, COEFSIZE))) {
fputs("Failed array allocation\n", stderr);
return -1;
}
/* Init with random data, sort by tuple sum */
srand48(111);
for (i = 0; i < len; i++) {
#if 1
x0 [i][0] = x0 [i][1] = x0 [i][2] = drand48();
#else
for (j = 0; j < 3; j++)
x0 [i][j] = drand48();
#endif
}
qsort(x0, len, COEFSIZE, tupleCmp);
memcpy(x1, x0, (size_t)len * COEFSIZE);
/* Transform and reconstruct data */
d4WaveletXform(x1, len);
d4WaveletInvXform(x1, len);
for (i = 0; i < len; i++)
#if 0
printf(
"%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n",
x0 [i][0], x0 [i][1], x0 [i][2],
x1 [i][0], x1 [i][1], x1 [i][2]
);
#else
printf("%.6f\t%.6f\n", x0 [i][0], x1 [i][0]);
#endif
free(x0);
free(x1);
return 0;
}
#endif
c4science · Help