clc;clear;close all;tic
T=imread('picture.png');T=rgb2gray(T);
T=double(T)/255;[m,n]=size(T);
subplot(151);imshow(T,[]);title('object');
T=reshape(T,[],1);%转换为一个列向量
iteration=10000; %用iteration个散斑照射物体
s=rng;mu=10;
O=zeros(iteration,n*n);Ri_sum=0;Ixy_sum=zeros(1,m*n);
for i=1:iteration
Ixy=dct(randn(m,n));
Ixy=reshape(Ixy,1,[]);%转换为一行向量
Ixy_sum=Ixy_sum+Ixy;
Ri=sum(Ixy(:));
Ri_sum=Ri_sum+Ri;
cc(i)=corr2(Ixy,T');
O(i,:)=O(i,:)+Ixy;%将所有的散斑组成和一个大矩阵
end
B=O*T;B_av=sum(B(:))/iteration;
Ri_av=Ri_sum/iteration;
Ixy_av=Ixy_sum/iteration;
%%
rng(s);G=zeros(1,m*n);test=zeros(1,m*n);Bi=zeros(iteration,1);
for i=1:iteration
Ixy=dct(randn(m,n));
Ixy=reshape(Ixy,1,[]);
Ri=sum(Ixy(:));
B_copy=B(i,1)-(B_av/Ri_av)*Ri;
Bi(i,1)=Bi(i,1)+B_copy;
G=G+B_copy*(O(i,:)-Ixy_av);
end
G=G/iteration;G=reshape(G,m,n);G=MatMap(G,0,1);
subplot(152);imshow(G,[]);title('all');
T=reshape(T,m,n);cc_ori=corr2(T,G);
toc