一、实验目的

1、实现卡方分析算法；

2、利用卡方分析能够区分经LSB隐写的载密图像和未经隐写算法修改的原始载体图像。

二、问题描述

1.输入的形式和输入值的范围；

84张bmp格式的标准图像；

2.输出的形式；

task1：载秘图像和原始图像的卡方分析对比图以及分割线；

task2：使用task1中的分割值T对图像类型进行测试，得到预测准确度；

3.程序所能达到的功能。

通过卡方分布检验或者特征值检验可以区分载密图像和原始图像。

task1：通过84张标准库图像嵌入秘密信息得到84张载秘图像，对共168张图像进行两种不同的卡方分析，得到两种不同卡方分析方法下的卡方值对比图和分割值T；

task2：利用84张标准库图像，将前42张嵌入秘密信息形成载密图像，后42不嵌入秘密信息作为原始图像，测试由task1生成的分割值T用于区分图像类型（原始图像或载秘图像）的预测准确度。

三、实验原理

简单LSB隐写算法像素值修改分析:

像素值2i改为2i+1或将2i+1改为2i，嵌入的数据d=0/1，其概率为50%。

记灰度值为j的像素数目为hj

隐写图像与原始图像相比：

隐写图像的h2i、h2i+1的值会更接近（直方图成对出现），图1为原始图像的直方图，图2为隐写图像的直方图。



找出最小的满足符合统计图像的，则可通过判断图像的与之间的大小关系得知该图像属于隐写图像或是原始图像。

四、实验设计

主函数task1.m和task2.m调用函数StgPrb.m和函数Diff.m

task1利用84张标准图像作为输入，通过嵌入秘密信息，得到84张载密图像，即利用共168张图像进行卡方分析，得到预测分割值T以及卡方值对比图。

卡方分析方式1的原始图像与载密图像卡方值对比图，如图3所示；



对卡方分析方式1的结果进行对数处理，使得其结果更明显，缩小卡方值的大小差异，入图4；



卡方分析方式2的原始图像与载密图像卡方值对比图，如图5所示；

task2利用84张标准图像作为输入，其中前42张图像嵌入秘密信息使得其图像嵌入率为1，后42张图像不嵌入任何信息，然后通过调用函数StgPrb.m和函数Diff.m得到图像的卡方值，将其与分割值T进行对比判断，预测该图像的具体类型。将预测类型与实际类型进行对比，即可以得到预测准确率。

其中使用函数StgPrb.m的预测精度为100%，使用函数Diff.m的预测精度为98.81%。

如下图6为检测方法1图像的实际类型与预测类型，其中1表示载密图像，0表示原始图像



如下图7为检测方法2图像的实际类型与预测类型，其中1表示载密图像，0表示原始图像。


 5.总结

在设计完Diff.m函数之后，发现不同图像的卡方值差距悬殊，通过log函数降低不同图像卡方值的数量级，使得不同类型的图像的卡方值结果更明显。

在网上发现一种更优的卡方分析方式，其通过chi2cdf函数巧妙的将散乱的数值变得更具规律，最终使得大多数载秘图像的卡方值达到1，原始图像的卡方值为0。测试以0.5为分割值判断图像类型时，预测精准度高达100%。


6.附录代码
task1.m

%% 清空环境
clc
clear
close all
 
%% 数据测试（84个训练图像）
cover=cell(84,1);
stego=cell(84,1);
sizes=zeros(84,2);
p1=zeros(84,2);
p2=zeros(84,2);
for k=1:84
    secret=rand(512)<0.5;
    %% 读取图像
    text=['original/',num2str(k),'.bmp'];
    cover{k}=imread(text);
    [sizes(k,1),sizes(k,2)]=size(cover{k});
   
    %% 求卡方值
    p1(k,1)=StgPrb(cover{k});
    %% 求特征值
    p2(k,1)=Diff(cover{k});
    
    %% 嵌入信息
    for i=1:sizes(k,1)
        for j=1:sizes(k,2)
            stego{k}(i,j)=bitset(cover{k}(i,j),1,secret(i,j));
        end
    end
    
    %% 求卡方值
    p1(k,2)=StgPrb(stego{k});
    %% 求特征值
    p2(k,2)=Diff(stego{k});
end
 
T=max(p2(:,2));
 
figure();
plot([1:84],p2(:,1),'r+',[1:84],p2(:,2),'b+',[1:84],ones(84,1).*T,'g-');
legend('原始图像','嵌入率为1的载密图像','分割线T');
xlabel('图像编号');
ylabel('特征值');
 
figure();
p2=log(p2);
T=log(T);
plot([1:84],p2(:,1),'r+',[1:84],p2(:,2),'b+',[1:84],ones(84,1).*T,'g-');
legend('原始图像','嵌入率为1的载密图像','分割线T');
xlabel('图像编号');
ylabel('log(特征值)');
 
figure();
plot([1:84],p1(:,1),'r+',[1:84],p1(:,2),'b+');
legend('原始图像','嵌入率为1的载密图像');
xlabel('图像编号');
ylabel('卡方值');
 
save T.mat T


task2.m

clc
clear
close all
 
load T.mat
 
%% 数据测试（84个测试图像）
cover=cell(84,1);
p1=zeros(84,1);
p2=zeros(84,1);
predict=zeros(84,2);
for k=1:42
    secret=rand(512)<0.5;
    %% 读取图像
    text=['original/',num2str(k),'.bmp'];
    cover{k}=imread(text);
    [sizes(k,1),sizes(k,2)]=size(cover{k});
    
    %% 嵌入信息
    for i=1:sizes(k,1)
        for j=1:sizes(k,2)
            cover{k}(i,j)=bitset(cover{k}(i,j),1,secret(i,j));
        end
    end
    
    %% 求卡方值
    p1(k)=StgPrb(cover{k});
    %% 求特征值
    p2(k)=Diff(cover{k});
    p2(k)=log(p2(k));
    
    if p1(k)>0.5
        predict(k,1)=1;
    else
        predict(k,1)=0;
    end
    
    if p2(k)>T
        predict(k,2)=0;
    else
        predict(k,2)=1;
    end
end
 
for k=43:84
    %% 读取图像
    text=['original/',num2str(k),'.bmp'];
    cover{k}=imread(text);
    [sizes(k,1),sizes(k,2)]=size(cover{k});
 
    %% 求卡方值
    p1(k)=StgPrb(cover{k});
    %% 求特征值
    p2(k)=Diff(cover{k});
    p2(k)=log(p2(k));
    
    if p1(k)>0.5
        predict(k,1)=1;
    else
        predict(k,1)=0;
    end
    
    if p2(k)>T
        predict(k,2)=0;
    else
        predict(k,2)=1;
    end
end
 
%% 求预测准确度
predict_rate1=0;
predict_rate2=0;
for k=1:42
    if predict(k,1)==1
        predict_rate1=predict_rate1+1;
    end
    if predict(k,2)==1
        predict_rate2=predict_rate2+1;
    end
end
 
for k=43:84
    if predict(k,1)==0
        predict_rate1=predict_rate1+1;
    end
    if predict(k,2)==0
        predict_rate2=predict_rate2+1;
    end
end
predict_rate1=predict_rate1/84
predict_rate2=predict_rate2/84
dx=[1:84];
real=[ones(1,42),zeros(1,42)];
figure();
plot(dx,real,'r*',dx,predict(:,1),'bo');
title("检验方法1");
legend('实际类型','预测类型');
figure();
plot(dx,real,'r*',dx,predict(:,2),'bo');
title("检验方法2");
legend('实际类型','预测类型');


StgPrb.m

function p=StgPrb(x)
n=sum(hist(x,[0:255]),2);
h2i=n([3:2:255]);
h2is=(h2i+n([4:2:256]))/2;
filter=(h2is~=0);
k=sum(filter);
idx=zeros(1,k);
for i=1:127
    if filter(i)==1
        idx(sum(filter(1:i)))=i;
    end
end
r=sum(((h2i(idx)-h2is(idx)).^2)./(h2is(idx)));
p=1-chi2cdf(r,k-1);


Diff.m

function p=Diff(x)
p=0;
[n,~]=imhist(x);
for i=0:127
    p=p+(n(2*i+1)-n(2*i+2))^2;
end