function [x_out, y_out, z_out] = surf2solid(x,y,z, S, delta)
%
% file:      	surf2solid.m, (c) Matthew Roughan, Mon Feb 14 2011
% directory:   /home/mroughan/src/matlab/STL/
% created: 	Mon Feb 14 2011 
% author:  	Matthew Roughan 
% email:   	matthew.roughan@adelaide.edu.au
% 
% Convert a surface into a solid object for printing on a 3D printer
%
%   This routine only works for surfaces defined on a rectangular
%   grid.  The matrices X and Y define the axis limits only.
%
% The assumption is that 
%     -- the surface is a function z(x,y)
%     -- the input are given in triples (x,y, z(x,y)) similar to the matrix input option
%        of surf, mesh or other similar commands
%           i.e., x, and y specify regular, rectangular lattice points
%        but note that this doesn't support the 1, or 2 argument versions of these commands
%        and assumes that x and y are matrices, not vectors
%     -- x values increase along rows
%        y values along columns, e.g., as if they are the outputs of 
%          [x,y] = meshgrid( 1:3, 1:4 )
%
% The output will  be a solid with height given by z at each point (x,y), and
%     -- the surface z will be shifted by z_0, such that z+z_0 >= delta
%        where delta is the minimum z height
%     -- the (x,y) co-ordinates will be shifted so that (x+x_0, y_y_o) is "centered"
%     -- the x, y, and z values will be scaled so that the largest value of
%                (max(x) - min(x), (max(y)-min(y)), and 2*(max(z)-min(z))
%        is such that the largest of these dimensions is size S where S in measured in mm
%          (NB: matlab co-ordinates come out roughly as mm in the makerbot)
%     -- the surface z(x,y)+z_0 will have vertical walls going down to the (x,y) plane,
%        and the bottom surface of the object will be the part of the (x,y) plane including
%        the (x,y) points of the input (suitable rescaled)
%                  (e.g. see what meshz does)
%
% NB: fair bit of code is inspired by meshz
%

err_id = 'surf2solid:InvalidInput';

% check inputs
if (nargin < 3)
  error(err_id,'Need to input at least 3 arguments.')
end
if ischar(x) || ischar(y) || ischar(z)
  error(err_id,'Input should not be characters.')
end
[m,n] = size(z);
[mx,nx] = size(x);
[my,ny] = size(y);
if (m == 1 || n == 1)
  error(err_id,'Arrays must have size >1 in each direction.')
end
if (~isequal(size(z),size(x)) || ~isequal(size(z),size(y)))
  error(err_id,'Input arrays (x,y,z) should all be the same size.')
end
if (~all(isfinite(z)))
  error(err_id,'z must have finite values.')
end
if (~all(isfinite(y)))
  error(err_id,'y must have finite values.')
end
if (~all(isfinite(x)))
  error(err_id,'x must have finite values.')
end
if (nargin < 4)
  S = 40; % make the default max length along an axis 40mm
end
if (S>100)
  error(err_id,'S must be <= 100 due to size of build platform')
end
if (nargin < 5)
  delta = 3; % make the default minimum height 3 mm 
end
if (delta < 0)
  error(err_id,'delta must be >=0 to build')
end


% scale co-ordinates
x_diff = max(max(x)) - min(min(x));
y_diff = max(max(y)) - min(min(y));
z_diff = 2*(max(max(z)) - min(min(z)));
max_diff = max([x_diff; y_diff; z_diff]);
scale_factor = S/max_diff
x_s = scale_factor * x;
y_s = scale_factor * y;
z_s = scale_factor * z;
mx_x = max(max(x_s));
mn_x = min(min(x_s));
mx_y = max(max(y_s));
mn_y = min(min(y_s));
mx_z = max(max(z_s));
mn_z = min(min(z_s));

% shift the min z co-ordinate
z = z_s + delta - mn_z;

% center the (x,y) co-ordinates
x = x_s  - (mx_x + mn_x)/2;
y = y_s  - (mx_y + mn_y)/2;
xmin = min(min(x));
ymin = min(min(y));
xmax = max(max(x));
ymax = max(max(y));


% create an extended Z matrix
% z_out = [[   0         0    zeros(1,n) 0       0     ];
% 	 [   0         0    z(1,:)     0       0     ];
% 	 [zeros(m,1) z(:,1) z       z(:,n) zeros(m,1)];
% 	 [   0         0    z(m,:)     0       0     ];
% 	 [   0         0    zeros(1,n) 0       0     ];
% 	];
z_out = [[z(1,1) z(1,:)  z(1,n)];
	 [z(:,1) z       z(:,n)];
	 [z(m,1) z(m,:)  z(m,n)];
	 ];
z_out = zero_pad(z_out, 2);


% extend the x and y arrays
d = [1 1];
mm = [m m];
nn = [n n];
x_out = [[ x(d,d)  x(d,:)  x(d,nn) ];
	 [ x(:,d)  x       x(:,nn) ];
	 [ x(mm,d) x(mm,:) x(mm,nn) ];
	];
x_out = zero_pad(x_out, 1);
y_out = [[ y(d,d)  y(d,:)  y(d,nn) ];
	 [ y(:,d)  y       y(:,nn) ];
	 [ y(mm,d) y(mm,:) y(mm,nn) ];
	];
y_out = zero_pad(y_out, 1);


% mesh(x_out, y_out, z_out);




function B = zero_pad(A, k)
% 
% pad zeros around entire outside of array to width k
%
%%%  B = padarray( padarray(A, k)', k)';
  
[mA,nA] = size(A);
B = [[zeros(k,k)  zeros(k,nA) zeros(k,k)];
     [zeros(mA,k) A           zeros(mA,k)];
     [zeros(k,k)  zeros(k,nA) zeros(k,k)];
    ];
