Matlab script to create a volume by rotating a 2D function of (r,z) about an axis at r=0

lathefunc.m

function V = lathefunc( F, rmax, zmax, outsize, z_mirror)
% V = lathefunc( F, rmax, zmax, out_size, zmirror )
%
% Create a volume by rotating a 2D function of (r,z) about an axis at r=0.
%
% F is a 2D grid dataset to be expanded
% Optional rmax and zmax is the function width and depth in given units,
%     default is size(F)
% Optional out_size defines an output volume size, default is to preserve
%     units and use [size(F,1)*2, size(F,2)]
% Optional z_mirror distance defines the distance between mirrored
%     applicators (use [-1] for default outsize)
%
% > lathefunc( peaks(20), 5, 5 )
%
% > lathefunc( peaks(20), 5, 5, -1, 10 )
%
% It is straightfoward to simply select and copy a large data set out of
% a spreadheet into a variable and use that as well.
% 
% Merck
% Fall, 2016

% if mod(size(F),2) == 1
%     F = F(:,1:end-1);
% end

if nargin < 5
    % No mirror
    z_mirror = -1;
end

if nargin < 4 || outsize(1) < 0
    outsize = [size(F,1)*2, size(F,1)*2, size(F,2)];
end

if nargin < 2
    rmax = size(F,1);
    zmax = size(F,2);
end

% Mirror and duplicate to make a full 2D cross section for reference
ref = [F(end:-1:1,:)', F'];

% Idea here is to create a grid of x,y,z values, convert them into polar
% coordinates and use that for a quick lookup into the original array.
x_ruler = linspace(-rmax,rmax,outsize(1) );
z_ruler = linspace(0,zmax,outsize(3) );

[cx,cy,cz] = meshgrid(x_ruler, x_ruler, z_ruler);
[th,r,z] = cart2pol( cx, cy, cz );

% Interpolate the polar coordinates out of the original coordinates of F
[fx,fy] = meshgrid(z_ruler, linspace(0,rmax,size(F,1)));
V = interp2(fx,fy,F,z,r);

% If there is a positive mirror_dist, assume dual applicators at that
% distance in z.
if z_mirror > 0
    % How big is a pixel in z
    z_spacing = z_ruler(2) - z_ruler(1);
    % How long does the new grid need to be in z
    new_zout = ceil(z_mirror/z_spacing)+1;
    old_zout = min(outsize(3), new_zout);

    Vnew = zeros(size(V,1), size(V,2), new_zout);
    Vnew(:,:,1:old_zout) = Vnew(:,:,1:old_zout) + V(:,:,1:old_zout);
    Vnew(:,:,end:-1:end-old_zout+1) = Vnew(:,:,end:-1:end-old_zout+1) + V(:,:,1:old_zout);
    V = Vnew;
    
    z_ruler = linspace(0,z_mirror,new_zout );
    [cx,cy,cz] = meshgrid(x_ruler, x_ruler, z_ruler);

end

% Print polarmap and reference
figure(1);
clf;
subplot(2,3,1);
imagesc( th(:,:,floor(end/2)) );
axis square;
title('Polar theta');
subplot(2,3,2);
imagesc( r(:,:,floor(end/2)) );
axis square;
title('Polar rho');
subplot(2,3,3);
imagesc( z(:,:,floor(end/2)) );
axis square;
imagesc(ref);
title('Reference cross-section');
axis equal;

% Print cross-sections of the volume
subplot(2,3,4);
imagesc( V(:,:,floor(end/2)) );
axis square;
title('Axial cross-section');
subplot(2,3,5);
imagesc( squeeze(V(:,floor(end/2),:))' );
title('Orthogonal cross-section');
axis equal;
subplot(2,3,6);
% imagesc( squeeze(out(end/2,:,:))' );
contourslice(cx, cy, cz, V, linspace(-rmax,rmax,5), [], linspace(0,zmax,5));
title('Volumetric contours');