% This example script doesn't do anything interesting
% it just shows you how to use the functions

im=double(dicomread('i5633321.MRDC.8'));
[nx ny] = size(im);

figure(278);
subplot(221);
imagesc(im);

% Uncomment appropriate line to get the desired undersampling pattern
U = zeros(nx,ny);
%U(:,2:2:end)=1;                         % every other line
U(:,ny/2+vdgrid(1024,512,300,10))=1;     % 1D Gaussian
%U=vdgrid2(1024,1024*512,400,10);        % 2D Gaussian
u = find(U);

k = ift(im);

imhat = ft(U.*k);
subplot(222);
imagesc(abs(imhat));
subplot(223);
imagesc(abs(imhat)-im);

lev = 2;  % number of levels of recursive wavelet decompositions
s = waveletcdf97(imhat,lev);

subplot(224);
imagesc(abs(s));

% How to enforce data consistency
khat = ift(imhat);
khat(u) = k(u);
imhat = abs(ft(khat)); %abs() effectively incorporates your knowledge that
                       % the image phase is zero as well