function score = DCTex(im_ref,im_dis)
% (c) ZHANG Fan, CUHK, June 17, 2011.
% F. Zhang, L. Ma, S. Li, & K. N. Ngan, "Practical image quality metric 
% applied to image coding," IEEE TMM, 2011, 13(5).

blocksize = 8;

%% Pre-calculate DCT bases and CSF
if ~exist('DCT.mat','file')
    winsize = blocksize;
    a = ones(1,winsize)/sqrt(winsize);
    a(2:end) = a(2:end)*sqrt(2);
    p = (0:winsize-1)*pi*0.5/winsize;
    W = zeros(winsize^2,winsize^2);
    for u = 0:winsize-1
        for v = 0:winsize-1
            for x = 0:winsize-1
                for y = 0:winsize-1
                    W(v*winsize+u+1,y*winsize+x+1) = ...
                        a(u+1)*a(v+1)*cos((2*x+1)*p(u+1))*cos((2*y+1)*p(v+1));
                end
            end
        end
    end
    [X,Y] = meshgrid(0:7,0:7);
    T = sqrt(X.^2 + Y.^2);
    t = (T(:)+10).*exp(-T(:));      % Luminance channel a4 = 10, a5 = 1
    t1 = (T(:)+10).*exp(-1.5*T(:)); % Chromatic channel a4 = 10, a5 = 1.5
    save DCT.mat W t t1;
else
    load DCT.mat W t t1;
end

%% Color space transform 
if ndims(im_ref) == 3
    im_ref = rgb2ycbcr(im_ref);
end
if ndims(im_dis) == 3
    im_dis = rgb2ycbcr(im_dis);
end
[row,col,clr] = size(im_ref);

%% Crop image to a size of 16n x 16m
row = floor(row/blocksize/2)*blocksize*2;
col = floor(col/blocksize/2)*blocksize*2;
im_ref = double(im_ref(1:row,1:col,:));
im_dis = double(im_dis(1:row,1:col,:));

%% Image blocking, vectorizing and DCT
[X,Y] = meshgrid(1:blocksize:row,0:blocksize*row:(col-1)*row);
Bindex = X +Y;
Bindex = (Bindex(:))';
[X,Y] = meshgrid(0:7,0:row:7*row);  
Boffset = X + Y;
sI = repmat(Bindex,[64,1]) + repmat(Boffset(:),[1,length(Bindex)]);
[X,Y] = meshgrid(1:blocksize:row/2,0:blocksize*row/2:(col/2-1)*row/2);
Bindex = X + Y;
Bindex = (Bindex(:))';
[X,Y] = meshgrid(0:7,0:row/2:7*row/2);  
Boffset = X + Y;
sI1 = repmat(Bindex,[64,1]) + repmat(Boffset(:),[1,length(Bindex)]);

S = W*(im_ref(sI)-im_dis(sI));
smooth = std(im_ref(sI),1); 
rough = var(mean(im_ref(sI)))/var(im_ref(1:row*col));
im1 = imresize(im_ref(:,:,1),0.5);
smooth1 = std(im1(sI1),1);
im1 = imresize(im_ref(:,:,2),0.5);
im2 = imresize(im_dis(:,:,2),0.5);
S1(:,:,1) = W*(im1(sI1)-im2(sI1));
im1 = imresize(im_ref(:,:,3),0.5);
im2 = imresize(im_dis(:,:,3),0.5);
S1(:,:,2) = W*(im1(sI1)-im2(sI1));

score = sum(sum(S.^2.*repmat(t,[1,size(S,2)]))./(smooth+20))...
+.8*sum(sum(sum(S1(:,:,:).^2.*repmat(t1,[1,size(S1,2),2]))./(repmat(smooth1+20,[1 1 2]))));
score = log(score*rough/numel(im_ref)+1);
end


