% grid parameters
d = 1; % grid in mm
fov = 160; % mm
N = fov/d;

% object parameters
sprad = 40;  % radius of sphere in mm
arrad = 10;   % radius of arrow shaft in mm
arlen = 120; % length of arrow in mm
headlen = 20; % length of arrowhead in mm
headang = 30; % angle of arrowhead cone in degrees

% rotflag=1 for object for lab project
% set rotflag=0 to display unrotated object (easier to see)
rotflag = 1;

% make rotated axes: magic angle from Z, 45° between X and Y
[X Y Z] = ndgrid(d*(-N/2:N/2-1)); % in mm
phi = rotflag*atan(sqrt(2)); % magic angle
cp=cos(phi); sp=sin(phi);
theta = rotflag*pi/4;
ct=cos(theta); st=sin(theta);
% 
rot = [cp 0 sp; 0 1 0; -sp 0 cp]*[ct st 0; -st ct 0; 0 0 1];
rotax = rot * [reshape(X,1,[]);reshape(Y,1,[]);reshape(Z,1,[])];
Xrot = reshape(rotax(1,:),size(X));
Yrot = reshape(rotax(2,:),size(Y));
Zrot = reshape(rotax(3,:),size(Z));

% make object
sphere = (X.^2 + Y.^2 + Z.^2) < sprad^2;  % no need to rotate sphere
shaft = (Xrot.^2 + Yrot.^2) < arrad^2 & abs(Zrot)<(arlen/2);
head = Zrot>=arlen/2 & Zrot<=(arlen/2+headlen) & (Xrot.^2+Yrot.^2) < (arlen/2+headlen-Zrot).^2;
arrow = max(shaft,head);
s = max(sphere,3*arrow); % arrow is twice as bright as ball

% smooth it
sig = 2; % sigma of gaussian in mm
fsize = ceil(2*sig/d);
[fx fy fz] = ndgrid(-fsize:fsize); % 5x5x5 pixel gaussian filter
h = exp(-(fx.^2+fy.^2+fz.^2)/(sig/d)^2);
ss = imfilter(s,h);

% uncomment this to look at FT of object
% log helps to visualize smaller components
%ss = log(abs(ft(ss)));

% draw it
figure; % every slice
imagesc(tile(permute(ss,[1 3 2])));
axis image;

figure; % 3 projections
subplot(131); imagesc(squeeze(sum(ss,1))); axis image;
subplot(132); imagesc(squeeze(sum(ss,2))); axis image;
subplot(133); imagesc(squeeze(sum(ss,3))); axis image;
