
function [smap_result] = com_smap(img)
 

    
    
for img_num = 1:1


% img = imread([testdatadir subdir imlist(img_num).name]);

% img = imread(img_file);

%% 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], 'bicubic');
new_g_img = imresize(g_img, [new_row new_col], 'bicubic');
new_b_img = imresize(b_img, [new_row new_col], 'bicubic');
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, 'bicubic');
cr_img = imresize(ycc_img(:, :, 3), 0.5, 'bicubic');


%% 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_ac_coeff = zeros(y_row_blk_num, y_col_blk_num, 2, 7);
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_ac_coeff(i, j, 1, :) = y_dct((i - 1) * 8 + 1, (j - 1) * 8 + 2:(j - 1) * 8 + 8);
%         y_ac_coeff(i, j, 2, :) = y_dct((i - 1) * 8 + 2:(i - 1) * 8 + 8,
%         dc_temp = y_dct((i - 1) * 8 + 1 : i * 8, (j - 1) * 8 + 1 : j *
%         8);
        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


first_min = min(min(ycc_dc_coeff(:, :, 1)));
first_max = max(max(ycc_dc_coeff(:, :, 1)));

second_min = min(min(ycc_dc_coeff(:, :, 2)));
second_max = max(max(ycc_dc_coeff(:, :, 2)));

thrid_min = min(min(ycc_dc_coeff(:, :, 3)));
thrid_max = max(max(ycc_dc_coeff(:, :, 3)));


y_ac_coeff = zeros(y_row_blk_num, y_col_blk_num, 3);
for i = 1:y_row_blk_num
    for j = 1:y_col_blk_num
        y_ac_coeff(i, j, 1) = y_dct_coeff(i, j, 1, 2) + y_dct_coeff(i, j, 1, 3) + y_dct_coeff(i, j, 2, 1) + y_dct_coeff(i, j, 2, 2) + y_dct_coeff(i, j, 3, 1);
        y_ac_coeff(i, j, 2) = sum(y_dct_coeff(i, j, 1, 4:7)) + sum(y_dct_coeff(i, j, 4:7, 1)) + y_dct_coeff(i, j, 2, 3) + y_dct_coeff(i, j, 3, 2) + y_dct_coeff(i, j, 3, 3) + y_dct_coeff(i, j, 4, 4);
        y_ac_coeff(i, j, 3) = y_dct_coeff(i, j, 1, 8) + sum(y_dct_coeff(i, j, 2, 4:8)) + sum(y_dct_coeff(i, j, 3, 4:8)) + sum(y_dct_coeff(i, j, 4, 2:8)) - y_dct_coeff(i, j, 4, 4) + sum(y_dct_coeff(i, j, 5, 2:8)) + sum(y_dct_coeff(i, j, 6, 2:8)) + sum(y_dct_coeff(i, j, 7, 2:8)) + sum(y_dct_coeff(i, j, 8, 1:8));
    end
end



fourth_min = min(min(min(y_ac_coeff)));
fourth_max = max(max(max(y_ac_coeff)));
 
%imshow(ycc_dc_coeff(:, :, 3), []);

%% calculate the RGB DCT DC coefficients, this matrix is used to calculate the RG, BY, I channels for saliency map
rgb_dc_coeff = ycbcr2rgb(uint8(ycc_dc_coeff/8));


%% calculate the RG, BY and I channels.
r_dc_coeff = rgb_dc_coeff(:, :, 1);
g_dc_coeff = rgb_dc_coeff(:, :, 2);
b_dc_coeff = rgb_dc_coeff(:, :, 3);

r_channel = double(r_dc_coeff) - double((g_dc_coeff + b_dc_coeff) / 2);
g_channel = double(g_dc_coeff) - double((r_dc_coeff + b_dc_coeff) / 2);
b_channel = double(b_dc_coeff) - double((r_dc_coeff + g_dc_coeff) / 2);
y_channel = double((r_dc_coeff + g_dc_coeff) / 2) - double(abs(r_dc_coeff - g_dc_coeff) / 2) - double(b_dc_coeff); % new yellow color channel

rg_channel = (double(g_channel) - double(r_channel));
by_channel = (double(b_channel) - double(y_channel));
i_channel = (r_dc_coeff + g_dc_coeff + b_dc_coeff) / 3;


%% rg is the Cb, by is the Cr.
% rg_channel = y_ac_coeff(:, :, 2);
% by_channel = y_ac_coeff(:, :, 3);
% i_channel = y_ac_coeff(:, :, 1);
% i_channel = mat2gray(i_channel);
% rg_channel = mat2gray(rg_channel);
% by_channel = mat2gray(by_channel);


%% 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('444444444444443');


% 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);
%         dist(i, j, :, :) = ((i - array_x).^2 + (j - array_y).^2);
        csf(i, j, :, :) = 64./(exp(256 * 0.106 * pi * (atan(8 * dist(i, j, :, :) / 1536) + 2.3) / (2.3 * 180)));
%         csf(i, j, :, :) = 64./(exp(256 * 0.106 * pi * (atan(8 * dist(i, j, :, :) / 900)) / (2.3 * 180)));
%         csf(i, j, :, :) = 1./(double(dist(i, j, :, :)) + 1);
%         rg_diff(i, j, :, :) = double(rg_channel(i, j)) - double(rg_channel);
%         by_diff(i, j, :, :) = double(by_channel(i, j)) - double(by_channel);
%         i_diff(i, j, :, :) = double(i_channel(i, j)) - double(i_channel);

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

        %% just for new test
        for m = i:y_row_blk_num
            for n = j:y_col_blk_num
%                 if(dist_param == 0)
%                     y_ac_diff(i, j, m, n) = HausdorffDist(y_ac_coeff(i, j, :), y_ac_coeff(m, n, :));
%                 else
% the difference between patches is the norm L2
                  y_ac_diff(i, j, m, n) = norm(reshape(y_ac_coeff(i, j, :) - y_ac_coeff(m, n, :), 3, 1), 2);
%                 end
            end
        end
    end
end

% disp('5555555555555');


%% just for new test
for i = 1:y_row_blk_num
    for j = 1:y_col_blk_num
        for m = 1:i
            for n = 1:j
                y_ac_diff(i, j, m, n) = y_ac_diff(m, n, i, j);
            end
        end        
    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

% new_rg_smap = mat2gray(rg_smap);
% new_by_smap = mat2gray(by_smap);
% new_i_smap = mat2gray(i_smap);
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));
min_ac = min(min(ac_smap));
max_ac = max(max(ac_smap));
% new_rg_smap = 1 - (rg_smap - min_rg)/(max_rg - min_rg);
% new_by_smap = 1 - (by_smap - min_by)/(max_by - min_by);
% new_i_smap = 1 - (i_smap - min_i)/(max_i - min_i);
% new_ac_smap = (ac_smap - min_ac)/(max_ac - min_ac);

%% 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);



% imwrite(new_rg_smap, [smapdir 'Cb_smap/' subdir imlist(img_num).name]);
% imwrite(new_by_smap, [smapdir 'Cr_smap/' subdir imlist(img_num).name]);
% imwrite(new_i_smap, [smapdir 'Y_smap/' subdir imlist(img_num).name]);
% imwrite(new_ac_smap, [smapdir 'ac_smap/' subdir imlist(img_num).name]);




%% FOR NEW TEST
% smap = (new_rg_smap + new_by_smap + new_i_smap + new_ac_smap)/4;
smap = (new_rg_smap + new_by_smap + new_i_smap + new_ac_smap)/4;
smap = mat2gray(smap);
% smap = 1 - smap;
% smap = mat2gray(smap);
% figure(5);
% imshow(smap);

% csf = 64/(exp(256 * 0.106 * pi * (atan(8 * d / 1536) + 2.3) / (2.3 * 180)));


final_smap = imresize(smap, [row, col]);



final_smap = mat2gray(final_smap);

[salient_row salient_col salient_value] = find(final_smap > 0.9);



%% csf for enhance
salient_distance = zeros(row, col);
salient_csf = zeros(row, col);
for i = 1:row
    for j = 1:col
        
%         distance(i, j, :, :) = sqrt((i - array_x).^2 + (j - array_y).^2);
%         csf(i, j, :, :) = 64./(exp(256 * 0.106 * pi * (atan(8 * distance(i, j, :, :) / 1536) + 2.3) / (2.3 * 180)));       
%         
        salient_distance(i, j) = min(sqrt((i - salient_row).^2 + (j - salient_col).^2));
        
    end
end
salient_distance = mat2gray(salient_distance);

salient_csf(:, :) = 64./(exp(256 * 0.106 * pi * (atan(8 * salient_distance./ 1536) + 2.3) / (2.3 * 180)));

salient_csf = mat2gray(salient_csf);
final_smap = final_smap .*  salient_csf;

%  final_final_smap = imfilter(final_final_smap, fspecial('disk', 3));
final_smap = mat2gray(final_smap);




smap_result = final_smap;


% imwrite(final_smap, [smapdir 'final_smap_csf_enhance/' subdir imlist(img_num).name]);
% imwrite(smap, sprintf('abs_new_out/%d.jpg', img_num));

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


% calculate three saliency maps for RG, BY and I channels respectively.

end
end




