
%% this code is a simple version of the following publication.
%% Yuming Fang, Zhenzhong Chen, Weisi Lin, Chia-Wen Lin, 'Saliency Detection in the Compressed Domain for Adaptive 
%% Image Retargeting. IEEE Transactions on Image Processing 21(9): 3888-3901 (2012).
%% please cite the above publication when using the code.

function [new_i_smap, new_rg_smap, new_by_smap] = DC_original_smap(img)


%% resize the image into the new size with 16x*16y
r_img = img(:, :, 1);
g_img = img(:, :, 2);
b_img = img(:, :, 3);
[row, col] = size(r_img);
new_row = ceil(row/16) * 16;
new_col = ceil(col/16) * 16;
new_r_img = imresize(r_img, [new_row new_col], 'bilinear');
new_g_img = imresize(g_img, [new_row new_col], 'bilinear');
new_b_img = imresize(b_img, [new_row new_col], 'bilinear');
new_img(:, :, 1) = new_r_img;
new_img(:, :, 2) = new_g_img;
new_img(:, :, 3) = new_b_img;

ycc_img = rgb2ycbcr(new_img);

%% the image is transfered into YCbCr with 4:2:0
y_img = ycc_img(:, :, 1);
cb_img = imresize(ycc_img(:, :, 2), 0.5, 'bilinear');
cr_img = imresize(ycc_img(:, :, 3), 0.5, 'bilinear');


%% divide the image into 8*8 block
[y_row y_col] = size(y_img);
y_row_blk_num = y_row/8;
y_col_blk_num = y_col/8;
y_dct = zeros(y_row, y_col);

%% calculate the dct coefficients for Y channel.
%% obtain the YCrCb DCT DC coefficients which are used to calculate the saliency value for patches (8*8).
ycc_dc_coeff = zeros(y_row_blk_num, y_col_blk_num, 3);

y_dct_coeff = zeros(y_row_blk_num, y_col_blk_num, 8, 8);
for i = 1:y_row_blk_num
    for j = 1:y_col_blk_num
        y_dct((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8) = dct2(y_img((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8));
        ycc_dc_coeff(i, j, 1) = y_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        y_dct_coeff(i, j, :, :) = y_dct((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8);

        if ycc_dc_coeff(i, j, 1) < 100
            display('Y the value of dc is less than 100');
            display(i);
            display(j);
        end
    end
end

disp('1111111111111');


%% calculate the dct coefficients for Cr and Cb channels.
[c_row c_col] = size(cb_img);
c_row_blk_num = c_row/8;
c_col_blk_num = c_col/8;
cb_dct = zeros(c_row, c_col);
cr_dct = zeros(c_row, c_col);
%small_ycc_dc_coeff = zeros(c_row_blk_num);
for i = 1:c_row_blk_num
    for j = 1:c_col_blk_num
        cb_dct((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8) = dct2(cb_img((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8));
        cr_dct((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8) = dct2(cr_img((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j * 8));
        ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 1, 2) = cb_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 2, 2) = cb_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 2, (j - 1) * 2 + 1, 2) = cb_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 2, (j - 1) * 2 + 2, 2) = cb_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        
        ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 1, 3) = cr_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 2, 3) = cr_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 2, (j - 1) * 2 + 1, 3) = cr_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        ycc_dc_coeff((i - 1) * 2 + 2, (j - 1) * 2 + 2, 3) = cr_dct((i - 1) * 8 + 1, (j - 1) * 8 + 1);
        
        if ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 1, 2) < 100
            display('cb the value of dc is less than 100');
            display(i);
            display(j);
        end
        if  ycc_dc_coeff((i - 1) * 2 + 1, (j - 1) * 2 + 1, 3) < 100
            display('cr the value of dc is less than 100');
            display(i);
            display(j);
        end
                
    end
end


rg_channel = ycc_dc_coeff(:, :, 2);
by_channel = ycc_dc_coeff(:, :, 3);
i_channel = ycc_dc_coeff(:, :, 1);


%% calculate the human visual sensitivity
array_x = zeros(y_row_blk_num, y_col_blk_num);
array_y = zeros(y_row_blk_num, y_col_blk_num);
for i = 1:y_row_blk_num
    for j = 1:y_col_blk_num
        array_x(i, j) = i;
        array_y(i, j) = j;
    end
end

disp('444444444444444');

% compute the csfs and differences between patches
dist = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
csf = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
rg_diff = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
by_diff = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
i_diff = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
% y_ac_diff = zeros(y_row_blk_num, y_col_blk_num, y_row_blk_num, y_col_blk_num);
for i = 1 : y_row_blk_num
    for j = 1 : y_col_blk_num
        dist(i, j, :, :) = sqrt((i - array_x).^2 + (j - array_y).^2);
        csf(i, j, :, :) = (1/(20*(sqrt(2*pi))))*exp(-(dist(i, j, :, :).^2/(2*20.^2)));

%% the new test
        rg_diff(i, j, :, :) = abs(double(rg_channel(i, j)) - double(rg_channel))./ (abs(rg_channel(i, j)) + abs(rg_channel));
        by_diff(i, j, :, :) = abs(double(by_channel(i, j)) - double(by_channel))./ (abs(by_channel(i, j)) + abs(by_channel));
        i_diff(i, j, :, :) = abs(double(i_channel(i, j)) - double(i_channel)) ./ (abs(i_channel(i, j)) + abs(i_channel));


    end
end



min_csf = min(min(min(min(csf))));
max_csf = max(max(max(max(csf))));
csf = (csf - min_csf)/(max_csf-min_csf);

rg_smap = zeros(y_row_blk_num, y_col_blk_num);
by_smap = zeros(y_row_blk_num, y_col_blk_num);
i_smap = zeros(y_row_blk_num, y_col_blk_num);
% ac_smap = zeros(y_row_blk_num, y_col_blk_num);
for i = 1 : y_row_blk_num
    for j = 1 : y_col_blk_num
        rg_smap(i, j) = sum(sum(rg_diff(i, j, :, :) .* csf(i, j, :, :)));
        by_smap(i, j) = sum(sum(by_diff(i, j, :, :) .* csf(i, j, :, :)));
        i_smap(i, j) = sum(sum(i_diff(i, j, :, :) .* csf(i, j, :, :)));
%         ac_smap(i, j) = sum(sum(y_ac_diff(i, j, :, :) .* csf(i, j, :, :)));
    end
end


min_rg = min(min(rg_smap));
max_rg = max(max(rg_smap));
min_by = min(min(by_smap));
max_by = max(max(by_smap));
min_i = min(min(i_smap));
max_i = max(max(i_smap));


%% FOR NEW TEST
new_rg_smap = (rg_smap - min_rg)/(max_rg - min_rg);
new_by_smap = (by_smap - min_by)/(max_by - min_by);
new_i_smap = (i_smap - min_i)/(max_i - min_i);
% new_ac_smap = (ac_smap - min_ac)/(max_ac - min_ac);

new_rg_smap = imresize(new_rg_smap, [row, col]);
new_by_smap = imresize(new_by_smap, [row, col]);
new_i_smap = imresize(new_i_smap, [row, col]);
% new_ac_smap = imresize(new_ac_smap, [row, col]);

new_rg_smap = mat2gray(new_rg_smap);
new_by_smap = mat2gray(new_by_smap);
new_i_smap = mat2gray(new_i_smap);



 clear new_img; 
 clear ycc_dc_coeff;
%  clear y_ac_coeff;
 clear y_dct_coeff;



end






