Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86836743
wavelet2.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Oct 8, 21:41
Size
27 KB
Mime Type
text/x-c
Expires
Thu, Oct 10, 21:41 (2 d)
Engine
blob
Format
Raw Data
Handle
21490163
Attached To
R10977 RADIANCE Photon Map
wavelet2.c
View Options
/*
=========================================================================
2-dimensional Daubechies D2 and D1 (Haar) wavelet transforms on
(2^l) x (2^l) sized arrays of 3-tuples, where l > 1. Array boundaries
are extended as periodic, which does not require padding coefficients.
Compile with -DWAVELET_TEST_2D to build standalone unit tests.
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>
/* 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 */
/* 2-tap Haar wavelet/scaling func coeffs */
const
WAVELET_COEFF
h1
=
1
/
SQRT2
;
/* 4-tap wavelet/scaling func coeffs; default is Daubechies D2 */
#ifdef WAVELET_BIOR
/* Biorthogonal 3.1 */
const
WAVELET_COEFF
h2
[
4
]
=
{
-
0.3535533905932738
,
1.0606601717798214
,
1.0606601717798214
,
-
0.3535533905932738
},
g2
[
4
]
=
{
-
0.1767766952966369
,
0.5303300858899107
,
-
0.5303300858899107
,
0.1767766952966369
},
/* Inverse coeffs */
h2i
[
4
]
=
{
0.1767766952966369
,
0.5303300858899107
,
0.5303300858899107
,
0.1767766952966369
},
/* = { -g2 [0], g2 [1], -g2 [2], g2 [3] } */
g2i
[
4
]
=
{
-
0.3535533905932738
,
-
1.0606601717798214
,
1.0606601717798214
,
0.3535533905932738
};
/* = { h2 [0], -h2 [1], h2 [2], -h2 [3] } */
#else
/* Daubechies D2 */
const
WAVELET_COEFF
h2
[
4
]
=
{
H2NORM
*
(
1
+
SQRT3
),
H2NORM
*
(
3
+
SQRT3
),
H2NORM
*
(
3
-
SQRT3
),
H2NORM
*
(
1
-
SQRT3
)
},
g2
[
4
]
=
{
H2NORM
*
(
1
-
SQRT3
),
-
H2NORM
*
(
3
-
SQRT3
),
H2NORM
*
(
3
+
SQRT3
),
-
H2NORM
*
(
1
+
SQRT3
)
},
/* = { h2 [3], -h2 [2], h2 [1], -h2 [0] } */
/* Inverse coeffs */
h2i
[
4
]
=
{
H2NORM
*
(
1
-
SQRT3
),
H2NORM
*
(
3
-
SQRT3
),
H2NORM
*
(
3
+
SQRT3
),
H2NORM
*
(
1
+
SQRT3
)
},
/* = { h2 [3], h2 [2], h2 [1], h2 [0] } */
g2i
[
4
]
=
{
-
H2NORM
*
(
1
+
SQRT3
),
H2NORM
*
(
3
+
SQRT3
),
-
H2NORM
*
(
3
-
SQRT3
),
H2NORM
*
(
1
-
SQRT3
)
};
/* = { g2 [3], g2 [2], g2 [1], g2 [0] } */
#endif
WaveletMatrix2
allocWaveletMatrix2
(
unsigned
len
)
/*
Allocate and return a 2D coefficient array of size (len x len). Returns
NULL if allocation failed.
*/
{
unsigned
i
;
WaveletMatrix2
y
=
NULL
;
if
(
len
>=
2
)
{
if
(
!
(
y
=
calloc
(
len
,
sizeof
(
WaveletCoeff3
*
))))
return
NULL
;
for
(
i
=
0
;
i
<
len
;
i
++
)
if
(
!
(
y
[
i
]
=
calloc
(
len
,
WAVELET_COEFFSIZE
)))
return
NULL
;
}
return
y
;
}
void
freeWaveletMatrix2
(
WaveletMatrix2
y
,
unsigned
len
)
/*
Free previously allocated 2D coefficient array y of size (len x len)
*/
{
unsigned
i
,
j
;
if
(
y
)
{
for
(
i
=
0
;
i
<
len
;
i
++
)
free
(
y
[
i
]);
free
(
y
);
}
}
void
zeroCoeffs2
(
WaveletMatrix2
y
,
unsigned
len
)
/* Set 2D array coefficients to zero */
{
unsigned
i
;
if
(
y
)
for
(
i
=
0
;
i
<
len
;
i
++
)
memset
(
y
[
i
],
0
,
len
*
WAVELET_COEFFSIZE
);
}
#if 1
void
clearCoeffs2
(
WaveletMatrix2
y
,
unsigned
len
)
/* Clear 2D array to facilitate debugging */
{
unsigned
i
,
j
,
k
;
if
(
y
)
for
(
i
=
0
;
i
<
len
;
i
++
)
for
(
j
=
0
;
j
<
len
;
j
++
)
for
(
k
=
0
;
k
<
3
;
k
++
)
y
[
i
]
[
j
]
[
k
]
=
NAN
;
}
static
char
*
coeffStr
(
const
WaveletCoeff3
coeff
)
/* Format coefficient in string for dump. Cleared coeffs are represented
* as dots to facilitate debugging. */
{
static
char
str
[
9
];
if
(
coeffIsEmpty
(
coeff
))
/* Coeff is blank */
return
" .. "
;
snprintf
(
str
,
sizeof
(
str
),
"% 5.2f"
,
coeffAvg
(
coeff
));
return
str
;
}
void
dumpCoeffs2
(
const
WaveletMatrix2
y1
,
const
WaveletMatrix2
y2
,
unsigned
y1Len
,
unsigned
y2Len
,
float
thresh
,
FILE
*
stream
)
/*
Multipurpose routine to dump 2D arrays y1 and y2 of dims y1Len resp
y2Len. These are presented side-by-side, separated by an arrow to
indicate the input (y1) and output (y2) of a transform step.
If either y1 or y2 is NULL, it will be omitted from the output,
enabling dumping a single array.
If thresh > 0, coefficients with absolute magnitude below this value
are marked with brackets ('[]') as thresholded.
If stream is not NULL, the coefficients will be dumped to the
specified file, else to stdout.
*/
{
unsigned
i
,
j
,
k
;
FILE
*
out
=
stream
?
stream
:
stdout
;
for
(
i
=
0
;
i
<
max
(
y1Len
,
y2Len
);
i
++
)
{
if
(
y1
)
{
for
(
j
=
0
;
j
<
y1Len
;
j
++
)
if
(
i
<
y1Len
)
fprintf
(
out
,
"%s
\t
"
,
coeffThresh
(
y1
,
i
,
j
,
thresh
)
?
"[ .. ]"
:
coeffStr
(
y1
[
i
][
j
])
);
else
fputc
(
'\t'
,
out
);
}
if
(
y2
)
{
fprintf
(
out
,
" -->
\t
"
);
for
(
j
=
0
;
j
<
y2Len
;
j
++
)
if
(
i
<
y2Len
)
fprintf
(
out
,
"%s
\t
"
,
coeffThresh
(
y2
,
i
,
j
,
thresh
)
?
"[ .. ]"
:
coeffStr
(
y2
[
i
][
j
])
);
else
fputc
(
'\t'
,
out
);
}
fputc
(
'\n'
,
out
);
}
}
float
rmseCoeffs2
(
const
WaveletMatrix2
y1
,
const
WaveletMatrix2
y2
,
unsigned
len
)
/* Calculate RMSE between 2D matrices y1 and y2 */
{
unsigned
i
,
j
;
float
d
,
rmse
=
0
;
for
(
i
=
0
;
i
<
len
;
i
++
)
for
(
j
=
0
;
j
<
len
;
j
++
)
{
#if 0
d = (coeffAvg(y1 [i][j]) - coeffAvg(y2 [i][j])) /
coeffAvg(y1 [i][j]);
#else
d
=
coeffAvg
(
y1
[
i
][
j
])
-
coeffAvg
(
y2
[
i
][
j
]);
#endif
rmse
+=
d
*
d
;
}
return
sqrt
(
rmse
/
((
float
)
len
*
len
));
}
#endif
static
int
d1Step2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Single step of forward 2D Daubechies D1 (Haar) wavelet transform on array
y of size (2^l) x (2^l) containing 3-tuples, where l >= 1. 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 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
,
hlen
=
1
<<
(
l
-
1
);
unsigned
h
,
i
,
j
,
k
;
#ifdef WAVELET_DBG
static
unsigned
axis
=
0
;
#endif
if
(
len
<
2
||
!
y
||
!
yt
)
/* Input shorter than wavelet support, or no input/output */
return
-
1
;
#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
;
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
[
hlen
+
h
]
[
i
]
[
k
]
=
h1
*
(
y
[
i
]
[
j
]
[
k
]
-
y
[
i
]
[
j
+
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'
);
axis
^=
1
;
#endif
return
0
;
}
static
int
d1InvStep2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Single step of inverse 2D Daubechies D1 (Haar) wavelet transform on
coefficient array y of size (2^l) x (2^l) containing 3-tuples, where
l>=1. 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 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
,
hlen
=
1
<<
(
l
-
1
);
unsigned
h
,
i
,
j
,
k
;
#ifdef WAVELET_DBG
static
unsigned
axis
=
1
;
#endif
if
(
len
<
2
||
!
y
||
!
yt
)
/* Too few coeffs for reconstruction, or no input/output */
return
-
1
;
#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
;
j
+=
2
)
{
h
=
j
>>
1
;
for
(
k
=
0
;
k
<
3
;
k
++
)
{
yt
[
i
]
[
j
]
[
k
]
=
h1
*
(
y
[
h
]
[
i
]
[
k
]
+
/* Avg */
y
[
hlen
+
h
]
[
i
]
[
k
]
/* Diff */
);
yt
[
i
]
[
j
+
1
]
[
k
]
=
h1
*
(
y
[
h
]
[
i
]
[
k
]
-
/* Avg */
y
[
hlen
+
h
]
[
i
]
[
k
]
/* Diff */
);
}
}
}
#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'
);
axis
^=
1
;
#endif
return
0
;
}
static
int
d4Step2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Single step of forward 2D Daubechies D2 wavelet transform on array y of
size (2^l) x (2^l) containing 3-tuples, where l >= 2. 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 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
,
hlen
=
1
<<
(
l
-
1
);
unsigned
h
,
i
,
j
,
k
;
#ifdef WAVELET_DBG
static
unsigned
axis
=
0
;
#endif
if
(
len
<
4
||
!
y
||
!
yt
)
/* Input shorter than wavelet support, or no input/output */
return
-
1
;
#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
++
)
{
/* Transform until upper boundary */
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
]
=
h2
[
0
]
*
y
[
i
]
[
j
]
[
k
]
+
h2
[
1
]
*
y
[
i
]
[
j
+
1
]
[
k
]
+
h2
[
2
]
*
y
[
i
]
[
j
+
2
]
[
k
]
+
h2
[
3
]
*
y
[
i
]
[
j
+
3
]
[
k
];
/* Detail/diff/highpass */
yt
[
hlen
+
h
]
[
i
]
[
k
]
=
g2
[
0
]
*
y
[
i
]
[
j
]
[
k
]
+
g2
[
1
]
*
y
[
i
]
[
j
+
1
]
[
k
]
+
g2
[
2
]
*
y
[
i
]
[
j
+
2
]
[
k
]
+
g2
[
3
]
*
y
[
i
]
[
j
+
3
]
[
k
];
}
}
/* Transform at upper boundary with wraparound.
Note j is set to last index from previous loop */
h
=
j
>>
1
;
for
(
k
=
0
;
k
<
3
;
k
++
)
{
/* Smooth/approx/avg/lowpass */
yt
[
h
]
[
i
]
[
k
]
=
h2
[
0
]
*
y
[
i
]
[
j
]
[
k
]
+
h2
[
1
]
*
y
[
i
]
[
j
+
1
]
[
k
]
+
h2
[
2
]
*
y
[
i
]
[
0
]
[
k
]
+
h2
[
3
]
*
y
[
i
]
[
1
]
[
k
];
/* Detail/diff/highpass */
yt
[
hlen
+
h
]
[
i
]
[
k
]
=
g2
[
0
]
*
y
[
i
]
[
j
]
[
k
]
+
g2
[
1
]
*
y
[
i
]
[
j
+
1
]
[
k
]
+
g2
[
2
]
*
y
[
i
]
[
0
]
[
k
]
+
g2
[
3
]
*
y
[
i
]
[
1
]
[
k
];
}
}
#ifdef WAVELET_DBG
printf
(
"%s FWD D2 (%d x %d)
\n
"
,
axis
?
"VERT"
:
"HORIZ"
,
len
,
len
);
dumpCoeffs2
(
y
,
yt
,
len
,
len
,
0
,
NULL
);
putchar
(
'\n'
);
axis
^=
1
;
#endif
return
0
;
}
static
int
d4InvStep2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Single step of inverse 2D Daubechies D2 wavelet transform on coefficient
array y of size (2^l) x (2^l) containing 3-tuples, where l >= 2. 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 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
,
hlen
=
1
<<
(
l
-
1
);
unsigned
h
,
i
,
j
,
k
;
#ifdef WAVELET_DBG
static
unsigned
axis
=
1
;
#endif
if
(
len
<
4
||
!
y
||
!
yt
)
/* Too few coeffs for reconstruction, or no input/output */
return
-
1
;
#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
++
)
{
/* Invert at lower boundary with wraparound */
for
(
k
=
0
;
k
<
3
;
k
++
)
{
yt
[
i
]
[
0
]
[
k
]
=
h2i
[
1
]
*
y
[
hlen
-
1
]
[
i
]
[
k
]
+
/* Last avg */
g2i
[
1
]
*
y
[
len
-
1
]
[
i
]
[
k
]
+
/* Last diff */
h2i
[
3
]
*
y
[
0
]
[
i
]
[
k
]
+
/* First avg */
g2i
[
3
]
*
y
[
hlen
]
[
i
]
[
k
];
/* First diff */
yt
[
i
]
[
1
]
[
k
]
=
h2i
[
0
]
*
y
[
hlen
-
1
]
[
i
]
[
k
]
+
g2i
[
0
]
*
y
[
len
-
1
]
[
i
]
[
k
]
+
h2i
[
2
]
*
y
[
0
]
[
i
]
[
k
]
+
g2i
[
2
]
*
y
[
hlen
]
[
i
]
[
k
];
}
/* Invert until upper boundary */
for
(
j
=
2
;
j
<
len
;
j
+=
2
)
{
h
=
(
j
>>
1
)
-
1
;
for
(
k
=
0
;
k
<
3
;
k
++
)
{
yt
[
i
]
[
j
]
[
k
]
=
h2i
[
1
]
*
y
[
h
]
[
i
]
[
k
]
+
/* Avg */
g2i
[
1
]
*
y
[
hlen
+
h
]
[
i
]
[
k
]
+
/* Diff */
h2i
[
3
]
*
y
[
h
+
1
]
[
i
]
[
k
]
+
/* Next avg */
g2i
[
3
]
*
y
[
hlen
+
h
+
1
]
[
i
]
[
k
];
/* Next diff */
yt
[
i
]
[
j
+
1
]
[
k
]
=
h2i
[
0
]
*
y
[
h
]
[
i
]
[
k
]
+
g2i
[
0
]
*
y
[
hlen
+
h
]
[
i
]
[
k
]
+
h2i
[
2
]
*
y
[
h
+
1
]
[
i
]
[
k
]
+
g2i
[
2
]
*
y
[
hlen
+
h
+
1
]
[
i
]
[
k
];
}
}
}
#ifdef WAVELET_DBG
printf
(
"%s INV D2 (%d x %d)
\n
"
,
axis
?
"VERT"
:
"HORIZ"
,
len
,
len
);
dumpCoeffs2
(
y
,
yt
,
len
,
len
,
0
,
NULL
);
putchar
(
'\n'
);
axis
^=
1
;
#endif
return
0
;
}
int
waveletXform2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Perform full 2D multiresolution forward wavelet transform on array y of
size (2^l) x (2^l) containing original signal as 3-tuples, where l >= 1.
Note no intra-tuple transform occurs.
The wavelet coefficients are returned in array y, containing the coarsest
approximation in y [0][0] followed by horizontal/vertical details in
order of increasing resolution/frequency.
A preallocated array yt of identical dimensions to y can be supplied as
buffer for intermediate results. If yt == NULL, a buffer is
automatically allocated and freed on demand, but this is inefficient for
frequent calls. It is recommended to preallocate yt to the maximum
expected size. The dimensions of yt are not checked; this is the
caller's responsibility.
Returns 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
;
unsigned
li
;
WaveletMatrix2
ytloc
=
NULL
;
/* Skip transform if input too short or missing */
if
(
l
<
1
||
!
y
)
return
-
1
;
if
(
!
yt
)
/* No buffer supplied; allocate one on demand */
if
(
!
(
yt
=
ytloc
=
allocWaveletMatrix2
(
len
)))
{
fprintf
(
stderr
,
"ERROR - Failed allocating %dx%d buffer array"
" in WaveletXform2()"
,
len
,
len
);
return
-
1
;
}
for
(
li
=
l
;
li
>
1
;
li
--
)
{
/* Apply horizontal & vertical Daubechies D2 transform, swapping input
and transposed output array */
if
(
d4Step2
(
y
,
yt
,
li
)
||
d4Step2
(
yt
,
y
,
li
))
return
-
1
;
}
#ifdef WAVELET_FINAL_HAAR
/* Apply horizontal & vertical Daubechies D1 (Haar) transform at coarsest
resolution (li==1) to obtain single approximation coefficient at
y [0][0]; all other coeffs are details. */
if
(
d1Step2
(
y
,
yt
,
li
)
||
d1Step2
(
yt
,
y
,
li
))
return
-
1
;
#endif
/* NOTE: All coefficients now in y */
if
(
ytloc
)
/* Free yt if allocated on demand */
freeWaveletMatrix2
(
ytloc
,
len
);
return
0
;
}
int
waveletInvXform2
(
WaveletMatrix2
y
,
WaveletMatrix2
yt
,
unsigned
l
)
/*
Perform full 2D multiresolution inverse wavelet transform on array y of
size (2^l) x (2^l) containing wavelet coefficients as 3-tuples, where
l >= 1. Note no intra-tuple transform occurs.
A preallocated array yt of identical dimensions to y can be supplied as
buffer for intermediate results. If yt == NULL, a buffer is
automatically allocated and freed on demand, but this is inefficient for
frequent calls. It is recommended to preallocate yt to the maximum
expected size. The dimensions of yt are not checked; this is the
caller's responsibility.
The reconstructed signal is returned in array y.
Returns 0 on success, else -1.
*/
{
const
unsigned
len
=
1
<<
l
;
unsigned
li
;
WaveletMatrix2
ytloc
=
NULL
;
/* Skip inverse transform if input too short or missing */
if
(
l
<
1
||
!
y
)
return
-
1
;
if
(
!
yt
)
/* No buffer supplied; allocate one on demand */
if
(
!
(
yt
=
ytloc
=
allocWaveletMatrix2
(
len
)))
{
fprintf
(
stderr
,
"ERROR - Failed allocating %dx%d buffer array"
" in WaveletInvXform2()"
,
len
,
len
);
return
-
1
;
}
#ifdef WAVELET_FINAL_HAAR
/* Invert horizontal & vertical Haar transform at coarsest level (li==1),
swapping input and transposed output array */
if
(
d1InvStep2
(
y
,
yt
,
1
)
||
d1InvStep2
(
yt
,
y
,
1
))
return
-
1
;
#endif
for
(
li
=
2
;
li
<=
l
;
li
++
)
{
/* Invert horizontal & vertical Daubechies D2 transform, swapping
input and transposed output arrays */
if
(
d4InvStep2
(
y
,
yt
,
li
)
||
d4InvStep2
(
yt
,
y
,
li
))
return
-
1
;
}
/* NOTE: Reconstructed signal now in y */
if
(
ytloc
)
/* Free yt if allocated on demand */
freeWaveletMatrix2
(
ytloc
,
len
);
return
0
;
}
#ifdef WAVELET_TEST_2D
#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 > 5)
#error Invalid wavelet unit test data initialiation option
#endif
int
main
(
int
argc
,
char
*
argv
[])
{
int
i
,
j
,
k
,
l
;
unsigned
len
,
numThresh
=
0
;
WaveletMatrix2
y0
=
NULL
,
y
=
NULL
;
FILE
*
dataFile
=
NULL
;
WAVELET_COEFF
inData
,
thresh
=
0
;
#ifdef WAVELET_TEST_mRGBE
#define TINY 1e-6
#define HUGE 1e10
mRGBE
mrgbeCoeff
;
mRGBERange
mrgbeRange
;
WaveletMatrix2
ymrgbe
=
NULL
;
#endif
if
(
argc
<
2
)
{
fprintf
(
stderr
,
"%s <l> [threshold] [dataFile]
\n
"
,
argv
[
0
]);
fputs
(
"Missing array resolution l > 1, "
"compression threshold >= 0
\n
"
,
stderr
);
return
-
1
;
}
if
(
!
(
l
=
atoi
(
argv
[
1
]))
||
l
<
1
)
{
fputs
(
"Invalid array resolution l
\n
"
,
stderr
);
return
-
1
;
}
else
len
=
1
<<
l
;
if
(
argc
>
2
&&
(
thresh
=
atof
(
argv
[
2
]))
<
0
)
{
fprintf
(
stderr
,
"Invalid threshold %.3f
\n
"
,
thresh
);
return
-
1
;
}
/* Allocate arrays for original and reconstruction */
if
(
!
(
y0
=
allocWaveletMatrix2
(
len
))
||
!
(
y
=
allocWaveletMatrix2
(
len
))
#ifdef WAVELET_TEST_mRGBE
||
!
(
ymrgbe
=
allocWaveletMatrix2
(
len
))
#endif
)
{
fprintf
(
stderr
,
"Failed allocating %dx%d array
\n
"
,
len
,
len
);
return
-
1
;
}
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
<
len
;
i
++
)
{
for
(
j
=
0
;
j
<
len
;
j
++
)
{
if
(
feof
(
dataFile
))
{
fprintf
(
stderr
,
"Premature end of file reading data from %s
\n
"
,
argv
[
2
]
);
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
<
len
;
i
++
)
{
for
(
j
=
0
;
j
<
len
;
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, indentical for all channels */
k
?
y
[
i
][
j
][
k
-
1
]
:
drand48
();
#elif WAVELET_TEST_INIT == 2
/* Monotonically increasing along axis 0 */
i
;
#elif WAVELET_TEST_INIT == 3
/* Monotonically increasing along axis 1 */
j
;
#elif WAVELET_TEST_INIT == 4
/* Monotonically increasing along both axes */
i
*
j
;
#elif WAVELET_TEST_INIT == 5
/* Monotonically increasing by linear index */
i
*
len
+
j
;
#
endif
}
}
}
}
#ifdef WAVELET_DBG
puts
(
"----------------------- FWD XFORM -------------------------
\n
"
);
#endif
/* Forward xform */
if
(
waveletXform2
(
y
,
NULL
,
l
))
{
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: y [0][0] is omitted it contains the coarsest
* approximation coefficient, which is essential for the
* reconstruction! */
for
(
i
=
0
;
i
<
len
;
i
++
)
for
(
j
=
0
;
j
<
len
;
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
,
len
,
0
,
thresh
,
NULL
);
if
(
numThresh
)
printf
(
"
\n
%d/%d coefficients thresholded = %.1f%% compression"
,
numThresh
,
len
*
len
,
100.
*
numThresh
/
(
len
*
len
)
);
#
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
;
/* Find min/max coeff and init mRGBE range */
for
(
i
=
0
;
i
<
len
;
i
++
)
for
(
j
=
0
;
j
<
len
;
j
++
)
for
(
k
=
0
;
k
<
3
;
k
++
)
{
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 approximation
* coefficient at y [0][0] */
for
(
i
=
0
;
i
<
len
;
i
++
)
for
(
j
=
0
;
j
<
len
;
j
++
)
if
(
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
(
y
,
ymrgbe
,
len
,
len
,
0
,
NULL
);
#endif
#endif
#ifdef WAVELET_DBG
puts
(
"
\n
----------------------- INV XFORM -------------------------
\n
"
);
#endif
/* Inverse xform (also using mRGBE coeffs if enabled) */
if
(
waveletInvXform2
(
y
,
NULL
,
l
)
#ifdef WAVELET_TEST_mRGBE
||
waveletInvXform2
(
ymrgbe
,
NULL
,
l
)
#endif
)
{
fputs
(
"
\n
Inverse xform failed
\n
"
,
stderr
);
return
-
1
;
}
#ifdef WAVELET_DBG
puts
(
"
\n
--------------------- ORIG vs. INV XFORM ------------------------
\n
"
);
dumpCoeffs2
(
y0
,
y
,
len
,
len
,
0
,
NULL
);
#endif
printf
(
"
\n
Avg RMSE = %.2f
\n
"
,
rmseCoeffs2
(
y0
,
y
,
len
));
#ifdef WAVELET_TEST_mRGBE
#ifdef WAVELET_DBG
puts
(
"
\n
------------------ ORIG vs. INV XFORM + mRGBE -------------------
\n
"
);
dumpCoeffs2
(
y0
,
ymrgbe
,
len
,
len
,
0
,
NULL
);
#endif
printf
(
"
\n
Avg RMSE with mRGBE enc = %.2f
\n
"
,
rmseCoeffs2
(
y0
,
ymrgbe
,
len
)
);
#endif
freeWaveletMatrix2
(
y0
,
len
);
freeWaveletMatrix2
(
y
,
len
);
#ifdef WAVELET_TEST_mRGBE
freeWaveletMatrix2
(
ymrgbe
,
len
);
#endif
return
0
;
}
#endif
/* WAVELET_TEST_2D */
Event Timeline
Log In to Comment