DCT 2D without FFT

Question

I'm trying to implement DCT (Discrete Cosine Transform) in Matlab but without using Fast Fourier Transform, just by using the next formula:

I know this can be inefficient but this way I will get how it works.

First I divided my grayscale image in 8x8 blocks, then I apply the formula to each block.

for i=1:8:h
    for j=1:8:w
        dctMatrix(i:(i-1)+block,j:(j-1)+block) = dctII(img(i:(i-1)+block,j:(j-1)+block), block);
    end
end

My dctII function looks like this:

   function [newB] = dctII(segmento, b)
    [h w] = size(segmento);
    segmento = double(segmento);
    newB = zeros(b,b);
    for u=0:h-1
        for v=0:w-1
            if u == 0
                Cu = 1/sqrt(2);
            else
                Cu = 1;
            end
            if v == 0
                Cv = 1/sqrt(2);
            else
                Cv = 1;
            end
            sumRes = summation(segmento,u,v,b);
            dct = (1/4)*Cu*Cv*sumRes;
            segmento(u+1,v+1) = dct;
        end
    end
    newB = segmento;
end

I also created a summation function to keep things more readable (at least for me).

function [sum] = summation(segmento,u,v,b)
    [h w] = size(segmento);
    sum = 0;
    for x=0:h-1
        for y=0:w-1
            sum = sum + (double(segmento(x+1,y+1))*cos((((2*x)+1)*u*pi)/(2*b))*cos((((2*y)+1)*v*pi)/(2*b)));
        end
    end
end

The problem is that the result of my algorithm differs by far the result of Matlab prebuilt function dct2. Maybe I didn't get DCT algorithm at all. Do you know what I'm doing wrong? I know all these nested loops kill performance serverely but I can't imagine how to solve this without using FFT.

Any help would be really appreciated, thanks.

Answer 1

Already solved!

My results differ from Matlab prebuilt function dct2 because that function doesn't consider an 8x8 block division like I do, so the matrices aren't equal.

I also made slight modifications to my code, I read that before working with pixels values you need to substract 128 to each value in order to work within the range [-128, 127]. This is my code:

function [newB] = dctII(segmento, b)
    [h w] = size(segmento);
    segmento = double(segmento) - 128;
    newB = zeros(b,b);
    for u=0:h-1
        for v=0:w-1
            if u == 0
                Cu = 1/sqrt(2);
            else
                Cu = 1;
            end
            if v == 0
                Cv = 1/sqrt(2);
            else
                Cv = 1;
            end
            sumRes = summation(segmento,u,v,b);
            dct = (1/4)*Cu*Cv*sumRes;
            newB(u+1,v+1) = dct;
        end
    end
end

function [sum] = summation(segmento,u,v,b)
    [h w] = size(segmento);
    sum = 0;
    for x=0:h-1
        for y=0:w-1
            sum = sum + (double(segmento(x+1,y+1))*cos(((2*x)+1)*u*pi/(2*b))*cos(((2*y)+1)*v*pi/(2*b)));
        end
    end
end

Is not very efficient but it is an implementation of the formula. Hope that helps