% NxN test target for tomographic reconstruction
function [object, x, y] = target1(N)

x = linspace(-1,1,N)' * ones(1,N);
y = ones(1,N)' * linspace(-1,1,N);
r = sqrt(x.^2 + y.^2);

object = zeros(N,N);

% Hyperbolic paraboloid confined to a circle.
object += 500*(r<0.8).*abs((x-0.2).*(y+0.1)+0.1);

% Smooth the image...
object = smooth2(object,32);

% Solid cardioid
object += 100*circshift( r < (1+x./r)./2, [-ceil(N/4), ceil(N/16)] );

% A little more smoothing...
object = smooth2(object,3);

% Hollow ellipse
ellipses = sqrt(1.4*(x-0.25).^2 + (y+0.2).^2);
object += 130*((ellipses<0.4) & (ellipses>0.35));

% Linear feature
object += 200*(-0.2>(x+2*y) & (x+2*y)>-0.3 ...
   & x<0.6 & y<-0.01);

% elliptical shading
object .*= (ellipses' ./ (1+(5*ellipses').^2));

% A little more smoothing, to minimize negatives on rotation...
object = smooth2(object,2);

% rotation
object = imrotate(object, 10, 'bicubic', 'crop');
object(find(object<0))=0; % threshold >0
object( find(isna(object)) )=0; % replace NA with 0.0

colormap(gray(256))
imagesc(x,y,object)
title('Test Target (target1.m)')
xlabel('y')
ylabel('x')
axis('equal')
