/***************************************************************** * tmn (TMN encoder) * Copyright (C) 1995 Telenor R&D * Karl Olav Lillevold * * These routines are translated from Gisle Bjøntegaards's FORTRAN * routines by Robert.Danielsen@nta.no * *****************************************************************/ #include "macros.h" #include "sim.h" #include #ifndef PI # ifdef M_PI # define PI M_PI # else # define PI 3.14159265358979323846 # endif #endif int zigzag[8][8] = { {0, 1, 5, 6,14,15,27,28}, {2, 4, 7,13,16,26,29,42}, {3, 8,12,17,25,30,41,43}, {9,11,18,24,31,40,44,53}, {10,19,23,32,39,45,52,54}, {20,22,33,38,46,51,55,60}, {21,34,37,47,50,56,59,61}, {35,36,48,49,57,58,62,63}, }; /********************************************************************** * * Name: Dct * Description: Does dct on an 8x8 block, does zigzag-scanning of * coefficients * * Input: 64 pixels in a 1D array * Returns: 64 coefficients in a 1D array * Side effects: * * Date: 930128 Author: Robert.Danielsen@nta.no * **********************************************************************/ int Dct( int *block, int *coeff) { int j1, i, j, k; float b[8]; float b1[8]; float d[8][8]; float f0=.7071068,f1=.4903926,f2=.4619398,f3=.4157348,f4=.3535534; float f5=.2777851,f6=.1913417,f7=.0975452; for (i = 0, k = 0; i < 8; i++, k += 8) { for (j = 0; j < 8; j++) { b[j] = block[k+j]; } /* Horizontal transform */ for (j = 0; j < 4; j++) { j1 = 7 - j; b1[j] = b[j] + b[j1]; b1[j1] = b[j] - b[j1]; } b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; b[4] = b1[4]; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[7] = b1[7]; d[i][0] = (b[0] + b[1]) * f4; d[i][4] = (b[0] - b[1]) * f4; d[i][2] = b[2] * f6 + b[3] * f2; d[i][6] = b[3] * f6 - b[2] * f2; b1[4] = b[4] + b[5]; b1[7] = b[7] + b[6]; b1[5] = b[4] - b[5]; b1[6] = b[7] - b[6]; d[i][1] = b1[4] * f7 + b1[7] * f1; d[i][5] = b1[5] * f3 + b1[6] * f5; d[i][7] = b1[7] * f7 - b1[4] * f1; d[i][3] = b1[6] * f3 - b1[5] * f5; } /* Vertical transform */ for (i = 0; i < 8; i++) { for (j = 0; j < 4; j++) { j1 = 7 - j; b1[j] = d[j][i] + d[j1][i]; b1[j1] = d[j][i] - d[j1][i]; } b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; b[4] = b1[4]; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[7] = b1[7]; d[0][i] = (b[0] + b[1]) * f4; d[4][i] = (b[0] - b[1]) * f4; d[2][i] = b[2] * f6 + b[3] * f2; d[6][i] = b[3] * f6 - b[2] * f2; b1[4] = b[4] + b[5]; b1[7] = b[7] + b[6]; b1[5] = b[4] - b[5]; b1[6] = b[7] - b[6]; d[1][i] = b1[4] * f7 + b1[7] * f1; d[5][i] = b1[5] * f3 + b1[6] * f5; d[7][i] = b1[7] * f7 - b1[4] * f1; d[3][i] = b1[6] * f3 - b1[5] * f5; } /* Zigzag - scanning */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { *(coeff + zigzag[i][j]) = (int)(d[i][j]); } } return 0; } #ifdef FASTIDCT /********************************************************************** * * Name: idct * Description: Descans zigzag-scanned coefficients and does * inverse dct on 64 coefficients * single precision floats * * Input: 64 coefficients, block for 64 pixels * Returns: 0 * Side effects: * * Date: 930128 Author: Robert.Danielsen@nta.no * **********************************************************************/ int idct(int *coeff,int *block) { int j1, i, j; double b[8], b1[8], d[8][8]; double f0=.7071068, f1=.4903926, f2=.4619398, f3=.4157348; double f4=.3535534; double f5=.2777851, f6=.1913417, f7=.0975452; double e, f, g, h; /* Horizontal */ /* Descan coefficients first */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { b[j] = *( coeff + zigzag[i][j]); } e = b[1] * f7 - b[7] * f1; h = b[7] * f7 + b[1] * f1; f = b[5] * f3 - b[3] * f5; g = b[3] * f3 + b[5] * f5; b1[0] = (b[0] + b[4]) * f4; b1[1] = (b[0] - b[4]) * f4; b1[2] = b[2] * f6 - b[6] * f2; b1[3] = b[6] * f6 + b[2] * f2; b[4] = e + f; b1[5] = e - f; b1[6] = h - g; b[7] = h + g; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; for (j = 0; j < 4; j++) { j1 = 7 - j; d[i][j] = b[j] + b[j1]; d[i][j1] = b[j] - b[j1]; } } /* Vertical */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { b[j] = d[j][i]; } e = b[1] * f7 - b[7] * f1; h = b[7] * f7 + b[1] * f1; f = b[5] * f3 - b[3] * f5; g = b[3] * f3 + b[5] * f5; b1[0] = (b[0] + b[4]) * f4; b1[1] = (b[0] - b[4]) * f4; b1[2] = b[2] * f6 - b[6] * f2; b1[3] = b[6] * f6 + b[2] * f2; b[4] = e + f; b1[5] = e - f; b1[6] = h - g; b[7] = h + g; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; for (j = 0; j < 4; j++) { j1 = 7 - j; d[j][i] = b[j] + b[j1]; d[j1][i] = b[j] - b[j1]; } } for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { *(block + i * 8 + j) = mnint(d[i][j]); } } return 0; } #else /* Perform IEEE 1180 reference (64-bit floating point, separable 8x1 * direct matrix multiply) Inverse Discrete Cosine Transform */ /* Here we use math.h to generate constants. Compiler results may vary a little */ /* private data */ /* cosine transform matrix for 8x1 IDCT */ static double c[8][8]; /* initialize DCT coefficient matrix */ void init_idctref() { int freq, time; double scale; for (freq=0; freq < 8; freq++) { scale = (freq == 0) ? sqrt(0.125) : 0.5; for (time=0; time<8; time++) c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5)); } } /* perform IDCT matrix multiply for 8x8 coefficient block */ void idctref(int *coeff, int *block) { int i, j, k, v; double partial_product; double tmp[64]; int tmp2[64]; extern int zigzag[8][8]; for (i=0; i<8; i++) for (j=0; j<8; j++) tmp2[j+i*8] = *(coeff + zigzag[i][j]); for (i=0; i<8; i++) for (j=0; j<8; j++) { partial_product = 0.0; for (k=0; k<8; k++) partial_product+= c[k][j]*tmp2[8*i+k]; tmp[8*i+j] = partial_product; } /* Transpose operation is integrated into address mapping by switching loop order of i and j */ for (j=0; j<8; j++) for (i=0; i<8; i++) { partial_product = 0.0; for (k=0; k<8; k++) partial_product+= c[k][i]*tmp[8*k+j]; v = floor(partial_product+0.5); block[8*i+j] = (v<-256) ? -256 : ((v>255) ? 255 : v); } } #endif