/*
* =====================================================================================
*
* Filename: dct.cpp
*
* Description:
*
* Version: 1.0
* Created: 28/04/2011 16:35:21
* Revision: none
* Compiler: gcc
*
* Author: John Ryland (jryland), jryland@xiaofrog.com
* Company: InvertedLogic
*
* =====================================================================================
*/
#include <stdio.h>
#include <math.h>
#include "dct.h"
// 98.07% quality
static const int quantizationMatrixA[64] =
{
1, 1, 1, 1, 1, 1, 1, 2,
1, 1, 1, 1, 1, 2, 2, 1,
1, 1, 1, 1, 1, 2, 2, 2,
1, 1, 1, 1, 1, 2, 2, 2,
1, 1, 1, 2, 2, 3, 3, 2,
1, 1, 1, 2, 2, 3, 3, 2,
1, 2, 2, 2, 3, 3, 3, 3,
2, 2, 3, 3, 3, 3, 3, 3
};
// 92.6% quality
static const int quantizationMatrixB[64] =
{
1, 1, 1, 2, 3, 6, 8, 10,
1, 1, 2, 3, 4, 8, 9, 8,
2, 2, 2, 3, 6, 8, 10, 8,
2, 2, 3, 4, 7, 12, 11, 9,
3, 3, 8, 11, 10, 16, 15, 11,
3, 5, 8, 10, 12, 15, 16, 13,
7, 10, 11, 12, 15, 17, 17, 14,
14, 13, 13, 15, 15, 14, 14, 14
};
// normal quality
static const int quantizationMatrixC[64] =
{
4, 1, 1, 2, 3, 5, 6, 7,
1, 1, 1, 2, 3, 7, 7, 7,
1, 1, 2, 3, 5, 7, 8, 7,
1, 2, 2, 3, 6, 11, 10, 7,
2, 2, 4, 7, 8, 14, 13, 9,
3, 4, 7, 8, 10, 13, 14, 11,
6, 8, 10, 11, 13, 15, 15, 13,
9, 11, 12, 12, 14, 12, 13, 12
};
// IJG reference matrix
static const int quantizationMatrixD[64] =
{
16, 11, 10, 16, 24, 40, 51, 61,
12, 12, 14, 19, 26, 58, 60, 55,
14, 13, 16, 24, 40, 57, 69, 56,
14, 17, 22, 29, 51, 87, 80, 62,
18, 22, 37, 56, 68, 109, 103, 77,
24, 35, 55, 64, 81, 104, 113, 92,
49, 64, 78, 87, 103, 121, 120, 101,
72, 92, 95, 98, 112, 100, 103, 99
};
const int *quantizationMatrix = &quantizationMatrixD[0];
static const int zigZagReorderTable[64] =
{
0, 1, 8, 16, 9, 2, 3, 10,
17, 24, 32, 25, 18, 11, 4, 5,
12, 19, 26, 33, 40, 48, 41, 34,
27, 20, 13, 6, 7, 14, 21, 28,
35, 42, 49, 56, 57, 50, 43, 36,
29, 22, 15, 23, 30, 37, 44, 51,
58, 59, 52, 45, 38, 31, 39, 46,
53, 60, 61, 54, 47, 55, 62, 63
};
// I think with a re-ordering of the coefficients
// it is possible to build in the zig-zag reordering without
// explicit extra steps.
// Also can simplify the decode steps with improved decode
// coefficients.
int inverseZigZagReorderTable[64];
//float DCTcoefficients[64*64];
long long intDCTcoefficients[64*64];
long long decodeDCTcoefficients[64*64];
//float DCToffsets[64];
long long intDCToffsets[64];
void initTables()
{
for (int v = 0; v < 8; v++) {
for (int u = 0; u < 8; u++) {
float offset = 0.0;
for (int y = 0; y < 8; y++) {
for (int x = 0; x < 8; x++) {
double u1 = ( u == 0 ) ? sqrt(1/8.0) : sqrt(2/8.0);
double v1 = ( v == 0 ) ? sqrt(1/8.0) : sqrt(2/8.0);
double coefficient = u1 * cos(M_PI/8.0 * (x + 0.5) * u) *
v1 * cos(M_PI/8.0 * (y + 0.5) * v) * 65536.0;
decodeDCTcoefficients[(y * 8 + x) * 64 + v * 8 + u] = round(coefficient);
coefficient /= quantizationMatrix[v*8+u];
//DCTcoefficients[(v * 8 + u) * 64 + y * 8 + x] = coefficient;// / quantizationMatrix[v*8+u];
intDCTcoefficients[(v * 8 + u) * 64 + y * 8 + x] = round(coefficient);
offset += -128 * coefficient;
//printf("%i,", (int)round(coefficient * 65536.0));
}
}
//DCToffsets[v*8+u] = offset;// / quantizationMatrix[v*8+u];
intDCToffsets[v*8+u] = round(offset) + 65536 / 2;// / quantizationMatrix[v*8+u];
//printf("%f,", offset);
}
}
for (int v = 0; v < 64; v++) {
inverseZigZagReorderTable[zigZagReorderTable[v]] = v;
}
}
void printBlock(int blockIn[64])
{
for (int y = 0; y < 8; y++) {
for (int x = 0; x < 8; x++) {
printf("%i,", blockIn[y*8+x]);
}
printf("\n");
}
}
void encodeDCT(int blockIn[64], int blockOut[64])
{
int coefficientIndex = 0;
for (int v = 0; v < 64; v++) {
int ival = 0;
for (int y = 0; y < 64; y++) {
ival += blockIn[y] * intDCTcoefficients[coefficientIndex++];
}
blockOut[inverseZigZagReorderTable[v]] = (ival + intDCToffsets[v]) >> 16;
}
}
void decodeDCT(int blockIn[64], int blockOut[64])
{
int dc[64];
for (int v = 0; v < 64; v++) {
dc[v] = blockIn[inverseZigZagReorderTable[v]] * quantizationMatrix[v];
}
int coefficientIndex = 0;
for (int y = 0; y < 64; y++) {
int ival = 0;
for (int v = 0; v < 64; v++) {
ival += dc[v] * decodeDCTcoefficients[coefficientIndex++];
}
blockOut[y] = ((ival + (65536/2)) >> 16) + 128;
}
}
void testDCT()
{
int exampleBlock[64] = {
52, 55, 61, 66, 70, 61, 64, 73,
63, 59, 55, 90,109, 85, 69, 72,
62, 59, 68,113,144,104, 66, 73,
63, 58, 71,122,154,106, 70, 69,
67, 61, 68,104,126, 88, 68, 70,
79, 65, 60, 70, 77, 68, 58, 75,
85, 71, 64, 59, 55, 61, 65, 83,
87, 79, 69, 68, 65, 76, 78, 94
};
// encode
printf("Encode:\n");
int exampleBlockDCT[64];
encodeDCT(exampleBlock, exampleBlockDCT);
printBlock(exampleBlockDCT);
// decode
printf("Decode:\n");
int decodeBlock[64];
decodeDCT(exampleBlockDCT, decodeBlock);
printBlock(decodeBlock);
// error
printf("Error:\n");
int terr = 0;
for (int y = 0; y < 8; y++) {
for (int x = 0; x < 8; x++) {
int error = exampleBlock[y*8+x] - decodeBlock[y*8+x];
printf("%i,", error);
terr += fabs(error);
}
printf("\n");
}
printf("error: %f\n", terr / 64.0);
}
