%?分段匹配追蹤（StOMP）算法。
%??讀文件
%?X=imread（‘C:\Users\zkdn\Desktop\MATLABdata\I30\I30.jpg‘）;
%?X=double（X）;
%?[ab]=size（X）;
load?C:\Users\zkdn\Desktop\MATLABdata\I1\I1.mat
X=I1;
background=imopen（Xstrel（‘disk‘15））;%獲取背景信息
X=imsubtract（Xbackground）;%利用函數去除背景

X=double（X）;
[ab]=size（X）;
%??小波變換矩陣生成
ww=dwt2（a‘haar‘）;
%??小波變換讓圖像稀疏化（注意該步驟會耗費時間，但是會增大稀疏度）
X1=ww*sparse（X）*ww‘;
%?X1=X;
X1=full（X1）;

%??隨機矩陣生成
M=200;
R=randn（Ma）;
%?R=randn（a）;
%?R=mapminmax（R0255）;
%?R=round（R）;

%??測量值
Y=R*X1;

%??OMP算法
%??恢復矩陣
X2=zeros（ab）;?
%??按列循環
for?i=1:b?
????%??通過OMP，返回每一列信號對應的恢復值（小波域）
????rec=CS_StOMP（Y（:i）Ra）;%[?theta?]?=?CS_StOMP（?yASts?）
????%??恢復值矩陣，用于反變換
????X2（:i）=rec;
end
%%%%
%?figure（4）;
%?X2=full（X2）;
%?imshow（uint8（X2））;
%?title（‘OMP計算后圖像‘）;

%??原始圖像
figure（1）;
imshow（uint8（X））;
title（‘原始圖像‘）;

%??變換圖像
%?figure（2）;
%?imshow（uint8（X1））;
%?title（‘小波變換后的圖像‘）;

%??壓縮傳感恢復的圖像
figure（3）;
%??小波反變換
X3=ww‘*sparse（X2）*ww;??
%?X3=X2;
X3=full（X3）;
imshow（uint8（X3））;
title（‘STOMP恢復的圖像‘）;

%??誤差（PSNR）
%??MSE誤差
errorx=sum（sum（abs（X3-X）.^2））;????????
%??PSNR
psnr=10*log10（255*255/（errorx/a/b））???



function?[?theta?]?=?CS_StOMP（?yAS?）??
%CS_StOMP?Summary?of?this?function?goes?here??
%Version:?1.0?written?by?jbb0523?@2015-04-29??
%???Detailed?explanation?goes?here??
%???y?=?Phi?*?x??
%???x?=?Psi?*?theta??
%???y?=?Phi*Psi?*?theta??
%???令?A?=?Phi*Psi?則y=A*theta??
%???S?is?the?maximum?number?of?StOMP?iterations?to?perform??
%???ts?is?the?threshold?parameter??
%???現在已知y和A，求theta??
%???Reference:Donoho?D?L，Tsaig?Y，Drori?I，Starck?J?L．Sparse?solution?of??
%???underdete