%***********************************************
% filename:visible_wmk.m
% function:在DCT閾中嵌入可視水印
% date:2006.03.29
% programmer:wt
%************************************************
function visible_wmk()
clear,clc
% 原始圖像
I=imread('E:\work\picture\wmk_test.jpg');
% 水印圖像
W=imread('E:\work\picture\wmk.bmp');
if isrgb(I)
    I=rgb2gray(I);
end
if isrgb(W)
    W=rgb2gray(W);
end
I=double(I);
W=double(W);
[m n]=size(I);
W=imresize(W,[m n]);
D=zeros(m,n);
% 原圖像和水印圖像都被分成大小為8*8的塊
block=8;
blockno_m=m/block;
blockno_n=n/block;
% 計(jì)算原圖像中每個(gè)塊的alpha_n和beta_n
% 利用Sobel邊緣算子,確定邊界塊,matlab里edge函數(shù)默認(rèn)使用sobel檢測(cè)器
[alpha,beta]=compute_dct(I,block,blockno_m,blockno_n);   %alpha應(yīng)該在0.95-0.98之間,beta應(yīng)該在0.07-0.17之間
alpha_max=max(max(alpha)),
alpha_min=min(min(alpha)),
beta_max=max(max(beta)),
beta_min=min(min(beta)),
[BW,t]=edge(I);
%t=3;
for i=1:blockno_m
    for j=1:blockno_n
        x=(i-1)*block+1;      y=(j-1)*block+1;
        block_8_8=I(x:x+block-1,y:y+block-1);               % 截取圖像塊 8x8 
        wmk_block=W(x:x+block-1,y:y+block-1);
        block_dct=dct2(block_8_8);                          % 對(duì)圖像進(jìn)行二維 DCT 變換
        wmk_block_dct=dct2(wmk_block);
        BW_8_8=BW(x:x+block-1,y:y+block-1);
        T=sum(sum(BW_8_8));                               % 對(duì)每個(gè)塊的邊緣信息求和
        if T>t                                           % 如果是邊緣塊,alpha取最大值,beta取最小值
            alpha(i,j)=alpha_max;
            beta(i,j)=beta_min;
        end
        % 利用公式1,修改原圖像第n個(gè)塊的DCT系數(shù)
        %block_dct=alpha(i,j)*block_dct+(beta(i,j)+0.07)*wmk_block_dct;  % beta加了一個(gè)系數(shù)0.07
        block_dct=alpha(i,j)*block_dct+beta(i,j)*wmk_block_dct;  
        
        % 逆DCT變換,得到加水印圖像
        block_dct=idct2(block_dct);                       
        D(x:x+block-1,y:y+block-1)=block_dct;              % 恢復(fù)圖像的像素


    end
end
subplot(2,2,1),imshow(uint8(I)),title('原始圖像');
subplot(2,2,2),imshow(uint8(W)),title('水印圖像');
subplot(2,2,3),imshow(uint8(D)),title('加入水印后的圖像');



return
%*************************************************
% sub1 計(jì)算原圖像中每個(gè)塊的alpha和beta
%*************************************************

function [alpha,beta]=compute_dct(I,block,blockno_m,blockno_n)
for i=1:blockno_m
    for j=1:blockno_n
        x=(i-1)*block+1;      y=(j-1)*block+1;
        block_8_8=I(x:x+block-1,y:y+block-1);               % 截取圖像塊 8x8 
        block_dct=dct2(block_8_8);                          % 對(duì)圖像進(jìn)行二維 DCT 變換
        c=block_dct;
        c_00(i,j)=c(1);                                     % 取圖像塊的直流分量組成矩陣c_00
        u_AC=(sum(sum(c))-c(1))/(block*block);          %AC系數(shù)的均值
        delta(i,j)=sum(sum((c-u_AC).*(c-u_AC)))/(block*block);   % 利用公式6計(jì)算delta
     end
end
c_00_max=max(max(c_00));                                % 直流分量里的最大值
u_n=c_00/c_00_max;                                      % 利用公式4,對(duì)原圖像的每個(gè)塊,計(jì)算規(guī)范化平均灰度值u_n,然后量化到0.1-1.0之間
 
%u_n_max=max(max(u_n)),
% u_n_min=min(min(u_n)),
u=sum(sum(c_00))/(blockno_m*blockno_n)/c_00_max;                            % 利用公式5計(jì)算圖像的平均灰度值u
%delta,
delta_log=log(delta);                                % 取delta的自然對(duì)數(shù)
delta_log_max=max(max(delta_log));
delta_n=delta_log/delta_log_max;                     % 利用公式7和ACDCT系數(shù),計(jì)算規(guī)范化方差delta_n,量化到0.1-1.0之間
% delta_n_max=max(max(delta_n)),
% delta_n_min=min(min(delta_n)),
% 利用公式2和3,計(jì)算alpha和beta
tt=exp(-(u_n-u).*(u_n-u));
alpha=delta_n.*tt;
beta=(1-tt)./delta_n;


return