/*************************************************************************

This software module was originally developed by 
	Stefan Rauthenberg (rauthenberg@HHI.DE), HHI
	(date: January, 1998)

in the course of development of the MPEG-4 Video (ISO/IEC 14496-2). 
This software module is an implementation of a part of one or more MPEG-4 Video tools 
as specified by the MPEG-4 Video. 
ISO/IEC gives users of the MPEG-4 Video free license to this software module or modifications 
thereof for use in hardware or software products claiming conformance to the MPEG-4 Video. 
Those intending to use this software module in hardware or software products are advised that its use may infringe existing patents. 
The original developer of this software module and his/her company, 
the subsequent editors and their companies, 
and ISO/IEC have no liability for use of this software module or modifications thereof in an implementation. 
Copyright is not released for non MPEG-4 Video conforming products. 
Microsoft retains full right to use the code for his/her own purpose, 
assign or donate the code to a third party and to inhibit third parties from using the code for non <MPEG standard> conforming products. 
This copyright notice must be included in all copies or derivative works. 

Copyright (c) 1998.

Module Name:

	sadct.cpp

Abstract:

	SADCT and inverse SADCT

Revision History:
  
*************************************************************************/


#include "typeapi.h"
#include "dct.hpp"
#include <math.h>
#if defined(__DEBUG_SADCT_) && !defined(NDEBUG)
#include <iostream.h>
#endif
#ifdef __MFC_
#ifdef _DEBUG
#undef THIS_FILE
static char BASED_CODE THIS_FILE[] = __FILE__;
#endif

#define new DEBUG_NEW				   
#endif // __MFC_

CSADCT::CSADCT():
	m_N(8),
	m_mat_tmp1(0),
	m_row_buf(0),
	m_ly(0),
	m_lx(0),
	m_mask(0)
{

	allocMatrix(&m_mat_tmp1, m_N, m_N);
	
	m_row_buf = new Float[m_N];
	m_ly = new Int[m_N];
	m_lx = new Int[m_N];
	
	allocMatrix(&m_mask, m_N, m_N);
	
	allocMatrix(&m_in, m_N, m_N);
	allocMatrix(&m_out, m_N, m_N);
	
#if defined(_MSC_VER)
	M_PI = 2.0 * asin(1.0);
	M_SQRT2 = sqrt(2.0);
#endif

// Schueuer HHI: added for fast_sadct
#ifdef _FAST_SADCT_
	c_buf = new Float[m_N];
	allocMatrix(&tmp_out, m_N, m_N);


	f0_2 = 0.707107;
	f0_3 = 0.577350;
	f1_3 = 0.707107;
	f2_3 = 0.408248;
	f3_3 = 0.816497;
	f0_4 = 0.500000;
	f1_4 = 0.653281;
	f2_4 = 0.270598;
	f0_5 = 0.447214;
	f1_5 = 0.601501;
	f2_5 = 0.371748;
	f3_5 = 0.511667;
	f4_5 = 0.195440;
	f5_5 = 0.632456;
	f0_6 = 0.408248;
	f1_6 = 0.557678;
	f2_6 = 0.500000; 
	f3_6 = 0.288675;
	f4_6 = 0.149429;
	f5_6 = 0.577350;
	f0_7 = 0.377964;
	f1_7 = 0.521121;
	f2_7 = 0.481588; 
	f3_7 = 0.417907;
	f4_7 = 0.333269;
	f5_7 = 0.231921;
	f6_7 = 0.118942;
	f7_7 = 0.534522;
	f0_8 = 0.7071068;
	f1_8 = 0.4903926;
	f2_8 = 0.4619398;
	f3_8 = 0.4157348;
	f4_8 = 0.3535534;
	f5_8 = 0.2777851;
	f6_8 = 0.1913417;
	f7_8 = 0.0975452;


	sq[0] = 0.00000000000;
	sq[1] = 1.00000000000;
	sq[2] =	1.41421356237;
	sq[3] =	1.73205080757;
	sq[4] =	2.00000000000;
	sq[5] =	2.23606797750;
	sq[6] =	2.44948974278;
	sq[7] =	2.64575131106;
	sq[8] =	2.82842712475;
#endif
}

CSADCT::~CSADCT()
{
	freeMatrix(m_mat_tmp1, m_N);

	delete [] m_row_buf;
	delete [] m_ly;
	delete [] m_lx;

	freeMatrix(m_mask, m_N);
	
	freeMatrix(m_in, m_N);
	freeMatrix(m_out, m_N);

// Schueuer HHI 
#ifdef _FAST_SADCT_
	delete [] c_buf;
	freeMatrix(tmp_out, m_N);
#endif

}

Void CSADCT::allocMatrix(Float ***mat, Int nr, Int nc) 
{
	Float **m = new Float* [nr];
	m[0] = new Float[nr*nc];
	
	for (Int i=1; i<nr; i++) {
		m[i] = m[i-1] + nc;
	}
	*mat = m;
}

Void CSADCT::allocMatrix(PixelC ***mat, Int nr, Int nc) 
{
	PixelC **m = new PixelC* [nr];
	m[0] = new PixelC[nr*nc];
	
	for (Int i=1; i<nr; i++) {
		m[i] = m[i-1] + nc;
	}
	*mat = m;
}

Float ***CSADCT::allocDctTable(Int n) 
{
	Float ***tbl = new Float** [n+1];
	tbl[0] = 0;		// transformation length of zero is impossible
	
	
	for (Int i=1; i<=n; i++) {
		allocMatrix(&tbl[i], n, n);
	}
	
	return (tbl);
}

Float ***CInvSADCT::allocReorderTable(Int n)
{
	Float ***tbl = new Float** [n];
	
	for (Int i=0; i<n; i++) {
		tbl[i] = new Float* [n];
		// to play it safe and get an exception if trying to access
		// an uninitialized poInter.
		memset(tbl[i], 0, sizeof(Float*)*n);
	}
	return tbl;
}

Void CInvSADCT::freeReorderTable(Float ***tbl, Int n)
{
	
	if (tbl) {
		for (Int i=0; i<n; i++)
			delete [] tbl[i];
		delete [] tbl;
	}
}



// parameter nr not used yet because the data elements are allocated in contiguous memory and 
// in separate rows.
Void CSADCT::freeMatrix(Float **mat, Int nr) 
{
	if (mat) {
		delete [] mat[0];
		delete [] mat;
	}
}

Void CSADCT::freeMatrix(PixelC **mat, Int nr) 
{
	if (mat) {
		delete [] mat[0];
		delete [] mat;
	}
}

Void CSADCT::freeDctTable(Float ***tbl, Int n) 
{
	if (tbl) {
		for (Int i=1; i<=n; i++) {
			freeMatrix(tbl[i], n);
		}
 		delete [] tbl;
	}
}

Void CSADCT::getRowLength(Int *lx, const PixelC* ppxlcMask, Int iStride)
{
	prepareMask(ppxlcMask, iStride);
	getRowLengthInternal(lx, m_mask, m_N, m_N);
}

/*
 *	returns in `lx' the number of active pels per line after shifting
 *	the pels towards the upper left corner.  This function may be
 *	called by a decoder in order to find out how the pels of a partial
 *	block are arranged for transmission.  
 */
Void CSADCT::getRowLengthInternal(Int *lx, PixelC **mask, Int bky, Int bkx)
{
	Int iy_out = 0;
	Int iy, ix, l;
	
	for (ix=0; ix<bkx; ix++) {
		l = 0;
		for (iy=0; iy<bky; iy++) {
			if ( mask[iy][ix] ) 
				l++;
		}
		if ( l ) 
			m_ly[iy_out++] = l;
	}
	
	for (ix=iy_out; ix<bkx; ix++) 
		m_ly[ix] = 0;
	
	for (iy=0; iy<bky; iy++) {
		l = 0;
		for (ix=0; ix<bkx; ix++) {
			if ( m_ly[ix] > iy )
				l++;
		}
		lx[iy] = l;
	}	
}

Void CSADCT::prepareMask(const PixelC* rgchMask, Int stride)
{
	for (Int iy=0; iy<m_N; iy++) {
    	PixelC* dstPtr = m_mask[iy];
        const PixelC* srcPtr = rgchMask;
        for (Int ix=0; ix<m_N; ix++) 
        	*dstPtr++ = (*srcPtr++ != transpValue) ? 1 : 0;
                    
		rgchMask += stride;
	}
}

Void CSADCT::prepareInputBlock(Float **matDst, const Int *rgiSrc, Int stride)
{
	Float *rowDstPtr;
	const Int *rowSrcPtr;
	
	for (Int iy=0; iy<m_N; iy++) {
		rowDstPtr = matDst[iy];
		rowSrcPtr = rgiSrc + iy*stride;
		for (Int ix=0; ix<m_N; ix++) {
			*rowDstPtr++ = *rowSrcPtr++;
		}
	}
}

Void CSADCT::prepareInputBlock(Float **matDst, const PixelC *rgiSrc, Int stride)
{
	Float *rowDstPtr;
	const PixelC *rowSrcPtr;

	for (Int iy=0; iy<m_N; iy++) {
		rowDstPtr = matDst[iy];
		rowSrcPtr = rgiSrc + iy*stride;
		for (Int ix=0; ix<m_N; ix++) {
			*rowDstPtr++ = *rowSrcPtr++;
		}
	}
}


CFwdSADCT::CFwdSADCT(UInt nBits) :
	CFwdBlockDCT(nBits)
{
	m_dct_matrix = allocDctTable(m_N);
	initTrfTables();
}

CFwdSADCT::~CFwdSADCT()
{
	freeDctTable(m_dct_matrix, m_N);
}

Void CFwdSADCT::initTrfTables(Float scale)
{
	Float **mat, a, factcos;
	Int u, x;
	Int n;
	
	for (n=1; n<=m_N; n++) {
		mat = m_dct_matrix[n];
		factcos = M_PI/(2*n);
		a = scale * sqrt(2.0 / n);   
		for (u=0; u<n; u++) {
			for (x=0; x<n; x++) {
				mat[u][x] = a * cos(factcos*u*(2*x+1));
				if ( u == 0 )
					mat[u][x] /= M_SQRT2;
			}
		}
	}
}

// shape adaptive dct for inter coded blocks.
Void CFwdSADCT::apply(const Int* rgiSrc, Int nColSrc, Int* rgiDst, Int nColDst, const PixelC* rgchMask, Int nColMask, Int *lx)
{
	if (rgchMask) {
		prepareMask(rgchMask, nColMask);
		prepareInputBlock(m_in, rgiSrc, nColSrc);

// Schueuer HHI: added for fast_sadct
#ifdef _FAST_SADCT_
		fast_transform(m_out, lx, m_in, m_mask, m_N, m_N);
#else
		transform(m_out, lx, m_in, m_mask, m_N, m_N);
#endif

		copyBack(rgiDst, nColDst, m_out, lx);
	}
	else
		CBlockDCT::apply(rgiSrc, nColSrc, rgiDst, nColDst, NULL, 0, NULL);
}

// inverse shape adaptive dct for intra coded blocks.
Void CFwdSADCT::apply (const PixelC* rgchSrc, Int nColSrc, Int* rgiDst, Int nColDst, const PixelC* rgchMask, Int nColMask, Int *lx)
{
	if (rgchMask) {
		prepareMask(rgchMask, nColMask);
		prepareInputBlock(m_in, rgchSrc, nColSrc);


// HHI Schueuer: inserted for fastsadct
#ifdef _FAST_SADCT_
		fast_deltaDCTransform(m_out, lx, m_in, m_mask, m_N, m_N);	
#else
		deltaDCTransform(m_out, lx, m_in, m_mask, m_N, m_N);
#endif
		// dirty hack for the AC/DC prediction: The transparent pixels of the
        // first row and columns must be cleared acc. to the sadct proposal.
        memset(rgiDst, 0, m_N*sizeof(PixelI));
        PixelI* rgiDstPtr = rgiDst+nColDst;
        for (int i = 1; i < m_N; i++ ) {
        	*rgiDstPtr = 0;
            rgiDstPtr += nColDst;
        }       

		copyBack(rgiDst, nColDst, m_out, lx);
	}
	else 
		CBlockDCT::apply(rgchSrc, nColSrc, rgiDst, nColDst, NULL, 0, NULL);
}

Void CFwdSADCT::copyBack(PixelI *rgiDst, Int nColDst, Float **in, Int *lx)
{
	Float *rowSrcPtr;
	PixelI *rowDstPtr;
	int i, j;
	
	for (i=0; i<m_N && lx[i]; i++) {
		rowSrcPtr = in[i];
		rowDstPtr = rgiDst;
		for (j = 0; j < lx[i]; j++) { 
			*rowDstPtr = (*rowSrcPtr < 0) ? (PixelI)(*rowSrcPtr - .5) :
			(PixelI)(*rowSrcPtr + .5);
			rowDstPtr++;
			rowSrcPtr++;
		}
		rgiDst += nColDst;
	}
}

Void CFwdSADCT::deltaDCTransform(Float **out, Int *lx, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, j;
	Float mean_value;
	Int active_pels;

	out[0][0] = 0.0;
	
	/* compute meanvalue */
	mean_value = 0.0;
	active_pels  = 0;
	for (i = 0; i < bky; i++) {
		for (j = 0; j < bkx; j++) {
			active_pels+= mask[i][j];
			mean_value += in[i][j] * mask[i][j];
		}
	}
	
	if (active_pels)
		mean_value = mean_value / (Float)active_pels;
	mean_value = mean_value + 0.5;
	mean_value = (int)mean_value;
	
	for (i = 0; i < bky; i++) {
		for (j = 0; j < bkx; j++) {
			in[i][j] -= mean_value;
		}
	}
	
 	transform(out, lx, in, mask, bky, bkx);
	
	/* copy meanvalue to DC-coefficient */
	out[0][0] = mean_value * 8.0;
	
}

// HHI Schueuer: inserted for fast sadct
#ifdef _FAST_SADCT_
Void CFwdSADCT::fast_deltaDCTransform(Float **out, Int *lx, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, j, jmax, k;
	Float mean_value;
	Float *row, *row_coeff;
	Int active_pels;

	out[0][0] = 0.0;
	
	/* compute meanvalue */
	// mean_value = 0.0;
	// active_pels  = 0;

	fastshiftupTranspose(m_mat_tmp1, m_ly, in, mask, &mean_value, &active_pels, bky, bkx);

	if (active_pels) {
		mean_value /= (Float) active_pels;
		if (mean_value > 0)
			mean_value = (Int) (mean_value + 0.5);
		else
			mean_value = (Int) (mean_value - 0.5);
	}

	memset(lx, 0, sizeof(Int)*bky);
	
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		row = m_mat_tmp1[i];
		for (j = 0; j < jmax; j++)
			row[j] -= mean_value;
		switch (jmax) {
		case 1:
			c_buf[0] = row[0];
			break;      
		case 2:
			dct_vec2 (row, c_buf);
			break;
		case 3:
			dct_vec3 (row, c_buf);
			break;
		case 4:
			dct_vec4 (row, c_buf);
			break;      
		case 5:
			dct_vec5 (row, c_buf);
			break;
		case 6:
			dct_vec6 (row, c_buf);
			break;
		case 7:
			dct_vec7 (row, c_buf);
			break;
		case 8:
			dct_vec8 (row, c_buf);
			break;
	}      
	for (k=0; k<jmax; k++) {
      tmp_out[k][lx[k]] = c_buf[k];
      lx[k]++;
    }
  }
	
	/* and finally the horizontal transformation */
	for (i=0; i<bky && lx[i]; i++) {
		jmax = lx[i];
		row = out[i];
		row_coeff = tmp_out[i];
		switch (jmax) {
		case 1:
			*row = row_coeff[0];
			break;      
		case 2:
			dct_vec2 (row_coeff, row);
		 break;
		case 3:
			dct_vec3 (row_coeff, row);
			break;
		case 4:
			dct_vec4 (row_coeff, row);
			break;      
		case 5:
			dct_vec5 (row_coeff, row);
			break;
		case 6:
			dct_vec6 (row_coeff, row);
			break;
		case 7:
			dct_vec7 (row_coeff, row);
			break;
		case 8:
			dct_vec8 (row_coeff, row);
			break;