opencv/3rdparty/openexr/IlmImf/ImfWav.cpp

///////////////////////////////////////////////////////////////////////////
//
// 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