function [ POINTS ] = ras2sld(Xc, UPSAMPLE, FILL, OUTPUTDIR)
% [ POINTS ] = ras2sld(Xc, UPSAMPLE, FILL, OUTPUTDIR)
%
% Interpolate closed curves to create a three-dimensional solid.
%
% Xc {1xNFILES} - a cell array of data matrices, as returned by otl2ras.
%                 Size of Xc{i} is at least [3xNDAT(i)].  
%                 If Xc{i} has more than 3 rows, extras are ignored.
%
% UPSAMPLE (scalar)   - factor by which X,Y,Z should be upsampled (default: 2)
%
% FILL (scalar) - if FILL==1, the object is "filled in," i.e. ras2sld returns coordinates of
%                 every point inside the object.  If FILL==0, FILL==[], or FILL is unspecified,
%                 ras2sld returns only the coordinates of points on the surface of the object.
%
% TYPE (string) - type of interpolation
%
% POINTS [3xprod(NPTS)] - a list of points in 3D which form the surface
%
% Optional Inputs:
% OUTPUTDIR (string)  - if specified, distance images are written to this directory
%                       (if OUTPUTDIR is '', distance images are not written out)
%
% Xc is treated as a stack of outlines.  Between each set of outlines,
% new outlines are created using Raya and Udupa's "shape-based
% interpolation" algorithm (IEEE Trans. on Medical Imaging, 9(1):32-42, 1990).
% The distance-map calculation is performed using the routine otl2dmap.m.
%
% The volume of a voxel in the interpolated image is roughly (1/UPSAMPLE)^3.
% X and Y are upsampled by UPSAMPLE before passing to otl2dmap, and interpolation
% is created at intervals of roughly (1/UPSAMPLE) along the Z axis.
%
% The resulting solid object is converted to a surface by finding all
% voxel faces with a negative distance on one side and a non-negative distance
% on the other side.  The (X,Y,Z) coordinates of the centerpoint of each such
% face are returned in the array POINTS.
%
% If OUTPUTDIR is specified, distance images will also be stored in
% OUTPUTDIR using the function dmapwr.m, in files of the form <OUTPUTDIR>/<Z>.mr.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1:
%  Check the arguments, and open various files and directories
%

% Define UPSAMPLE
if nargin<2, UPSAMPLE=2; end
if nargin<3, FILL=0; end
if length(FILL)==0, FILL=0; end
if nargin<4, OUTPUTDIR=[]; end

% Initialize the output POINTS array to have the correct number of rows
POINTS = zeros(3,0);

% Find out whether or not we should write out DMAP images
write_dmaps = 0;
if length(OUTPUTDIR)>0,
  if length(dir(OUTPUTDIR)) > 0, 
    write_dmaps = 1; 
  else,
    disp(sprintf('interp3d will not write DMAPs to a nonexistant directory: %s',OUTPUTDIR));
  end;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2:
% Figure out which of the cardinal directions is Z (between image slices),
% and which directions are X and Y (within each image slice).
% First, find the sum and sum-of-squares matrices:
SS = zeros(3,3);
for i=1:length(Xc),
  [ foo, N ] = size(Xc{i});
  if N>1,
    SS = SS + (N-1)*cov(Xc{i}');
  end
end

% ras2sld assumes that Z is one of the cardinal directions, so
% we can just assign Z to be the direction with minimum within-group variance.
[ minvar, zdir ] = min(diag(SS));
xydir = [1 2 3];
xydir = xydir(xydir ~= zdir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3:
%   Work our way up through Z; for each pair of outlines,
%   create all interpolations between iplane-1 and iplane.
%
for iplane=2:length(Xc),
  
  % Check to make sure that both of the maps have non-zero outline length
  if all( [size(Xc{iplane-1}) size(Xc{iplane})]>0 ),
    STARTX = Xc{iplane-1}(xydir,:);
    ENDX = Xc{iplane}(xydir,:);

    % Create the STARTMAP and ENDMAP
    % The mapsize is determined by the range of outline points in STARTX and ENDX
    maxX = ceil(max([STARTX(1:2,:)'; ENDX(1:2,:)'])) + 1;
    minX = floor(min([STARTX(1:2,:)'; ENDX(1:2,:)'])) - 1;
    mapsize = [ maxX(2)-minX(2), maxX(1)-minX(1) ];
    STARTMAP = otl2dmap(UPSAMPLE*(STARTX(1,:)-minX(1)), UPSAMPLE*(STARTX(2,:)-minX(2)), UPSAMPLE*mapsize, UPSAMPLE*2 );
    ENDMAP = otl2dmap(UPSAMPLE*(ENDX(1,:)-minX(1)), UPSAMPLE*(ENDX(2,:)-minX(2)), UPSAMPLE*mapsize, UPSAMPLE*2 );
    
    % Lowest value of Z, and spacing between planes
    STARTZ = mean(Xc{iplane-1}(zdir,:)')
    ENDZ = mean(Xc{iplane}(zdir,:)');
    SPACING = abs(ENDZ - STARTZ);

    % OLDMAP starts out equal to the STARTMAP
    OLDMAP = STARTMAP;
    oldZ = STARTZ;
    
    % How many interpolations from iplane-1 to iplane?
    % Figure out where to put the interpolated planes.
    % Include STARTMAP only when iplane==2 (first time through).
    HOWMANYPLANES = ceil(UPSAMPLE*SPACING);
    fracs = [ 0:HOWMANYPLANES ] / HOWMANYPLANES;
    if iplane>2 & all(size(Xc{iplane-2})>0), 
      fracs = fracs(2:length(fracs));
    end
    
    % Create each interpolated layer
    for frac=fracs,
      % Create the interpolated image, and save it if desired
      Z = (1-frac) * STARTZ + frac * ENDZ;
      INTERMAP = (1-frac) * STARTMAP + frac * ENDMAP;
      if write_dmaps, ras2sld_dmapwr(INTERMAP, Z, OUTPUTDIR, minX, maxX);  end;
      
      % Find all zero-crossings between OLDMAP and INTERMAP, and add them to POINTS
      if frac > 0 & FILL == 0,
	[Y, X] = find(xor(INTERMAP<0, OLDMAP<0));
	if length(X)>0,
	  newpts = zeros(3,length(X));
	  newpts(xydir,:) = [ X'/UPSAMPLE + minX(1); Y'/UPSAMPLE + minX(2)];
	  newpts(zdir,:) = (Z+oldZ)/2;
	  POINTS = [ POINTS newpts ];
	end
      end
      
      % If FILL==0, find all zero-crossings (all edges) inside INTERMAP
      % If FILL==1, find all places where INTERMAP > 0 (points inside the object)
      if FILL > 0,
	[ Y, X ] = find(INTERMAP > 0);
	
	% FILL==1: Dither the boundary
	% Find places where INTERMAP==0, and choose half of them at random
	% to append on the end of Y, X.
	[ Y0, X0 ] = find(INTERMAP == 0);
	which = find(rand(size(Y0)) > 0.5);
	Y = [ Y; Y0(which) ];
	X = [ X; X0(which) ];
	
      else,
	[ X, Y ] = zerocrossings(INTERMAP);
      end
      
      % Add points to the output array
      if length(X)>0,
	newpts = zeros(3,length(X));
	newpts(xydir,:) = [ X'/UPSAMPLE + minX(1); Y'/UPSAMPLE + minX(2)];
	newpts(zdir,:) = Z;
	POINTS = [ POINTS newpts ];
      end

      % Create OLDMAP and oldZ for the next loop
      OLDMAP = INTERMAP;  
      oldZ = Z;
    end;
  end
end  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   That's all, folks!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Internal function zerocrossings
% Find all edges within a plane where the intensity changes from below zero to above zero
function [X, Y] = zerocrossings(Im)

[maxY,maxX] = size(Im);

% BG is one where Im >0, BL is Im<0, B0 is Im==0
BG = (Im > 0);
BL = (Im < 0);
B0 = (B==0);

% find all points at which Im==0
[Y0,X0] = find( B0 );

% The rest of this function finds any zero-crossing points which may
% not have been marked with a zero pixel in the DMAP.  I'm not even sure
% if this code is necessary (it may be that DMAP marks all zero-crossings 
% with a zero-pixel), but this code is here to make sure.

% find all points such that Im(Y1,X1) > 0, but Im(Y1,X1+1) < 0
% .. the outline point is at Y1, X1+0.5
[Y1,X1] = find( BG & [BL(:, 2:maxX) zeros(maxY,1)]);
X1 = X1+0.5;

% find all points such that Im(Y2,X2) > 0, but Im(Y2, X2-1) < 0
% .. the outline point is at Y2, X2-0.5
[Y2,X2] = find( BG & [zeros(maxY,1) BL(:, 1:(maxX-1))]);
X2 = X2-0.5;

% find all points such that Im(Y3,X3) > 0, but Im(Y3+1,X3) < 0
% .. the outline point is at Y3+0.5, X3
[Y3,X3] = find( BG & [BL(2:maxY, :); zeros(1,maxX)]);
Y3 = Y3+0.5;		 

% find all points such that Im(Y3,X3) > 0, but Im(Y3-1,X3) < 0
% .. the outline point is at Y4-0.5, X4
[Y4,X4] = find( BG & [zeros(1,maxX); BL(1:(maxY-1), :)] );
Y4 = Y4-0.5;

% Sort together all of these points, and remove duplications
XY = unique([ X0 Y0; X1 Y1; X2 Y2; X3 Y3; X4 Y4 ], 'rows' );
if prod(size(XY))==0, XY=zeros(0,2); end
X = XY(:,1);
Y = XY(:,2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Internal function ras2sld_dmapwr
% Write out specified DMAP
function ras2sld_dmapwr(DMAP, Z, OUTPUTDIR, minX, maxX);

gehead = struct('mri', struct('corners', zeros(4,3), 'user', zeros(1,27)));
gehead.mri.corners(:,1:2) = [ minX([2 1]); maxX(2) minX(1); minX(2) maxX(1);  maxX([2 1]) ];
gehead.mri.corners(:,3) = repmat(Z, [4 1]);

gefilename = fullfile(OUTPUTDIR, sprintf('%3.3g.mr', Z));
bytes = dmapwr(DMAP, gehead, gefilename);

