tempestpr.m

傅立叶采样数据生成程序

clear all;
x=load('d:\zhang\sparseMRI\datamat.dat');
figure(80)
imshow(x,[0 255]);
im = x>100;
TVWeight = 0.05; 	% Weight for TV penalty
xfmWeight = 0.0;	% Weight for Transform L1 penalty
Itnlim = 15;
XFM = 1;
% generate variable density random sampling
DN = size(x); 	% size of k-space
N = DN;
pctg = [1/15];
P = 1000;
pdf = genPDF(DN,P,pctg , 2 ,0,0);
mask = genSampling(pdf,10,2);
data = x;
%generate Fourier sampling operator and data
phmask = zpad(hamming(6)*hamming(6)',N(1),N(2)); %mask to grab center frequency
phmask = phmask/max(phmask(:));			 %for low-order phase estimation and correction
ph = exp(i*angle((ifft2c(data.*phmask)))); % estimate phase for phase correction

%generate Fourier sampling operator
FT = p2DFT(mask, N, ph, 2);

% scale data
im_dc = FT'*(data.*mask./pdf);
data = data/max(abs(im_dc(:)));
im_dc = im_dc/max(abs(im_dc(:)));

% initialize Parameters for reconstruction
param = init;
param.FT = FT;
param.XFM = XFM;
param.TV = TVOP;
param.data = data;
param.TVWeight =TVWeight;     % TV penalty 
param.xfmWeight = xfmWeight;  % L1 wavelet penalty
param.Itnlim = Itnlim;

res = data;
% do iterations
for n=1:5
	res = fnlCg(res,param);
	im_res = XFM'*res;
	figure(100), imshow(abs(im_res),[0 255]), drawnow
end