391 lines
8.4 KiB
C++
391 lines
8.4 KiB
C++
///////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
|
|
// Digital Ltd. LLC
|
|
//
|
|
// All rights reserved.
|
|
//
|
|
// Redistribution and use in source and binary forms, with or without
|
|
// modification, are permitted provided that the following conditions are
|
|
// met:
|
|
// * Redistributions of source code must retain the above copyright
|
|
// notice, this list of conditions and the following disclaimer.
|
|
// * Redistributions in binary form must reproduce the above
|
|
// copyright notice, this list of conditions and the following disclaimer
|
|
// in the documentation and/or other materials provided with the
|
|
// distribution.
|
|
// * Neither the name of Industrial Light & Magic nor the names of
|
|
// its contributors may be used to endorse or promote products derived
|
|
// from this software without specific prior written permission.
|
|
//
|
|
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
|
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
|
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
|
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
|
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
|
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
|
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
//
|
|
///////////////////////////////////////////////////////////////////////////
|
|
|
|
|
|
|
|
//-----------------------------------------------------------------------------
|
|
//
|
|
// 16-bit Haar Wavelet encoding and decoding
|
|
//
|
|
// The source code in this file is derived from the encoding
|
|
// and decoding routines written by Christian Rouet for his
|
|
// PIZ image file format.
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
|
|
#include <ImfWav.h>
|
|
|
|
namespace Imf {
|
|
namespace {
|
|
|
|
|
|
//
|
|
// Wavelet basis functions without modulo arithmetic; they produce
|
|
// the best compression ratios when the wavelet-transformed data are
|
|
// Huffman-encoded, but the wavelet transform works only for 14-bit
|
|
// data (untransformed data values must be less than (1 << 14)).
|
|
//
|
|
|
|
inline void
|
|
wenc14 (unsigned short a, unsigned short b,
|
|
unsigned short &l, unsigned short &h)
|
|
{
|
|
short as = a;
|
|
short bs = b;
|
|
|
|
short ms = (as + bs) >> 1;
|
|
short ds = as - bs;
|
|
|
|
l = ms;
|
|
h = ds;
|
|
}
|
|
|
|
|
|
inline void
|
|
wdec14 (unsigned short l, unsigned short h,
|
|
unsigned short &a, unsigned short &b)
|
|
{
|
|
short ls = l;
|
|
short hs = h;
|
|
|
|
int hi = hs;
|
|
int ai = ls + (hi & 1) + (hi >> 1);
|
|
|
|
short as = ai;
|
|
short bs = ai - hi;
|
|
|
|
a = as;
|
|
b = bs;
|
|
}
|
|
|
|
|
|
//
|
|
// Wavelet basis functions with modulo arithmetic; they work with full
|
|
// 16-bit data, but Huffman-encoding the wavelet-transformed data doesn't
|
|
// compress the data quite as well.
|
|
//
|
|
|
|
const int NBITS = 16;
|
|
const int A_OFFSET = 1 << (NBITS - 1);
|
|
const int M_OFFSET = 1 << (NBITS - 1);
|
|
const int MOD_MASK = (1 << NBITS) - 1;
|
|
|
|
|
|
inline void
|
|
wenc16 (unsigned short a, unsigned short b,
|
|
unsigned short &l, unsigned short &h)
|
|
{
|
|
int ao = (a + A_OFFSET) & MOD_MASK;
|
|
int m = ((ao + b) >> 1);
|
|
int d = ao - b;
|
|
|
|
if (d < 0)
|
|
m = (m + M_OFFSET) & MOD_MASK;
|
|
|
|
d &= MOD_MASK;
|
|
|
|
l = m;
|
|
h = d;
|
|
}
|
|
|
|
|
|
inline void
|
|
wdec16 (unsigned short l, unsigned short h,
|
|
unsigned short &a, unsigned short &b)
|
|
{
|
|
int m = l;
|
|
int d = h;
|
|
int bb = (m - (d >> 1)) & MOD_MASK;
|
|
int aa = (d + bb - A_OFFSET) & MOD_MASK;
|
|
b = bb;
|
|
a = aa;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
|
|
//
|
|
// 2D Wavelet encoding:
|
|
//
|
|
|
|
void
|
|
wav2Encode
|
|
(unsigned short* in, // io: values are transformed in place
|
|
int nx, // i : x size
|
|
int ox, // i : x offset
|
|
int ny, // i : y size
|
|
int oy, // i : y offset
|
|
unsigned short mx) // i : maximum in[x][y] value
|
|
{
|
|
bool w14 = (mx < (1 << 14));
|
|
int n = (nx > ny)? ny: nx;
|
|
int p = 1; // == 1 << level
|
|
int p2 = 2; // == 1 << (level+1)
|
|
|
|
//
|
|
// Hierachical loop on smaller dimension n
|
|
//
|
|
|
|
while (p2 <= n)
|
|
{
|
|
unsigned short *py = in;
|
|
unsigned short *ey = in + oy * (ny - p2);
|
|
int oy1 = oy * p;
|
|
int oy2 = oy * p2;
|
|
int ox1 = ox * p;
|
|
int ox2 = ox * p2;
|
|
unsigned short i00,i01,i10,i11;
|
|
|
|
//
|
|
// Y loop
|
|
//
|
|
|
|
for (; py <= ey; py += oy2)
|
|
{
|
|
unsigned short *px = py;
|
|
unsigned short *ex = py + ox * (nx - p2);
|
|
|
|
//
|
|
// X loop
|
|
//
|
|
|
|
for (; px <= ex; px += ox2)
|
|
{
|
|
unsigned short *p01 = px + ox1;
|
|
unsigned short *p10 = px + oy1;
|
|
unsigned short *p11 = p10 + ox1;
|
|
|
|
//
|
|
// 2D wavelet encoding
|
|
//
|
|
|
|
if (w14)
|
|
{
|
|
wenc14 (*px, *p01, i00, i01);
|
|
wenc14 (*p10, *p11, i10, i11);
|
|
wenc14 (i00, i10, *px, *p10);
|
|
wenc14 (i01, i11, *p01, *p11);
|
|
}
|
|
else
|
|
{
|
|
wenc16 (*px, *p01, i00, i01);
|
|
wenc16 (*p10, *p11, i10, i11);
|
|
wenc16 (i00, i10, *px, *p10);
|
|
wenc16 (i01, i11, *p01, *p11);
|
|
}
|
|
}
|
|
|
|
//
|
|
// Encode (1D) odd column (still in Y loop)
|
|
//
|
|
|
|
if (nx & p)
|
|
{
|
|
unsigned short *p10 = px + oy1;
|
|
|
|
if (w14)
|
|
wenc14 (*px, *p10, i00, *p10);
|
|
else
|
|
wenc16 (*px, *p10, i00, *p10);
|
|
|
|
*px= i00;
|
|
}
|
|
}
|
|
|
|
//
|
|
// Encode (1D) odd line (must loop in X)
|
|
//
|
|
|
|
if (ny & p)
|
|
{
|
|
unsigned short *px = py;
|
|
unsigned short *ex = py + ox * (nx - p2);
|
|
|
|
for (; px <= ex; px += ox2)
|
|
{
|
|
unsigned short *p01 = px + ox1;
|
|
|
|
if (w14)
|
|
wenc14 (*px, *p01, i00, *p01);
|
|
else
|
|
wenc16 (*px, *p01, i00, *p01);
|
|
|
|
*px= i00;
|
|
}
|
|
}
|
|
|
|
//
|
|
// Next level
|
|
//
|
|
|
|
p = p2;
|
|
p2 <<= 1;
|
|
}
|
|
}
|
|
|
|
|
|
//
|
|
// 2D Wavelet decoding:
|
|
//
|
|
|
|
void
|
|
wav2Decode
|
|
(unsigned short* in, // io: values are transformed in place
|
|
int nx, // i : x size
|
|
int ox, // i : x offset
|
|
int ny, // i : y size
|
|
int oy, // i : y offset
|
|
unsigned short mx) // i : maximum in[x][y] value
|
|
{
|
|
bool w14 = (mx < (1 << 14));
|
|
int n = (nx > ny)? ny: nx;
|
|
int p = 1;
|
|
int p2;
|
|
|
|
//
|
|
// Search max level
|
|
//
|
|
|
|
while (p <= n)
|
|
p <<= 1;
|
|
|
|
p >>= 1;
|
|
p2 = p;
|
|
p >>= 1;
|
|
|
|
//
|
|
// Hierarchical loop on smaller dimension n
|
|
//
|
|
|
|
while (p >= 1)
|
|
{
|
|
unsigned short *py = in;
|
|
unsigned short *ey = in + oy * (ny - p2);
|
|
int oy1 = oy * p;
|
|
int oy2 = oy * p2;
|
|
int ox1 = ox * p;
|
|
int ox2 = ox * p2;
|
|
unsigned short i00,i01,i10,i11;
|
|
|
|
//
|
|
// Y loop
|
|
//
|
|
|
|
for (; py <= ey; py += oy2)
|
|
{
|
|
unsigned short *px = py;
|
|
unsigned short *ex = py + ox * (nx - p2);
|
|
|
|
//
|
|
// X loop
|
|
//
|
|
|
|
for (; px <= ex; px += ox2)
|
|
{
|
|
unsigned short *p01 = px + ox1;
|
|
unsigned short *p10 = px + oy1;
|
|
unsigned short *p11 = p10 + ox1;
|
|
|
|
//
|
|
// 2D wavelet decoding
|
|
//
|
|
|
|
if (w14)
|
|
{
|
|
wdec14 (*px, *p10, i00, i10);
|
|
wdec14 (*p01, *p11, i01, i11);
|
|
wdec14 (i00, i01, *px, *p01);
|
|
wdec14 (i10, i11, *p10, *p11);
|
|
}
|
|
else
|
|
{
|
|
wdec16 (*px, *p10, i00, i10);
|
|
wdec16 (*p01, *p11, i01, i11);
|
|
wdec16 (i00, i01, *px, *p01);
|
|
wdec16 (i10, i11, *p10, *p11);
|
|
}
|
|
}
|
|
|
|
//
|
|
// Decode (1D) odd column (still in Y loop)
|
|
//
|
|
|
|
if (nx & p)
|
|
{
|
|
unsigned short *p10 = px + oy1;
|
|
|
|
if (w14)
|
|
wdec14 (*px, *p10, i00, *p10);
|
|
else
|
|
wdec16 (*px, *p10, i00, *p10);
|
|
|
|
*px= i00;
|
|
}
|
|
}
|
|
|
|
//
|
|
// Decode (1D) odd line (must loop in X)
|
|
//
|
|
|
|
if (ny & p)
|
|
{
|
|
unsigned short *px = py;
|
|
unsigned short *ex = py + ox * (nx - p2);
|
|
|
|
for (; px <= ex; px += ox2)
|
|
{
|
|
unsigned short *p01 = px + ox1;
|
|
|
|
if (w14)
|
|
wdec14 (*px, *p01, i00, *p01);
|
|
else
|
|
wdec16 (*px, *p01, i00, *p01);
|
|
|
|
*px= i00;
|
|
}
|
|
}
|
|
|
|
//
|
|
// Next level
|
|
//
|
|
|
|
p2 = p;
|
|
p >>= 1;
|
|
}
|
|
}
|
|
|
|
|
|
} // namespace Imf
|