ITU - Telecommunications Standardization Sector
STUDY GROUP 16
Video Coding Experts Group (Question 15)
_________________
Eighth Meeting: Berlin, Germany, 3-6 August, 1999	Document  Q15-H-29
Filename: q15h29.rtf
Generated: 27 July ’98	

Question:	Q.15/SG16	
Source:	Wolfgang Sörgel and 
Thomas Wiegand
Telecommunications Laboratory
University Erlangen-Nuremberg
Erlangen, Germany	Tel:
Fax:
Email: 	+49 9131 8527781
+49 9131 8528849
wsoergel@nt.e-technik.uni-erlangen.de	
Title:	Verification of IDCT specification	
Purpose:	Information	
_____________________________
Summary
This document reports the results of the verification of a fixed length IDCT specification as proposed in Q15-G-34. Following the procedure outlined in Annex A "Inverse Transform Accuracy Specification" we show that the proposed algorithm is compliant with all accuracy requirements.

Results

Random integer pixel data were generated using the routine specified in A.9 / H.263, DCT transformed according to A.2 / H.263 and rounded and clipped to the range -2048 to 2047. For each 8 by 8 block of test data a reference IDCT was computed according to section 6.2.4 of H.263 and rounded to the range -256 to 255. For each transformed block of test data, these results were compared to the results obtained from the IDCT to be tested. 
Test data sets of 10000 blocks of size 8 by 8 were generated for the range -256 to 255, -5 to 5 and -300 to 300. For each test block and each test block with changed sign
•	the peak error was at most 1 in amplitude.
•	the mean square error per pixel did not exceed 0.06.
•	the overall mean square error did not exceed 0.02.
•	the mean error per pixel did not exceed 0.015 in magnitude.
•	the overall mean error did not exceed 0.0015 in magnitude.
All zeros in produces all zeros out.

Table 1 shows the resulting maximal mean and mean square errors per pixel and the overall mean and mean square errors. Table 2 shows the results for the same random data as used before but with changed sign for each pixel.
 
Table 1. Errors Between Reference IDCT and Proposed IDCT.

Data Range	max(MSE) / pixel	overall MSE	max(ME) / pixel	overall ME	
- 256 to 255	0.0125	0.009373	0.003	0.000011	
-5 to 5	0.0135	0.010091	0.0086	0.000034	
-300 to 300	0.0102	0.008188	0.0025	0.000050	

Table 2. Errors Between Reference IDCT and Proposed IDCT, Sign of Test Data Changed

Data Range	max(MSE) / pixel	overall MSE	max(ME) / pixel	overall ME	
- 256 to 255	0.0116	0.009430	0.0022	0.000045	
-5 to 5	0.0129	0.010078	0.0095	0.000047	
-300 to 300	0.0105	0.00823	0.0022	0.000030	

Conclusions

A test to verify the compliance of the IDCT proposed in Q 15-G-34 with the "Inverse Transform Accuracy Specification" as specified in Annex A / H.263 has been conducted. According to the results, the proposed IDCT is found to be compliant with the specified error tolerances.

Appendix

The following appendix contains the complete source code in "C" for the program used to conduct the experiments:

/* IDCT Compliance test c.f. ITU-T Recommendation H.263 Annex A */
/* Compilation: cc -o test_idct test_idct.c idct.c -lm */ 
/* Usage: test_idct n_datasets [negate]
   n_datasets: Number of tests (10000 for Annex A)
   negate: if given, change sign of each data pixel (A.9) */
/* Wolfgang Soergel <wsoergel@lnt.de> July 1999*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*Type definitions and prototypes of idct in question: */
typedef short int REGISTER;
void idct(REGISTER block[64]);

/* Type definition for tests: */
#define MAX(a,b) (a>b ? a:b);

/* Reference dct */
void r_dct(REGISTER block[64],REGISTER coeffs[64]);

/* Reference idct: */
void r_idct(REGISTER coeffs[64], REGISTER block[64]);

/* round and clip to integer range */
int round_clip(double x, int l, int h);

/* random number generator c.f. Annex A: */
long itu_rand(long l, long h);

long itu_rand(long L, long H)
{
  static long randx=1;
  static double z=(double)0x7fffffff;

  long i,j;
  double x;
  
  randx= (randx*1103515245) + 12345;
  i=randx & 0x7ffffffe;
  x=((double)i)/z;
  x *= (double)(L+H+1);
  j=x;
  
  return(j-L);
}

/* Reference IDCT, naive implementation: */void r_idct(REGISTER coeffs[64], REGISTER block[64]){  int x,y,u,v, ind;  double tmp_u, tmp_v;  double C_u,C_v;    double M_PI_16 = M_PI/16.0;  for(x=0;x<8;x++){    for(y=0;y<8;y++){      tmp_u=0.0;
      for(u=0;u<8;u++){
	C_u= (u > 0 ? 1: M_SQRT1_2 );
	tmp_v=0.0;
	for(v=0;v<8;v++){
	  C_v=(v > 0 ? 1:M_SQRT1_2);
	  tmp_v += C_v*coeffs[8*u+v]*cos(M_PI_16*(2*y+1)*v);
	}
	tmp_u += tmp_v*C_u*cos(M_PI_16*(2*x+1)*u);
      }
      block[8*x+y]=round_clip(tmp_u/4.0,-256,255);
     }
  }
}

void r_dct(REGISTER block[64], REGISTER coeffs[64]){  int x,y,u,v, ind;  double C_u, C_v;  double tmp_x, tmp_y;  double M_PI_16 = M_PI/16.0;  for(u=0;u<8;u++){    C_u= (u > 0 ? 1 : M_SQRT1_2 );
    for(v=0;v<8;v++){
      C_v= (v > 0 ? 1 : M_SQRT1_2 );
      tmp_x=0.0;
      for(x=0;x<8;x++){
	tmp_y=0.0;
	for(y=0;y<8;y++){
	  tmp_y+=block[8*x+y]*cos(M_PI_16*(2*y+1)*v);
	}
	tmp_x+=tmp_y*cos(M_PI_16*(2*x+1)*u);
      }
      coeffs[8*u+v]=round_clip(tmp_x/4.0 * C_u * C_v,-2048,2047);       
    }   
  }
}
int round_clip(double x, int l, int h){  int i;    i = x>0 ?  (int) (x+0.5) : (int) (x-0.5) ;    return (i > l ? (i < h ? i : h) : l);  }int main(int argc, char ** argv)
{
  REGISTER block[64], coeffs[64], ref_block[64], test_block[64];
  double MSE, mse[64], ME, me[64];
  REGISTER peak_error, err, sign=1;
  int u,v,i,n, k,n_datasets=1;  
  int L[3],H[3];
  
  L[0]=256;L[1]=5;L[2]=300;
  H[0]=255;H[1]=5;H[2]=300;
  
  if(argc>1)
    n_datasets=atoi(argv[1]);
  if(argc>2)
    sign=-1;;

  fprintf(stderr,"idct accuracy check c.f. H.263 Annex A, sign=%i\n\n",sign);
  for(k=0;k<3;k++){

    /* reset Error meassures: */	
    for(i=0;i<64;i++){
      mse[i]=0.0;
      me[i]=0.0;
    }
    MSE=0.0;
    ME=0.0;
    peak_error=0;
    
    /* Do the experiments */
    for(n=0;n<n_datasets;n++){
      
      /*Generate random block */
      for(u=0;u<8;u++)
	for(v=0;v<8;v++)
	  block[8*u+v]= sign * itu_rand(L[k],H[k]);

      /* Compute the reference dct: */
      r_dct(block,coeffs);

      /* Compute the reference idct: */ 
      r_idct(coeffs, ref_block);

      /* Compute the idct */
      idct(coeffs);

      /* Meassurements */
      for(i=0;i<64;i++){
	err=coeffs[i]-ref_block[i];
	  /* Check for peak error */
	peak_error= MAX(peak_error,abs(err));
	if(peak_error>1)
	  fprintf(stderr,"Failure:[%i,%i],n=%i: peak_error[%i]=%i\n",
		  -L[k],H[k],n,i,coeffs[i]-ref_block[i]);
	
	/* Mean Square Error */
	mse[i]+=(err)*(err);
	/* Mean Error */
	me[i]+=err;
      }      
    }
    
    fprintf(stderr,"\n\nIntervall: [%i,%i]\n\n",-L[k],H[k]);
    /* Assesment of errors */
    for(i=0;i<64;i++){
      me[i]/=(double)n_datasets;
      fprintf(stderr,"i=%i, me=%f ",i,me[i]);
      if(fabs(me[i])>0.015)
	fprintf(stderr,"Failure:[%i,%i] mean error[%i]=%f\n", 
		-L[k],H[k],i,me[i]);  
      
      mse[i]/=(double)n_datasets;
      fprintf(stderr,"mse=%f\n",mse[i]);
      if(fabs(mse[i])>0.06)
	fprintf(stderr,"Failure:[%i,%i] mean square error[%i]=%f\n",
		-L[k],H[k],i,mse[i]); 
      
      ME+=me[i];
      MSE+=mse[i];      
    }
    ME/=64.0;
    MSE/=64.0;
    fprintf(stderr,"Overall Mean Error[%i,%i]=%f , Mean Squared Error = %f\n",
	    -L[k],H[k],ME,MSE);    
    if(fabs(ME>0.0015))
      fprintf(stderr,"Failure:[%i,%i] overall mean error=%f\n",
	      -L[k],H[k],mse[i]); 
    if(fabs(MSE>0.02))
      fprintf(stderr,"Failure:[%i,%i] overall mean square error=%f\n",
	      -L[k],H[k],mse[i]); 
    
  }

  return 0;
}