DCT 2D without FFT

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.

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