www.pudn.com > dm_dct_fenkuai.rar > dm_dct_fenkuai.m, change:2007-09-17,size:5566b


function [psnr,err]=dm_dct_fenkuai(step,attack_style,attack_strength) 
%%% attack_style 用于设置选择何种攻击 
%%% step 用于设置量化步长 
%%% err 返回误码率 
%%% psnr 返回峰值信噪比 
%%% attack_strength 用于控制攻击的强弱 
%%%%%%%%%%%%%%%% 读入原始载体图像 %%%%%%%%%%%%%%%%%%%%%%%%%%% 
[x,map]=imread('lena.bmp'); 
[row,col]=size(x); 
M=row; 
N=col; 
MN=col*row; 
figure(1),imshow(uint8(x));%%; image(x), colormap(map); 
title('原始图像','Fontsize',16,'color','blue'); 
x_source=x; 
fun1=@dct2; 
J1=blkproc(x,[8 8],fun1); %%%分割图像进行二维DCT变换 
%%%%%%%%%%%%%%%% 读入原始水印 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[m,mapm]=imread('laugh-single.bmp'); 
figure(2),imshow(m); 
title('嵌入的水印图像','Fontsize',8,'color','blue'); 
[rowm,colm]=size(m); 
MNm=rowm*colm; 
Mm=rowm; 
Nm=colm; 
mm=m(1:MNm); %%只有两个值,0和255 
Rm=MNm/MN; 
%%%%%%%%%%%%%%% 抖动量化实现水印信息的嵌入 %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% 计算(i,k)处的频率,用来判断(i,k)处是否进行抖动量化 
for i=1:M 
    for k=1:N 
        if mod(i,8)==0||mod(k,8)==0 
            fre(i,k)=1; 
        else 
            fre(i,k)=(((mod(i,8)-1)/16)^2+((mod(k,8)-1)/16)^2)^(1/2); 
        end 
    end 
end 
%% 计算(i,k)属于哪个块,用来判断(i,k)如何进行抖动量化 
for i=1:M 
    for k=1:N 
        block(i,k)=(M/8)*(ceil(i/8)-1)+ceil(k/8); 
    end 
end 
b1=block(M,N); %总块数 
L=fix(b1/MNm); %多少块嵌入1bit,为了解码方便,确保L为奇数 
%进行抖动量化 
%不进行纠错编码的情况下,产生抖动序列 
z=mm; 
num=1; 
num1=1; 
num2=1; 
for i=1:fix(MN*Rm) 
    d0(i)=-step/4; 
    d1(i)=d0(i)+step/2; 
end 
for i=1:M 
    for k=1:N 
        if fre(i,k)>0.25 
            J11(i,k)=J1(i,k); 
        else 
            if block(i,k)>L*MNm 
                J11(i,k)=J1(i,k); 
            else 
                if double(z(ceil(block(i,k)/L)))-255==0 %如果灰度为255 
                    J11(i,k)=round((J1(i,k)+d1(ceil(block(i,k)/L)))/step)*step-d1(ceil(block(i,k)/L)); 
                else 
                    J11(i,k)=round((J1(i,k)+d0(ceil(block(i,k)/L)))/step)*step-d0(ceil(block(i,k)/L)); 
                end 
            end 
        end 
    end 
end 
fun2=@idct2; 
J2=blkproc(J11,[8 8],fun2); %%%二维DCT逆变换 
figure(3);imshow(uint8(round(J2))); 
imwrite(uint8(round(J2)),'dm_watermarked.bmp'); 
%%%%%%%%%%%%%%%% 计算峰值信噪比PSNR %%%%%%%%%%%%%%%%%%%%%%%% 
x_temp1=J2-double(x_source); 
figure(4),imshow(uint8(round(100*x_temp1))); 
imwrite((uint8(round(100*x_temp1))),'dm_chazhi.bmp'); 
x_temp2=x_temp1( : ); 
x_temp3=abs(x_temp2); 
x_temp4=x_temp3'*x_temp3; 
d_embed=x_temp4/(M*M); 
SDR1=255*255/d_embed; 
psnr=10*log10(SDR1); 
%%%%%%%%%%%%%%%% 攻击 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
J2=imread('dm_watermarked.bmp'); 
if attack_style==1 
    %%放大两倍操作(提取前先缩小两倍) 
    xxx1=imresize(J2,2,'bicubic'); 
    xxx2=imresize(xxx1,1/2,'bicubic'); 
    yy=double(xxx2); 
end 
if attack_style==2 
    %%放大4倍操作(提取前先缩小4倍) 
    xxx1=imresize(J2,4,'bicubic'); 
    xxx2=imresize(xxx1,1/4,'bicubic'); 
    yy=double(xxx2); 
end 
if attack_style==3 
    %%缩小1/4操作 
    xxx1=imresize(J2,3/4,'bicubic'); 
    xxx2=imresize(xxx1,4/3,'bicubic'); 
    yy=double(xxx2); 
end 
if attack_style==4 
    %% 
    xxx1=imresize(J2,2/4,'bicubic'); 
    xxx2=imresize(xxx1,4/2,'bicubic'); 
    yy=double(xxx2); 
end 
%%3*3空域低通滤波 
if attack_style==5 
    B=(1/9)*ones(3,3); 
    xxx2=filter2(B,J2); 
    yy=double(xxx2); 
end 
%%4领域平均 
if attack_style==6 
    B=[0 1 0;1 0 1;0 1 0]*(1/4); 
    xxx2=filter2(B,J2); 
    yy=double(xxx2); 
end 
%%8领域平均 
if attack_style==7 
    B=[1 1 1;1 0 1;1 1 1]*(1/8); 
    xxx2=filter2(B,J2); 
    yy=double(xxx2); 
end 
%%窗口中值滤波 
if attack_style==8 
    xxx2=medfilt2(J2); %默认3*3 
    yy=double(xxx2); 
end 
if attack_style==9 
    % a1=input('Please input length of window a1:'); 
    % b1=input('Please input length of window b1:'); 
    a1=1; 
    b1=3; 
    xxx2=medfilt2(J2,[a1 b1]); 
    save al al; 
    save b1 b2; 
    yy=double(xxx2); 
end 
%% 裁减 
if attack_style==10 
    for i=128-44:128+45 
        for j=128-45:128+44 
            J2(i,j)=0; 
        end 
    end 
    yy=double(J2); 
end 
if attack_style==11 
    for i=128-64:128+63 
        for j=128-64:128+63 
            J2(i,j)=0; 
        end 
    end 
    yy=double(J2); 
end 
if attack_style==12 
    yy=imnoise(uint8(round(J2)),'gaussian',0,attack_strength); %高斯噪声 
end 
if attack_style==13 
    imwrite(uint8(round(J2)),'jpeg_n.jpg','jpg','Quality',attack_strength); 
    [yy,map]=imread('jpeg_n.jpg','jpg'); 
end 
if attack_style>13 
    yy=J2; 
end 
J3=blkproc(yy,[8 8],fun1); 
%%%%%%%%%%%%%%%%%%%%% 水印提取 %%%%%%%%%%%%%%%%% 
t1=zeros(L*MNm,1); 
t0=zeros(L*MNm,1); 
for i=1:M 
    for k=1:N 
        if fre(i,k)<=0.25&&block(i,k)<=L*MNm 
            sy0(i,k)=round((J3(i,k)+d0(ceil(block(i,k)/L)))/step)*step-d0(ceil(block(i,k)/L)); 
            sy1(i,k)=round((J3(i,k)+d1(ceil(block(i,k)/L)))/step)*step-d1(ceil(block(i,k)/L)); 
            j=block(i,k); 
            t1(j)=t1(j)+(J3(i,k)-sy1(i,k))^2; 
            t0(j)=t0(j)+(J3(i,k)-sy0(i,k))^2; 
        end 
    end 
end 
%每块分别解码 
for j=1:L*MNm 
    if t0(j)-t1(j)>0 
        t(j)=255; 
    else 
        t(j)=0; 
    end 
end 
%%每5块联合判决(观察即可) 
for j=1:MNm 
    tt(j)=sum(t((j-1)*L+1:j*L)); 
end 
for j=1:MNm 
    if tt(j)<(L/2)*255 
        m2(j)=0; 
    else 
        m2(j)=255; 
    end 
end 
%%%%%%%%%%%%%% 计算误码率 %%%%%%%%%%%%%%%%%%%] 
err=sum(xor(m2,mm)); 
err=err/1024; 
mdd=reshape(m2,rowm,colm); 
figure(5),imshow(uint8(round(mdd))); 
title('提取的水印图像','Fontsize',8,'color','blue');