quicktime/encore50/text_dct.c

   1
   2 /**************************************************************************
   3  *                                                                        *
   4  * This code is developed by Adam Li.  This software is an                *
   5  * implementation of a part of one or more MPEG-4 Video tools as          *
   6  * specified in ISO/IEC 14496-2 standard.  Those intending to use this    *
   7  * software module in hardware or software products are advised that its  *
   8  * use may infringe existing patents or copyrights, and any such use      *
   9  * would be at such party's own risk.  The original developer of this     *
  10  * software module and his/her company, and subsequent editors and their  *
  11  * companies (including Project Mayo), will have no liability for use of  *
  12  * this software or modifications or derivatives thereof.                 *
  13  *                                                                        *
  14  * Project Mayo gives users of the Codec a license to this software       *
  15  * module or modifications thereof for use in hardware or software        *
  16  * products claiming conformance to the MPEG-4 Video Standard as          *
  17  * described in the Open DivX license.                                    *
  18  *                                                                        *
  19  * The complete Open DivX license can be found at                         *
  20  * http://www.projectmayo.com/opendivx/license.php .                      *
  21  *                                                                        *
  22  **************************************************************************/
  23
  24 /**************************************************************************
  25  *
  26  *  text_dct.c
  27  *
  28  *  Copyright (C) 2001  Project Mayo
  29  *
  30  *  Adam Li
  31  *
  32  *  DivX Advance Research Center <darc@projectmayo.com>
  33  *
  34  **************************************************************************/
  35
  36 /* This file contains some functions for fDCT/iDCT transformation.        */
  37 /* Some codes of this project come from SSG MPEG-2 implementation.        */
  38 /* Please see seperate acknowledgement file for a list of contributors.   */
  39
  40 #include <math.h>
  41
  42 /* The first part of it is for the forward DCT */
  43
  44 #ifndef PI
  45 # ifdef M_PI
  46 #  define PI M_PI
  47 # else
  48 #  define PI 3.14159265358979323846
  49 # endif
  50 #endif
  51
  52 /* private data */
  53 static double c[8][8]; /* transform coefficients */
  54
  55 void fdct_enc(block)
  56 short *block;
  57 {
  58   int i, j, k;
  59   double s;
  60   double tmp[64];
  61
  62   for (i=0; i<8; i++)
  63     for (j=0; j<8; j++)
  64     {
  65       s = 0.0;
  66
  67       for (k=0; k<8; k++)
  68         s += c[j][k] * block[8*i+k];
  69
  70       tmp[8*i+j] = s;
  71     }
  72
  73   for (j=0; j<8; j++)
  74     for (i=0; i<8; i++)
  75     {
  76       s = 0.0;
  77
  78       for (k=0; k<8; k++)
  79         s += c[i][k] * tmp[8*k+j];
  80
  81       block[8*i+j] = (int)floor(s+0.499999);
  82       /*
  83        * reason for adding 0.499999 instead of 0.5:
  84        * s is quite often x.5 (at least for i and/or j = 0 or 4)
  85        * and setting the rounding threshold exactly to 0.5 leads to an
  86        * extremely high arithmetic implementation dependency of the result;
  87        * s being between x.5 and x.500001 (which is now incorrectly rounded
  88        * downwards instead of upwards) is assumed to occur less often
  89        * (if at all)
  90        */
  91     }
  92 }
  93
  94 void init_fdct_enc()
  95 {
  96   int i, j;
  97   double s;
  98
  99   for (i=0; i<8; i++)
 100   {
 101     s = (i==0) ? sqrt(0.125) : 0.5;
 102
 103     for (j=0; j<8; j++)
 104       c[i][j] = s * cos((PI/8.0)*i*(j+0.5));
 105   }
 106 }
 107
 108
 109
 110
 111 /* the second part of it is for the inverse DCT */
 112
 113 #define W1 2841 /* 2048*sqrt(2)*cos(1*pi/16) */
 114 #define W2 2676 /* 2048*sqrt(2)*cos(2*pi/16) */
 115 #define W3 2408 /* 2048*sqrt(2)*cos(3*pi/16) */
 116 #define W5 1609 /* 2048*sqrt(2)*cos(5*pi/16) */
 117 #define W6 1108 /* 2048*sqrt(2)*cos(6*pi/16) */
 118 #define W7 565  /* 2048*sqrt(2)*cos(7*pi/16) */
 119
 120 /* private data */
 121 static short iclip[1024]; /* clipping table */
 122 static short *iclp;
 123
 124 /* private prototypes */
 125 static void idctrow_enc (short *blk);
 126 static void idctcol_enc (short *blk);
 127
 128 /* two dimensional inverse discrete cosine transform */
 129 void idct_enc(block)
 130 short *block;
 131 {
 132   int i;
 133
 134   for (i=0; i<8; i++)
 135     idctrow_enc(block+8*i);
 136
 137   for (i=0; i<8; i++)
 138     idctcol_enc(block+i);
 139 }
 140
 141 void init_idct_enc()
 142 {
 143   int i;
 144
 145   iclp = iclip+512;
 146   for (i= -512; i<512; i++)
 147     iclp[i] = (i<-256) ? -256 : ((i>255) ? 255 : i);
 148 }
 149
 150 /* row (horizontal) IDCT
 151  *
 152  *           7                       pi         1
 153  * dst[k] = sum c[l] * src[l] * cos( -- * ( k + - ) * l )
 154  *          l=0                      8          2
 155  *
 156  * where: c[0]    = 128
 157  *        c[1..7] = 128*sqrt(2)
 158  */
 159
 160 static void idctrow_enc(blk)
 161 short *blk;
 162 {
 163   int x0, x1, x2, x3, x4, x5, x6, x7, x8;
 164
 165   /* shortcut */
 166   if (!((x1 = blk[4]<<11) | (x2 = blk[6]) | (x3 = blk[2]) |
 167         (x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3])))
 168   {
 169     blk[0]=blk[1]=blk[2]=blk[3]=blk[4]=blk[5]=blk[6]=blk[7]=blk[0]<<3;
 170     return;
 171   }
 172
 173   x0 = (blk[0]<<11) + 128; /* for proper rounding in the fourth stage */
 174
 175   /* first stage */
 176   x8 = W7*(x4+x5);
 177   x4 = x8 + (W1-W7)*x4;
 178   x5 = x8 - (W1+W7)*x5;
 179   x8 = W3*(x6+x7);
 180   x6 = x8 - (W3-W5)*x6;
 181   x7 = x8 - (W3+W5)*x7;
 182
 183   /* second stage */
 184   x8 = x0 + x1;
 185   x0 -= x1;
 186   x1 = W6*(x3+x2);
 187   x2 = x1 - (W2+W6)*x2;
 188   x3 = x1 + (W2-W6)*x3;
 189   x1 = x4 + x6;
 190   x4 -= x6;
 191   x6 = x5 + x7;
 192   x5 -= x7;
 193
 194   /* third stage */
 195   x7 = x8 + x3;
 196   x8 -= x3;
 197   x3 = x0 + x2;
 198   x0 -= x2;
 199   x2 = (181*(x4+x5)+128)>>8;
 200   x4 = (181*(x4-x5)+128)>>8;
 201
 202   /* fourth stage */
 203   blk[0] = (x7+x1)>>8;
 204   blk[1] = (x3+x2)>>8;
 205   blk[2] = (x0+x4)>>8;
 206   blk[3] = (x8+x6)>>8;
 207   blk[4] = (x8-x6)>>8;
 208   blk[5] = (x0-x4)>>8;
 209   blk[6] = (x3-x2)>>8;
 210   blk[7] = (x7-x1)>>8;
 211 }
 212
 213 /* column (vertical) IDCT
 214  *
 215  *             7                         pi         1
 216  * dst[8*k] = sum c[l] * src[8*l] * cos( -- * ( k + - ) * l )
 217  *            l=0                        8          2
 218  *
 219  * where: c[0]    = 1/1024
 220  *        c[1..7] = (1/1024)*sqrt(2)
 221  */
 222 static void idctcol_enc(blk)
 223 short *blk;
 224 {
 225   int x0, x1, x2, x3, x4, x5, x6, x7, x8;
 226
 227   /* shortcut */
 228   if (!((x1 = (blk[8*4]<<8)) | (x2 = blk[8*6]) | (x3 = blk[8*2]) |
 229         (x4 = blk[8*1]) | (x5 = blk[8*7]) | (x6 = blk[8*5]) | (x7 = blk[8*3])))
 230   {
 231     blk[8*0]=blk[8*1]=blk[8*2]=blk[8*3]=blk[8*4]=blk[8*5]=blk[8*6]=blk[8*7]=
 232       iclp[(blk[8*0]+32)>>6];
 233     return;
 234   }
 235
 236   x0 = (blk[8*0]<<8) + 8192;
 237
 238   /* first stage */
 239   x8 = W7*(x4+x5) + 4;
 240   x4 = (x8+(W1-W7)*x4)>>3;
 241   x5 = (x8-(W1+W7)*x5)>>3;
 242   x8 = W3*(x6+x7) + 4;
 243   x6 = (x8-(W3-W5)*x6)>>3;
 244   x7 = (x8-(W3+W5)*x7)>>3;
 245
 246   /* second stage */
 247   x8 = x0 + x1;
 248   x0 -= x1;
 249   x1 = W6*(x3+x2) + 4;
 250   x2 = (x1-(W2+W6)*x2)>>3;
 251   x3 = (x1+(W2-W6)*x3)>>3;
 252   x1 = x4 + x6;
 253   x4 -= x6;
 254   x6 = x5 + x7;
 255   x5 -= x7;
 256
 257   /* third stage */
 258   x7 = x8 + x3;
 259   x8 -= x3;
 260   x3 = x0 + x2;
 261   x0 -= x2;
 262   x2 = (181*(x4+x5)+128)>>8;
 263   x4 = (181*(x4-x5)+128)>>8;
 264
 265   /* fourth stage */
 266   blk[8*0] = iclp[(x7+x1)>>14];
 267   blk[8*1] = iclp[(x3+x2)>>14];
 268   blk[8*2] = iclp[(x0+x4)>>14];
 269   blk[8*3] = iclp[(x8+x6)>>14];
 270   blk[8*4] = iclp[(x8-x6)>>14];
 271   blk[8*5] = iclp[(x0-x4)>>14];
 272   blk[8*6] = iclp[(x3-x2)>>14];
 273   blk[8*7] = iclp[(x7-x1)>>14];
 274 }
 275
 276
 277