lib/h.263/dct.c

   1 /* dct.c
   2
   3    DCT implementations and test routines.
   4
   5    Copyright (c) 2003 Matthias Kramm <kramm@quiss.org>
   6
   7    This program is free software; you can redistribute it and/or modify
   8    it under the terms of the GNU General Public License as published by
   9    the Free Software Foundation; either version 2 of the License, or
  10    (at your option) any later version.
  11
  12    This program is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15    GNU General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with this program; if not, write to the Free Software
  19    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */
  20
  21 #include <memory.h>
  22
  23 int zigzagtable[64] = {
  24     0, 1, 5, 6, 14, 15, 27, 28,
  25     2, 4, 7, 13, 16, 26, 29, 42,
  26     3, 8, 12, 17, 25, 30, 41, 43,
  27     9, 11, 18, 24, 31, 40, 44, 53,
  28     10, 19, 23, 32, 39, 45, 52, 54,
  29     20, 22, 33, 38, 46, 51, 55, 60,
  30     21, 34, 37, 47, 50, 56, 59, 61,
  31     35, 36, 48, 49, 57, 58, 62, 63};
  32
  33 static double table[8][8] =
  34 {
  35 {0.707106781186548,0.707106781186548,0.707106781186548,0.707106781186548,0.707106781186548,0.707106781186548,0.707106781186548,0.707106781186548},
  36 {0.980785280403230,0.831469612302545,0.555570233019602,0.195090322016128,-0.195090322016128,-0.555570233019602,-0.831469612302545,-0.980785280403230},
  37 {0.923879532511287,0.382683432365090,-0.382683432365090,-0.923879532511287,-0.923879532511287,-0.382683432365090,0.382683432365090,0.923879532511287},
  38 {0.831469612302545,-0.195090322016128,-0.980785280403230,-0.555570233019602,0.555570233019602,0.980785280403230,0.195090322016129,-0.831469612302545},
  39 {0.707106781186548,-0.707106781186547,-0.707106781186548,0.707106781186547,0.707106781186548,-0.707106781186547,-0.707106781186547,0.707106781186547},
  40 {0.555570233019602,-0.980785280403230,0.195090322016128,0.831469612302545,-0.831469612302545,-0.195090322016128,0.980785280403231,-0.555570233019602},
  41 {0.382683432365090,-0.923879532511287,0.923879532511287,-0.382683432365090,-0.382683432365091,0.923879532511287,-0.923879532511286,0.382683432365090},
  42 {0.195090322016128,-0.555570233019602,0.831469612302545,-0.980785280403231,0.980785280403230,-0.831469612302545,0.555570233019602,-0.195090322016129}
  43 };
  44
  45 void dct(int*src)
  46 {
  47     double tmp[64];
  48     int x,y,u,v,t;
  49
  50     for(v=0;v<8;v++)
  51     for(u=0;u<8;u++)
  52     {
  53         double c = 0;
  54         for(x=0;x<8;x++)
  55         {
  56             c+=table[u][x]*src[v*8+x];
  57         }
  58         tmp[v*8+u] = c;
  59     }
  60     for(u=0;u<8;u++)
  61     for(v=0;v<8;v++)
  62     {
  63         double c = 0;
  64         for(y=0;y<8;y++)
  65         {
  66             c+=table[v][y]*tmp[y*8+u];
  67         }
  68         src[v*8+u] = (int)(c*0.25+0.5);
  69     }
  70 }
  71
  72 void idct(int*src)
  73 {
  74     double tmp[64];
  75     int x,y,u,v;
  76     for(y=0;y<8;y++)
  77     for(x=0;x<8;x++)
  78     {
  79         double c = 0;
  80         for(u=0;u<8;u++)
  81         {
  82             c+=table[u][x]*src[y*8+u];
  83         }
  84         tmp[y*8+x] = c;
  85     }
  86     for(y=0;y<8;y++)
  87     for(x=0;x<8;x++)
  88     {
  89         double c = 0;
  90         for(v=0;v<8;v++)
  91         {
  92             c+=table[v][y]*tmp[v*8+x];
  93         }
  94         src[y*8+x] = (int)(c*0.25+0.5);
  95     }
  96 }
  97
  98 static double c[8] = {1.0,
  99 0.980785280403230, // cos(Pi*1/16), sin(Pi*7/16)
 100 0.923879532511287, // cos(Pi*2/16), sin(Pi*6/16)
 101 0.831469612302545, // cos(Pi*3/16), sin(Pi*5/16)
 102 0.707106781186548, // cos(Pi*4/16), sin(Pi*4/16), 1/sqrt(2)
 103 0.555570233019602, // cos(Pi*5/16), sin(Pi*3/16)
 104 0.382683432365090, // cos(Pi*6/16), sin(Pi*2/16)
 105 0.195090322016128 // cos(Pi*7/16), sin(Pi*1/16)
 106 };
 107
 108 static double cc[8];
 109 static int ccquant = -1;
 110
 111 void preparequant(int quant)
 112 {
 113     if(ccquant == quant)
 114         return;
 115     cc[0] = c[0]/(quant*2*4);
 116     cc[1] = c[1]/(quant*2*4);
 117     cc[2] = c[2]/(quant*2*4);
 118     cc[3] = c[3]/(quant*2*4);
 119     cc[4] = c[4]/(quant*2*4);
 120     cc[5] = c[5]/(quant*2*4);
 121     cc[6] = c[6]/(quant*2*4);
 122     cc[7] = c[7]/(quant*2*4);
 123     ccquant = quant;
 124 }
 125
 126 inline static void innerdct(const double*a,double*b, const double*c)
 127 {
 128     // c1*c7*2 = c6
 129     // c2*c6*2 = c4
 130     // c3*c5*2 = c2
 131     // c4*c4*2 = 1
 132
 133      //{  1,  3,  5,  7, -7, -5, -3, -1},
 134      //{  3, -7, -1, -5,  5,  1,  7, -3},
 135      //{  5, -1,  7,  3, -3, -7,  1, -5},
 136      //{  7, -5,  3, -1,  1, -3,  5, -7}
 137     double b0,b1,b2,b3,b4,b5;
 138     b2 = (a[0]+a[7]);
 139     b3 = (a[1]+a[6]);
 140     b4 = (a[2]+a[5]);
 141     b5 = (a[3]+a[4]);
 142
 143     b0 = (b2+b5)*c[4];
 144     b1 = (b3+b4)*c[4];
 145     b[0*8] = b0 + b1;
 146     b[4*8] = b0 - b1;
 147     b[2*8] = (b2-b5)*c[2] + (b3-b4)*c[6];
 148     b[6*8] = (b2-b5)*c[6] + (b4-b3)*c[2];
 149
 150     b0 = (a[0]-a[7]);
 151     b1 = (a[1]-a[6]);
 152     b2 = (a[2]-a[5]);
 153     b3 = (a[3]-a[4]);
 154
 155     b[1*8] = b0*c[1] + b1*c[3] + b2*c[5] + b3*c[7];
 156     b[3*8] = b0*c[3] - b1*c[7] - b2*c[1] - b3*c[5];
 157     b[5*8] = b0*c[5] - b1*c[1] + b2*c[7] + b3*c[3];
 158     b[7*8] = b0*c[7] - b1*c[5] + b2*c[3] - b3*c[1];
 159 }
 160
 161 void dct2(int*src, int*dest)
 162 {
 163     double tmp[64], tmp2[64];
 164     double*p;
 165     int u,x,v,t;
 166
 167     for(t=0;t<64;t++)
 168         tmp2[t] = src[t];
 169
 170     for(v=0;v<8;v++)
 171     {
 172         double* a=&tmp2[v*8];
 173         double* b=&tmp[v];
 174         innerdct(a,b,c);
 175     }
 176     for(v=0;v<8;v++)
 177     {
 178         double* a=&tmp[v*8];
 179         double* b=&tmp2[v];
 180         innerdct(a,b,cc);
 181     }
 182     for(t=0;t<64;t++) {
 183         int v = (int)(tmp2[t]);
 184         dest[zigzagtable[t]] = v;
 185     }
 186 }
 187
 188
 189 void zigzag(int*src)
 190 {
 191     int tmp[64];
 192     int t;
 193     for(t=0;t<64;t++) {
 194         tmp[zigzagtable[t]] = src[t];
 195     }
 196     memcpy(src, tmp, sizeof(int)*64);
 197 }
 198
 199