Homec4science
Diffusion RADIANCE Photon Map (precompContrib) Browse

RadiancePmap/wavelet.c
70bd99e96f17precompContrib

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 d4LiftStep (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;
/* Haar wavelet for 1st iteration to handle boundary effects */
for (j = 0; j < 3; j++) {
/* Predict (= detail = high pass) */
y [len >> 1][j] = pred [j] = x [1][j] - x [0][j];
/* Update (= average = low pass) */
y [0][j] = updt [j] = x [0][j] + 0.5 * pred [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];
if (i > 0) {
/* D4 wavelet for subsequent iterations */
/* 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];
}
}
if (i < len) {
/* Haar wavelet for last iteration to handle boundary effects */
for (j = 0; j < 3; j++) {
/* Predict (= detail = high pass) */
y [i + len >> 1][j] = pred [j] = x [i + 1][j] - x [i][j];
/* Update (= average = low pass) */
y [i >> 1][j] = updt [j] = x [i][j] + 0.5 * pred [j];
}
}
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
static void d4InvLiftStep (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;
/* Haar wavelet for 1st iteration to handle boundary effects */
for (j = 0; j < 3; j++) {
/* Invert predict/update from upper/lower half of x */
pred [j] = x [len >> 1][j];
updt [j] = x [0][j];
y [0][j] = updt [j] - 0.5 * pred [j];
y [1][j] = pred [j] + y [0][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];
if (i > 0) {
/* D4 wavelet for subsequent iterations */
/* 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];
}
}
if (i < len) {
/* Haar wavelet for last iteration to handle boundary effects */
for (j = 0; j < 3; j++) {
/* Invert predict/update from upper/lower half of x */
pred [j] = x [i + len >> 1][j];
updt [j] = x [i >> 1][j];
y [i][j] = updt [j] - 0.5 * pred [j];
y [i + 1][j] = pred [j] + y [i][j];
}
}
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
static void haarConvStep (double x [][3], int len)
/*
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.
*/
{
/* Aux array, coeffs */
double y [len][3];
const double c0 = 1 / sqrt(2);
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 smoothing 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("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
static void haarInvConvStep (double x [][3], int len)
/*
Single step (at fixed resolution) of an inverse Haar 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, coeffs */
double y [len][3];
const double c0 = 1 / sqrt(2);
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++) {
y [i][j] = c0 * (x [hi][j] + x [hl + hi][j]);
y [i + 1][j] = c0 * (x [hi][j] - x [hl + hi][j]);
}
}
if (len & 1) {
/* Odd length; pad additional value from last smoothing term */
for (j = 0; j < 3; j++)
y [len - 1][j] = x [hl - 1][j];
}
#ifdef _WAVELET_DBG
printf("INV HAAR XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
static void qHaarConvStep (double x [][3], int len)
/*
Single step (at fixed resolution) of a forward quadratic 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.
*/
{
/* Aux array, coeffs */
double y [len][3];
const double h0 = 0.5, g0 = sqrt(3) * h0;
const unsigned hl = -(-len >> 1);
int i, hi, j;
if (len >= 3) {
/* First and intermediate iterations */
for (i = 0; i <= len - 3; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
/* Smooth */
y [hi][j] =
h0 * (x [i][j] + x [i + 1][j] + x [i + 2][j]);
/* Detail */
y [hl + hi][j] =
g0 * x [i + 1][j] - h0 * (x [i][j] + x [i + 2][j]);
}
}
if (i <= len - 2) {
/* Last iteration with boundary extension */
hi = i >> 1;
for (j = 0; j < 3; j++) {
/* Smooth */
y [hi][j] =
h0 * (x [i][j] + x [i + 1][j] + x [0][j]);
/* Detail */
y [hl + hi][j] =
g0 * x [i + 1][j] - h0 * (x [i][j] + x [0][j]);
}
}
if (len & 1) {
/* Odd length; pad with last value as additional smoothing term */
for (j = 0; j < 3; j++)
y [hl - 1][j] = x [len - 1][j];
}
#ifdef _WAVELET_DBG
printf("FWD QHAAR XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
static void qHaarInvConvStep (double x [][3], int len)
/*
Single step (at fixed resolution) of an inverse quadratic Haar 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, coeffs */
double y [len][3];
const double h0 = 0.5, g0 = sqrt(3) * h0;
const unsigned hl = -(-len >> 1);
int i, hi, j;
if (len >= 3) {
/* First iteration with boundary extension */
hi = len >> 1;
for (j = 0; j < 3; j++) {
y [0][j] = h0 * (
x [0][j] - x [hl][j] +
x [hi - 1][j] - x [len - 1][j]
);
y [1][j] = h0 * x [0][j] + g0 * x [hl][j];
}
/* Intermediate and last iterations */
for (i = 0; i <= len - 4; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
y [i + 2][j] = h0 * (
x [hi][j] - x [hl + hi][j] +
x [hi + 1][j] - x [hl + hi + 1][j]
);
y [i + 3][j] = h0 * x [hi][j] + g0 * x [hl + hi][j];
}
}
if (len & 1) {
/* Odd length; pad additional value from last smoothing term */
for (j = 0; j < 3; j++)
y [len - 1][j] = x [hl - 1][j];
}
#ifdef _WAVELET_DBG
printf("INV QHAAR XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
static void d4ConvStep (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.
*/
{
/* Aux array, coeffs */
double y [len][3];
const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3),
h0 = norm * (1 + sqrt3), h1 = norm * (3 + sqrt3),
h2 = norm * (3 - sqrt3), h3 = norm * (1 - sqrt3),
g0 = h3, g1 = -h2, g2 = h1, g3 = -h0;
int i, hi, j;
const unsigned hl = -(-len >> 1);
if (len >= 4) {
/* First and intermediate iterations */
for (i = 0; i <= len - 4; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
/* Smooth */
y [hi][j] =
h0 * x [i][j] + h1 * x [i + 1][j] +
h2 * x [i + 2][j] + h3 * x [i + 3][j];
/* Detail */
y [hl + hi][j] =
g0 * x [i][j] + g1 * x [i + 1][j] +
g2 * x [i + 2][j] + g3 * x [i + 3][j];
}
}
if (i <= len - 2) {
/* Last iteration with wrapped boundary extension */
hi = i >> 1;
for (j = 0; j < 3; j++) {
/* Smooth */
y [hi][j] =
h0 * x [i][j] + h1 * x [i + 1][j] +
h2 * x [0][j] + h3 * x [1][j];
/* Detail */
y [hl + hi][j] =
g0 * x [i][j] + g1 * x [i + 1][j] +
g2 * x [0][j] + g3 * x [1][j];
}
}
if (len & 1) {
/* Odd length; pad with last value as additional smoothing term */
for (j = 0; j < 3; j++)
y [hl - 1][j] = x [len - 1][j];
}
#ifdef _WAVELET_DBG
printf("FWD D4 XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
static void d4InvConvStep (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, coeffs */
double y [len][3];
const double norm = 0.25 / sqrt(2), sqrt3 = sqrt(3),
h0 = norm * (1 + sqrt3), h1 = norm * (3 + sqrt3),
h2 = norm * (3 - sqrt3), h3 = norm * (1 - sqrt3),
g0 = h3, g1 = -h2, g2 = h1, g3 = -h0;
int i, hi, j;
const unsigned hl = -(-len >> 1);
if (len >= 4) {
/* First iteration with wrapped boundary extension */
hi = len >> 1;
for (j = 0; j < 3; j++) {
y [0][j] =
h0 * x [0][j] + g0 * x [hl][j] +
h2 * x [hi - 1][j] + g2 * x [len - 1][j];
y [1][j] =
h1 * x [0][j] + g1 * x [hl][j] +
h3 * x [hi - 1][j] + g3 * x [len - 1][j];
}
/* Intermediate and last iterations */
for (i = 0; i <= len - 4; i += 2) {
hi = i >> 1;
for (j = 0; j < 3; j++) {
y [i + 2][j] =
h2 * x [hi][j] + g2 * x [hl + hi][j] +
h0 * x [hi + 1][j] + g0 * x [hl + hi + 1][j];
y [i + 3][j] =
h3 * x [hi][j] + g3 * x [hl + hi][j] +
h1 * x [hi + 1][j] + g1 * x [hl + hi + 1][j];
}
}
if (len & 1) {
/* Odd length; pad additional value from last smoothing term */
for (j = 0; j < 3; j++)
y [len - 1][j] = x [hl - 1][j];
}
#ifdef _WAVELET_DBG
printf("INV D4 XFORM (len = %d)\n", len);
for (i = 0; i < len; i++)
printf("%.6f\t-->\t%.6f\n", x [i][0], y [i][0]);
putchar('\n');
#endif
/* Return coeffs in x */
memcpy(x, y, len * COEFSIZE);
}
}
void waveletXform (double x [][3], unsigned len)
/* Perform full multiresolution forward 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. */
{
int l, j;
/* Shifting negated length leverages the fact that shift always rounds
* down; as a result, odd lengths are rounded *up* to account for extra
* low frequency coefficient (cf. Stefan Moebius' "TurboWavelets") */
for (l = len; l >= 2; l = -(-l >> 1)) {
if (l >= 4)
d4ConvStep(x, l);
else
/* Apply Haar xform in final (and penultimate, if is len odd)
iteration to a obtain single smoothing an detail term */
haarConvStep(x, l);
}
}
void waveletInvXform (double x [][3], unsigned len)
/* Perform full multiresolution inverse 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. */
{
int l, i, j, log2len;
/* Get log2(len); reconstruct lengths used in forward xform, preserving
* odd lengths by signed shifting (which always rounds down) */
for (log2len = 0, l = len; l > 1; log2len++, l = -(-l >> 1));
for (i = 1; i <= log2len; i++) {
l = -(-(int)len >> (log2len - i));
if (l >= 4)
d4InvConvStep(x, l);
else
/* Inverse Haar xform in first (and second, if len is odd)
iteration to expand single smoothing and detail term */
haarInvConvStep(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 */
waveletXform(x1, len);
#ifdef _WAVELET_DBG
puts("-----------------------------------------------------------\n");
#endif
waveletInvXform(x1, len);
#ifdef _WAVELET_DBG
puts("-----------------------------------------------------------\n");
#endif
puts("ORIG\t vs.\tINV XFORM");
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