Trying to implement The inverse of a dct 8*8 matrix

Question

I have managed to calculate the dct of an 8*8 matrix and I am having trouble doing the inverse. Can anyone have a look at this code and tell me what I am doing now. I should be getting the exact same values as before but im getting different values. I am reading in the input from a csv file and out putting it to anther csv file. Its programmed in c

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

double idct,
Cu,
sum,
Cv;

int i,
j,
u,
v;


float idctMatrix[8][8],
greyLevel;

FILE * fp = fopen("mydata.csv", "r");
FILE * wp = fopen("idct.csv", "w");
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);
            // Level around 0
            greyLevel = idctMatrix[u][v];
            idct = (greyLevel * cos((2 * i + 1) * u * M_PI / 16.0) *
                    cos((2 * j + 1) * v * M_PI / 16.0));
            sum += idct;
            }               
        }

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

the original matrix is:

{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}};

the dct is :

2040   0  -0   0   0   0  -0  -0
   0   0   0   0  -0   0  -0   0
  -0   0  -0   0   0   0   0   0
   0  -0  -0  -0   0  -0  -0   0
   0   0  -0   0  -0  -0  -0   0
   0  -0  -0  -0  -0   0  -0  -0
  -0  -0  -0   0   0   0   0  -0
  -0   0   0   0  -0   0  -0   0

the idct calculated should be the same as the original

Please give an example of an input that it fails on, along with expected output, and the actual output. Also, what have you learnt by running this in the debugger? — Oliver Charlesworth, Dec 18 '11 at 18:18
Also, there are many definitions of the DCT (http://en.wikipedia.org/wiki/Discrete_cosine_transform#Formal_definition). Which one do you want? — Oliver Charlesworth, Dec 18 '11 at 18:22
You should get something deeper into the code. We won't go through it step-by-step. Are intermediate results correct? Are the output values approximately equal to the values you are expecting? — Sebastian, Dec 18 '11 at 18:30
ive updated the what the original matrix is and what the dct i have calculated looks like — user427641, Dec 18 '11 at 18:35
@user427641: You should really try debugging this yourself (e.g. stepping through it in the debugger, or printing intermediate values) to learn at what point its behaviour diverges from what you were expecting. Only then should you ask a Stack Overflow question. Otherwise, you're asking us to do your debugging for you. — Oliver Charlesworth, Dec 18 '11 at 18:37
Are you sure this is the code you're debugging? Apart from being incomplete, I'm surprised that fprinting to a file opened in "r" mode works. — , Dec 18 '11 at 18:38
Please actually read then link provided by Oli Charlesworth and check the implementation of your DCT (i assume you want DCT-II). Applying IDCT is straight forward, since you only have to calculate the DCT again and apply a scaling factor (which is listed on the particular Wikipedia page). Don't forget that your data must be within a given range, so you may need to normalize it. Your code is somewhat lengthy for my personal taste, computing the DCT is relatively short within two loops and you also don't need any if-else-constructs. — Sebastian, Dec 18 '11 at 18:50
didn't I answer that one already with working source code? [here](http://stackoverflow.com/questions/8310749/discrete-cosine-transform-dct-implementation-c/8311564#8311564), so why are you still stuck with it? Even the same wiki links were given.... — Bort, Dec 19 '11 at 15:09
Your question title talks about matrix inverse but the content appears to be about dct inverse. — David Heffernan, Dec 19 '11 at 19:50
i am having a lot of trouble pin pointing where I am going wrong, hoping someone else can spot something — user427641, Dec 19 '11 at 20:03
Also, as I can see, you have a function `idct_func` which is commented out, so you don't use it, but you have the implementation there. That just makes the code look big and unlikely for the people to read. — Shahbaz, Dec 19 '11 at 20:17
You should also try to show your efforts on the problem. For example, use a debugger and sau what happens etc. — Shahbaz, Dec 19 '11 at 20:18
@shahbaz yes I am calling my idct_func..its at the end of my dct_func. Please read it carefully and make an attempted at a constructive answer rather than moaning. I am just looking for someone who knows more than me to point out a simple mistake i am making — user427641, Dec 19 '11 at 20:24

Alexey Frunze · Accepted Answer · 2011-12-21T04:32:44.143

You're trying to calculate the IDCT in-place, using idctMatrix[][] as both input and output and hence the input is being modified in the process. That's wrong. You need to have a separate 2-dimensional array for input (or output).

It also appears that the u- and v-dependent scaling factors Cu and Cv are applied outside of the u- and v- loops. That would be wrong too.

EDIT: Here's the code if you can't understand the words:

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

#ifndef M_PI
#define M_PI 3.14159265358979324
#endif

void Compute8x8Dct(const double in[8][8], double out[8][8])
{
  int i, j, u, v;
  double s;

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

      for (u = 0; u < 8; u++)
        for (v = 0; v < 8; v++)
          s += in[u][v] * cos((2 * u + 1) * i * M_PI / 16) *
                          cos((2 * v + 1) * j * M_PI / 16) *
               ((i == 0) ? 1 / sqrt(2) : 1) *
               ((j == 0) ? 1 / sqrt(2) : 1);

      out[i][j] = s / 4;
    }
}

void Compute8x8Idct(const double in[8][8], double out[8][8])
{
  int i, j, u, v;
  double s;

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

      for (u = 0; u < 8; u++)
        for (v = 0; v < 8; v++)
          s += in[u][v] * cos((2 * i + 1) * u * M_PI / 16) *
                          cos((2 * j + 1) * v * M_PI / 16) *
               ((u == 0) ? 1 / sqrt(2) : 1.) *
               ((v == 0) ? 1 / sqrt(2) : 1.);

      out[i][j] = s / 4;
    }
}

void Print8x8(const char* title, const double in[8][8])
{
  int i, j;

  printf("%s\n", title);

  for (i = 0; i < 8; i++)
  {
    for (j = 0; j < 8; j++)
      printf("%8.3f ", in[i][j]);
    printf("\n");
  }
}

int main(void)
{
  double pic1[8][8], dct[8][8], pic2[8][8];
  int i, j;

  for (i = 0; i < 8; i++)
    for (j = 0; j < 8; j++)
#if 01
      pic1[i][j] = 255;
#else
      pic1[i][j] = (i ^ j) & 1;
#endif
  Print8x8("pic1:", pic1);
  Compute8x8Dct(pic1, dct);
  Print8x8("dct:", dct);
  Compute8x8Idct(dct, pic2);
  Print8x8("pic2:", pic2);

  return 0;
}

Here's its output:

pic1:
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
dct:
2040.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
   0.000    0.000    0.000    0.000    0.000    0.000    0.000    0.000 
pic2:
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000 
 255.000  255.000  255.000  255.000  255.000  255.000  255.000  255.000

apologies alex i missed one line in my code, it takes in the inMatrix[8][8] and the output matrix is the idctMatrix[8][8], thanks for the help anyway — user427641, Dec 18 '11 at 19:48
It appears that you've neglected my answer completely, at least I'm not seeing any changes in the code towards fixing the indicated problems. Good job, keep it up. — Alexey Frunze, Dec 20 '11 at 09:03
I have not neglected your answer, it doesnt seem to be correct. I thought it was — user427641, Dec 20 '11 at 12:10
See the code in the updated answer. The code reflects the words of the initial answer exactly. — Alexey Frunze, Dec 21 '11 at 04:34
much appreciated alex thanks, i understand ur your code much more — user427641, Dec 21 '11 at 10:48

Trying to implement The inverse of a dct 8*8 matrix

1 Answers1

Linked