Homec4science
Diffusion RADIANCE Photon Map (master) Browse

RadiancePmap/wavelet3.c
1f0e61fbf860master

Last Change
Raw File

wavelet3.c
View Options

/*
=========================================================================
2-dimensional Daubechies D2 and D1 (Haar) wavelet transforms on arbitrary
sized arrays of 3-tuples. Array boundaries are extended as defined by
WAVELET_EXTEND_MODE, which introduces padding coefficients.
Compile with -DWAVELET_TEST_2DPAD 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 "wavelet2.h"
#include <math.h>
#include <string.h>
#include <stdlib.h>
static unsigned padD1Step2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned len
)
/*
Single iteration of forward 2D Daubechies D1 (Haar) wavelet transform
with padding on array y of arbitrary size (len x len) containing
3-tuples. This variant uses constant boundary extension. 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 array dimension len on success, else 0.
*/
{
static unsigned axis = 0;
const unsigned hLen = -(-(signed)len >> 1);
int h, i, j, k;
if (len < 2 || !y || !yt)
/* Input shorter than wavelet support, or no input/output */
return 0;
#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 - 2; 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 [h + hLen] [i] [k] = h1 * (
y [i] [j ] [k] -
y [i] [j + 1] [k]
);
}
}
if (len & 1) {
/* Odd length; pad with last value as additional approx term */
for (k = 0; k < 3; k++)
/* Smooth/approx/avg/lowpass */
yt [hLen - 1] [i] [k] = y [i] [len - 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');
#endif
axis ^= 1;
return len;
}
static unsigned padD1InvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned len
)
/*
Single iteration of inverse 2D Daubechies D1 (Haar) wavelet transform
with padding on array y of arbitrary size (len x len) containing
3-tuples. 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 array dimension len on success, else 0.
*/
{
static unsigned axis = 1;
const unsigned hLen = -(-(signed)len >> 1);
int h, i, j, k;
if (len < 2 || !y || !yt)
/* Too few coeffs for reconstruction, or no input/output */
return 0;
#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 - 2; j += 2) {
h = j >> 1;
for (k = 0; k < 3; k++) {
yt [i] [j ] [k] = h1 * (
y [h ] [i] [k] + /* Approx */
y [h + hLen] [i] [k] /* Detail */
);
yt [i] [j + 1] [k] = h1 * (
y [h ] [i] [k] - /* Approx */
y [h + hLen] [i] [k] /* Detail */
);
}
}
if (len & 1) {
/* Odd length; pad additional value from last approximation term */
for (k = 0; k < 3; k++)
yt [i] [len - 1] [k] = y [hLen - 1] [i] [k];
}
}
#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');
#endif
axis ^= 1;
return len;
}
static unsigned padD2Step2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned yLen, unsigned ytOffset
)
/*
Single iteration of forward 2D Daubechies D2 wavelet transform with
boundary padding on array y containing (yLen x yLen) 3-tuples. This
variant extends the input array beyond its boundary as defined by
WAVELET_EXTEND_MODE, which results in padding coefficients. 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 output wavelet coefficients are returned in the *TRANSPOSED* array
yt. 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.
Both the input array y and output array yt must be large enough to also
accommodate the additional padding coefficients introduced at the
boundaries. The number of output coefficient pairs ytLen returned in yt
is obtained as return value, and can be interrogated beforehand without
performing a transform by passing either y or yt as NULL.
The approximation coefficients are returned in yt [0..ytLen-1] [i], while
the details are returned in yt [ytOffset..ytOffset+ytLen-1] [i]. The
offset ytOffset must be supplied by the caller to accommodate the detail
coefficients (plus padding!) without clobbering the approximation
coefficients (again, plus padding!). Consequently, there is a gap
between the approximation coefficients in the upper left of the matrix,
and the detail coefficients in the lower right. This gap shrinks and
shifts toward the upper left with each iteration.
The following schema indicates the matrix contents during the transform
(see the unit test output when compiled with WAVELET_TEST_2DPAD):
----> j
+----------+--------------+ +-----------+-------------+
| | y_ij | | | S(y_ji) | |
| +----------+ | (H_k xform)^T +-----------+ |
v | | ----------> | |
i | | +-----------+ |
| | | D(y_ji) | |
+-------------------------+ +-------------------------+
----> i
+----------+--------------+ +----------+---+----------+
| | S(y_ji) | | |S(S(y_ij))| |S(D(y_ij))|
| +----------+ | (V_k xform)^T +----------+ +----------+
v | | ----------> | |
j +----------+ | +----------+ +----------+
| D(y_ji) | | |D(S(y_ij))| |D(D(y_ij))|
+-------------------------+ +----------+---+----------+
... ... ...
... ... ...
+------+------+--------+--------+----------------+
|SS_n-1|SD_n-1| | | |
|------+------+ SD_n-2 | | |
+DS_n-1|DD_n-1+--------+ . . . | SD_0 |
+------+----+-+--------+ | |
| DS_n-2 | | DD_n-2 | +----------------+
(V_n-1 xform)^T +-----------+ +--------+ . . . |
----------> | . . . . . |
| . . . . . |
+----------------+ +----------------+
| | | |
| DS_0 | . . . . . . | DD_0 |
| | | |
| | | |
+----------------+--------------+----------------+
where i, j are the array indices of the outer resp. inner loop, S(y_ij)
D(y_ij) are approximation resp. details coefficients of input y_ij
(itself either the original signal or approximation/detail coefficients
from previous iterations), and H_k, V_k are the horizontal resp.
vertical transforms of the k-th (out of the maximum n) iteration.
*/
{
static unsigned axis = 0;
unsigned h, i, j, k, l, i0 = 0, iN, wIdx, dExt;
/* Number of output coeff pairs ytLen, incl. boundary padding, and
max dimension of entire array */
const unsigned ytLen = (yLen + 3) >> 1, ytTotalLen = ytOffset + ytLen;
WAVELET_COEFF yExt;
#if (WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD1 || \
WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD2 \
)
/* Gradient boundary extension stuff */
WAVELET_COEFF dy1_0, dy1_1, dy2_0;
#endif
if (y && yt && yLen >= 2) {
#ifdef WAVELET_DBG
/* Clear output array yt to facilitate debugging */
clearCoeffs2(yt, ytTotalLen);
#else
/* Zero output (sub)array yt; note this will not affect accumulated
* coefficients from previous iterations, as they lie outside the
* array bounds */
zeroCoeffs2(yt, ytTotalLen);
#endif
/* FWD TRANSFORM GOTCHAS:
*
* Each full iteration of the forward transform consists of a
* horizontal step followed by a vertical one. 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.
*
* In the initial horizontal transform, the input array only contains
* a set of approximation coefficients (or the original signal) at
* y [0..yLen-1] [0..yLen-1].
*
* In the subsequent vertical transform, the input array then contains
* two sets of coefficients from the preceding horizontal transform:
* approximation coeffs y [0..ytLen-1] [0..ytLen-1] and detail coeffs
* y [ytOffset..ytTotalLen-1] [ytOffset..ytTotalLen-1]. In this case,
* the loop below is iterated twice with an offset yOffset to select
* both sets of coefficients to transform. Note the detail coeffs are
* selected first to facilitate the looping logic; if i0==0 at the end
* of the loop (i.e. approx coeffs were selected), it terminates.
*/
do {
i0 ^= ytOffset * axis;
iN = i0 ? ytTotalLen : yLen;
for (i = i0; i < iN; i++) {
/* First iteration (j = -2); pad coefficients at left edge and
* extend input beyond boundary according to extension mode */
for (k = 0; k < 3; k++) {
#if (WAVELET_EXTEND_MODE == WAVELET_EXTEND_CONST || \
!defined(WAVELET_EXTEND_MODE) \
)
/* Constant boundary extension (=default) */
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
(h2 [0] + h2 [1] + h2 [2]) * y [i] [0] [k] +
(h2 [3] ) * y [i] [1] [k];
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
(g2 [0] + g2 [1] + g2 [2]) * y [i] [0] [k] +
(g2 [3] ) * y [i] [1] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD1
/* 1st order gradient boundary extension */
dy1_0 = y [i] [0] [k] - y [i] [1] [k];
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
h2 [0] * (y [i] [0] [k] + 2 * dy1_0) +
h2 [1] * (y [i] [0] [k] + dy1_0) +
h2 [2] * (y [i] [0] [k] ) +
h2 [3] * (y [i] [1] [k] );
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
g2 [0] * (y [i] [0] [k] + 2 * dy1_0) +
g2 [1] * (y [i] [0] [k] + dy1_0) +
g2 [2] * (y [i] [0] [k] ) +
g2 [3] * (y [i] [1] [k] );
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD2
/* 2nd order gradient boundary extension */
dy1_0 = y [i] [0] [k] - y [i] [1] [k];
dy1_1 = y [i] [1] [k] - y [i] [2] [k];
dy2_0 = dy1_0 - dy1_1;
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
/* h2 [0] * (y [i] [0] [k] + 2 * (dy1_0 + 1 * dy2_0)) +
h2 [1] * (y [i] [0] [k] + 1 * (dy1_0 + 0 * dy2_0)) +
h2 [2] * (y [i] [0] [k] + 0 * (dy1_0 )) +
h2 [3] * (y [i] [1] [k] ); */
h2 [0] * (y [i] [0] [k] + 2 * (dy1_0 + dy2_0)) +
h2 [1] * (y [i] [0] [k] + dy1_0 ) +
h2 [2] * (y [i] [0] [k] ) +
h2 [3] * (y [i] [1] [k] );
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
/* g2 [0] * (y [i] [0] [k] + 2 * (dy1_0 + 1 * dy2_0)) +
g2 [1] * (y [i] [0] [k] + 1 * (dy1_0 + 0 * dy2_0)) +
g2 [2] * (y [i] [0] [k] + 0 * (dy1_0 )) +
g2 [3] * (y [i] [1] [k] ); */
g2 [0] * (y [i] [0] [k] + 2 * (dy1_0 + dy2_0)) +
g2 [1] * (y [i] [0] [k] + dy1_0 ) +
g2 [2] * (y [i] [0] [k] ) +
g2 [3] * (y [i] [1] [k] );
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_CYCL
/* Cyclic (perdiodic / wraparound) boundary extension */
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
h2 [0] * y [i] [yLen - 2] [k] +
h2 [1] * y [i] [yLen - 1] [k] +
h2 [2] * y [i] [0 ] [k] +
h2 [3] * y [i] [1 ] [k];
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
g2 [0] * y [i] [yLen - 2] [k] +
g2 [1] * y [i] [yLen - 1] [k] +
g2 [0] * y [i] [0 ] [k] +
g2 [3] * y [i] [1 ] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_SYMM
/* Symmetric (reflected) boundary extension */
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
h2 [0] * y [i] [1] [k] +
h2 [1] * y [i] [0] [k] +
h2 [2] * y [i] [0] [k] +
h2 [3] * y [i] [1] [k];
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
g2 [0] * y [i] [1] [k] +
g2 [1] * y [i] [0] [k] +
g2 [2] * y [i] [0] [k] +
g2 [3] * y [i] [1] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_ZERO
/* Constant zero boundary extension */
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
h2 [2] * y [i] [0] [k] + h2 [3] * y [i] [1] [k];
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
g2 [2] * y [i] [0] [k] + g2 [3] * y [i] [1] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_INV
/* Inverted (negative reflected) boundary extension */
/* Smooth/approx/avg/lowpass */
yt [0 ] [i] [k] =
-h2 [0] * y [i] [1] [k] +
-h2 [1] * y [i] [0] [k] +
h2 [2] * y [i] [0] [k] +
h2 [3] * y [i] [1] [k];
/* Detail/diff/highpass */
yt [ytOffset] [i] [k] =
-g2 [0] * y [i] [1] [k] +
-g2 [1] * y [i] [0] [k] +
g2 [2] * y [i] [0] [k] +
g2 [3] * y [i] [1] [k];
#endif
}
/* Intermediate and last iterations (j >= 0) */
for (j = 0; j < yLen; j += 2) {
h = (j >> 1) + 1;
for (k = 0; k < 3; k++) {
#ifdef WAVELET_DBG
/* Only necessary in debugging mode, since otherwise
already cleared by zeroCoeffs2() */
yt [h] [i] [k] = yt [h + ytOffset] [i] [k] = 0;
#endif
/* Convolve up to boundary of wavelet or input signal,
* whichever comes first */
for (l = j; l < min(j + 4, yLen); l++) {
/* Smooth/approx/avg/lowpass */
yt [h ] [i] [k] += h2 [l - j] * y [i] [l] [k];
/* Detail/diff/highpass */
yt [h + ytOffset] [i] [k] += g2 [l - j] * y [i] [l] [k];
}
/* If necessary, convolve beyond input signal boundary
* with extension. This is skipped if the wavelet boundary
* was already reached */
for (; l < j + 4; l++) {
wIdx = l - j;
#if (WAVELET_EXTEND_MODE == WAVELET_EXTEND_CONST || \
!defined(WAVELET_EXTEND_MODE) \
)
/* Constant boundary extension (=default) */
yExt = y [i] [yLen - 1] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD1
/* 1st order gradient boundary extension */
dy1_0 = y [i] [yLen - 1] [k] - y [i] [yLen - 2] [k];
dExt = l - yLen + 1;
yExt = y [i] [yLen - 1] [k] + dExt * dy1_0;
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_GRAD2
/* 2nd order gradient boundary extension */
dy1_0 = y [i] [yLen - 1] [k] - y [i] [yLen - 2] [k];
dy1_1 = y [i] [yLen - 2] [k] - y [i] [yLen - 3] [k];
dy2_0 = dy1_0 - dy1_1;
dExt = l - yLen + 1;
yExt = y [i] [yLen - 1] [k] +
dExt * (dy1_0 + (dExt ? dExt - 1 : 0) * dy2_0);
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_CYCL
/* Cyclic (perdiodic/wraparound) boundary extension */
yExt = y [i] [l - yLen] [k];
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_SYMM
/* Symmetric (reflected) boundary extension */
yExt = y [i] [((yLen - 1) << 1) - l + 1] [k];
/* = (yLen - 1) - (l - (yLen - 1)) + 1 */
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_ZERO
/* Constant zero boundary extension */
yExt = 0;
#elif WAVELET_EXTEND_MODE == WAVELET_EXTEND_INV
/* Inverted (negative reflected) boundary extension */
yExt = -y [i] [((yLen - 1) << 1) - l + 1] [k];
#endif
/* Smooth/approx/avg/lowpass */
yt [h ] [i] [k] += h2 [wIdx] * yExt;
/* Detail/diff/highpass */
yt [h + ytOffset] [i] [k] += g2 [wIdx] * yExt;
}
}
}
}
} while (i0);
#ifdef WAVELET_DBG
if (axis)
printf("VERT FWD D2: (%dx%d) H-approx, (%dx%d) H-detail\n"
"-> (%dx%d) V-approx of H-approx, (%dx%d) V-approx of H-detail,\n"
" (%dx%d) V-detail of H-approx, (%dx%d) V-detail of V-detail\n",
ytLen, yLen, ytLen, yLen,
ytLen, ytLen, ytLen, ytLen, ytLen, ytLen, ytLen, ytLen
);
else
printf(
"HORIZ FWD D2: (%dx%d) -> (%dx%d) H-approx, (%dx%d) H-detail\n",
yLen, yLen,
ytLen, yLen, ytLen, yLen
);
dumpCoeffs2(y, yt, ytTotalLen, ytTotalLen, 0, NULL);
putchar('\n');
#endif
axis ^= 1;
}
return ytLen;
}
static unsigned padD2InvStep2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned yLen, unsigned yOffset
)
/*
Single iteration of inverse 2D Daubechies D2 wavelet transform with
boundary padding on array y containing pairs of (yLen x yLen)
approximation and detail coefficients in 3-tuples. 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. 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.
The number of output coefficients is smaller than the input as padding
coefficients are collapsed at the boundaries during the inverse
transform. The number of output coefficients per dimension ytLen is
obtained as returned value, and can be interrogated beforehand without
performing a transform by passing either y or yt as NULL.
The inversion expects approximation coeffs y [0..yLen-1][0..yLen-1]
and y [0..yLen-1] [yOffset..yOffset+yLen-1], and detail coefficients
y [yOffset..yOffset+yLen-1] [0..yLen-1] and
y [yOffset..yOffset+yLen-1] [yOffset..yOffset+yLen-1].
The offset yOffset must be supplied by the caller to account for the fact
that the number of output coefficients ytLen is smaller than yLen due to
the collapsed padding coefficients. Consequently there is a gap between
the inverted coefficients located in the corners of the matrix. This gap
increases with each iteration, until all coefficients have been consumed
by the inversion.
*/
{
static unsigned axis = 1;
unsigned h, i, j, k, l, i0 = 0, iN;
/* Number of output coeff (pairs) ytLen, max dimension of entire array */
unsigned ytLen = (yLen - 1) << 1;
const unsigned yTotalLen = yOffset + yLen;
if (y && yt && ytLen >= 2) {
#ifdef WAVELET_DBG
/* Clear output array yt to facilitate debugging */
clearCoeffs2(yt, yTotalLen);
#else
/* Zero output (sub)array yt; note this will not affect accumulated
* coefficients from previous iterations, as they lie outside the
* array bounds */
zeroCoeffs2(yt, yTotalLen);
#endif
/* INV TRANSFORM GOTCHAS:
*
* i, j are swapped relative to the forward transform, as axis order
* is now reversed (vertical before horizontal xform). 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.
*
* In the inverse vertical transform, the input array contains two
* sets of coefficients:
* (a) approx of approx coeffs y [0..yLen-1] [0..yLen-1] in the
* upper left, and details of approximation coefficients
* y [yOffset..yTotalLen-1] [0..yLen-1] in the lower left.
* (b) approx of detail coeffs y [0..yLen-1] [yOffset..yTotalLen-1]
* in the upper right, and details of detail coefficients
* y [yOffset..yTotalLen-1] [yOffset..yTotalLen-1] in lower right.
*
* The loop below is iterated twice with an offset to select both sets
* of coefficients to inverse the transform. Note set (b) is selected
* first to facilitate the looping logic; if i0==0 at the end of the
* loop (i.e. set (a) was selected), it terminates.
*/
do {
i0 ^= yOffset * (axis);
iN = i0 ? yTotalLen : ytLen;
for (i = i0; i < iN; i++) {
for (j = 0; j < ytLen; j += 2) {
h = j >> 1;
for (k = 0; k < 3; k++) {
yt [i] [j ] [k] =
h2i [1] * y [h ] [i] [k] + /* Approx */
g2i [1] * y [h + yOffset ] [i] [k] + /* Detail */
h2i [3] * y [h + 1 ] [i] [k] + /* Approx */
g2i [3] * y [h + 1 + yOffset] [i] [k]; /* Detail */
yt [i] [j + 1] [k] =
h2i [0] * y [h ] [i] [k] + /* Approx */
g2i [0] * y [h + yOffset ] [i] [k] + /* Detail */
h2i [2] * y [h + 1 ] [i] [k] + /* Approx */
g2i [2] * y [h + 1 + yOffset] [i] [k]; /* Detail */
}
}
}
} while (i0);
#ifdef WAVELET_DBG
if (axis)
printf("VERT INV D2: "
"(%dx%d) V-approx of H-approx, (%dx%d) V-approx of H-detail,\n"
"(%dx%d) V-detail of H-approx, (%dx%d) V-detail of V-detail\n"
"-> (%dx%d) H-approx, (%dx%d) H-detail\n",
yLen, yLen, yLen, yLen, yLen, yLen, yLen, yLen,
yLen, ytLen, yLen, ytLen
);
else
printf(
"HORIZ INV D2: (%dx%d) H-approx, (%dx%d) H-detail -> (%dx%d)\n",
yLen, ytLen, yLen, ytLen,
ytLen, ytLen
);
dumpCoeffs2(y, yt, yTotalLen, yTotalLen, 0, NULL);
putchar('\n');
#endif
axis ^= 1;
}
return ytLen;
}
unsigned padWaveletXform2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned yLen0, unsigned *numNonZero
)
/* Perform full 2D multiresolution forward wavelet transform with padded
boundary coefficients on array y containing (yLen0 x yLen0) samples of
3-tuples. Note no intra-tuple transform occurs. An additional array yt
of identical dimension to y is required as buffer for intermediate
results.
The wavelet coefficients are returned in array y, containing the coarsest
approximation coefficients in y [0][0] followed by horizontal/vertical
detail coefficients in order of increasing resolution/frequency.
NOTE: Both y and yt must be large enough to accommodate the additional
padding coefficients introduced at the array boundaries during the
transform. It is the caller's responsibility to preallocate and
dimension these arrays accordingly with allocWaveletMatrix2(). The
dimension of the output array (which generally exceeds the input
dimension yLen0) must be interrogated beforehand by calling the function
with y or yt set to NULL.
The returned value is the output dimension, or 0 if the transform failed.
In addition, the number of nonzero coefficients is returned in
numNonZero (if not NULL), which is a subset of the number of output
coefficients, i.e. (output dimension)^2.
*/
{
/* Output size (number of coeffs, including boundary padding coeffs),
output array offset and pointers */
unsigned yLen, ytLen, ytOffset, totalLen = 0, totalNZ = 0;
/* Transform if enough samples, else just return output size.
NOTE: the D2 transform will be skipped for input sizes yLen0 < 4. */
if (yLen0 >= 2) {
/* Figure out total number of coefficients (= maximum size of output
array including boundary padding) by summing the number of detail
coeffs returned per dimension by each xform iteration. Also sum
number of nonzero detail coefficients. */
for (yLen = yLen0, totalLen = 0;
yLen > WAVELET_PADD2_MINLEN;
totalLen += yLen = padD2Step2(NULL, NULL, yLen, 0),
totalNZ += yLen * yLen
);
/* Finally add number of approx coeffs at coarsest resolution */
totalLen += yLen;
totalNZ = 3 * totalNZ + yLen * yLen;
if (y && yt) {
/* Apply Daubechies D2 transform, accounting for boundary padding */
for (yLen = yLen0, ytOffset = totalLen;
yLen > WAVELET_PADD2_MINLEN;
yLen = ytLen
) {
/* Get num of coeff pairs returned by the next xform iteration,
and calculate the offset of detail coeffs in output array */
ytOffset -= ytLen = padD2Step2(NULL, NULL, yLen, 0);
/* Apply horizontal & vertical Daubechies D2 transform, swapping
input and transposed output array. The resulting detail
coefficients are prepended to those from the previous
iteration from position ytOffset; consequently, the returned
number of coefficient pairs is the input size for next
iteration */
if (!padD2Step2(y, yt, yLen, ytOffset) ||
!padD2Step2(yt, y, yLen, ytOffset)
)
return 0;
}
#ifdef WAVELET_FINAL_HAAR
/* Apply horizontal & vertical Daubechies D1 (Haar) transform in
penultimate and final iterations (since yLen is odd) to
decompose to a single approximation coefficient y [0] [0] at the
coarsest resolution; all other coeffs are details. Obviously no
boundary padding is necessary here. */
for (; yLen >= 2; yLen--) {
if (!padD1Step2(y, yt, yLen) || !padD1Step2(yt, y, yLen))
return 0;
}
#endif
}
}
/* Optionally return number of nonzero coeffs */
if (numNonZero)
*numNonZero = totalNZ;
/* All coeffs now in y; return output size */
return totalLen;
}
unsigned padWaveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 yt,
unsigned yLen0, unsigned yLenOrg
)
/* Perform full 2D multiresolution inverse wavelet transform on array y
containing (yLen0 x yLen0) boundary padded wavelet coefficients as
3-tuples. The original number of samples (yLenOrg x yLenOrg) is required
for the reconstruction. An additional array yt of identical dimension to
y is required as buffer for intermediate results.
The reconstructed samples are returned in y. This is smaller than (yLen0
x yLen0) since the boundary coefficients from the forward transform are
collapsed.
The returned value is the reconstructed number of samples per dimension
(yLenOrg), or 0 if the inverse transform failed.
*/
{
unsigned yLen = 0, nextyLen, ytLen, ytOffset, i, j, d2Iter;
if (yLenOrg < 2)
/* Too few coefficients for reconstruction */
return 0;
/* Figure out num Daubechies iterations based on original num coeffs */
for (d2Iter = 0, yLen = yLenOrg;
yLen > WAVELET_PADD2_MINLEN;
d2Iter++, yLen = padD2Step2(NULL, NULL, yLen, 0)
);
#ifdef WAVELET_FINAL_HAAR
/* Invert horizontal & vertical Daubechies D1 (Haar) transform at the
two coarsest levels */
for (yLen = 2; yLen <= min(3, yLen0); yLen++)
if (!padD1InvStep2(y, yt, yLen) || !padD1InvStep2(yt, y, yLen))
return 0;
#endif
/* Invert horiz & vert Daubechies D2 transform on subsequent levels */
for (i = d2Iter; i; i--) {
/* Get num detail coeffs and their offset for current iteration */
for (j = 0, ytOffset = yLen0, yLen = nextyLen = yLenOrg; j < i; j++) {
nextyLen = yLen;
ytOffset -= yLen = padD2Step2(NULL, NULL, yLen, 0);
}
if (!(ytLen = padD2InvStep2(y, yt, yLen, ytOffset)) ||
!(ytLen = padD2InvStep2(yt, y, yLen, ytOffset))
)
return 0;
/* XXX: Drop excess reconstructed approximation coeffs to match number
of detail coeffs in next iteration. */
yLen = min(ytLen, nextyLen);
}
return yLenOrg;
}
#ifdef WAVELET_TEST_2DPAD
#define TINY 1e-6
#define HUGE 1e10
#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 > 7)
#error Invalid wavelet unit test data initialiation option
#endif
/* Define to save coeffs and inv. xformed data to file */
#define WAVELET_TEST_COEFFOUT "wavelet3-test-coeff.dat"
#define WAVELET_TEST_COEFFOUT_mRGBE "wavelet3-test-coeff-mrgbe.dat"
#define WAVELET_TEST_OUT "wavelet3-test-xform.dat"
#define WAVELET_TEST_OUT_mRGBE "wavelet3-test-xform-mrgbe.dat"
#ifdef WAVELET_TEST_LOG
#define WAVELET_TEST_VAL(c) ((c) > TINY ? log(c) : 0)
#define WAVELET_TEST_IVAL(c) (fabs(c) > TINY ? exp(c) : 0)
#else
#define WAVELET_TEST_VAL(c) (c)
#define WAVELET_TEST_IVAL(c) (c)
#endif
static void avgDetailCoeffs (const WaveletMatrix2 y, unsigned yLen,
WaveletCoeff3 avg
)
{
unsigned i, j, k, numCoeffs = 0;
avg [0] = avg [1] = avg [2] = 0;
for (i = 0; i < yLen; i++)
for (j = 0; j < yLen; j++)
/* Skip zero and cleared coeffs */
if ((i > 2 || j > 2) && !coeffIsZero(y [i][j]) &&
!coeffIsEmpty(y [i][j])
) {
for (k = 0; k < 3; k++)
avg [k] += y [i][j][k];
numCoeffs++;
}
for (k = 0; k < 3; k++)
avg [k] /= numCoeffs;
}
int main (int argc, char *argv [])
{
int i, j, k;
unsigned y0Len, yLen, numThresh = 0, numNonzero = 0;
WaveletMatrix2 y0 = NULL, y = NULL, yt = NULL;
FILE *dataFile = NULL;
WAVELET_COEFF inData, thresh = 0;
WaveletCoeff3 avgDetail;
#ifdef WAVELET_TEST_mRGBE
mRGBE mrgbeCoeff;
mRGBERange mrgbeRange;
WaveletMatrix2 ymrgbe = NULL;
#endif
if (argc < 2) {
fprintf(stderr, "%s <dim> [threshold] [dataFile]\n", argv [0]);
fputs("Missing array dimension dim >= 2\n", stderr);
return -1;
}
if (!(y0Len = atoi(argv [1])) || y0Len < 2) {
fprintf(stderr, "Invalid array dimension %d\n", y0Len);
return -1;
}
if (argc > 2 && (thresh = atof(argv [2])) < 0) {
fprintf(stderr, "Invalid threshold %.3f\n", thresh);
return -1;
}
/* Get size of (padded) output array */
if (!(yLen = padWaveletXform2(NULL, NULL, y0Len, &numNonzero))) {
fputs("Can't determine number of coefficients "
"(input too short?)\n", stderr
);
return -1;
}
/* Allocate & clear arrays for original and reconstruction */
if (!(y0 = allocWaveletMatrix2(y0Len)) ||
!(y = allocWaveletMatrix2(yLen)) ||
!(yt = allocWaveletMatrix2(yLen))
#ifdef WAVELET_TEST_mRGBE
|| !(ymrgbe = allocWaveletMatrix2(yLen))
#endif
) {
fprintf(stderr, "Failed allocating %dx%d array\n", yLen, yLen);
return -1;
}
clearCoeffs2(y0, y0Len);
clearCoeffs2(y, yLen);
clearCoeffs2(yt, yLen);
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 < y0Len; i++) {
for (j = 0; j < y0Len; j++) {
if (feof(dataFile)) {
fprintf(stderr,
"Premature end of file reading data from %s\n",
argv [3]
);
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 < y0Len; i++) {
for (j = 0; j < y0Len; 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, but correlated channels */
k ? (0.1 + 0.9 * drand48()) * y [i] [j] [k - 1]
: drand48();
#elif WAVELET_TEST_INIT == 2
/* Random data, but identical channels */
k ? y [i] [j] [k - 1] : drand48();
#elif WAVELET_TEST_INIT == 3
/* Monotonically increasing along axis 0 (avoiding
* zero in case of log encoding) */
i + 1;
#elif WAVELET_TEST_INIT == 4
/* Monotonically increasing along axis 1 (avoiding zero
* in case of log encoding) */
j + 1;
#elif WAVELET_TEST_INIT == 5
/* Monotonically increasing along both axes (avoiding
* zero in case of log encoding) */
i * j + 1;
#elif WAVELET_TEST_INIT == 6
/* Monotonically increasing by linear index (avoiding
* zero in case of log encoding) */
i * y0Len + j + 1;
#elif WAVELET_TEST_INIT == 7
/* Symmetric "bump" function (avoiding zero in case of
* log encoding) */
(1.1 - fabs(j - y0Len/2. + 0.5) / (y0Len/2. - 0.5)) *
(1.1 - fabs(i - y0Len/2. + 0.5) / (y0Len/2. - 0.5));
#endif
}
}
}
}
#ifdef WAVELET_DBG
puts("----------------------- FWD XFORM -------------------------\n");
#endif
#ifdef WAVELET_TEST_LOG
/* Apply logarithm to input before xform */
for (i = 0; i < y0Len; i++)
for (j = 0; j < y0Len; j++)
for (k = 0; k < 3; k++)
y [i] [j] [k] = WAVELET_TEST_VAL(y [i] [j] [k]);
#endif
/* Forward xform */
if (!padWaveletXform2(y, yt, y0Len, NULL)) {
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: the coarsest approximation coefficients in the upper left of
* y are omitted as they are essential for reconstruction! */
for (i = numThresh = 0; i < yLen; i++)
for (j = 0; j < yLen; 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, yLen, 0, thresh, NULL);
printf("\n%d/%d nonzero coefficients", numNonzero, yLen * yLen);
avgDetailCoeffs(y, yLen, avgDetail);
printf("\nAvg = [%.3g, %.3g, %.3g]",
avgDetail [0], avgDetail [1], avgDetail [2]
);
if (numThresh)
printf("\n%d/%d coefficients thresholded = %.1f%% compression",
numThresh, yLen * yLen, 100. * numThresh / (yLen * yLen)
);
#ifdef WAVELET_TEST_COEFFOUT
if (!(dataFile = fopen(WAVELET_TEST_COEFFOUT, "w"))) {
fprintf(stderr, "Cannot open output file %s\n",
WAVELET_TEST_COEFFOUT
);
return -1;
}
dumpCoeffs2(y, NULL, yLen, 0, 0, dataFile);
fclose(dataFile);
fprintf(stderr, "\nDumped to %s\n", WAVELET_TEST_COEFFOUT);
#endif
#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;
clearCoeffs2(ymrgbe, yLen);
/* Find min/max coeff and init mRGBE range, omitting upper left 3x3
* approximation coefficients, since these aren't encoded */
for (i = 0; i < yLen; i++)
for (j = 0; j < yLen; j++)
if ((i >= WAVELET_PADD2_APPROXDIM ||
j >= WAVELET_PADD2_APPROXDIM
) && !coeffIsEmpty(y [i][j]) && !coeffIsZero(y [i][j])
)
for (k = 0; k < 3; k++) {
/* Skip empty (NaN) coeffs */
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 upper 3x3
* approximation coefficients */
for (i = 0; i < yLen; i++)
for (j = 0; j < yLen; j++)
if ((i >= WAVELET_PADD2_APPROXDIM ||
j >= WAVELET_PADD2_APPROXDIM
) && !coeffIsEmpty(y [i][j]) && !coeffIsZero(y [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(ymrgbe, NULL, yLen, yLen, 0, NULL);
avgDetailCoeffs(y, yLen, avgDetail);
printf("\nAvg = [%.3g, %.3g, %.3g]",
avgDetail [0], avgDetail [1], avgDetail [2]
);
#ifdef WAVELET_TEST_COEFFOUT
if (!(dataFile = fopen(WAVELET_TEST_COEFFOUT_mRGBE, "w"))) {
fprintf(stderr, "Cannot open output file %s\n",
WAVELET_TEST_COEFFOUT_mRGBE
);
return -1;
}
dumpCoeffs2(ymrgbe, NULL, yLen, 0, 0, dataFile);
fclose(dataFile);
fprintf(stderr, "\nDumped to %s\n", WAVELET_TEST_COEFFOUT_mRGBE);
#endif
#endif
#endif
#ifdef WAVELET_DBG
puts("\n----------------------- INV XFORM -------------------------\n");
#endif
/* Inverse xform */
if (!padWaveletInvXform2(y, yt, yLen, y0Len)) {
fputs("\nInverse xform failed\n", stderr);
return -1;
}
#ifdef WAVELET_TEST_mRGBE
#ifdef WAVELET_DBG
puts("\n------------------- INV XFORM (mRGBE) ---------------------\n");
#endif
/* Inverse xform using mRGBE coeffs if enabled */
if (!padWaveletInvXform2(ymrgbe, yt, yLen, y0Len)) {
fputs("\nInverse xform failed\n", stderr);
return -1;
}
#endif
#ifdef WAVELET_TEST_LOG
/* Invert logarithm after xform */
for (i = 0; i < y0Len; i++)
for (j = 0; j < y0Len; j++)
for (k = 0; k < 3; k++) {
y [i] [j] [k] = WAVELET_TEST_IVAL(y [i] [j] [k]);
#ifdef WAVELET_TEST_mRGBE
ymrgbe [i] [j] [k] = WAVELET_TEST_IVAL(ymrgbe [i] [j] [k]);
#endif
}
#endif
#ifdef WAVELET_DBG
puts("\n--------------------- ORIG vs. INV XFORM ------------------------\n");
dumpCoeffs2(y0, y, y0Len, y0Len, 0, NULL);
#ifdef WAVELET_TEST_OUT
if (!(dataFile = fopen(WAVELET_TEST_OUT, "w"))) {
fprintf(stderr, "Cannot open output file %s\n",
WAVELET_TEST_OUT
);
return -1;
}
dumpCoeffs2(y, NULL, y0Len, 0, 0, dataFile);
fclose(dataFile);
fprintf(stderr, "\nDumped to %s\n", WAVELET_TEST_OUT);
#endif
#endif
printf("\nAvg RMSE = %.2f\n", rmseCoeffs2(y0, y, y0Len));
#ifdef WAVELET_TEST_mRGBE
#ifdef WAVELET_DBG
puts("\n------------------ ORIG vs. INV XFORM + mRGBE -------------------\n");
dumpCoeffs2(y0, ymrgbe, y0Len, y0Len, 0, NULL);
#ifdef WAVELET_TEST_OUT_mRGBE
if (!(dataFile = fopen(WAVELET_TEST_OUT_mRGBE, "w"))) {
fprintf(stderr, "Cannot open output file %s\n",
WAVELET_TEST_OUT_mRGBE
);
return -1;
}
dumpCoeffs2(ymrgbe, NULL, y0Len, 0, 0, dataFile);
fclose(dataFile);
fprintf(stderr, "\nDumped to %s\n", WAVELET_TEST_OUT_mRGBE);
#endif
#endif
printf("\nAvg RMSE with mRGBE enc = %.2f\n",
rmseCoeffs2(y0, ymrgbe, y0Len)
);
#endif
freeWaveletMatrix2(y0, y0Len);
freeWaveletMatrix2(y, yLen);
freeWaveletMatrix2(yt, yLen);
#ifdef WAVELET_TEST_mRGBE
freeWaveletMatrix2(ymrgbe, yLen);
#endif
return 0;
}
#endif /* WAVELET_TEST_2DPAD */
c4science · Help