% Our own phantom defined in the square [0,1]^2.
% It consists of a few rectangles.
%
% Arguments:
% M     the final image will have size MxM
%
% Returns:
% res   pixel image, a discrete model of a nonnegative X-ray attenuation coefficient
%
% Samuli Siltanen Feb 2017

function res = SqPhantom(M)

% Initialize as zero
res = zeros(M,M);

% Construct MxM evaluation points inside the square [0,1]^2
pixelwidth = 1/M;
tmp   = pixelwidth/2 + [0:(M-1)]*pixelwidth;
[X,Y] = meshgrid(tmp);
Y     = flipud(Y);

% First rectangle
ind1 = (X>1/6)&(X<1/pi)&(Y>1/4)&(Y<.44);
res(ind1) = 1;

% Second rectangle
ind2 = (X>.65)&(X<.8)&(Y>.2)&(Y<.7);
res(ind2) = 1.2;

% Third rectangle
ind3 = (X>.11)&(X<.44)&(Y>.67)&(Y<.78);
res(ind3) = 0.7;


