灰度图像的中值滤波取证

实验报告
实验题目：灰度图像的中值滤波取证
课程题目：多媒体安全技术
学号：5120180885
1、实验目的和要求
2.1 h0特征散点图
2.2 rho=h0/h1特征散点图
2.3 加权rho 特征散点图
2、实验设备（环境）
Matlab
3、实验步骤
（1）h0特征散点图，使用matlab编写如下代码

（2）rho=h0/h1特征散点图，使用matlab编写如下代码

（3）加权rho 特征散点图，使用matlab编写如下代码

4、实验结果
（1）h0特征散点图

（2）rho=h0/h1特征散点图

（3）加权rho 特征散点图

5、讨论和分析（分析实验原理，为什么会得到这样的结果以及在试验中应注意的问题）
实验原理：中值滤波是一种非线性数字滤波器技术，经常用于去除图像或者其它信号中的噪声。这个设计思想就是检查输入信号中的采样并判断它是否代表了信号，使用奇数个采样组成的观察窗实现这项功能。观察窗口中的数值进行排序，位于观察窗中间的中值作为输出。然后，丢弃最早的值，取得新的采样，重复上面的计算过程。中值滤波是图像处理中的一个常用步骤，它对于斑点噪声和椒盐噪声来说尤其有用。保存边缘的特性使它在不希望出现边缘模糊的场合也很有用。
原始图像，经过中值滤波处理，图像更加趋于平滑，趋于一个平稳的值，去除噪声的同时保护信号的边缘。所以看出图像经过不断处理，最终趋于稳定，接近于一条直线。
实验中应注意的问题：整个实验自己写起来代码还是有些困难的，在学长和老师帮助下解决一个一个小问题，最终成功实现。需要注意其中很多小问题，比如分号，小写L和1的区别等等。复习了Plot画图的过程，不仅可以画直线，也可以用点或者*来代替。

代码

mat1

file_path = 'D:\schoolstudy\大三下\隐写课\实验2\100\';
img_path_list = dir(strcat(file_path,'*.tif'));
img_num = length(img_path_list);
images = cell(1,img_num);

h01=[];
h02=[];
% h03=[];
% h04=[];
% h05=[];
% h06=[];

if img_num>0
    for i = 1:img_num
    img_name = img_path_list(i).name;
    
fprintf('正在处理%s\n',img_name);

[k,l] = deal(0,1);
I = imread(img_name);
I = double(I);
I_med = medfilt2(I,[3,3],'symmetric');

h01(i) = med_dif(I,k,l);
h02(i) = med_dif(I_med,k,l);
% h03(i) = med_dif2(I,k,l);
% h04(i) = med_dif2(I_med,k,l);
% h05(i) = h01(i)/h03(i);
% h06(i) = h02(i)/h04(i);

    end 
end

plot([1:100],h01(1:100),'k*',[1:100],h02(1:100),'b');
% plot([1:100],h05(1:100),'r*',[1:100],h06(1:100),'b.');
xlabel('图片序号');
ylabel('h0');

function h0_hat = med_dif(I,k,l)
    I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
    s = size(I);
    y = s(1)*s(2);
    h0 = sum(I2(:)==0);
    h0_hat = h0/y;
end
% 
% function h0_hat = med_dif2(I,k,l)
%     I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
%     s = size(I);
%     y = s(1)*s(2);
%     h0 = sum(I2(:)==1);
%     h0_hat = h0/y;
% end

mat2

file_path = 'D:\schoolstudy\大三下\隐写课\实验2\100\';
img_path_list = dir(strcat(file_path,'*.tif'));
img_num = length(img_path_list);
images = cell(1,img_num);

h01=[];
h02=[];
h03=[];
h04=[];
h05=[];
h06=[];

if img_num>0
    for i = 1:img_num
    img_name = img_path_list(i).name;
    
fprintf('正在处理%s\n',img_name);

[k,l] = deal(0,1);
 I = imread(img_name);
I = double(I);
I_med = medfilt2(I,[3,3],'symmetric');

h01(i) = med_dif(I,k,l);
h02(i) = med_dif(I_med,k,l);
h03(i) = med_dif2(I,k,l);
h04(i) = med_dif2(I_med,k,l);
h05(i) = h01(i)/h03(i);
h06(i) = h02(i)/h04(i);

    end 
end

% plot([1:100],h01(1:100),'k.',[1:100],h02(1:100),'b.');
plot([1:100],h05(1:100),'r*',[1:100],h06(1:100),'b.');
xlabel('图片序号');
ylabel('h0/h1');

function h0_hat = med_dif(I,k,l)
    I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
    s = size(I);
    y = s(1)*s(2);
    h0 = sum(I2(:)==0);
    h0_hat = h0/y;
end

function h0_hat = med_dif2(I,k,l)
    I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
    s = size(I);
    y = s(1)*s(2);
    h0 = sum(I2(:)==1);
    h0_hat = h0/y;
end

mat3

file_path = 'D:\schoolstudy\大三下\隐写课\实验2\100\';
img_path_list = dir(strcat(file_path,'*.tif'));
img_num = length(img_path_list);
images = cell(1,img_num);

h01=[];
h02=[];
h03=[];
h04=[];
h05=[];
h06=[];
q = [];
q1=[];

if img_num>0
    for i = 1:img_num
    img_name = img_path_list(i).name;
    
fprintf('正在处理%s\n',img_name);

[k,l] = deal(0,1);
I = imread(img_name);
I = double(I);
I_med = medfilt2(I,[3,3],'symmetric');

% h01(i) = med_dif(I,k,l);
% h02(i) = med_dif(I_med,k,l);
% h03(i) = med_dif2(I,k,l);
% h04(i) = med_dif2(I_med,k,l);
% h05(i) = h01(i)/h03(i);
% h06(i) = h02(i)/h04(i);

q(i) = get_med_qb(I,32,k,l);
q1(i) = get_med_qb(I_med,32,k,l);
    end 
end

% plot([1:100],h01(1:100),'k.',[1:100],h02(1:100),'b.');
% plot([1:100],h05(1:100),'r*',[1:100],h06(1:100),'b.');
% xlabel('图片序号');
% ylabel('h0/h1');
% scatter([1:100],q(1:100));
% hold on
% scatter([1:100],q1(1:100));
plot([1:100],q(1:100),'r*',[1:100],q1(1:100),'b.');

function h0_hat = med_dif(I,k,l)
    I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
    s = size(I);
    y = s(1)*s(2);
    h0 = sum(I2(:)==0);
    h0_hat = h0/y;
end

function h0_hat = med_dif2(I,k,l)
    I2 = I(1:end-k,1:end-1)-I(1+k:end,1+l:end);
    s = size(I);
    y = s(1)*s(2);
    h0 = sum(I2(:)==1);
    h0_hat = h0/y;
end

function q = get_med_qb(img,B,k,l)
    s = size(img);
    tt =1;
    qb_hat = zeros(floor((s(1)/B))*floor((s(2)/B)),1);
    for i = 1:B:s(1)
        for j = 1:B:s(2)
            temp = img(i:min(i+B-1,s(1)),j:min(j+B-1,s(2)));
            I_d = temp(1:end-k,1:end-1)-temp(1+k:end,1+l:end);
            h0 = sum(I_d(:)==0);
            h1 = sum(I_d(:)==1);
                if h1 ==0
                    h1=1;
                end
                qb = h0/h1;
                wb = 1-(h0/(B*B-B));
                qb_hat(tt) = qb*wb;
                tt = tt+1;
        end
    end
    q = median(qb_hat(:));
end

rs_program

file_path = 'D:\schoolstudy\大三下\隐写课\实验一\100\';
img_path_list = dir(strcat(file_path,'*.tif'));
img_num = length(img_path_list);
images = cell(1,img_num);
rs_final=zeros(10,4);

if img_num>0
    for i = 1:img_num
    img_name = img_path_list(i).name;
    
fprintf('正在处理%s\n',img_name);
t=imread(img_name);
[p,q]=size(t);
I=t(1:p,1:q);
sz=size(I);
m=floor(sz(1)/8);
n=floor(sz(2)/8);
rs=zeros(10,4);
cor=zeros(1,3);
M=randsrc(8,8,[0 1]);
tmp=zeros(8,8);
for i=1 : 10
    r=floor(sz(1)*i/10);
    c=floor(sz(2)*i/10);
    msg=randsrc(r,c,[0 1]);
    s=I;
    s(1:r,1:c)=bitset(s(1:r,1:c),1,msg);

for j=1 : m
    rv=[(j-1)*8+1:j*8];
    for k=1 : n
        cv=[(k-1)*8+1:k*8];
        tmp=s(rv,cv);
        cor(1)=SpaceCor(tmp);
        cor(2)=SpaceCor(f1(tmp,M));
        cor(3)=SpaceCor(f_1(tmp,M));
        if cor(2)>cor(1)
            rs(i,1)=rs(i,1)+1;
        else
            if cor(2)<cor(1)
                rs(i,2)=rs(i,2)+1;
            end
        end
        if cor(3)>cor(1)
            rs(i,3)=rs(i,3)+1;
        else
            if cor(3)<cor(1)
                rs(i,4)=rs(i,4)+1;
            end
        end
    end
end
end
rs=rs/ (m * n);
for i=1 : 10
    rs_final(i,1) = rs_final(i,1)+rs(i,1);
    rs_final(i,2) = rs_final(i,2)+rs(i,2);
    rs_final(i,3) = rs_final(i,3)+rs(i,3);
    rs_final(i,4) = rs_final(i,4)+rs(i,4);
end
    end
end
rs_final = rs_final/100;
figure;
plot([0.1:0.1:1.0],rs_final(:,1),[0.1:0.1:1.0],rs_final(:,2),[0.1:0.1:1.0],rs_final(:,3),[0.1:0.1:1.0],rs_final(:,4));
legend('Rm','Sm','R-m','S-m');

function y=SpaceCor(x)
idx = [1,1; 1,2; 2,1; 3,1; 2,2; 1,3; 1,4; 2,3; 3,2; 4,1; 5,1; 4,2; 3,3; 2,4; 1,5;...
    1,6; 2,5; 3,4; 4,3; 5,2; 6,1; 7,1; 6,2; 5,3; 4,4; 3,5; 2,6; 1,7;...
    1,8; 2,7; 3,6; 4,5; 5,4; 6,3; 7,2; 8,1; 8,2; 7,3; 6,4; 5,5; 4,6; 3,7; 2,8;...
    3,8; 4,7; 5,6; 6,5; 7,4; 8,3; 8,4; 7,5; 6,6; 5,7; 4,8; 5,8; 6,7; 7,6; 8,5;...
    8,6; 7,7; 6,8; 7,8; 8,7; 8,8;];
n=64;
y=0;
for i=1 : n-1
    r1=idx(i,1);
    c1=idx(i,2);
    r2=idx(i+1,1);
    c2=idx(i+1,2);
    y=y + abs(x(r1,c1)- x(r2,c2));
end
end

function y=f_1(x,M)
szx=size(x);
szm=size(M);
if szx ~= szm
    fprintf(1,'the size of x and M must be the same!\n');
end
y=x;
for i=1:szx(1)
    for j=1:szx(2)
        if M(i,j)~=0
            odd=mod(x(i,j),2);
            if odd == 0
                y(i,j)=x(i,j) - 1;
            else
                y(i,j)=x(i,j) + 1;
            end
        end
    end
end
end

function y = f1(x,M)
szx=size(x);
szm=size(M);
if szx~=szm
    fprintf(1,'the size of x and M must be the same!\n');
end
y=x;
for i=1 :szx(1)
    for j=1: szx(2)
        if M(i,j)~=0
            odd = mod(x(i,j),2);
            if odd ==0
                y(i,j)=x(i,j)+1;
            else
                y(i,j)=x(i,j)-1;
            end
        end
    end
end
end