Main Page | Class List | File List | Class Members | File Members

vid_dct.c

Go to the documentation of this file.
00001 /************************************************************************
00002  *
00003  *  vid_dct.c, part of tmn (TMN encoder)
00004  *  Copyright (C) 1995, 1996  Telenor R&D, Norway
00005  *        Karl Olav Lillevold <Karl.Lillevold@nta.no>
00006  *  
00007  *  Contacts: 
00008  *  Karl Olav Lillevold               <Karl.Lillevold@nta.no>, or
00009  *  Robert Danielsen                  <Robert.Danielsen@nta.no>
00010  *
00011  *  Telenor Research and Development  http://www.nta.no/brukere/DVC/
00012  *  P.O.Box 83                        tel.:   +47 63 84 84 00
00013  *  N-2007 Kjeller, Norway            fax.:   +47 63 81 00 76
00014  *  
00015  ************************************************************************/
00016 
00017 /*
00018  * Disclaimer of Warranty
00019  *
00020  * These software programs are available to the user without any
00021  * license fee or royalty on an "as is" basis.  Telenor Research and
00022  * Development disclaims any and all warranties, whether express,
00023  * implied, or statuary, including any implied warranties or
00024  * merchantability or of fitness for a particular purpose.  In no
00025  * event shall the copyright-holder be liable for any incidental,
00026  * punitive, or consequential damages of any kind whatsoever arising
00027  * from the use of these programs.
00028  *
00029  * This disclaimer of warranty extends to the user of these programs
00030  * and user's customers, employees, agents, transferees, successors,
00031  * and assigns.
00032  *
00033  * Telenor Research and Development does not represent or warrant that
00034  * the programs furnished hereunder are free of infringement of any
00035  * third-party patents.
00036  *
00037  * Commercial implementations of H.263, including shareware, are
00038  * subject to royalty fees to patent holders.  Many of these patents
00039  * are general enough such that they are unavoidable regardless of
00040  * implementation design.
00041  * */
00042 
00043 
00044 /*****************************************************************
00045  * 
00046  * These routines are translated from Gisle Bjøntegaards's FORTRAN
00047  * routines by Robert.Danielsen@nta.no
00048  *
00049  *****************************************************************/
00050 
00051 #include "vid_macros.h"
00052 #include "vid_sim.h"
00053 
00054 #include <math.h>
00055 
00056 #ifndef PI
00057 # ifdef M_PI
00058 #  define PI M_PI
00059 # else
00060 #  define PI 3.14159265358979323846
00061 # endif
00062 #endif
00063 
00064 int     zigzag[8][8] = {
00065   {0, 1, 5, 6,14,15,27,28},
00066   {2, 4, 7,13,16,26,29,42},
00067   {3, 8,12,17,25,30,41,43},
00068   {9,11,18,24,31,40,44,53},
00069   {10,19,23,32,39,45,52,54},
00070   {20,22,33,38,46,51,55,60},
00071   {21,34,37,47,50,56,59,61},
00072   {35,36,48,49,57,58,62,63},
00073 };
00074 
00075 /**********************************************************************
00076  *
00077  *      Name:        Dct
00078  *      Description:    Does dct on an 8x8 block, does zigzag-scanning of
00079  *        coefficients
00080  *
00081  *      Input:        64 pixels in a 1D array
00082  *      Returns:        64 coefficients in a 1D array
00083  *      Side effects:   
00084  *
00085  *      Date: 930128    Author: Robert.Danielsen@nta.no
00086  *
00087  **********************************************************************/
00088 
00089 int Dct( int *block, int *coeff)
00090 {
00091   int        j1, i, j, k;
00092   float b[8];
00093   float        b1[8];
00094   float        d[8][8];
00095   float f0=(float).7071068;
00096   float f1=(float).4903926;
00097   float f2=(float).4619398;
00098   float f3=(float).4157348;
00099   float f4=(float).3535534;
00100   float f5=(float).2777851;
00101   float f6=(float).1913417;
00102   float f7=(float).0975452;
00103 
00104   for (i = 0, k = 0; i < 8; i++, k += 8) {
00105     for (j = 0; j < 8; j++) {
00106       b[j] = (float)block[k+j];
00107     }
00108     /* Horizontal transform */
00109     for (j = 0; j < 4; j++) {
00110       j1 = 7 - j;
00111       b1[j] = b[j] + b[j1];
00112       b1[j1] = b[j] - b[j1];
00113     }
00114     b[0] = b1[0] + b1[3];
00115     b[1] = b1[1] + b1[2];
00116     b[2] = b1[1] - b1[2];
00117     b[3] = b1[0] - b1[3];
00118     b[4] = b1[4];
00119     b[5] = (b1[6] - b1[5]) * f0;
00120     b[6] = (b1[6] + b1[5]) * f0;
00121     b[7] = b1[7];
00122     d[i][0] = (b[0] + b[1]) * f4;
00123     d[i][4] = (b[0] - b[1]) * f4;
00124     d[i][2] = b[2] * f6 + b[3] * f2;
00125     d[i][6] = b[3] * f6 - b[2] * f2;
00126     b1[4] = b[4] + b[5];
00127     b1[7] = b[7] + b[6];
00128     b1[5] = b[4] - b[5];
00129     b1[6] = b[7] - b[6];
00130     d[i][1] = b1[4] * f7 + b1[7] * f1;
00131     d[i][5] = b1[5] * f3 + b1[6] * f5;
00132     d[i][7] = b1[7] * f7 - b1[4] * f1;
00133     d[i][3] = b1[6] * f3 - b1[5] * f5;
00134   }
00135   /* Vertical transform */
00136   for (i = 0; i < 8; i++) {
00137     for (j = 0; j < 4; j++) {
00138       j1 = 7 - j;
00139       b1[j] = d[j][i] + d[j1][i];
00140       b1[j1] = d[j][i] - d[j1][i];
00141     }
00142     b[0] = b1[0] + b1[3];
00143     b[1] = b1[1] + b1[2];
00144     b[2] = b1[1] - b1[2];
00145     b[3] = b1[0] - b1[3];
00146     b[4] = b1[4];
00147     b[5] = (b1[6] - b1[5]) * f0;
00148     b[6] = (b1[6] + b1[5]) * f0;
00149     b[7] = b1[7];
00150     d[0][i] = (b[0] + b[1]) * f4;
00151     d[4][i] = (b[0] - b[1]) * f4;
00152     d[2][i] = b[2] * f6 + b[3] * f2;
00153     d[6][i] = b[3] * f6 - b[2] * f2;
00154     b1[4] = b[4] + b[5];
00155     b1[7] = b[7] + b[6];
00156     b1[5] = b[4] - b[5];
00157     b1[6] = b[7] - b[6];
00158     d[1][i] = b1[4] * f7 + b1[7] * f1;
00159     d[5][i] = b1[5] * f3 + b1[6] * f5;
00160     d[7][i] = b1[7] * f7 - b1[4] * f1;
00161     d[3][i] = b1[6] * f3 - b1[5] * f5;
00162   }
00163   /* Zigzag - scanning */
00164   for (i = 0; i < 8; i++) {
00165     for (j = 0; j < 8; j++) {
00166       *(coeff + zigzag[i][j]) = (int)(d[i][j]);
00167     }
00168   }
00169   return 0;
00170 }
00171 
00172 #ifdef FASTIDCT
00173 
00174 /**********************************************************************
00175  *
00176  *      Name:        idct
00177  *      Description:    Descans zigzag-scanned coefficients and does
00178  *        inverse dct on 64 coefficients
00179  *                      single precision floats
00180  *
00181  *      Input:        64 coefficients, block for 64 pixels
00182  *      Returns:        0
00183  *      Side effects:   
00184  *
00185  *      Date: 930128    Author: Robert.Danielsen@nta.no
00186  *
00187  **********************************************************************/
00188 
00189 int idct(int *coeff,int *block)
00190 {
00191   int                j1, i, j;
00192   double b[8], b1[8], d[8][8];
00193   double f0=.7071068;
00194   double f1=.4903926;
00195   double f2=.4619398;
00196   double f3=.4157348;
00197   double f4=.3535534;
00198   double f5=.2777851;
00199   double f6=.1913417;
00200   double f7=.0975452;
00201   double e, f, g, h;
00202 
00203   /* Horizontal */
00204 
00205   /* Descan coefficients first */
00206 
00207   for (i = 0; i < 8; i++) {
00208     for (j = 0; j < 8; j++) {
00209       b[j] = *( coeff + zigzag[i][j]);
00210     }
00211     e = b[1] * f7 - b[7] * f1;
00212     h = b[7] * f7 + b[1] * f1;
00213     f = b[5] * f3 - b[3] * f5;
00214     g = b[3] * f3 + b[5] * f5;
00215 
00216     b1[0] = (b[0] + b[4]) * f4;
00217     b1[1] = (b[0] - b[4]) * f4;
00218     b1[2] = b[2] * f6 - b[6] * f2;
00219     b1[3] = b[6] * f6 + b[2] * f2;
00220     b[4] = e + f;
00221     b1[5] = e - f;
00222     b1[6] = h - g;
00223     b[7] = h + g;
00224     
00225     b[5] = (b1[6] - b1[5]) * f0;
00226     b[6] = (b1[6] + b1[5]) * f0;
00227     b[0] = b1[0] + b1[3];
00228     b[1] = b1[1] + b1[2];
00229     b[2] = b1[1] - b1[2];
00230     b[3] = b1[0] - b1[3];
00231 
00232     for (j = 0; j < 4; j++) {
00233       j1 = 7 - j;
00234       d[i][j] = b[j] + b[j1];
00235       d[i][j1] = b[j] - b[j1];
00236     }
00237   }
00238 
00239   /* Vertical */
00240   
00241   for (i = 0; i < 8; i++) {
00242     for (j = 0; j < 8; j++) {
00243       b[j] = d[j][i];
00244     }
00245     e = b[1] * f7 - b[7] * f1;
00246     h = b[7] * f7 + b[1] * f1;
00247     f = b[5] * f3 - b[3] * f5;
00248     g = b[3] * f3 + b[5] * f5;
00249 
00250     b1[0] = (b[0] + b[4]) * f4;
00251     b1[1] = (b[0] - b[4]) * f4;
00252     b1[2] = b[2] * f6 - b[6] * f2;
00253     b1[3] = b[6] * f6 + b[2] * f2;
00254     b[4] = e + f;
00255     b1[5] = e - f;
00256     b1[6] = h - g;
00257     b[7] = h + g;
00258 
00259     b[5] = (b1[6] - b1[5]) * f0;
00260     b[6] = (b1[6] + b1[5]) * f0;
00261     b[0] = b1[0] + b1[3];
00262     b[1] = b1[1] + b1[2];
00263     b[2] = b1[1] - b1[2];
00264     b[3] = b1[0] - b1[3];
00265 
00266     for (j = 0; j < 4; j++) {
00267       j1 = 7 - j;
00268       d[j][i] = b[j] + b[j1];
00269       d[j1][i] = b[j] - b[j1];
00270     }
00271   }
00272 
00273   for (i = 0; i < 8; i++) {
00274     for (j = 0; j < 8; j++) {
00275       *(block + i * 8 + j) = mnint(d[i][j]);
00276     }
00277   }
00278   return 0;
00279 }
00280 
00281 #else
00282 /*  Perform IEEE 1180 reference (64-bit floating point, separable 8x1
00283  *  direct matrix multiply) Inverse Discrete Cosine Transform
00284 */
00285 
00286 
00287 /* Here we use math.h to generate constants.  Compiler results may
00288    vary a little */
00289 
00290 
00291 /* private data */
00292 
00293 /* cosine transform matrix for 8x1 IDCT */
00294 static double c[8][8];
00295 
00296 /* initialize DCT coefficient matrix */
00297 
00298 void init_idctref()
00299 {
00300   int freq, time;
00301   double scale;
00302 
00303   for (freq=0; freq < 8; freq++)
00304   {
00305     scale = (freq == 0) ? sqrt(0.125) : 0.5;
00306     for (time=0; time<8; time++)
00307       c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5));
00308   }
00309 }
00310 
00311 /* perform IDCT matrix multiply for 8x8 coefficient block */
00312 
00313 void idctref(int *coeff, int *block)
00314 {
00315   int i, j, k, v;
00316   double partial_product;
00317   double tmp[64];
00318   int tmp2[64];
00319   extern int zigzag[8][8];
00320 
00321   for (i=0; i<8; i++)
00322     for (j=0; j<8; j++)
00323       tmp2[j+i*8] = *(coeff + zigzag[i][j]);
00324 
00325   for (i=0; i<8; i++)
00326     for (j=0; j<8; j++)
00327     {
00328       partial_product = 0.0;
00329 
00330       for (k=0; k<8; k++)
00331         partial_product+= c[k][j]*tmp2[8*i+k];
00332 
00333       tmp[8*i+j] = partial_product;
00334     }
00335 
00336   /* Transpose operation is integrated into address mapping by switching 
00337      loop order of i and j */
00338 
00339   for (j=0; j<8; j++)
00340     for (i=0; i<8; i++)
00341     {
00342       partial_product = 0.0;
00343 
00344       for (k=0; k<8; k++)
00345         partial_product+= c[k][i]*tmp[8*k+j];
00346 
00347       v = (int)floor(partial_product+0.5);
00348       block[8*i+j] = (v<-256) ? -256 : ((v>255) ? 255 : v);
00349     }
00350 
00351 
00352 }
00353 #endif

Generated on Sun Jul 16 16:27:45 2006 by  doxygen 1.3.9.1