tubeSurf( data, r, noCPoints, mapData );

function [X, Y, Z, C] = tubeSurf( data, r, noCPoints, mapData );
%
% Purpose : to generate a tubular surface given a sequence of points along
% a 3d curve.
%
% Use (syntax): 
%   [X, Y, Z] = tubeSurf( data, r );
%   returns the surface with a default number of points along the
%   circumfrence.
%   or  
%   [X, Y, Z, C] = tubeSurf( data, r, noCPoints, mapData );
%   also returns a colormap suitable for use with the SURF command. It maps
%   the mapData linearly along the surface
%
% Input Parameters :
%   data:       3D data points along the curve as a matrix of column vectors
%   r:          the radius of the tube
%   noCPoints:  numbe of point along the circumfrene of the circle forming
%               the tube. Optional, default = 20
%   mapData:    Data which will be mapped on the the color map matrix, it
%               must be the same length as the numbe of points supplied
%
% Return Parameters :
%   X, Y, Z:    The surface coordinates as required for the SURF command
%   C:          A colormap matrix which corresponds to mapplin the map data
%               linearly along the tube.
%
% Requires:
%
% Description and algorithms:
%
% References :
% 
% Author :  Paul O'Leary
% Date :    22. May 2005
% Version : 1.0
%---------------------------------------------------------------------------
% (c) 2005, O'Leary, Harker, University of Leoben, Leoben, Austria
% email: [e-mail], url: automation.unileoben.ac.at
%---------------------------------------------------------------------------
% History:
%   Date:           Comment:
%   22. May 2005    Original Version 1.0
%--------------------------------------------------------------------------

% Test the correctness of the input data
%
[n, noPoints] = size( data );
if ~( (n==4) | (n==3) )
    error('The data must be a matrix of homogeneous column vectors, i.e. 4 x n');
end;
%
if (nargin ==4)&(~(length(mapData)==noPoints))
    error('The mapData must be the same length as the number of points');
end;   
%
% setup default values
%
if nargin < 3
    noCPoints = 20;
    mapData = ones(1,noPoints);
end;
%
if nargin < 4
    mapData = ones(1,noPoints);
end;
%
epsilon = 1e-10;
%
% normalize the data for color mapping
%
if nargin == 4
    minMap = min( mapData );
    maxMap = max( mapData );
    scale = maxMap - minMap;
    if scale == 0
        error('Invalid map data, it has zero scale');
    end;
    mapData = 2*(mapData - minMap)/scale - 1;
end;
%
% Pick up the data and calculate the gradients
%
x = data(1,:);
y = data(2,:);
z = data(3,:);
%
dx = gradient( x );
dy = gradient( y );
dz = gradient( z );
%
% Cacculate the lengths of the vernicular segments
%
Ls = sqrt( dx.^2 + dy.^2 + dz.^2);
%
% Eliminate the zero length elements
%
indices = find( abs(Ls) > epsilon );
dx = dx(indices);
dy = dy(indices);
dz = dz(indices);
Ls = Ls(indices);
%
% Give a worning if zero elements are found.
%
if length(indices) ~= length(Ls)
    warning('zero length segments have been ignored');
end;
%
% Calculate the required angles of rotation.
%
sinY = - dz./ Ls;
cosY = sqrt(1 - sinY.^2);
%
phiZ = atan2( dy, dx);
sinZ = sin( phiZ );
cosZ = cos( phiZ );
%
% Generate the circle on the x z plane
%
phiC = linspace(0,2*pi,noCPoints);
xc = zeros(size(phiC));
yc = r*cos( phiC );
zc = r*sin( phiC );
cc = ones(size(xc));
circ = [xc;yc;zc];
%
% Generate the surface
%
X = [];
Y = [];
Z = [];
IC = [];
for k=1:noPoints
    Ry = [cosY(k),0,sinY(k); 0,1,0; -sinY(k),0,cosY(k)];
    Rz = [cosZ(k),-sinZ(k),0; sinZ(k),cosZ(k),0; 0,0,1];
    %
    R = Rz * Ry;
    %
    temp =  R*circ;
    %
    X = [X; temp(1,:)+x(k)];
    Y = [Y; temp(2,:)+y(k)];
    Z = [Z; temp(3,:)+z(k)];
    IC = [IC; mapData(k)*cc];
end;
%
% return the color mapped matrix if required.
%
if nargout == 4
    C = IC;
end;