sadct.cpp

《Visual C++小波变换技术与工程实践》靳济芳编著的光盘程序。

		}
	}
 
	

	
	/* copy meanvalue to DC-coefficient */
	out[0][0] = mean_value * 8.0;
	
}
#endif

Void CFwdSADCT::transform(Float **out, Int *lx, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, j, jmax, k;
	Float **trf_mat, *row;
	Float c;
  
	shiftupTranspose(m_mat_tmp1, m_ly, in, mask, bky, bkx);
	
	memset(lx, 0, sizeof(Int)*bky);
	
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		trf_mat = m_dct_matrix[jmax];
		row = m_mat_tmp1[i];
		for (k=0; k<jmax; k++) {
			for (c=0,j=0; j<jmax; j++) 
				c += trf_mat[k][j] * row[j];
			out[k][lx[k]] = c;
			lx[k]++;
		}
	}
	
	/* and finally the horizontal transformation */
	for (i=0; i<bky && lx[i]; i++) {
		jmax = lx[i];
		trf_mat = m_dct_matrix[jmax];
		memcpy(m_row_buf, out[i], jmax*sizeof(Float));
		row = out[i];
		for (k=0; k<jmax; k++) {
			for (c=0,j=0; j<jmax; j++) 
				c += trf_mat[k][j] * m_row_buf[j];
			*row++ = c;
		}    
	}

}

// Schueuer HHI : added for fast_sadct
#ifdef _FAST_SADCT_
Void CFwdSADCT::fast_transform(Float **out, Int *lx, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, jmax, k;
	Double *row, *row_coeff;
  
	shiftupTranspose(m_mat_tmp1, m_ly, in, mask, bky, bkx);
	
	memset(lx, 0, sizeof(Int)*bky);
	
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		row = m_mat_tmp1[i];
	 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;
		}
	}

}
#endif

//  HHI Schueuer: inserted for fast sadct
#ifdef _FAST_SADCT_
Void CFwdSADCT::fastshiftupTranspose(Float **out, Int *ly, Float **in, PixelC **mask, 
							Float *mean, Int *active_pels, Int bky, Int bkx)
{
  Int iy_out = 0, ix_out;
  Int iy, ix, l;

  *mean = 0.0;
  *active_pels = 0;

  for (ix = 0; ix < bkx; ix++) {
    ix_out = l = 0;
    for (iy = 0; iy < bky; iy++) {
      if ( mask[iy][ix] ) {
		out[iy_out][ix_out++] = in[iy][ix];
		*mean += in[iy][ix];
		l++;
      }
    }
    if ( l ) { 
      ly[iy_out++] = l;
      *active_pels += l;
    }
  }	
  /* initialize the length of the unoccupied columns to zero. The term 
     column refers to the pel positions in `in'.  In `out' columns are
     saved as rows (transposition) to speed up calculation. */
  for (ix=iy_out; ix<bkx; ix++) 
    ly[ix] = 0;

}
#endif

/*
 *	`sadct_shiftup_transpose' shifts upwards pels marked in `mask' towards
 *	the top margin of the rectangular block.  If a column consists of more
 *	than one stretch of pels, the gaps between these stretches will be
 *	eleminated.  The resulting columns are transposed and type converted
 *	into `out' and are occupying the upper part of the block without
 *	any gaps in vertical direction.
 *
 */
Void CFwdSADCT::shiftupTranspose(Float **out, Int *ly, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int iy_out = 0, ix_out;
	Int iy, ix, l;
	
	for (ix=0; ix<bkx; ix++) {
		ix_out = l = 0;
		for (iy=0; iy<bky; iy++) {
			if ( mask[iy][ix] ) {
				out[iy_out][ix_out++] = in[iy][ix];
				l++;
			}
		}
		if ( l ) 
			ly[iy_out++] = l;
	}
	
	/* initialize the length of the unoccupied columns to zero. The term 
	column refers to the pel positions in `in'.  In `out' columns are
	saved as rows (transposition) to speed up calculation. */
	for (ix=iy_out; ix<bkx; ix++) 
		ly[ix] = 0;
	
}

// Schueuer HHI : added for fast_sadct
#ifdef _FAST_SADCT_
Int CFwdSADCT::dct_vec2 (Double *vec, Double *coeff)
{

  coeff[0] = (vec[0] + vec[1]) * f0_2;
  coeff[1] = (vec[0] - vec[1]) * f0_2;

  return 0;
}

Int CFwdSADCT::dct_vec3 (Double *vec, Double *coeff)
{
  Double b;

  b = vec[0] + vec[2];

  coeff[0] = (vec[1] + b) * f0_3;
  coeff[1] = (vec[0] - vec[2]) * f1_3;
  coeff[2] = b * f2_3 -  vec[1] * f3_3;

  return 0;
}

Int CFwdSADCT::dct_vec4 (Double *vec, Double *coeff)
{
  Double b[4];

  /* stage 1 */
  b[0] = vec[0] + vec[3];
  b[1] = vec[1] + vec[2];
  b[2] = vec[1] - vec[2];
  b[3] = vec[0] - vec[3];

  /* stage 2 */
  coeff[0] = (b[0] + b[1]) * f0_4;
  coeff[2] = (b[0] - b[1]) * f0_4;
  coeff[1] = b[2] * f2_4 + b[3] * f1_4;
  coeff[3] = b[3] * f2_4 - b[2] * f1_4;

  return 0;
}

Int CFwdSADCT::dct_vec5 (Double *vec, Double *coeff)
{
  Double b[5];

  b[0] = vec[0] + vec[4];
  b[1] = vec[0] - vec[4];
  b[2] = vec[1] + vec[3];
  b[3] = vec[1] - vec[3];
  b[4] = vec[2] * f5_5;

  coeff[0] = (b[0] + b[2] + vec[2]) * f0_5;
  coeff[1] = b[1] * f1_5 + b[3] * f2_5;
  coeff[2] = b[0] * f3_5 - b[2] * f4_5 - b[4];
  coeff[3] = b[1] * f2_5 - b[3] * f1_5;
  coeff[4] = b[0] * f4_5 - b[2] * f3_5 + b[4];

  return 0;
}

Int CFwdSADCT::dct_vec6 (Double *vec, Double *coeff)
{
  Double b[8];

  /* stage 1 */
  b[0] = vec[0] + vec[5];
  b[1] = vec[0] - vec[5];
  b[2] = vec[1] + vec[4];
  b[3] = vec[1] - vec[4];
  b[4] = vec[2] + vec[3];
  b[5] = vec[2] - vec[3];
  b[6] = b[3] * f0_6;
  b[7] = b[0] + b[4];

  /* stage 2 */
  coeff[0] = (b[7] + b[2]) * f0_6;
  coeff[1] = b[1] * f1_6 + b[5] * f4_6 + b[6];
  coeff[2] = (b[0] - b[4]) * f2_6;
  coeff[3] = (b[1] - b[3] -b[5]) * f0_6;
  coeff[4] = b[7] * f3_6 - b[2] * f5_6;
  coeff[5] = b[1] * f4_6 + b[5] * f1_6 - b[6];

  return 0;
}

Int CFwdSADCT::dct_vec7 (Double *vec, Double *coeff)
{
  Double b[7];

  b[0] = vec[0] + vec[6];
  b[1] = vec[0] - vec[6];
  b[2] = vec[1] + vec[5];
  b[3] = vec[1] - vec[5];
  b[4] = vec[2] + vec[4];
  b[5] = vec[2] - vec[4];
  b[6] = vec[3] * f7_7;

  coeff[0] = (b[0] + b[2] + b[4] + vec[3]) * f0_7;
  coeff[1] = b[1] * f1_7 + b[3] * f3_7 + b[5] * f5_7;
  coeff[2] = b[0] * f2_7 + b[2] * f6_7 - b[4] * f4_7 - b[6];
  coeff[3] = b[1] * f3_7 - b[3] * f5_7 - b[5] * f1_7;
  coeff[4] = b[0] * f4_7 - b[2] * f2_7 - b[4] * f6_7 + b[6];
  coeff[5] = b[1] * f5_7 - b[3] * f1_7 + b[5] * f3_7;
  coeff[6] = b[0] * f6_7 - b[2] * f4_7 + b[4] * f2_7 - b[6];
  
  return 0;
}

Int CFwdSADCT::dct_vec8 (Double *vec, Double *coeff)
{
  Int   j1, j;
  Double b[8];
  Double b1[8];

  /* stage 1 */
  for (j = 0; j < 4; j++)
    {
      j1 = 7 - j;
      b1[j] = vec[j] + vec[j1];
      b1[j1] = vec[j] - vec[j1];
    }

  /* stage 2 */
  b[0] = b1[0] + b1[3];
  b[1] = b1[1] + b1[2];
  b[2] = b1[1] - b1[2];
  b[3] = b1[0] - b1[3];
  b[4] = b1[4];
  b[5] = (b1[6] - b1[5]) * f0_8;
  b[6] = (b1[6] + b1[5]) * f0_8;
  b[7] = b1[7];

  /* stage 3/4 for the coeff. 0,2,4,6 */
  coeff[0] = (b[0] + b[1]) * f4_8;
  coeff[4] = (b[0] - b[1]) * f4_8;
  coeff[2] = b[2] * f6_8 + b[3] * f2_8;
  coeff[6] = b[3] * f6_8 - b[2] * f2_8;

  /* stage 3 */
  b1[4] = b[4] + b[5];
  b1[7] = b[7] + b[6];
  b1[5] = b[4] - b[5];
  b1[6] = b[7] - b[6];

  /* stage 4 for coeff. 1,3,5,7 */
  coeff[1] = b1[4] * f7_8 + b1[7] * f1_8;
  coeff[5] = b1[5] * f3_8 + b1[6] * f5_8;
  coeff[7] = b1[7] * f7_8 - b1[4] * f1_8;
  coeff[3] = b1[6] * f3_8 - b1[5] * f5_8;
    
  return 0;
}
#endif


CInvSADCT::CInvSADCT(UInt nBits) :
	CInvBlockDCT(nBits)
{
	m_reorder_h = allocReorderTable(m_N);
	m_reorder_v = allocReorderTable(m_N);
	m_idct_matrix = allocDctTable(m_N);

	initTrfTables();
}

CInvSADCT::~CInvSADCT()
{
	freeDctTable(m_idct_matrix, m_N);
	freeReorderTable(m_reorder_h, m_N);
	freeReorderTable(m_reorder_v, m_N);
}

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

// inverse shape adaptive dct for Intra coded blocks.
Void CInvSADCT::apply (const Int* rgiSrc, Int nColSrc, PixelC* rgchDst, Int nColDst, const PixelC* rgchMask, Int nColMask)
{
	if (rgchMask) { 
		prepareMask(rgchMask, nColMask);
		prepareInputBlock(m_in, rgiSrc, nColSrc);


//  HHI Schueuer: inserted for sadct
#ifdef _FAST_SADCT_
		fast_deltaDCTransform(m_out, m_in, m_mask, m_N, m_N);
#else
		deltaDCTransform(m_out, 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(rgchDst, 0, m_N*sizeof(PixelC));
        PixelC* rgchDstPtr = rgchDst+nColDst;
#ifdef __TRACE_DECODING_
		// This complete clean is for tracing/debugging only, because I don't want to see memory garbage in the trace file.
        for (int i=1; i<m_N; i++ ) {
        	memset(rgchDstPtr, 0, m_N*sizeof(PixelC));
            rgchDstPtr += nColDst;
        }
#else
        for (int i=1; i<m_N; i++ ) {
        	*rgchDstPtr = 0;
            rgchDstPtr += nColDst;
        }
#endif        
		copyBack(rgchDst, nColDst, m_out, m_mask);	
	}
	else 
		CBlockDCT::apply(rgiSrc, nColSrc, rgchDst, nColDst);
}

// inverse shape adative dct for inter coded blocks.  
Void CInvSADCT::apply (const PixelI* rgiSrc, Int nColSrc, PixelI* rgiDst, Int nColDst, const PixelC* rgchMask, Int nColMask)
{
	if (rgchMask) {
		// boundary block assumed, for details refer to the comment
		// inside the body of the other apply method.
		prepareMask(rgchMask, nColMask);
		prepareInputBlock(m_in, rgiSrc, nColSrc);

// Schueuer HHI: added for fast_sadct
#ifdef _FAST_SADCT_
		fast_transform(m_out, m_in, m_mask, m_N, m_N);
#else
 		transform(m_out, 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;
#ifdef __TRACE_DECODING_
		// This complete clean out is for tracing/debugging only, because I don't want to see memory garbage in the trace file.
        for (int i=1; i<m_N; i++ ) {
        	memset(rgiDstPtr, 0, m_N*sizeof(PixelI));
            rgiDstPtr += nColDst;
        }
#else
        for (int i=1; i<m_N; i++ ) {
        	*rgiDstPtr = 0;
            rgiDstPtr += nColDst;
        }
#endif       
		copyBack(rgiDst, nColDst, m_out, m_mask);
	}
	else 
		CBlockDCT::apply(rgiSrc, nColSrc, rgiDst, nColDst);
}


/*
 *	Inverse shape adaptive transformation of block `in'.  The spatial positions
 *	of valid pels are marked in `mask' by 1.  Please note that the
 *	dct coefficients encoding those pels are expected to be found in 
 *	the upper left corner of block `in'.
 *	  
 *
 *	The following drawing explains the relation between `in', `out'
 *	and `mask':
 *
 *	in -> I I I - - - - - 
 *	 	  I I - - - - - -
 *		  I - - - - - - -
 *		  - - ...
 *				                 out ->    - - - - O - - -         
 *	mask ->   - - - - 1 - - -         	   - - O O - - - -
 *		      - - 1 1 - - - -		       - - O O - - - -
 *		      - - 1 1 - - - -		       - - - O - - - -
 *		      - - - 1 - - - -		       - - - - - - - -
 *		      - - - - - - - -		       - - ...
 *		      - - ...
 *
 */
Void CInvSADCT::deltaDCTransform(Float **out, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, j, k, l;

	Float mean_value, check_sum = 0.0;
	Float ly_sum = 0.0, ddx;

#if 1
	/* reconstruction of meanvalue and zero setting of ddc */
	mean_value = (Int) ((in[0][0] / 8.0) + 0.5);
	in[0][0] = 0.0;
#endif	

	// Klaas Schueuer
	for (i = 0; i < 8; i++) 
		for (j = 0; j < 8; j++) 
			out [i][j] = 0.0;
	// end 

#if defined(__DEBUG_SADCT_) && !defined(NDEBUG)
	cout << "deltaDCTransform(): mean_value=" << mean_value << endl;
#endif

 	transform(out, in, mask, bky, bkx);
	
	/* computing of checksum and ddc correction */
	check_sum = 0;
	for (i = 0; i < bky; i++) {
		for (j = 0; j < bkx; j++)
			if (mask[i][j])
				check_sum += out[i][j];
	}	
	
	for (i=0; i<bkx; i++) {
		if (m_ly[i])
			ly_sum += sqrt ((Float)m_ly[i]);
	}
#if defined(__DEBUG_SADCT_) && !defined(NDEBUG)
	cout << "chksum=" << check_sum << "  ly_sum=" << ly_sum << '\n';
#endif
#if 1	
	k = 0;
	for (i = 0; i < bkx; i++) {
		l = 0;
		for (j = 0; j < bky; j++) {
			if (mask[j][i]) {
				if (l==0) {
					k++;
					l++;
					if (check_sum>0)  
						ddx = (Int) ((1.0 / (sqrt ((Float) m_ly[k-1]) * ly_sum) * check_sum ) + 0.5);
					else
						ddx = (Int) ((1.0 / (sqrt ((Float) m_ly[k-1]) * ly_sum) * check_sum ) - 0.5);
				}
#if defined(__DEBUG_SADCT_) && !defined(NDEBUG)
				cout << "delta correction: out[" << j << "][" << i << "]=" 
					 << out[j][i] << " ddx=" << ddx << '\n';
#endif
				out[j][i] += mean_value - ddx;