///////////////////////////////////////////////////////////////////////////

//

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

#include "ImfNamespace.h"

OPENEXR_IMF_INTERNAL_NAMESPACE_SOURCE_ENTER

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;

    }

}

OPENEXR_IMF_INTERNAL_NAMESPACE_SOURCE_EXIT