2D Discrete Cosine Transform using 1D DCT

Question

I am trying to implement the 2D Discrete Cosine Transform to an image by using 1D DCT operations. My output is incorrect if I compare it to the dct2 MATLAB function. I don't understand what went wrong in my code and where it's happening.

If someone can point out the mistake or any other advice, that would be really helpful.

Here is my code written in MATLAB

% main function
signal=rand(100);
signal_dct=myDCT(signal);
figure; imshow((signal_dct));

% function to calculate 2D DCT of an image
function res=myDCT(signal)

    signal=double(signal);

    l=size(signal,1);

    res=zeros(l);  %initialize the final result matrix

    for k=1:l     %calculate 1D DCT of each row of image
        res(k,:)=mdct(signal(k,:));  
    end

    for k=1:l   %calculate 1D DCT of each column of image
        res(:,k)=mdct(res(:,k));
    end
end

%% function to calculate 1D DFT of a 1D signal
function res=mdct(signal)

    l=size(signal,1);

    for i=1:l

        if i==1   %for signal index of 1, alpha is 1/sqrt(l)
            alpha=sqrt(1/l);
        else   %for signal index of greater than 1
            alpha=sqrt(2/l);
        end

        j=[1:l];
        % summation calculates single entry of res by applying the  
        % formula of DCT on the signal
        summation=sum(sum(signal(j)*cos((pi*(2*(j-1)+1)*(i-1))/(2*l))));
        res(i)=alpha*summation;        
    end
end

Answer 1

You are correct in that the 2D DCT is separable. You would just apply the 1D DCT to every row first, then take the intermediate result and apply it to the columns. However, you have two fundamental mistakes. Let's go through them.

Mistake #1 - Size of DCT isn't correct

Specifically, look at this statement in your mdct function:

l=size(signal,1);

Because you are applying the DCT for each row, then each column, the above would only work if you are applying the DCT to the columns. size(signal,1) would certainly give you the length of the input vector if the input was a column. However, if your input was a row, then the output of size(signal,1) would be 1. Therefore, you should replace size(signal,1) with numel so that you are for sure going to get the total number of elements - regardless if the input is a row or column.

Also, if you want to make the code compatible to do the summation in the DCT loop, you should make sure that the input is a row vector regardless. As such, do this instead:

l = numel(signal);
signal = signal(:).';

The first line determines how many elements we have for our input signal, and the second line ensures that we have a row vector. This is done by (:) to unroll the elements into a column vector, then doing .' after to ensure we transpose the result to get a row vector.

Mistake #2 - Summation statement isn't correct

Next, you're going to have to do element-wise multiplications in your summation to get what you're looking for. You also don't need the extra sum call there. It's superfluous. Therefore, modify your summation statement to this:

summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l)));

There's no need to do signal(j) because j spans the entire length of the vector, and you can just do that with signal.

---

Once I made these changes, and I did this on a smaller size matrix to ensure that we get the same results:

rng(123123);
signal=rand(7);
signal_dct=myDCT(signal);
signal_dct2 = dct2(signal);

The last line of code calls dct2 so that we can compare the results from your custom function and what dct2 gives us.

We get:

>> signal_dct

signal_dct =

    3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621
   -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581
   -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282
   -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427
   -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767
   -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789
   -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710

>> signal_dct2

signal_dct2 =

    3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621
   -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581
   -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282
   -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427
   -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767
   -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789
   -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710

---

As you can see, both results match up. Looks good to me!

---

Just to be sure we are consistent, this is the full code listing for both your functions, with the modifications I made:

% function to calculate 2D DCT of an image
function res=myDCT(signal)

    signal=double(signal);

    l=size(signal,1);

    res = zeros(l);

    for k=1:l     %calculate 1D DCT of each row of image
        res(k,:)=mdct(signal(k,:));  
    end

    for k=1:l   %calculate 1D DCT of each column of image
        res(:,k)=mdct(res(:,k));
    end
end

%% function to calculate 1D DFT of a 1D signal
function res=mdct(signal)

    %// Change
    l = numel(signal);
    signal = signal(:).';

    for i=1:l

        if i==1   %for signal index of 1, alpha is 1/sqrt(l)
            alpha=sqrt(1/l);
        else   %for signal index of greater than 1
            alpha=sqrt(2/l);
        end

        j=[1:l];
        % summation calculates single entry of res by applying the  
        % formula of DCT on the signal
        %// Change
        summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l)));
        res(i)=alpha*summation;        
    end
end