function I = csphantom (n, bandwidth, signalToNoiseRatio, b1scale, img_real)
%%CSPHANTOM Test image for compressive sensing reconstructions.
%
% IMG = CSPHANTOM(N) returns an NxN phantom with various features designed
% to fail a Fourier-based compressed sensing reconstruction.
%
% IMG = CSPHANTOM(N,BW) returns an NxN phantom with sampling bandwidth
% equal to BW*N. This more accurately simulates what occurs when a
% physical phantom is sampled in an MRI scanner. The default value of BW
% is Inf. BW can be thought of as the intrinsic phantom pixel size divided
% by the Fourier measurement resolution. This parameter affects the
% sharpness of small features and the level of Gibbs ringing visible.
%
% Copyright 2012 David S. Smith <david.smith@vanderbilt.edu>.
%
% This file is part of CSPHANTOM.
if nargin < 1; n = 256; end % default size if none given
if nargin < 2; bandwidth = Inf; end
if nargin < 3; signalToNoiseRatio = Inf; end
if nargin < 4; b1scale = 1; end
if nargin < 5; img_real = true; end
if ~isinf(bandwidth)
n = n * bandwidth; % temporarily upsample the phantom creation
end
if n < 8, error('Smallest sensible phantom is 8 x 8.'); end
% quadrant sizes - slightly less than one fourth of the area, to give room
% for a border
quadrant_size = round((n - round(0.1*n))/2);
% Set up background shading: low amplitude sinc, a la B1 inhomogeneities
bglevel = 0.5; % central background gray level
[X, Y] = meshgrid(1:(2*quadrant_size), 1:(2*quadrant_size));
R = sqrt((X-quadrant_size).^2 + (Y-quadrant_size).^2);
BG = bglevel * cos(b1scale*R/(0.4*pi*quadrant_size));
% go back to quadrant coordinates
[X, Y] = meshgrid(1:quadrant_size, 1:quadrant_size);
% QUADRANT II: low-contrast circles
nc = 4; % # rows and columns of circles
Q2 = BG(1:quadrant_size,1:quadrant_size);
inc = quadrant_size / (nc+1);
for j = 1:nc
rad = quadrant_size / (11 + j^2);
for k = 1:nc
val = 0.9*bglevel - 0.1/k;
rctr = round(j*inc);
cctr = round(k*inc);
Z2 = (X-cctr).^2 + (Y-rctr).^2 <= rad^2;
Q2(Z2) = val;
end
end
% QUADRANT I: diagonal ramp
Q1 = BG(1:quadrant_size,(quadrant_size+1):(2*quadrant_size));
scale = quadrant_size/5;
ctr = round(quadrant_size/2);
Z = abs(X-ctr) + abs(Y-ctr);
Z1 = Z <= scale;
Z2 = Z > scale & Z <= 2*scale;
Q1(Z1) = Q1(Z1) * 2 .* (1 - Z(Z1) / scale);
Q1(Z2) = Q1(Z2) .* (Z(Z2) / scale - 1);
d = round(3*quadrant_size/16);
r = [d quadrant_size-d];
c = [d quadrant_size-d];
w = round(d/4);
for j = 1:2
for k = 1:2
Q1(r(j)-w:r(j)+w,c(k)-w:c(k)+w) = Q1(r(j)-w:r(j)+w,c(k)-w:c(k)+w)*(1+k*j/64);
end
end
% QUADRANT III: Gaussian bumps and quadratic hole
ctr = round(quadrant_size/2);
scale = quadrant_size / 4;
R2 = (X-ctr).^2 + (Y-ctr).^2;
Q3 = BG((quadrant_size+1):(2*quadrant_size),1:quadrant_size);
Q3withHole = Q3 .* R2 / scale^2;
Q3 = min(Q3,Q3withHole);
ctr = round([quadrant_size/5 4*quadrant_size/5]);
scale = scale / 6;
for r = 1:2
for c = 1:2
Q3 = Q3 + (3-r)*c*0.1*bglevel * exp(-((X-ctr(r)).^2+(Y-ctr(c)).^2)/scale^2);
end
end
% QUADRANT IV: line pairs
Q4 = BG((quadrant_size+1):(2*quadrant_size),(quadrant_size+1):(2*quadrant_size));
border_thick = round(quadrant_size/8); % define a border around the line pair region
line_length = max(1,border_thick);
Q4floor = 0.1;
%
% HORIZONTAL LINES
%
line_sep = 1;
line_thick = 1;
nlines = 3;
sep_incr = 1;
if ~isinf(bandwidth)
line_sep = line_sep*round(bandwidth);
line_thick = line_thick*round(bandwidth);
sep_incr = sep_incr*round(bandwidth);
end
col = 1;
row = 1 + round(0.5*border_thick) + line_length;
while (row + line_thick) <= quadrant_size - border_thick
for t = 1:nlines
Q4(row:row+line_thick-1,col:col+line_length-1) = Q4floor;
row = row + line_thick + line_sep;
if (row + line_thick) > quadrant_size - border_thick, break; end
end
row = row + line_sep + line_thick;
line_sep = line_sep + sep_incr;
line_thick = line_thick + sep_incr;
end
%
% VERTICAL LINES
%
line_sep = 1;
line_thick = 1;
sep_incr = 1;
if ~isinf(bandwidth)
line_sep = line_sep*round(bandwidth);
line_thick = line_thick*round(bandwidth);
sep_incr = sep_incr*round(bandwidth);
end
col = 1 + round(0.5*border_thick) + line_length;
row = 1;
while (col + line_thick) <= quadrant_size - border_thick % vertical lines
%Q4(row:row+line_length-1,col:col+line_thick-1) = Q4floor;
for t = 1:nlines
Q4(row:row+line_length-1,col:col+line_thick-1) = Q4floor;
col = col + line_thick + line_sep;
if (col + line_thick) > quadrant_size - border_thick, break; end
end
col = col + line_sep + line_thick;
line_sep = line_sep + sep_incr;
line_thick = line_thick + sep_incr;
end
%
% CIRCLES
%
nlines = 3;
ctr = round(quadrant_size/2); % concentric circles
R = sqrt((X-ctr).^2 + (Y-ctr).^2);
circ_thick = 1;
circ_gap = 3;
circ_incr = 1;
if ~isinf(bandwidth)
circ_thick = circ_thick*round(bandwidth);
circ_gap = circ_gap*round(bandwidth);
circ_incr = circ_incr*round(bandwidth);
end
r = circ_thick + circ_gap;
white = false;
while (ctr + r) < (quadrant_size - 1.5*border_thick)
for t = 1:nlines
rim = abs(R - r) <= 2*circ_thick;
if white
Q4(rim) = Q4(rim).*(1+Q4floor*exp(-2*(R(rim)-r).^2/circ_thick^2));
else
Q4(rim) = Q4(rim).*(1-exp(-2*(R(rim)-r).^2/circ_thick^2));
end
r = r + circ_thick + circ_gap;
if (ctr + r) > (quadrant_size - 1.2*border_thick), break; end
end
if white
white = false;
else
white = true;
end
r = r + circ_gap;
end
% pad area of signal with 'air'
npadf = (n - 2*quadrant_size) / 2;
npad = round(npadf);
I = zeros(n);
I(npad+1:npad+2*quadrant_size,npad+1:npad+2*quadrant_size) = [Q2 Q1; Q3 Q4];
% downsample to specified sampling bandwidth
if ~isinf(bandwidth)
F = fft2(I);
n = n / bandwidth; % truncate Fourier data
F = F([1:ceil(n/2) end-floor(n/2)+1:end],[1:ceil(n/2) end-floor(n/2)+1:end]);
% ringing filter
W = ifftshift(tukey(n)*tukey(n)');
F = F .* W; % apply window
I = abs(ifft2(F)) / bandwidth^2;
end
% finally add some Gaussian complex noise, if SNR .ne. Inf
if ~isinf(signalToNoiseRatio)
sigma = sqrt(0.125 / (signalToNoiseRatio^2 + 1));
I = I + sigma*randn(size(I)) + 1i*sigma*randn((size(I)));
end
if img_real, I = abs(I); end
function w = tukey (N)
%%WINDOW Tukey window
alpha = 0.8;
w = ones(N,1);
n = -(N-1)/2 : -alpha*N/2;
L = length(n);
w(1:L) = 0.5*(1+cos(pi*(abs(n)-alpha*N/2)/((1-alpha)*N/2)));
w(N : -1 : N-L+1) = w(1:L);
