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 /.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);