Discrete Cosine Transform DCT implementation C

Question

Im trying to implement a forward and inverse Discrete Cosine Transform (DCT) in C. The code is to transorm a single input block of pixels to the transformation matrix via the dct() function and then back to the original pixel values via the idct() function. Please see the attached code. My output form the idct are consecutive values of 244, 116, 244, 116 etc. From the look of the idct values, it doesnt look like my program is working.. Could someone help me out and give me an idea of what results i should be expecting after each function? Obviously after the idct, i should be getting prety close to the original input Matrix.

Thanks

 # include <stdio.h>
 # define PI 3.14

void dct(float [][]);       // Function prototypes
void idct(float [][]);     // Function prototypes

void dct(float inMatrix[8][8]){

    double dct,
    Cu,
    sum,
    Cv;

    int i,
    j,
    u,
    h = 0,
    v;

    FILE * fp = fopen("mydata.csv", "w");

    float dctMatrix[8][8],
    greyLevel;                       

    for (u = 0; u < 8; ++u) {
        for (v = 0; v < 8; ++v) {

            if (u == 0) {
                Cu = 1.0 / sqrt(2.0);
            } else {
                Cu = 1.0;
            }

            if (v == 0) {
                Cv = 1.0 / sqrt(2.0);
            } else {
                Cu = (1.0);
            }   

            sum = 0.0;  

            for (i = 0; i < 8; i++) {
                for (j = 0; j < 8; j++) {

                    // Level around 0
                    greyLevel = inMatrix[i][j];

                    dct = greyLevel * cos((2 * i + 1) * u * PI / 16.0) *
                        cos((2 * j + 1) * v * PI / 16.0);

                    sum += dct;

                }               
            }
            dctMatrix[u][v] = 0.25 * Cu * Cv * sum;
            fprintf(fp, "\n %f", dctMatrix[u][v]);          
        }
        fprintf(fp, "\n");
    }  
    idct(dctMatrix);  
 }

void idct(float dctMatrix[8][8]){

    double idct,
    Cu,
    sum,
    Cv;

    int i,
    j,
    u,
    v;

    float idctMatrix[8][8],
    greyLevel;

    FILE * fp = fopen("mydata.csv", "a");

    fprintf(fp, "\n Inverse DCT");                     

    for (i = 0; i < 8; ++i) {
        for (j = 0; j < 8; ++j) { 

            sum = 0.0;  

        for (u = 0; u < 8; u++) {
            for (v = 0; v < 8; v++) {

            if (u == 0) {
                Cu = 1.0 / sqrt(2.0);
            } else {
                Cu = 1.0;
              }

            if (v == 0) {
                Cv = 1.0 / sqrt(2.0);
            } else {
                Cu = (1.0);
              }   

                    // Level around 0
                greyLevel = dctMatrix[u][v];

                idct = (greyLevel * cos((2 * i + 1) * u * PI / 16.0) *
                        cos((2 * j + 1) * v * PI / 16.0));

                sum += idct;

                }               
            }
            idctMatrix[i][j] = 0.25 * Cu * Cv * sum;
            fprintf(fp, "\n %f", idctMatrix[i][j]);         
        }
        fprintf(fp, "\n");
    }    
 }


int main() {

   float    
    testBlockA[8][8] = { {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255} },

    testBlockB[8][8] = {{255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255} };

    dct(testBlockB);
}

Hello paul. According to SO guidelines, to insert code "indent by four spaces" or to have a "preformatted block within a list, indent by 8 spaces"..I followed these instructions for both ways yet the ouput was what you seen..Perhaps you would care to enlighten me? — user915071, Nov 29 '11 at 13:38
Just highlight the code and click on the code formatting button `{}` - this automatically indents by 4 spaces to get the required formatting. I did this for you initially but somehow you managed to undo it. — Paul R, Nov 29 '11 at 14:34
perhaps the developers would care to edit the instructions because that is not the same as the instructions given. i had to make an edit to the question a minute after i entered it, but in the meantime you had formatted the code. — user915071, Nov 29 '11 at 18:32

Bort · Accepted Answer · 2011-11-29T15:30:47.473

There are at least two typos in the Cv constant assignment at the if statements:

    if (v == 0) {
        Cv = 1.0 / sqrt(2.0);
    } else {
        Cu = (1.0); // << this should be Cv = 1.0
    }

Didn't check too properly though. Using the german wikipedia about cosine transform, following code works... I didn't want to spent time on figuring out how you defined what conversion constant. I guess you need to make sure that you use the correct constants and inverse functions:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

void dct(float **DCTMatrix, float **Matrix, int N, int M);
void write_mat(FILE *fp, float **testRes, int N, int M);
void idct(float **Matrix, float **DCTMatrix, int N, int M);
float **calloc_mat(int dimX, int dimY);
void free_mat(float **p);


float **calloc_mat(int dimX, int dimY){
    float **m = calloc(dimX, sizeof(float*));
    float *p = calloc(dimX*dimY, sizeof(float));
    int i;
    for(i=0; i <dimX;i++){
    m[i] = &p[i*dimY];

    }
   return m;
}

void free_mat(float **m){
  free(m[0]);
  free(m);
}

void write_mat(FILE *fp, float **m, int N, int M){

   int i, j;
   for(i =0; i< N; i++){
    fprintf(fp, "%f", m[i][0]);
    for(j = 1; j < M; j++){
       fprintf(fp, "\t%f", m[i][j]);
        }   
    fprintf(fp, "\n");
   }
   fprintf(fp, "\n");
}

void dct(float **DCTMatrix, float **Matrix, int N, int M){

    int i, j, u, v;
    for (u = 0; u < N; ++u) {
        for (v = 0; v < M; ++v) {
        DCTMatrix[u][v] = 0;
            for (i = 0; i < N; i++) {
                for (j = 0; j < M; j++) {
                    DCTMatrix[u][v] += Matrix[i][j] * cos(M_PI/((float)N)*(i+1./2.)*u)*cos(M_PI/((float)M)*(j+1./2.)*v);
                }               
            }
        }
    }  
 }

void idct(float **Matrix, float **DCTMatrix, int N, int M){
    int i, j, u, v;

    for (u = 0; u < N; ++u) {
        for (v = 0; v < M; ++v) {
          Matrix[u][v] = 1/4.*DCTMatrix[0][0];
          for(i = 1; i < N; i++){
          Matrix[u][v] += 1/2.*DCTMatrix[i][0];
           }
           for(j = 1; j < M; j++){
          Matrix[u][v] += 1/2.*DCTMatrix[0][j];
           }

           for (i = 1; i < N; i++) {
                for (j = 1; j < M; j++) {
                    Matrix[u][v] += DCTMatrix[i][j] * cos(M_PI/((float)N)*(u+1./2.)*i)*cos(M_PI/((float)M)*(v+1./2.)*j);
                }               
            }
        Matrix[u][v] *= 2./((float)N)*2./((float)M);
        }
    }  
 }



int main() {

   float    
    testBlockA[8][8] = { {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255},
                         {255, 255, 255, 255, 255, 255, 255, 255} },

    testBlockB[8][8] = {{255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255},
                        {255, 0, 255, 0, 255, 0, 255, 0},
                        {0, 255, 0, 255, 0, 255, 0, 255} };

    FILE * fp = fopen("mydata.csv", "w");
    int dimX = 8, dimY = 8;
    int i, j;

    float **testBlock = calloc_mat(dimX, dimY);
    float **testDCT = calloc_mat(dimX, dimY);
    float **testiDCT = calloc_mat(dimX, dimY);

    for(i = 0; i<dimX; i++){
      for(j = 0; j<dimY; j++){
        testBlock[i][j] = testBlockB[i][j];
      }
    }

    dct(testDCT, testBlock, dimX, dimY);
    write_mat(fp, testDCT, dimX, dimY);

    idct(testiDCT, testDCT, dimX, dimY);
    write_mat(fp, testiDCT, dimX, dimY);

    fclose(fp);
    free_mat(testBlock);
    free_mat(testDCT);
    free_mat(testiDCT);

    return 0;
}

Edit the dct is based on the crossproduct of formula DCT-II in the wiki. the idct is based on the crossproduct of formula DCT-III with the normalization factor 2/N per dimension (since this is the inverse to DCT-II as mentioned in the text). Edit I am pretty sure that the factor in the inverse dct should sqrt(2) and not 1/sqrt(2) in your version.

Hi Bort, I tried compiling that in DEV-cpp but i get the following errors: "In function `float** calloc_mat(int, int)': error: invalid conversion from `void*' to `float**' error: invalid conversion from `void*' to `float**'"..i.e its comlaining about lines 13 and 14 in above code. any ideas? Why are you using pointers to pointers by the way? — user915071, Nov 29 '11 at 20:58
I am using pointers, because they allow me to reuse the code for larger and smaller images than 8 by 8. I assume you can fix the error by a cast of the returned void pointer of calloc, although this shouldn't be necessary any more with modern compilers (I was using gcc 4.5). Add the cast as `float **m = (float**) calloc(dimX, sizeof(float*));` and `float *p = (float *) calloc(dimX*dimY, sizeof(float));` — Bort, Dec 01 '11 at 11:56
It works. Compile it and run it with these commands: gcc dct.c -o dct -lm; ./dct — Juan Carlos Kuri Pinto, Mar 29 '16 at 18:32

unwind · Answer 2 · 2011-11-29T12:59:15.630

2

You don't

#include <math.h>

That might mean that the compiler is assuming things about the math functions that are not true, for instance that they all returnint. Please note that all functions you call need to be declared somewhere, C doesn't have a "built-in" sin() anymore than it has a built-in printf() (for the latter you correctly include stdin.h, of course).

Also, you can use M_PI once you've included <math.h>.

edited Nov 29 '11 at 12:59

answered Nov 29 '11 at 12:47

unwind

391,730
64
469
606

Hi unwind, I have made that adjustment but i am getting the same output..According to GNU C : "Function: double cos (double x) This function returns the cosine of x, where x is given in radians. The return value is in the range -1 to 1." Should i cast all int values in the dct equation to double? Also, the function returns a vlue in radians, so i should be converting back to float?? – user915071 Nov 29 '11 at 13:50

score 1 · Answer 3 · answered Sep 01 '16 at 09:37

In addition to previous answer about typos in Cv constant (in both dct() and idct() functions), you used inverse DCT formula (2nd) incorrectly. You had to multiply by Cv and Cu every time in the loop. So, the correct code of idct() should be:

void idct(float dctMatrix[8][8]){

    double idct,
    Cu,
    sum,
    Cv;

    int i,
    j,
    u,
    v;

    float idctMatrix[8][8],
    greyLevel;

    FILE * fp = fopen("mydata.csv", "a");
    fprintf(fp, "\n Inverse DCT");      

    for (i = 0; i < 8; ++i) {
        for (j = 0; j < 8; ++j) { 

            sum = 0.0;  

        for (u = 0; u < 8; u++) {
            for (v = 0; v < 8; v++) {

            if (u == 0) {
                Cu = 1.0 / sqrt(2.0);
            } else {
                Cu = 1.0;
              }

            if (v == 0) {
                Cv = 1.0 / sqrt(2.0);
            } else {
                Cv = (1.0); //mistake was here - the same is in dct()
              }   
                greyLevel = dctMatrix[u][v];

                 // Multiply by Cv and Cu here!
                idct = (greyLevel * Cu * Cv *          
                        cos((2 * i + 1) * u * PI / 16.0) *
                        cos((2 * j + 1) * v * PI / 16.0));

                sum += idct;
                }               
            }
            // not "* Cv * Cu" here!
            idctMatrix[i][j] = 0.25 * sum;           
            fprintf(fp, "\n %f", idctMatrix[i][j]);  
        }
        fprintf(fp, "\n");      
    }    
 }

In this case output values are close to 255, 0, 255, 0, etc.

Discrete Cosine Transform DCT implementation C

3 Answers3

Linked