I want to plot phantom image's sampling pattern given here.
1 view (last 30 days)
Show older comments
this is the image:

and this is the code:
path(path, './Optimization');
path(path, './Measurements');
path(path, './Data');
% Phantom
n = 256;
N = n*n;
X = phantom(n);
x = X(:);
% number of radial lines in the Fourier domain
L = 22;
% Fourier samples we are given
[M,Mh,mh,mhi] = LineMask(L,n);
OMEGA = mhi;
A = @(z) A_fhp(z, OMEGA);
At = @(z) At_fhp(z, OMEGA, n);
% measurements
y = A(x);
path(path, './Optimization');
path(path, './Measurements');
path(path, './Data');
% Phantom
n = 256;
N = n*n;
X = phantom(n);
x = X(:);
% number of radial lines in the Fourier domain
L = 22;
% Fourier samples we are given
[M,Mh,mh,mhi] = LineMask(L,n);
OMEGA = mhi;
A = @(z) A_fhp(z, OMEGA);
At = @(z) At_fhp(z, OMEGA, n);
% measurements
y = A(x);
% min l2 reconstruction (backprojection)
xbp = At(y);
Xbp = reshape(xbp,n,n);
% recovery
tic
tvI = sum(sum(sqrt([diff(X,1,2) zeros(n,1)].^2 + [diff(X,1,1); zeros(1,n)].^2 )));
fprintf('Original TV = %8.3f\n', tvI);
xp = tveq_logbarrier(xbp, A, At, y, 1e-1, 2, 1e-8, 600);
Xtv = reshape(xp, n, n);
toc
% A_fhp.m
%
% Takes measurements in the upper half-plane of the 2D Fourier transform.
%
% Usage: b = A_fhp(x, OMEGA)
%
% x - N vector
%
% b - K vector = [mean; real part(OMEGA); imag part(OMEGA)]
%
% OMEGA - K/2-1 vector denoting which Fourier coefficients to use
% (the real and imag parts of each freq are kept).
%
% Written by: Justin Romberg, Caltech
% Created: October 2005
% Email: jrom@acm.caltech.edu
%
function [M,Mh,mi,mhi] = LineMask(L,N)
thc = linspace(0, pi-pi/L, L);
%thc = linspace(pi/(2*L), pi-pi/(2*L), L);
M = zeros(N);
% full mask
for ll = 1:L
if ((thc(ll) <= pi/4) | (thc(ll) > 3*pi/4))
yr = round(tan(thc(ll))*(-N/2+1:N/2-1))+N/2+1;
for nn = 1:N-1
M(yr(nn),nn+1) = 1;
end
else
xc = round(cot(thc(ll))*(-N/2+1:N/2-1))+N/2+1;
for nn = 1:N-1
M(nn+1,xc(nn)) = 1;
end
end
end
% upper half plane mask (not including origin)
Mh = zeros(N);
Mh = M;
Mh(N/2+2:N,:) = 0;
Mh(N/2+1,N/2+1:N) = 0;
M = ifftshift(M);
mi = find(M);
Mh = ifftshift(Mh);
mhi = find(Mh);
end
It seems we must plot the output of "[M,Mh,mi,mhi] = LineMask(L,N)" function.
0 Comments
Accepted Answer
More Answers (0)
See Also
Categories
Find more on Image Processing Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!