I want to plot phantom image's sampling pattern given here.

1 view (last 30 days)
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.

Accepted Answer

Walter Roberson
Walter Roberson on 7 Feb 2017
imshow(M) perhaps
  2 Comments
sn at
sn at on 7 Feb 2017
Edited: Walter Roberson on 3 Jun 2018
I tried
>> [M,Mh,mi,mhi] = Line2Mask(22,256);
>> imshow(M)
and I got the following figure:
Not the same.
and for
>> [M,Mh,mi,mhi] = Line2Mask(22,256);
>> OMEGA=mhi;
>> imshow(mhi)
Warning: Image is too big to fit on screen; displaying at 25%
no picture is shown.

Sign in to comment.

More Answers (0)

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!