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

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: [e-mail]
%

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.

Answer 1

[Walter Roberson]

imshow(M) perhaps

[sn at]

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.

[Lauren Hirt]

That definitely looks like an fftshift issue.