/************************************************************************ * * dct.c, part of tmn (TMN encoder) * * Copyright (C) 1997 University of BC, Canada * * Contacts: * Michael Gallant * Guy Cote * Berna Erol * * UBC Image Processing Laboratory http://www.ee.ubc.ca/image * 2356 Main Mall tel.: +1 604 822 4051 * Vancouver BC Canada V6T1Z4 fax.: +1 604 822 5949 * * Copyright (C) 1995, 1996 Telenor R&D, Norway * * Contacts: * Robert Danielsen * * Telenor Research and Development http://www.nta.no/brukere/DVC/ * P.O.Box 83 tel.: +47 63 84 84 00 * N-2007 Kjeller, Norway fax.: +47 63 81 00 76 * ************************************************************************/ /* Disclaimer of Warranty * * These software programs are available to the user without any license fee * or royalty on an "as is" basis. The University of British Columbia * disclaims any and all warranties, whether express, implied, or * statuary, including any implied warranties or merchantability or of * fitness for a particular purpose. In no event shall the * copyright-holder be liable for any incidental, punitive, or * consequential damages of any kind whatsoever arising from the use of * these programs. * * This disclaimer of warranty extends to the user of these programs and * user's customers, employees, agents, transferees, successors, and * assigns. * * The University of British Columbia does not represent or warrant that the * programs furnished hereunder are free of infringement of any * third-party patents. * * Commercial implementations of H.263, including shareware, are subject to * royalty fees to patent holders. Many of these patents are general * enough such that they are unavoidable regardless of implementation * design. * */ /***************************************************************** * * These routines are translated from Gisle Bjøntegaards's FORTRAN * routines by Robert.Danielsen@nta.no * * 970715 Modified by Guy Cote to include new * scanning for (advanced intra coding mode of H.263+ * *****************************************************************/ #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}, }; int alternate_horizontal[8][8] = { {0, 1, 2, 3, 10, 11, 12, 13}, {4, 5, 8, 9, 17, 16, 15, 14}, {6, 7, 19, 18, 26, 27, 28, 29}, {20, 21, 24, 25, 30, 31, 32, 33}, {22, 23, 34, 35, 42, 43, 44, 45}, {36, 37, 40, 41, 46, 47, 48, 49}, {38, 39, 50, 51, 56, 57, 58, 59}, {52, 53, 54, 55, 60, 61, 62, 63}, }; int alternate_vertical[8][8] = { {0, 4, 6, 20, 22, 36, 38, 52}, {1, 5, 7, 21, 23, 37, 39, 53}, {2, 8, 19, 24, 34, 40, 50, 54}, {3, 9, 18, 25, 35, 41, 51, 55}, {10, 17, 26, 30, 42, 46, 56, 60}, {11, 16, 27, 31, 43, 47, 57, 61}, {12, 15, 28, 32, 44, 48, 58, 62}, {13, 14, 29, 33, 45, 49, 59, 63}, }; /********************************************************************** * * Name: Dct * Description: Does dct on an 8x8 block * * 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 = (float) .7071068; float f1 = (float) .4903926; float f2 = (float) .4619398; float f3 = (float) .4157348; float f4 = (float) .3535534; float f5 = (float) .2777851; float f6 = (float) .1913417; float f7 = (float) .0975452; for (i = 0, k = 0; i < 8; i++, k += 8) { for (j = 0; j < 8; j++) { b[j] = (float) 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; } for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { *(coeff + j + i * 8) = (int) (d[i][j]); } } return 0; } /********************************************************************** * * Name: Scan and Descan * Description: Does scanning of the an 8x8 block quantized dct * coefficients. * (zigzag, or alernate horizontal, alternate * vertical when advanced intra coding is used) * * Input: 64 quantized coefficients in a 1D array * Returns: 0 * Side effects: scans the dct qcoeff in the desired order * (zigzag if advanced intra coding is not used) * * Date: 970716 * Author: Guy Cote * **********************************************************************/ void Scan (int *qcoeff, int Intra_Mode) { int i, j, k; int store_qcoeff[384]; /* Zigzag - scanning of a Macroblock */ /* Choose Scanning if Advanced Intra Coding Mode is used */ for (i = 0; i < 384; i++) { store_qcoeff[i] = *(qcoeff + i); } for (k = 0; k < 6; k++) { /* Do for all blocks */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { if (advanced_intra_coding) { switch (Intra_Mode) { case INTRA_MODE_DC: *(qcoeff + zigzag[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]); break; case INTRA_MODE_VERT_AC: *(qcoeff + alternate_horizontal[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]); break; case INTRA_MODE_HORI_AC: *(qcoeff + alternate_vertical[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]); break; default: printf ("Invalid Intra_Mode in Advanced Intra Coding"); exit (-1); break; } } else { *(qcoeff + zigzag[i][j] + k * 64) = (int) (store_qcoeff[i * 8 + j + k * 64]); } } } } return; } void DeScan (int *qcoeff, int Intra_Mode) { /* Descan coefficients first */ /* Choose Scanning if Advanced Intra Coding Mode is used */ int store_qcoeff[384]; int i, j, k; for (i = 0; i < 384; i++) { store_qcoeff[i] = *(qcoeff + i); } for (k = 0; k < 6; k++) { /* Do for all blocks */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { if (advanced_intra_coding) { /* Descan with proper scanning matrix */ switch (Intra_Mode) { case INTRA_MODE_DC: *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + zigzag[i][j]); break; case INTRA_MODE_VERT_AC: *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + alternate_horizontal[i][j]); break; case INTRA_MODE_HORI_AC: *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + alternate_vertical[i][j]); break; default: printf ("Invalid Intra_Mode in Advanced Intra Coding\n"); exit (-1); break; } } else { *(qcoeff + k * 64 + i * 8 + j) = *(store_qcoeff + k * 64 + zigzag[i][j]); } } } } return; } #ifdef FASTIDCT /********************************************************************** * * Name: idct * Description: 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; double f1 = .4903926; double f2 = .4619398; double f3 = .4157348; double f4 = .3535534; double f5 = .2777851; double f6 = .1913417; double f7 = .0975452; double e, f, g, h; /* Horizontal */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) b[j] = coeff[j + i * 8]; 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]; for (i = 0; i < 8; i++) for (j = 0; j < 8; j++) tmp2[j + i * 8] = coeff[j + i * 8]; 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 = (int) floor (partial_product + 0.5); block[8 * i + j] = (v < -256) ? -256 : ((v > 255) ? 255 : v); } } #endif