sadct.cpp

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

			}
		}
	}
#endif	

}


// HHI Schueuer: inserted for fast sadct
#ifdef _FAST_SADCT_
Void CInvSADCT::fast_deltaDCTransform(Float **out, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, k, jmax;
	Float **dest;
	Float *in_ptr;


	Float mean_value;
	Float e1 = 0.0, e2 = 0.0, e_dc = 0.0;

	mean_value = in[0][0] / 8.0;
	in[0][0] = 0.0;

	build_v_reorder_tbl(m_ly, out, mask, bky, bkx);
	
	build_h_reorder_tbl(m_lx, m_ly, m_mat_tmp1, bky, bkx);

	/* inverse horizontal transformation */
	for (i=0; i<bky && m_lx[i]; i++) {
		jmax = m_lx[i];
		in_ptr = in[i];
		
		dest = m_reorder_h[i];
		switch (jmax) {
		case 1:
			c_buf[0] = in_ptr[0];
			break;
		case 2:
			 idct_vec2 (in_ptr, c_buf);
			break;
		case 3:
			idct_vec3 (in_ptr, c_buf);
			break;      
		case 4:
			idct_vec4 (in_ptr, c_buf);
			break;
		case 5:
			idct_vec5 (in_ptr, c_buf);
			break;
		case 6:
			idct_vec6 (in_ptr, c_buf);
			break;      
		case 7:
			idct_vec7 (in_ptr, c_buf);
			break;
		case 8:
			idct_vec8 (in_ptr, c_buf);
			break;
		}
		if ( i == 0) {
			for (k = 0; k < jmax; k++) {
				e1 += sq[m_ly[k]] * c_buf[k];
				e2 += sq[m_ly[k]];
			}
			e_dc = e1 / e2;
			for (k=0; k<jmax; k++) 
				*dest[k] = c_buf[k] - e_dc;
		}
		else {
			for (k=0; k<jmax; k++)
				*dest[k] = c_buf[k];
		}
	}
 
	/* inverse vertical transformation */
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		in_ptr = m_mat_tmp1[i];
		
		dest = m_reorder_v[i];
		switch (jmax) {
		case 1:
			c_buf[0] = in_ptr[0];
			break;
		case 2:
			idct_vec2 (in_ptr, c_buf);
			break;
		case 3:
			idct_vec3 (in_ptr, c_buf);
			break;      
		case 4:
			idct_vec4 (in_ptr, c_buf);
			break;
		case 5:
			idct_vec5 (in_ptr, c_buf);
			break;
		case 6:
			idct_vec6 (in_ptr, c_buf);
			break;      
		case 7:
			idct_vec7 (in_ptr, c_buf);
			break;
		case 8:
			idct_vec8 (in_ptr, c_buf);
			break;
		}
		for (k=0; k<jmax; k++) 
			*dest[k] = c_buf[k] + mean_value;   
	}

}
#endif


/*
 *	inverse sadct 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::transform(Float **out, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, j, k, jmax;
	Float **trf_mat, **dest;
	Float c;
	Float *in_ptr;

	build_v_reorder_tbl(m_ly, out, mask, bky, bkx);
	
	build_h_reorder_tbl(m_lx, m_ly, m_mat_tmp1, bky, bkx);

	/* inverse horizontal transformation */
	for (i=0; i<bky && m_lx[i]; i++) {
		jmax = m_lx[i];
		trf_mat = m_idct_matrix[jmax];
		in_ptr = in[i];
		
		dest = m_reorder_h[i];
		for (k=0; k<jmax; k++) {
			for (c=0,j=0; j<jmax; j++) 
				c += trf_mat[k][j] * in_ptr[j];
			*dest[k] = c;
		}    
	}
	
	/* inverse vertical transformation */
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		trf_mat = m_idct_matrix[jmax];
		in_ptr = m_mat_tmp1[i];
		
		dest = m_reorder_v[i];
		for (k=0; k<jmax; k++) {
			for (c=0,j=0; j<jmax; j++) 
				c += trf_mat[k][j] * in_ptr[j];
			*dest[k] = c;
		}    
	}
}

// HHI Schueuer: added for fast_sadct
#ifdef _FAST_SADCT_
Void CInvSADCT::fast_transform(Float **out, Float **in, PixelC **mask, Int bky, Int bkx)
{
	Int i, k, jmax;
	Float **dest;
	Float *in_ptr;

	build_v_reorder_tbl(m_ly, out, mask, bky, bkx);
	
	build_h_reorder_tbl(m_lx, m_ly, m_mat_tmp1, bky, bkx);

	/* inverse horizontal transformation */
	for (i=0; i<bky && m_lx[i]; i++) {
		jmax = m_lx[i];
		in_ptr = in[i];
		
		dest = m_reorder_h[i];
		switch (jmax) {
		case 1:
			c_buf[0] = in_ptr[0];
			break;
		case 2:
			 idct_vec2 (in_ptr, c_buf);
			break;
		case 3:
			idct_vec3 (in_ptr, c_buf);
			break;      
		case 4:
			idct_vec4 (in_ptr, c_buf);
			break;
		case 5:
			idct_vec5 (in_ptr, c_buf);
			break;
		case 6:
			idct_vec6 (in_ptr, c_buf);
			break;      
		case 7:
			idct_vec7 (in_ptr, c_buf);
			break;
		case 8:
			idct_vec8 (in_ptr, c_buf);
			break;
		}
		for (k=0; k<jmax; k++)
			*dest[k] = c_buf[k];    
	}
 
	/* inverse vertical transformation */
	for (i=0; i<bkx && m_ly[i]; i++) {
		jmax = m_ly[i];
		in_ptr = m_mat_tmp1[i];
		
		dest = m_reorder_v[i];
		switch (jmax) {
		case 1:
			c_buf[0] = in_ptr[0];
			break;
		case 2:
			idct_vec2 (in_ptr, c_buf);
			break;
		case 3:
			idct_vec3 (in_ptr, c_buf);
			break;      
		case 4:
			idct_vec4 (in_ptr, c_buf);
			break;
		case 5:
			idct_vec5 (in_ptr, c_buf);
			break;
		case 6:
			idct_vec6 (in_ptr, c_buf);
			break;      
		case 7:
			idct_vec7 (in_ptr, c_buf);
			break;
		case 8:
			idct_vec8 (in_ptr, c_buf);
			break;
		}
		for (k=0; k<jmax; k++)
			*dest[k] = c_buf[k];     
	}
}
#endif

// Schueuer HHI: added fast_sadct
#ifdef _FAST_SADCT_
Int CInvSADCT::idct_vec2 (Double *in, Double *out)
{

  out[0] = (in[0] + in[1]) * f0_2;
  out[1] = (in[0] - in[1]) * f0_2;

  return 0;
}

Int CInvSADCT::idct_vec3 (Double *in, Double *out)
{
  Double b[3];

  b[0] = in[0] * f0_3;
  b[1] = b[0] + in[2] * f2_3;
  b[2] = in[1] * f1_3;

  out[0] = b[1] + b[2];
  out[2] = b[1] - b[2];
  out[1] = b[0] - in[2] * f3_3;

  return 0;
}

Int CInvSADCT::idct_vec4 (Double *in, Double *out)
{
  Double b[4];

  /* stage 1 */

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

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

  return 0;
}

Int CInvSADCT::idct_vec5 (Double *in, Double *out)
{
  Double b[5];

  b[0] = in[0] * f0_5;
  b[1] = b[0] + in[2] * f3_5 + in[4] * f4_5;
  b[2] = in[1] * f1_5 + in[3] * f2_5;
  b[3] = b[0] - in[4] * f3_5 - in[2] * f4_5;
  b[4] = in[1] * f2_5 - in[3] * f1_5;

  out[0] = b[1] + b[2];
  out[1] = b[3] + b[4];
  out[2] = b[0] + (in[4] - in[2]) * f5_5;
  out[3] = b[3] - b[4];
  out[4] = b[1] - b[2];

  return 0;
}

Int CInvSADCT::idct_vec6 (Double *in, Double *out)
{
  Double b[10];
  Double b1[7];

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

  /* stage 2 */
  b1[0] = b[0] + b[4];
  b1[1] = b1[0] + b[2];
  b1[2] = b[1] + b[3] + b[5];
  b1[3] = b[0] - b[8];
  b1[4] = b[6] - b[3];
  b1[5] = b1[0] - b[2];
  b1[6] = b[7] - b[3] + b[9];

  /* stage 3 */
  out[0] = b1[1] + b1[2];
  out[1] = b1[3] + b1[4];
  out[2] = b1[5] + b1[6];
  out[3] = b1[5] - b1[6];
  out[4] = b1[3] - b1[4];
  out[5] = b1[1] - b1[2];

  return 0;
}

Int CInvSADCT::idct_vec7 (Double *in, Double *out)
{
  Double b[20];
  Double b1[6];

  /* stage 1 */
  b[0] = in[0] * f0_7;
  b[1] = in[1] * f1_7;
  b[2] = in[1] * f3_7;
  b[3] = in[1] * f5_7;
  b[4] = in[2] * f2_7;
  b[5] = in[2] * f6_7;
  b[6] = in[2] * f4_7;
  b[7] = in[3] * f3_7;
  b[8] = in[3] * f5_7;
  b[9] = in[3] * f1_7;
  b[10] = in[4] * f4_7;
  b[11] = in[4] * f2_7;
  b[12] = in[4] * f6_7;
  b[13] = in[5] * f5_7;
  b[14] = in[5] * f1_7;
  b[15] = in[5] * f3_7;
  b[16] = in[6] * f6_7;
  b[17] = in[6] * f4_7;
  b[18] = in[6] * f2_7;
  b[19] = (in[4] - in[2] - in[6]) * f7_7;
  
  /* stage 2 */
  b1[0] = b[0] + b[4] + b[10] + b[16];
  b1[1] = b[1] + b[7] + b[13];
  b1[2] = b[0] + b[5] - b[11] - b[17];
  b1[3] = b[2] - b[8] - b[14];
  b1[4] = b[0] - b[6] - b[12] + b[18];
  b1[5] = b[3] - b[9] + b[15];

  /* stage 3 */

  out[0] = b1[0] + b1[1];
  out[1] = b1[2] + b1[3];
  out[2] = b1[4] + b1[5];
  out[3] = b[0] + b[19];
  out[4] = b1[4] - b1[5];
  out[5] = b1[2] - b1[3];
  out[6] = b1[0] - b1[1];  
  
  return 0;
}

Int CInvSADCT::idct_vec8 (Double *in, Double *out)
{

  Int j1, j;
  Double tmp[8], tmp1[8];
  Double e, f, g, h;

  /* stage 1 for k = 4,5,6,7 */
  e = in[1] * f7_8 - in[7] * f1_8;
  h = in[7] * f7_8 + in[1] * f1_8;
  f = in[5] * f3_8 - in[3] * f5_8;
  g = in[3] * f3_8 + in[5] * f5_8;

  /* stage 1+2 for k = 0,1,2,3 */
  tmp1[0] = (in[0] + in[4]) * f4_8;
  tmp1[1] = (in[0] - in[4]) * f4_8;
  tmp1[2] = in[2] * f6_8 - in[6] * f2_8;
  tmp1[3] = in[6] * f6_8 + in[2] * f2_8;

  tmp[4] = e + f;
  tmp1[5] = e - f;
  tmp1[6] = h - g;
  tmp[7] = h + g;
  
  tmp[5] = (tmp1[6] - tmp1[5]) * f0_8;
  tmp[6] = (tmp1[6] + tmp1[5]) * f0_8;
  tmp[0] = tmp1[0] + tmp1[3];
  tmp[1] = tmp1[1] + tmp1[2];
  tmp[2] = tmp1[1] - tmp1[2];
  tmp[3] = tmp1[0] - tmp1[3];

  /* stage 4 */
  for (j = 0; j < 4; j++) {
    j1 = 7 - j;
    out[j] = tmp[j] + tmp[j1];
    out[j1] = tmp[j] - tmp[j1];
  }
  
  return 0;
}
#endif

	
Void CInvSADCT::build_v_reorder_tbl(Int *l_y, 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] ) {
				m_reorder_v[iy_out][ix_out++] = in[iy]+ix;
				l++;
			}
		}
		if ( l ) 
			l_y[iy_out++] = l;
	}
	
	/* initialize the length of the unoccupied columns resp. rows to zero. */
	for (ix=iy_out; ix<bkx; ix++) 
		l_y[ix] = 0;
	
}
 
Void CInvSADCT::build_h_reorder_tbl(Int *l_x, const Int *l_y, Float **in, Int bky, Int bkx)
{
	Int i, k, jmax;
	Int col_ind;
	Float *row;
	
	memset(l_x, 0, sizeof(Int)*bky);
	for (i=0; i<bkx && l_y[i]; i++) {
		jmax = l_y[i];
		row = in[i];
		for (k=0; k<jmax; k++) {
			col_ind = l_x[k];
			m_reorder_h[k][col_ind] = row + k;
			l_x[k]++;
		}
	}
}

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

Void CInvSADCT::copyBack(PixelC *rgchDst, Int nColDst, Float **in, PixelC **mask)
{
	Float *rowSrcPtr;
	PixelC *rowMaskPtr;
	PixelC *rowDstPtr;
	PixelI iValue;
	int i, j;
	
	for (i=0; i<m_N; i++) {
		rowSrcPtr = in[i];
		rowMaskPtr = mask[i];
		rowDstPtr = rgchDst;
		for (j = 0; j < m_N; j++) { 
			if (rowMaskPtr) {
				iValue = (*rowSrcPtr < 0) ? (PixelI)(*rowSrcPtr - .5) :
			(PixelI)(*rowSrcPtr + .5);
				*rowDstPtr = m_rgchClipTbl [iValue];
			}
			rowMaskPtr++;
			rowDstPtr++;
			rowSrcPtr++;
		}
		rgchDst += nColDst;
	}
}