function [m,ok] = meshFromCells( vvxs, vcells, centroids )
%m = meshFromCells( vvxs, vcells )
%   Create a mesh from a set of biological cells.
%
%   vvxs is an N*3 array of N vertexes, and vcells is a C*1 cell array
%   listing the vertexes of each cell.  centroids, if supplied, is a C*3
%   array giving the position of a point within each cell.  If not
%   supplied, it defaults to the average of all the vertexes belonging to
%   that cell.

    m = [];
    numedgevxs = size(vvxs,1);
    numcells = length(vcells);
    if (nargin < 3) || isempty(centroids)
        centroids = zeros( numcells, size(vvxs,2) );
        for i=1:numcells
            centroids(i,:) = sum( vvxs( vcells{i}, : ), 1 )/length(vcells{i});
        end
    end
    numstaredges = 0;
    for i=1:numcells
        numstaredges = numstaredges + length(vcells{i});
    end
    m.nodes = [ vvxs; centroids ];
    
    m.tricellvxs = zeros( numstaredges, 3 );
    m.edgeends = zeros( numstaredges*2, 2 );
    numFEs = 0;
    cellnumsegments = zeros( numcells, 1 );
    for i=1:numcells
        vvi = vcells{i}(:);
        cellnumsegments(i) = length(vvi);
        cellFEs = [ repmat( numedgevxs+i, length(vvi), 1 ), vvi(:), vvi([2:end,1]') ];
        m.tricellvxs( (numFEs+1):(numFEs+length(vvi)), : ) = cellFEs;
        numFEs = numFEs+length(vvi);
    end
    [m,ok] = setmeshfromnodes( m, [] );
    
    m.fourthlayer.iscellcentre = [ false( numedgevxs, 1 ); true( numcells, 1 ) ];
    m.fourthlayer.cellcentres = find(m.fourthlayer.iscellcentre);

    m.fourthlayer.iscellwall = ~any( m.fourthlayer.iscellcentre( m.edgeends ), 2 );
    m.fourthlayer.cellwallindexes = find(m.fourthlayer.iscellwall);

    numwalls = length(m.fourthlayer.cellwallindexes);
    m.fourthlayer.cellproperties = zeros( numcells, 0 );
    m.fourthlayer.wallproperties = zeros( numwalls, 0 );
    
    % To build the edge vertex properties, we have to, for each edge,
    % divide it into segments of equal length (approximately the same
    % segment length for all edges), and represent the properties of each
    % segment as three values per segment-end, excluding the farther end of
    % the final segment.
    
    wallends = reshape( m.nodes( m.edgeends( m.fourthlayer.cellwallindexes, : )', : ), ...,
                        2, [], 3 );
    walllengths = sqrt( sum( (wallends( 2, :, : ) - wallends( 1, :, : )).^2, 3 ) )';
    maxwl = max( walllengths );
    minwl = min( walllengths );
    MAXWALLSEGMENTS = 10;
    MINWALLSEGMENTS = 3;
    r = max( min( maxwl/minwl, MAXWALLSEGMENTS ), MINWALLSEGMENTS );
    m.fourthlayer.targetwallseglength = maxwl/r;
    
    
    m.fourthlayer.edgevertexproperties = cell( numwalls, 1 );
    for i=1:numwalls
        numsegs = max( 1, round( walllengths(i)/m.fourthlayer.targetwallseglength ) );
        m.fourthlayer.edgevertexproperties{i} = zeros( 3, numsegs, 0 );
    end
end
