图像隐写分析-WSVD小波奇异值分解特征编程实现

Farid等人通过对大量的载体图像分析发现，在小波状分解中存在强大的高阶统计规律【1】。嵌入信息之后会明显的改变这些统计特征，可以被用来作为分析的工具，并且可以结合支持向量机（线性和非线性）来检测这些统计偏差。小波特征集的方法依赖于构建载体图像的高阶统计模型并分析隐密图像与这些模型的偏差。

将图像进行三级小波分解，保留每次分解产生的三个方向的子带系数和低频分量的绝对值，一共有12个子带。特征的计算：使用小波对给定图像进行n级分解，统计模型由尺度i = 1，…，n - 1处的子带系数的平均值，方差，偏度和峰度组成，这些系数跟它们的位置，方向等都有一定相关性。第二组统计是基于系数幅度的最优线性预测器中的误差。为了理解算法的流程，以尺度i的垂直带Vi（x,y）为例子。

正是从这个错误，收集了额外的统计数据，即平均值，方差，偏度和峰度。 对于尺度为i = 1，…，n-1处的每个垂直子带重复该过程。对于水平和对角子带重复类似的过程，其中对H与D子带的计算中的最后两个加权量分别是 。

对于一个尺寸为M×N的矩阵F(i,j)，其均值，方差，偏度，峰度的计算方式如下：

为了检测隐藏信息的存在，他们采用了从多尺度分解得到的统计量的模型。该模型包括基本系数统计以及来自系数幅度的最优线性预测器的误差统计。这些高级的统计特征记录了载体图像的属性。更重要的是，当图像被嵌入秘密信息之后，这些特征会发生显著的改变，这对检测分析是有好处的。在测试实验中，使用Jsteg算法和两种版本的OutGuess算法作为测试样本，实验结果表明，该特征对OutGuess0.2算法的检测效果较差，总体效果一般，但是有一定的分类能力。

WSVD特征

WSVD是小波奇异值分解（Wavelet Singular Value Decomposition）的简称，它是在小波变换的基础上，对小波子带矩阵进行奇异值分解。奇异值分解可以作为矩阵化简的工具，而且奇异值分解对秘密信息的嵌入十分敏感。使用奇异值分解提取特征的方式同直接使用小波系数提取特征相比，可以大大的减少计算复杂度和特征维度。
在对一幅图像进行小波变换时会产生四个矩阵，尺寸为原矩阵大小的一半，分别是变换出来的低频分量和水平、垂直、对角线高频分量，这几个分量经过逆变换之后可以恢复原始图像。学者们发现大多数的嵌入秘密信息的过程都可以在小波域被捕获，所以小波变换也成了隐密分析的重要工具
奇异值分解作为一种矩阵分解的方法，可以用来作为分析特征，因为当图像嵌入秘密信息之后，矩阵的SVD的秩将会有一个升序的变化。为了分析嵌入过程，假定隐藏的信息是一种额外添加的、零均值的高斯白噪声z，嵌入在原始信号s中，来建立隐密信号x。这可以表示为：

一个矩阵（图像）的奇异值分解成分对嵌入的数据相当敏感，可以作为图像隐密分析中隐藏内容的标志。

本文采用的方案如下：首先将图像通过三级Haar小波变换分解，并且分解第一尺度对角线子带以提高系统的性能，这样一共得到16个矢量。然后，计算各子带的奇异值，从中提取五个典型的特征，这些功能包括前三阶中心矩、对数均值和条件数即最大奇异值的和最小奇异值的比。因此，在小波奇异值分解这部分一共得到80维特征。

经过上述的特征提取过程，DCT融合特征一共有193+81+80=274维特征用于隐密分析。

实现代码

依旧是matlab：

function y=wsvd80feature(I)  
%I为输入的灰度图像 

%y输出为16个子带，每个子带计算奇异值向量 取五个统计特征 共80个
%每个奇异值向量去三个特征矩 对数均值 最大奇异值的和最小奇异值的比 
%个人观察统计矩为2，3，4阶中心矩
%每次变换有一个低频成分 三个高频垂直 水平 对角线成分


%第一级尺度变换
%e=imread('d:\1.jpg');
%I=rgb2gray(e);
%I=I(1:512,1:512);
result=zeros(1,80);
temp=zeros(1,5);
[c1,l1]=wavedec2(I,1,'Haar');  %二维一层小波分解
ca1=wcodemat(appcoef2(c1,l1,'Haar',1),255);  %低频分量
ch1=wcodemat(detcoef2('h',c1,l1,1),255);  %高频水平分量
cv1=wcodemat(detcoef2('v',c1,l1,1),255);  %高频垂直分量
cd1=wcodemat(detcoef2('d',c1,l1,1),255);  %高频对角分量
%第一级尺度对角线分解
[c4,l4]=wavedec2(cd1,1,'Haar');  %二维一层小波分解
cad=wcodemat(appcoef2(c4,l4,'Haar',1),255);  %低频分量
chd=wcodemat(detcoef2('h',c4,l4,1),255);  %高频水平分量
cvd=wcodemat(detcoef2('v',c4,l4,1),255);  %高频垂直分量
cdd=wcodemat(detcoef2('d',c4,l4,1),255);  %高频对角分量
%第二级尺度变换
[c2,l2]=wavedec2(ca1,1,'Haar');
ca2=wcodemat(appcoef2(c2,l2,'Haar',1),255);
ch2=wcodemat(detcoef2('h',c2,l2,1),255);
cv2=wcodemat(detcoef2('v',c2,l2,1),255);
cd2=wcodemat(detcoef2('d',c2,l2,1),255);
%第三级尺度变换
[c3,l3]=wavedec2(ca2,1,'Haar');
ca3=wcodemat(appcoef2(c3,l3,'Haar',1),255);
ch3=wcodemat(detcoef2('h',c3,l3,1),255);
cv3=wcodemat(detcoef2('v',c3,l3,1),255);
cd3=wcodemat(detcoef2('d',c3,l3,1),255);


svd1=svd(ca1);
svd2=svd(ch1);
svd3=svd(cv1);
svd4=svd(cd1);
svd5=svd(cad);
svd6=svd(chd);
svd7=svd(cvd);
svd8=svd(cdd);
svd9=svd(ca2);
svd10=svd(ch2);
svd11=svd(cv2);
svd12=svd(cd2);
svd13=svd(ca3);
svd14=svd(ch3);
svd15=svd(cv3);
svd16=svd(cd3);

for ik=1:16
    switch ik
            case 1
           temp=feature5get(svd1);
             case 2
           temp=feature5get(svd2);
             case 3
           temp=feature5get(svd3);
             case 4
           temp=feature5get(svd4);
             case 5
           temp=feature5get(svd5);
             case 6
           temp=feature5get(svd6);
             case 7
           temp=feature5get(svd7);
             case 8
           temp=feature5get(svd8);
             case 9
           temp=feature5get(svd9);
             case 10
           temp=feature5get(svd10);
             case 11
           temp=feature5get(svd11);
             case 12
           temp=feature5get(svd12);
             case 13
           temp=feature5get(svd13);
             case 14
           temp=feature5get(svd14);
             case 15
           temp=feature5get(svd15);
             case 16
           temp=feature5get(svd16);       
    end 
     result((ik-1)*5+1:ik*5)=temp(1,:);
end

%c1=uint8([ca2,ch2;cv2,cd2]);
%subplot(1,2,1);
%imshow(I);
%subplot(1,2,2);
%imshow(c1);

y=result;
end

其中五个统计量的函数如下：

function y=feature5get(x)  
%子带特征值提取 五个
%三个统计矩 对数均值 最大最小比
result=[0,0,0,0,0];
% result(1)=log10(sum((x-mean(x)).^2)/length(x)); %二阶中心矩
% result(2)=log10(sum((x-mean(x)).^3)/length(x)); %三阶中心矩
% result(3)=log10(sum((x-mean(x)).^4)/length(x)); %四阶中心矩
% result(4)=sum(log(x))/length(x);
% result(5)=max(x)/min(x);

result(1)=sum((x-mean(x)).^2)/length(x); %二阶中心矩
result(2)=sum(((x-mean(x))/sqrt(result(1))).^3)/length(x); %三阶中心矩
result(3)=sum(((x-mean(x))/sqrt(result(1))).^4)/length(x); %四阶中心矩

empty=length(x(x==0));  %除去0值，防止无穷大
x(x==0)=1;

result(5)=log10(max(x)/min(x)); %最大值除以最小值 

result(4)=sum(log10(x))/(length(x)-empty);  %数据中有0 除掉吧 要不没办法计算 inf  减去0的数目
y=result;
end

这个特征用处还是有的，就是感觉原理有点不太理解

【1】 LYU S, FARID H. Steganalysis using higher-order image statistics[J]. IEEE Transactions on Information Forensics & Security, 2006, 1(1):111-119.