Commit initial reference implementation of IDCT.

diff --git a/Makefile b/Makefile
index a30f39098a04a9bbaaee13a6b376d68351fe0ad5..926fa57317ff1cb35e75d80356f0f95244be5387 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 CC=gcc
 CFLAGS=-std=gnu99 -O2 -msse4.1 -g -Wall -Wextra
-LDFLAGS=
+LDFLAGS=-lm
 
 all: tests
 
@@ -20,15 +20,21 @@ BYTESOURCE_TEST_OBJS=bytesource.o choice.o unstuff.o bytesource_test.o
 bytesource_test: $(BYTESOURCE_TEST_OBJS)
        $(CC) $(LDFLAGS) -o $@ $(BYTESOURCE_TEST_OBJS)
 
-tests: unstuff_test bitsource_test dehuff_test bytesource_test
+IDCT_TEST_OBJS=idct.o idct_test.o
+idct_test: $(IDCT_TEST_OBJS)
+       $(CC) $(LDFLAGS) -o $@ $(IDCT_TEST_OBJS)
+
+tests: unstuff_test bitsource_test dehuff_test bytesource_test idct_test
 clean:
        $(RM) $(UNSTUFF_TEST_OBJS) unstuff_test
        $(RM) $(BITSOURCE_TEST_OBJS) bitsource_test
        $(RM) $(DEHUFF_TEST_OBJS) dehuff_test
        $(RM) $(BYTESOURCE_TEST_OBJS) bytesource_test
+       $(RM) $(IDCT_TEST_OBJS) idct_test
 
 test: tests
        ./unstuff_test
        ./bitsource_test
        ./dehuff_test
        ./bytesource_test
+       ./idct_test

diff --git a/idct.c b/idct.c
new file mode 100644 (file)
index 0000000..54a9dc9
--- /dev/null
+++ b/idct.c
@@ -0,0 +1,39 @@
+#include <math.h>
+
+#include "idct.h"
+
+void idct_reference(const int16_t* input, const uint32_t* quant_table, uint8_t* output)
+{
+       double temp[DCTSIZE2];
+
+       for (unsigned y = 0; y < 8; ++y) {
+               for (unsigned x = 0; x < 8; ++x) {
+                       double acc = 0.0;
+                       for (unsigned u = 0; u < 8; ++u) {
+                               for (unsigned v = 0; v < 8; ++v) {
+                                       double c_u = (u == 0) ? 1/sqrt(2.0) : 1.0;
+                                       double c_v = (v == 0) ? 1/sqrt(2.0) : 1.0;
+                                       acc += c_u * c_v
+                                               * input[u * DCTSIZE + v] * quant_table[u * DCTSIZE + v]
+                                               * cos((2 * x + 1) * v * M_PI / 16.0)
+                                               * cos((2 * y + 1) * u * M_PI / 16.0);
+                               }
+                       }
+                       temp[y * DCTSIZE + x] = 0.25 * acc;
+               }
+       }
+
+       for (unsigned y = 0; y < 8; ++y) {
+               for (unsigned x = 0; x < 8; ++x) {
+                       double val = temp[y * DCTSIZE + x];
+                       if (val < 0.0) {
+                               output[y * DCTSIZE + x] = 0;
+                       } else if (val >= 255.0) {
+                               output[y * DCTSIZE + x] = 255;
+                       } else {
+                               output[y * DCTSIZE + x] = (uint8_t)(val + 0.5);
+                       }
+               }
+       }
+}
+

diff --git a/idct.h b/idct.h
new file mode 100644 (file)
index 0000000..ec41993
--- /dev/null
+++ b/idct.h
@@ -0,0 +1,14 @@
+#ifndef _IDCT_H
+#define _IDCT_H
+
+#include <stdint.h>
+
+#define DCTSIZE 8
+#define DCTSIZE2 64
+
+typedef void (idct_func_t)(const int16_t*, const uint32_t*, uint8_t*);
+
+// Non-factorized reference version (section A.3.3 of the JPEG standard).
+void idct_reference(const int16_t* input, const uint32_t* quant_table, uint8_t* output);
+
+#endif /* !defined(_IDCT_H) */

diff --git a/idct_test.c b/idct_test.c
new file mode 100644 (file)
index 0000000..7dfac43
--- /dev/null
+++ b/idct_test.c
@@ -0,0 +1,130 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+#include <time.h>
+#include <sys/time.h>
+
+#include "idct.h"
+
+// Generate random coefficients in the range [-15..15].
+void gen_random_coeffs(int16_t* dst, size_t len)
+{
+       // Standard NR LCG (we avoid rand() to get consistent behavior across platforms).
+       static uint32_t seed = 1234;
+       for (unsigned i = 0; i < len; ++i) {
+               seed = seed * 1664525U + 1013904223U;
+               if (seed >> 31) {
+                       dst[i] = (uint8_t)(seed >> 27) & 0x7;
+               } else {
+                       dst[i] = -((uint8_t)(seed >> 27) & 0x7);
+               }
+       }
+}
+
+// Test that the input is pretty close to the reference for random inputs. 
+// (If the reference funtion is given in, this becomes a simple test of its
+// determinism.)
+void test_random_inputs(idct_func_t* idct)
+{
+       int16_t coeff[DCTSIZE2]; 
+       uint32_t quant[DCTSIZE2];
+       uint8_t output[DCTSIZE2];
+       uint8_t reference[DCTSIZE2];
+               
+       // Unit quantization (ie., no scaling).
+       for (unsigned i = 0; i < DCTSIZE2; ++i) {
+               quant[i] = 1;
+       }
+
+       for (unsigned i = 0; i < 1000; ++i) {   
+               gen_random_coeffs(coeff, DCTSIZE2);
+
+               (*idct)(coeff, quant, output);
+               (idct_reference)(coeff, quant, reference);
+
+               // Find the RMS difference.
+               int diff_squared = 0;
+               for (unsigned i = 0; i < DCTSIZE2; ++i) {
+                       diff_squared += (output[i] - reference[i]) * (output[i] - reference[i]);
+               }
+
+               assert(diff_squared <= 5);
+       }
+}
+
+// Test that a single DC coefficient becomes spread out to all blocks.
+void test_dc_becomes_spread_out(idct_func_t* idct)
+{
+       int16_t coeff[DCTSIZE2] = { 0 }; 
+       uint32_t quant[DCTSIZE2];
+       uint8_t output[DCTSIZE2];
+
+       // Unit quantization (ie., no scaling).
+       for (unsigned i = 0; i < DCTSIZE2; ++i) {
+               quant[i] = 1;
+       }
+
+       for (unsigned i = 0; i < 255*8; ++i) {  
+               uint32_t reference_value = i / 8;
+               coeff[0] = i;
+
+               (*idct)(coeff, quant, output);
+
+               for (unsigned i = 0; i < DCTSIZE2; ++i) {
+                       assert(abs(output[i] - reference_value) <= 1);
+               }
+       }
+}
+
+double timediff(const struct timeval* a, const struct timeval* b)
+{
+       return (double)(b->tv_sec - a->tv_sec) +
+               (double)(b->tv_usec - a->tv_usec) * 1e-6;
+}
+
+void test_performance(idct_func_t* idct)
+{
+       const unsigned num_runs = (idct == idct_reference) ? 5000 : 5000000;
+
+       int16_t coeff[DCTSIZE2]; 
+       uint32_t quant[DCTSIZE2];
+       uint8_t output[DCTSIZE2];
+               
+       gen_random_coeffs(coeff, DCTSIZE2);
+
+       struct timeval start, now;
+       gettimeofday(&start, NULL);
+
+       for (unsigned i = 0; i < num_runs; ++i) {
+               (*idct)(coeff, quant, output);
+       }
+       
+       gettimeofday(&now, NULL);
+
+       double diff = timediff(&start, &now);
+       printf("%u runs in %.2f seconds = %.2f DCTs/sec\n",
+               num_runs, diff, num_runs / diff);
+}
+
+void test_all_idct(idct_func_t* idct)
+{
+       printf("  test_dc_becomes_spread_out()\n");
+       test_dc_becomes_spread_out(idct);       
+
+       printf("  test_random_inputs()\n");
+       test_random_inputs(idct);       
+
+       printf("  performance test: ");
+       test_performance(idct);
+}
+
+int main(void)
+{
+       printf("idct_reference:\n");
+       test_all_idct(idct_reference);
+
+       printf("All tests pass.\n");
+       return 0;
+}