visiblewmk.m

在图像中加入可视水印

%***********************************************
% 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;
% 计算原图像中每个块的alpha_n和beta_n
% 利用Sobel边缘算子,确定边界块,matlab里edge函数默认使用sobel检测器
[alpha,beta]=compute_dct(I,block,blockno_m,blockno_n);   %alpha应该在0.95-0.98之间,beta应该在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);                          % 对图像进行二维 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));                               % 对每个块的边缘信息求和
        if T>t                                           % 如果是边缘块,alpha取最大值,beta取最小值
            alpha(i,j)=alpha_max;
            beta(i,j)=beta_min;
        end
        % 利用公式1,修改原图像第n个块的DCT系数
        %block_dct=alpha(i,j)*block_dct+(beta(i,j)+0.07)*wmk_block_dct;  % beta加了一个系数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;              % 恢复图像的像素


    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 计算原图像中每个块的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);                          % 对图像进行二维 DCT 变换
        c=block_dct;
        c_00(i,j)=c(1);                                     % 取图像块的直流分量组成矩阵c_00
        u_AC=(sum(sum(c))-c(1))/(block*block);          %AC系数的均值
        delta(i,j)=sum(sum((c-u_AC).*(c-u_AC)))/(block*block);   % 利用公式6计算delta
     end
end
c_00_max=max(max(c_00));                                % 直流分量里的最大值
u_n=c_00/c_00_max;                                      % 利用公式4,对原图像的每个块,计算规范化平均灰度值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计算图像的平均灰度值u
%delta,
delta_log=log(delta);                                % 取delta的自然对数
delta_log_max=max(max(delta_log));
delta_n=delta_log/delta_log_max;                     % 利用公式7和ACDCT系数,计算规范化方差delta_n,量化到0.1-1.0之间
% delta_n_max=max(max(delta_n)),
% delta_n_min=min(min(delta_n)),
% 利用公式2和3,计算alpha和beta
tt=exp(-(u_n-u).*(u_n-u));
alpha=delta_n.*tt;
beta=(1-tt)./delta_n;


return