partial_fourier.m 2.21 KB
Newer Older
1
%% Load Data
2
% get matfile object
3
h = matfile('data.mat');
4
5

% get img variable from the mat-file
6
7
8
img = h.img;

%% 6/8 Partial Fourier sampling
9
% generate normally-distributed complex noise
10
n = randn(96) + 1j*randn(96);
11
12
13
14
15
16

% Fourier transform the image and add noise
y = fftshift(fft2(img),1) + n;

% set bottom 24/96 lines to 0
y(73:end,:) = 0;
17

18
19
20
21
% show sampling
figure();
imshow(log(abs(fftshift(y,2))), [], 'colormap', jet)

22
%% Estimate phase
23
24
25
26
27
28
29
30
% create zero-padded hanning filter for ky-filtering
filt = padarray(hann(48),24);

% generate low-res image with inverse Fourier transform
low = ifft2(ifftshift(y.*filt,1));

% get phase image
phs = exp(1j*angle(low));
31

32
33
34
35
36
37
38
39
40
41
% show phase estimate alongside true phase
figure();
subplot(1,2,1);
imshow(angle(img), [-pi,pi], 'colormap', hsv);
title('True image phase');

subplot(1,2,2);
imshow(angle(phs), [-pi,pi], 'colormap', hsv)
title('Estimated phase');

42
%% POCS reconstruction
43
% initialise image estimate to be zeros
44
est = zeros(96);
45
46

% set number of iterations
47
iters = 10;
48
49

% each iteration cycles between projections
50
for i = 1:iters
51
52
53
54
55
56
57
58
% projection onto data-consistent set:
    % use a-priori phase to get complex image
    est = est.*phs;

    % Fourier transform to get k-space
    est = fftshift(fft2(est), 1);

    % replace data with measured lines
59
    est(1:72,:) = y(1:72,:);
60
61
62
63
64
65
66
67
68
69
70
71
72

    % inverse Fourier transform to get back to image space
    est = ifft2(ifftshift(est, 1));

% projection onto non-negative reals
    % remove a-priori phase
    est = est.*conj(phs);

    % get real part
    est = real(est);

    % ensure output is non-negative
    est = max(est, 0);
73
74
end

75
76
77
78
79
80
81
82
83
84
85
86
87
%% Display error and plot reconstruction
% compute zero-filled recon
zf = ifft2(ifftshift(y, 1));

% compute rmse for zero-filled and POCS recon
err_zf = norm(zf(:) - img(:));
err_pocs = norm(est(:).*phs(:) - img(:));

% print errors
fprintf(1, 'RMSE for zero-filled recon: %f\n', err_zf);
fprintf(1, 'RMSE for POCS recon: %f\n', err_pocs);

% plot both recons side-by-side
88
figure();
89
90

% plot zero-filled
91
subplot(2,2,1);
92
imshow(abs(zf), [0 1]);
93
title('Zero-Filled');
94
95
subplot(2,2,3);
plot(abs(zf(:,48)), 'linewidth', 2);
96
97

% plot POCS
98
subplot(2,2,2);
99
100
imshow(est, [0 1]);
title('POCS recon');
101
102
subplot(2,2,4);
plot(abs(est(:,48)), 'linewidth', 2);