基于位平面分解的超混沌图像加密

clear all; clc;clear; P=imread('1.png'); P=rgb2gray(P); figure;imhist(P);title('原始图片直方图'); iptsetpref('imshowborder','tight'); figure;imshow(P);title('明文图像'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 原始图像相邻像素相关性分析 [M,N]=size(P); I1=P; NN=5000; %随机取5000对像素点 x1=ceil(rand(1,NN)*(M-1)); %生成5000个1~M-1的随机整数作为行 y1=ceil(rand(1,NN)*(N-1)); %生成5000个1~N-1的随机整数作为列 %预分配内存 XX_R_SP=zeros(1,NN);YY_R_SP=zeros(1,NN); %水平 XX_R_CZ=zeros(1,NN);YY_R_CZ=zeros(1,NN); %垂直 XX_R_DJX=zeros(1,NN);YY_R_DJX=zeros(1,NN); %对角线 for i=1:NN %水平 XX_R_SP(i)=I1(x1(i),y1(i)); YY_R_SP(i)=I1(x1(i)+1,y1(i)); %垂直 XX_R_CZ(i)=I1(x1(i),y1(i)); YY_R_CZ(i)=I1(x1(i),y1(i)+1); %对角线 XX_R_DJX(i)=I1(x1(i),y1(i)); YY_R_DJX(i)=I1(x1(i)+1,y1(i)+1); end %水平 figure;scatter(XX_R_SP,YY_R_SP,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻水平方向像素灰度值');title('原始图像水平相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); %垂直 figure;scatter(XX_R_CZ,YY_R_CZ,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻垂直方向像素灰度值');title('原始图像垂直相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); %对角线 figure;scatter(XX_R_DJX,YY_R_DJX,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻对角线方向像素灰度值');title('原始图像对角线相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); EX1_R=0;EY1_SP_R=0;DX1_R=0;DY1_SP_R=0;COVXY1_SP_R=0; %计算水平相关性时需要的变量 EY1_CZ_R=0;DY1_CZ_R=0;COVXY1_CZ_R=0; %垂直 EY1_DJX_R=0;DY1_DJX_R=0;COVXY1_DJX_R=0; %对角线 I1=double(I1); for i=1:NN %第一个像素点的E，水平、垂直、对角线时计算得出的第一个像素点的E相同，统一用EX1表示 EX1_R=EX1_R+I1(x1(i),y1(i)); %第二个像素点的E，水平、垂直、对角线的E分别对应EY1_SP、EY1_CZ、EY1_DJX EY1_SP_R=EY1_SP_R+I1(x1(i),y1(i)+1); EY1_CZ_R=EY1_CZ_R+I1(x1(i)+1,y1(i)); EY1_DJX_R=EY1_DJX_R+I1(x1(i)+1,y1(i)+1); end %统一在循环外除以像素点对数1000，可减少运算次数 EX1_R=EX1_R/NN; EY1_SP_R=EY1_SP_R/NN; EY1_CZ_R=EY1_CZ_R/NN; EY1_DJX_R=EY1_DJX_R/NN; for i=1:NN %第一个像素点的D，水平、垂直、对角线时计算得出第一个像素点的D相同，统一用DX表示 DX1_R=DX1_R+(I1(x1(i),y1(i))-EX1_R)^2; %第二个像素点的D，水平、垂直、对角线的D分别对应DY1_SP、DY1_CZ、DY1_DJX DY1_SP_R=DY1_SP_R+(I1(x1(i),y1(i)+1)-EY1_SP_R)^2; DY1_CZ_R=DY1_CZ_R+(I1(x1(i)+1,y1(i))-EY1_CZ_R)^2; DY1_DJX_R=DY1_DJX_R+(I1(x1(i)+1,y1(i)+1)-EY1_DJX_R)^2; %两个相邻像素点相关函数的计算，水平、垂直、对角线 COVXY1_SP_R=COVXY1_SP_R+(I1(x1(i),y1(i))-EX1_R)*(I1(x1(i),y1(i)+1)-EY1_SP_R); COVXY1_CZ_R=COVXY1_CZ_R+(I1(x1(i),y1(i))-EX1_R)*(I1(x1(i)+1,y1(i))-EY1_CZ_R); COVXY1_DJX_R=COVXY1_DJX_R+(I1(x1(i),y1(i))-EX1_R)*(I1(x1(i)+1,y1(i)+1)-EY1_DJX_R); end %统一在循环外除以像素点对数1000，可减少运算次数 DX1_R=DX1_R/NN; DY1_SP_R=DY1_SP_R/NN; DY1_CZ_R=DY1_CZ_R/NN; DY1_DJX_R=DY1_DJX_R/NN; COVXY1_SP_R=COVXY1_SP_R/NN; COVXY1_CZ_R=COVXY1_CZ_R/NN; COVXY1_DJX_R=COVXY1_DJX_R/NN; %水平、垂直、对角线的相关性 RXY1_SP_R=COVXY1_SP_R/sqrt(DX1_R*DY1_SP_R); RXY1_CZ_R=COVXY1_CZ_R/sqrt(DX1_R*DY1_CZ_R); RXY1_DJX_R=COVXY1_DJX_R/sqrt(DX1_R*DY1_DJX_R); %% 原始图片信息熵 T1_R=imhist(I1); %统计图像R通道灰度值从0~255的分布情况，存至T1 S1_R=sum(T1_R); %计算整幅图像R通道的灰度值 xxs1_R=0; %原始图片R通道相关性 for i=1:256 pp1_R=T1_R(i)/S1_R; if pp1_R~=0 xxs1_R=xxs1_R-pp1_R*log2(pp1_R); end end %%%%%%%%%%%%%%%%%%% [M,N]=size(P);P=double(P); % [h,w,n]=size(P); imageBit = myImageBitPlaneSlicing(P); % for循环显示比特平面 for i=1:1:8 % 得到第i比特面 temp = imageBit(:,:,i); % 得到非0值的坐标 index = find(temp>0); % 值为1的用255表示，突出显示效果，0值依旧为0，不处理 temp(index) = 255; % 矩阵类型转化为uint8并显示比特图像 temp = uint8(temp); subplot(3,3,i+1); imshow(temp); % 带变化参数的标题是用一个字符矩阵表示 str = ['第',num2str(i),'比特平面']; title(str); end B8=imageBit(:,:,8); %置乱与复原的共同参数 n=10; a=3;b=5; N=M; %置乱 Pn=zeros(M,N); for i=1:n for y=1:M for x=1:N xx=mod((x-1)+b*(y-1),N)+1; yy=mod(a*(x-1)+(a*b+1)*(y-1),N)+1; imgn(yy,xx)=B8(y,x); end end P=imgn; end figure; imshow(imgn,[]); imageBit(:,:,8)=imgn; result1 = image_construct(imageBit,4); figure; imshow(result1,[]); %%%%%%%%%%%%%%%%%%% n=2*M*N; h=0.002;t=800; a=10;b=8/3;c=28;r=-1; x0=1.1;y0=2.2;z0=3.3;w0=4.4; s=zeros(1,n); for i=1:n+t K11=a*(y0-x0)+w0;K12=a*(y0-(x0+K11*h/2))+w0; K13=a*(y0-(x0+K12*h/2))+w0;K14=a*(y0-(x0+h*K13))+w0; x1=x0+(K11+K12+K13+K14)*h/6; K21=c*x1-y0-x1*z0;K22=c*x1-(y0+K21*h/2)-x1*z0; K23=c*x1-(y0+K22*h/2)-x1*z0;K24=c*x1-(y0+h*K23)-x1*z0; y1=y0+(K21+K22+K23+K24)*h/6; K31=x1*y1-b*z0;K32=x1*y1-b*(z0+K31*h/2); K33=x1*y1-b*(z0+K32*h/2);K34=x1*y1-b*(z0+h*K33); z1=z0+(K31+K32+K33+K34)*h/6; K41=-y1*z1+r*w0;K42=-y1*z1+r*(w0+K41*h/2); K43=-y1*z1+r*(w0+K42*h/2);K44=-y1*z1+r*(w0+h*K43); w1=w0+(K41+K42+K43+K44)*h/6; x0=x1;y0=y1;z0=z1;w0=w1; if i>t s(i-t)=x1; if mod((i-t),3000)==0 x0=x0+h*sin(y0); end end end S=mod(floor(s*pow2(16)),256); S1=S(1:M*N);S2=S(M*N+1:2*M*N);B=zeros(M,N);C=zeros(M,N); S1=reshape(S1,M,N); S2=reshape(S2,M,N); S1=uint8(S1); S2=uint8(S2); R1=uint8(zeros(M,N)); k=1; for i=1:M for j=1:N P1(i,j)=bitxor(S1(k),result1(i,j));%异或操作，结果放到C矩阵中 k=k+1; end end figure; imshow(P1);%加密图像 figure;imhist(P1);title('加密图片直方图'); %% 加密图像相邻图像相关性分析 %相关性曲线 %水平 x1=ceil(rand(1,NN)*(M-1)); %生成5000个1~M-1的随机整数作为行 y1=ceil(rand(1,NN)*(N-1)); %生成5000个1~N-1的随机整数作为列 Q_R=uint8(P1); XX_R_SP=zeros(1,NN);YY_R_SP=zeros(1,NN); %预分配内存 %垂直 XX_R_CZ=zeros(1,NN);YY_R_CZ=zeros(1,NN); %预分配内存 %对角线 XX_R_DJX=zeros(1,NN);YY_R_DJX=zeros(1,NN); %预分配内存 for i=1:NN %水平 XX_R_SP(i)=Q_R(x1(i),y1(i)); YY_R_SP(i)=Q_R(x1(i)+1,y1(i)); %垂直 XX_R_CZ(i)=Q_R(x1(i),y1(i)); YY_R_CZ(i)=Q_R(x1(i),y1(i)+1); %对角线 XX_R_DJX(i)=Q_R(x1(i),y1(i)); YY_R_DJX(i)=Q_R(x1(i)+1,y1(i)+1); end %水平 figure;scatter(XX_R_SP,YY_R_SP,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻水平方向像素灰度值');title('加密图像水平相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); %垂直 figure;scatter(XX_R_CZ,YY_R_CZ,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻垂直方向像素灰度值');title('加密图像垂直相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); %对角线 figure;scatter(XX_R_DJX,YY_R_DJX,18,'filled');xlabel('随机点像素灰度值');ylabel('与该点相邻对角线方向像素灰度值');title('加密图像对角线相邻元素相关性点图');axis([0 255,0 255]);set(gca,'XTick',0:15:255);set(gca,'YTick',0:15:255); %% 加密后信息熵 T2_R=imhist(Q_R); S2_R=sum(T2_R); xxs2_R=0; for i=1:256 pp2_R=T2_R(i)/S2_R; if pp2_R~=0 xxs2_R=xxs2_R-pp2_R*log2(pp2_R); end end Q_R=double(Q_R); EX2_R=0;EY2_SP_R=0;DX2_R=0;DY2_SP_R=0;COVXY2_SP_R=0; %水平 EY2_CZ_R=0;DY2_CZ_R=0;COVXY2_CZ_R=0; %垂直 EY2_DJX_R=0;DY2_DJX_R=0;COVXY2_DJX_R=0; %对角线 for i=1:NN %第一个像素点的E，水平、垂直、对角线时计算得出的第一个像素点的E相同，统一用EX2表示 EX2_R=EX2_R+Q_R(x1(i),y1(i)); %第二个像素点的E，水平、垂直、对角线的E分别对应EY2_SP、EY2_CZ、EY2_DJX EY2_SP_R=EY2_SP_R+Q_R(x1(i),y1(i)+1); EY2_CZ_R=EY2_CZ_R+Q_R(x1(i)+1,y1(i)); EY2_DJX_R=EY2_DJX_R+Q_R(x1(i)+1,y1(i)+1); end %统一�