%??此函數用FISTA算法解決壓縮感知（傅里葉變換做測量，磁共振應用）
%??參考文獻：Amir?Beck?and?Marc?Teboulle?“A?Fast?Iterative?Shrinkage-Thresholding?Algorithm?for?
%??Linear?Inverse?Problems“?SIAM?J.?IMAGING?SCIENCES?Vol.?2?No.?1?pp.183-202?2009.
%??編程人：沙威?香港大學
%??Email:?wsha@eee.hku.hk;?dr.weisha@gmail.com
%??轉載時請保留上面注解

function?FISTA_FT;

clc;clear

%??圖像庫
A={‘phantom256.bmp‘?‘fruits256.bmp‘?‘cameraman256.bmp‘?‘lena256.bmp‘...
????‘peppers256.bmp‘?‘boat256.bmp‘?‘baboon256.bmp‘};

I=imread（A{3}）;???????????%??1表示最經典的phantom圖像，可以嘗試其它的圖像
X=double（I）/255;??????????%??原始矩陣歸一化
load?ww???????????????????%??小波變換矩陣
W_i=WL_for（wwX）;?????????%??小波變換
max_e=max（max（full（abs（W_i））））;?%??最大系數
threshold=0.005;??????????%??稀疏化閾值（小波變換）
th_v=（max_e*threshold）;???%??閾值

N=length（X）;??????????????%??原始圖像大小

%?空域隨機矩陣（Y=A*X測量，A病態，需要用共軛梯度技術，需要修正）
%?A=rand（N/4N）;??????????%??隨機測量矩陣（25%測量）
%?A=A/norm（A）;????????????%??測量矩陣歸一化

load?mask_radial;?????????%??mask矩陣（30%測量）
Y=FT_for（mask_matrixX）;??%??傅里葉變換系數作為測量

X_r=zeros（NN）;??%??恢復矩陣
N_iter=100;??????%??迭代次數
c=1;?????????????%??代理函數參數
gamma=1/c;???????%??迭代系數（可以優化）

t=1;?????????????%??自適應迭代系數

for?m=1:N_iter
????N_iter-m
????
????t_p=t;????????????????????????????%??記錄之前迭代系數
????X_p=X_r;??????????????????????????%??記錄之前的小波域系數

????X_t=WL_back（wwX_r）;??????????????%??小波反變換，得到時域圖像
????X_m=FT_for（mask_matrixX_t）;??????%??傅里葉變換域測量
????R_m=FT_back（mask_matrix（Y-X_m））;?%??測量殘差
????X_s=X_r+gamma*WL_for（wwR_m）;?????%??小波正變換，更新小波域系數
????X_r=Threshold_F（X_sth_v/c）;??????%??小波域系數軟閾值

????t=（1+sqrt（1+4*t^2））/2;????????????%??新的迭代系數
????X_r=X_r+（t_p-1）/t*（X_r-X_p）;??????%??更新小波域系數方法1
????%X_r=X_r+（t-1）/（t+1）*（X_r-X_p）;???%??更新小波域系數方法2
????
????th_v=th_v*0.9;????????????????????%??軟閾值松弛
????
????if?（th_v<1e-3）????????????????????%??閾值很小時候截斷迭代
????????break;
????end
????
????%??和原始圖像相比的峰值信噪比PSNR
????errorx=sum（sum（abs（WL_back（wwX_r）-X）.^2））;???%??MSE誤差
????psnr=10*log10（1*1/（errorx/N/N））;??%??PSNR
????disp（‘峰值信噪比：‘）
????disp（psnr）
?
end

figure（1）
subplot（221）
imshow（uint8（X*255））
title（‘原始圖像‘）
subplot（222）
imshow（uint8（real（Y）*255））
title（‘傅里葉域測量‘）
subplot（223）
imshow（uint8（full（WL_back（wwX_r））*255））
ss=sprintf（‘恢復圖像（PSNR:?%0.2f）‘psnr）;
title（ss）
subplot（224）
imshow（uint8（full（WL_back（wwX_r）-X）*255））
title（‘誤差圖像‘）

%??小波變換和反變換
function?MM=WL_for（wwM）
MM=（ww*（M）*ww‘）;?????%??2d?image
%MM=（ww*sparse（M））;??%??1d?signal

function?MM=WL_back（wwM）
MM=（ww‘*（M）*ww）;??????%??2d?image
%MM=（ww‘*sparse（M））;??%??1d?signal

%??閾值函數
function?X_th=Threshold_F（Xthreshold）

%??軟閾值
flag1=X>（threshold）;
flag2=X<-（threshold）;
flag3=flag1?|?flag2;
X_th=X-double（flag1）.*threshold;
X_th=X_th+double（flag2）.*threshold;
X_th=X_th.*flag3;

%?傅里葉變換和反變換
function?MM=FT_for（maskM）
M=real（M）;
MM=mask.*（fftshift（fft2（M）））;

function?MM=FT_back（maskM）
MM=real（ifft2（ifftshift（mask.*M）））;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件??????66614??2010-04-17?22:01??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\baboon256.bmp

?????文件??????66614??2010-04-17?22:02??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\boat256.bmp

?????文件??????66358??2010-04-17?22:22??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\brain256.bmp

?????文件??????66614??2010-04-17?22:07??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\cameraman256.bmp

?????文件???????3298??2013-12-14?22:23??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\FISTA_FT.m

?????文件???????3007??2013-12-14?23:27??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\FISTA_Radon.m

?????文件??????66614??2010-04-17?22:03??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\fruits256.bmp

?????文件??????66614??2010-04-16?17:30??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\lena256.bmp

?????文件???????6031??2013-12-08?20:17??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\mask_radial.mat

?????文件??????66614??2010-04-17?22:03??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\peppers256.bmp

?????文件???????3824??2010-04-17?22:21??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\phantom256.bmp

?????文件??????80168??2010-04-18?11:34??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA\ww.mat

?????目錄??????????0??2013-12-14?23:29??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）\FISTA

?????目錄??????????0??2014-03-03?10:16??FISTA?壓縮感知的圖像使用快速迭代shrinkage-thresholding算法（小波sparsifying轉換;傅氏域抽樣磁共振成像和氡轉換斷層測量）

-----------?---------??----------?-----??----

???????????????562370????????????????????14