% function [ SURF, POINTS ] = ras2surf(RAS, PPMM, TYPE, DTHRESH)
%
% Interpolate open curves to obtain uniformly-sampled points on a surface in 3D space.
%
% RAS {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.
% PPMM (scalar) - the average number of pixels per millimeter we try to return.
% TYPE (string) - the type of interpolation to use between planes -- may be
% any of the types understood by imresamp. Default is 'linear'
% DHTHRESH (scalar) - distance threshold for multiline concatenation
%
% POINTS [3 x prod(N1,N2)] - a list of points in 3D which form the surface
% SURF [N2 x N1 x 3] - a matrix of the points, organized as the surface is organized
%
% This function operates in the following stages
% 1. Each curve in Xc is partitioned into multiline objects, using ras2ml with
% distance threshold DTHRESH. The longest multiline in each plane is kept,
% while any shorter multilines are discarded.
% 2. The remaining multilines (one per Xc element) are resampled so that they each
% have N1 points, where N1 is chosen so that the longest multiline object will
% have PPMM points per millimeter. Resampling is done using imresamp, with type 'sinc'.
% 3. The resampled multilines are treated as N1 different curves. These N1 curves
% are resampled using spline interpolation, so that each one has N2 points,
% where N2 is chosen so that the longest such curve has about PPMM points per millimeter.
%
% In the future, it might be nice to use splines to resample contours with grossly
% non-equidistant sample points. For now, spline interpolation is not implemented;
% the extra code would be mostly wasted, since spline interp is so much slower than
% sinc interpolation.
%
function [ SURF, POINTS ] = ras2surf(RAS, PPMM, TYPE, DTHRESH)
% Default value of arguments
if nargin<2, PPMM=1; end
if nargin<3, TYPE='linear'; end
if nargin<4, DTHRESH=2; end
% Sort Xc into multi-line objects
[ MLINES, LENGTHS ] = ras2ml(RAS, DTHRESH);
% Initialize variables for storing important information
SS = [0 0 0];
npXc = [];
lXc = zeros(1,length(MLINES));
Xc = {};
nXc = 0;
% Go through, picking out the longest multiline (total Euclidean distance) from each non-empty plane
for iML=1:length(MLINES),
if all(size(MLINES{iML}) > 0),
% Increment nXc
nXc = nXc + 1;
% Pick out the longest multiline in this slice
[ maxl, maxi ] = max(LENGTHS{iML});
Xc{nXc} = MLINES{iML}{maxi};
lXc(nXc) = maxl;
% Figure out how many points are in each of the MLINES objects
npts = zeros(1, length(MLINES{iML}));
for j=1:length(MLINES{iML}), npts(j) = size(MLINES{iML}{j}, 2); end
npXc(nXc) = npts(maxi);
% If any points were thrown away, warn the user
if sum(npts) > npXc(nXc),
disp(sprintf('ras2surf: kept only %d contiguous points out of %d total in slice number %d.', ...
npXc(nXc), sum(npts), iML));
end
% Calculate SS, the average within-group sum-of-squares
SS = SS + (npXc(nXc) - 1) * var(Xc{nXc}');
end
end
% Sort the mlines from - to + in the Maximum Within-Group Variance Dimension (MWGVD)
[ maxSS, MWGVD ] = max(SS);
for iXc=1:nXc,
ave1 = mean(Xc{iXc}(:,1:floor(npXc(iXc)/2))');
ave2 = mean(Xc{iXc}(:,(floor(npXc(iXc)/2)+1):npXc(iXc))');
if ave1(MWGVD) > ave2(MWGVD),
Xc{iXc} = fliplr(Xc{iXc});
end
end
% N1 is the number of points required such that the longest multiline has PPMM points/mm
N1 = ceil(max(lXc) * PPMM);
ROWS = zeros([nXc, N1, 3]);
for ifile=1:nXc,
% Endpoints of resampled contour have same abscissa as endpoints of original
LX = (N1-1)/(size(Xc{ifile},2)-1);
% Resample the contour, and store in ROWS
for dim=1:3,
ROWS(ifile, :, dim) = imresamp(Xc{ifile}(dim,:), 1, LX, 'sinc');
end
end
% Calculate the cumulative distance from beginning of each curve to each point
cdist = [ zeros(1, N1); cumsum( sqrt( sum( diff(ROWS).^2, 3))) ];
% N2 is the number of points required such that the longest curve has PPMM points/mm
N2 = ceil( max(cdist(nXc,:)) * PPMM );
SURF = zeros([N2, N1, 3]);
if nargout>1, POINTS = zeros([3, N2*N1]); end
if strcmp(TYPE,'spline'),
% Spline interpolate each of the 'curves', using cdist as an abscissa
for idim=1:3,
for ic=1:N1,
SURF(:,ic,idim) = spline( cdist(:,ic), ROWS(:,ic,idim), (cdist(nXc,ic)/(N2-1))*[0:(N2-1)]' );
if nargout > 1,
POINTS(idim,((ic-1)*N2+1):(ic*N2)) = SURF(:,ic,idim)';
end
end
end
else
% Calculate the resampling ratio in the Y dimension
LY = (N2-1)/(nXc-1);
% Resample the 3 ROWS matrices to get the SURF matrix
for idim=1:3,
surf = imresamp(ROWS(:,:,idim), LY, 1, TYPE);
SURF(:,:,idim) = surf;
POINTS(idim,:) = surf(:)';
end
end