function m = leaf_makesecondlayer( m, varargin )
%m = leaf_makesecondlayer( m, ... )
%   Make a new Bio-A layer, either adding to or discarding any existing one.
%
%   Options:
%       mode:   One of the following strings:
%           'full': cover the entire surface with a continuous sheet of
%               cells.  There will be one cell per FE cell plus one cell
%               per FE vertex.
%           'grid': cover the entire surface with a continuous sheet of
%               square cells.
%           'voronoi': cover the entire surface with a continuous sheet of
%               cells.  The 'numcells' option specifies how many.  This is
%               only valid for flat or nearly flat meshes.
%           'single': make a single second layer cell in a random position.
%           'few': make a specified number of second layer cells, each in a
%               different randomly chosen FEM cell.
%           'each': make one second layer cell within each FEM cell.
%       absdiam:   Not valid for mode=full.  In all other cases, this is a real
%               number, being the diameter of a single cell.
%       absarea:   Not valid for mode=full.  In all other cases, this is a real
%               number, being the area of a single cell.
%       reldiam:   Not valid for mode=full.  In all other cases, this is a real
%               number, being the diameter of a single cell as a proportion
%               of the diameter of the current mesh (estimated as the
%               square root of the mesh area).
%       relarea:   Not valid for mode=full.  In all other cases, this is a real
%               number, being the area of a single cell as a proportion
%               of the area of the current mesh.
%       relinitarea:   Not valid for mode=full.  In all other cases, this is a real
%               number, being the area of a single cell as a proportion
%               of the area of the initial mesh.
%       abssize:   OBSOLETE, REMOVED.  Use 'absdiam' instead.
%       relsize:   OBSOLETE, REMOVED.  Use 'reldiam' instead.
%       relinitarea:   OBSOLETE, REMOVED.
%       relFEsize:   OBSOLETE, REMOVED.
%       axisratio:   The ratio of major to minor axis of the cells.  If
%                    this is not 1, the major axis will be parallel to the
%                    polariser gradient, if any; if there is no polariser
%                    gradient, or at the cell location it is beneath the
%                    threshold for it to be  effective, the cell will be
%                    created circular.  The default value is 1.
%       numcells:    Valid only for mode='voronoi', mode='each', or
%               mode='few'.  An integer specifying the number of cells to
%               create.
%       sides:  Each cell is created as a regular polygon having this many
%               sides.  The default is 12, but if you are using a very
%               large number of cells, rendering will be significantly
%               faster if you use fewer sides.
%       refinement:  NOT USED.  IGNORE.
%       add:    Boolean.  If true, existing cells are retained and new
%               cells are added to them.  If false, any existing biological
%               layer is discarded.  If 'mode' is 'full' or 'voronoi', the
%               old layer is always discarded, ignoring the 'add' argument.
%       allowoverlap:  Boolean.  If true, cells are allowed to overlap each
%               other.  If false, cells will not overlap each other.  You
%               may as a result end up with fewer cells that you asked for.
%       colors:  The colour of the new cells, as an RGB value.
%       colorvariation:  The amount of variation in the colour of the new
%               cells. Each component of the colour value will be randomly
%               chosen within this ratio of the value set by the 'color'
%               argument.  That is, a value of 0.1 will set each component
%               to between 0.9 and 1.1 times the corresponding component of
%               the specified colour.  (The variation is actually done in
%               HSV rather than RGB space, but the difference is slight.)
%       probperFE:  For 'each' mode, this specifies for each FE the
%               probability per unit area for creating a cell there.
%               Any value less than 0 is equivalent
%               to 0.  The values will be scaled so as to add to 1.
%       probpervx:  Like probperFE, but the values are given per vertex.
%               Alternatively, the value can be a morphogen name or index,
%               in which case the values of that morphogen at each vertex
%               will be used.
%       The options 'probperFE', and 'probpervx' are mutually
%       exclusive.
%       vertexdata:  Instead of generating the cells according to one of
%               the methods described above, take this N*3 array to be the
%               positions of all of the cell vertexes.
%       celldata:  To be used with vertexdata.  This is a cell array of
%               integer arrays, each integer array specifying the vertexes
%               of one cell in order around the cell.
%       vertexindexing: This is either 1 or 0.  If celldata is provided, then
%               this specifies whether the indexes in celldata are 1-based
%               or 0-based.  1-based is the default, since that is Matlab's
%               convention for all indexing of arrays, but external data
%               may often be zero-based.  If vertexindexing is not supplied
%               or is empty, then if zero occurs anywhere in celldata, it
%               will be assumed to be zero-based, otherwise 1-based.
%       If either 'vertexdata' or 'celldata' is supplied, the other must
%       also be, and all other options that specify the manner of creating
%       cells will be ignored.
%
%   At most one of ABSDIAM, ABSAREA, RELDIAM, and RELAREA should be given.
%
%   Equivalent GUI operation: clicking the "Make cells" button on the Bio-A
%   panel.  This is equivalent to m = leaf_makesecondlayer(m,'mode','full').
%
%   Topics: Bio layer.

% The second layer contains the following information:
% For each cell ci:
%*       cells(ci).vxs(:)       A list of all its vertexes, in clockwise order.
%*       cells(ci).edges(:)     A list of all its edges, in clockwise order.
%           These cannot be 2D arrays, since different cells may have
%           different numbers of vertexes or edges.
%*       cellcolor(ci,1:3):     Its colour.
%       celllabel(ci,1):       Its label (an arbitrary integer). (OBSOLETE.)
%       celltargetarea(ci)     The cells' target areas.  Initially equal to
%                              their current areas.
%       cellarea(ci)           The cells' current areas.
%       areamultiple(ci)       A morphogen, initially 1.  The effective
%                              target area is areamultiple*celltargetarea.
%       cloneindex(ci)         An integer, used to distinguish clones.
%                              Inherited by descendants.
%
% For each clone vertex vi:
%*       vxFEMcell(vi)          Its FEM cell index.
%*       vxBaryCoords(vi,1:3)   Its FEM cell barycentric coordinates.
%*       cell3dcoords(vi,1:3)   Its 3D coordinates (which can be calculated
%                               from the other data).
% For each clone edge ei:
%*       edges(ei,1:4)          The indexes of the clone vertexes at its ends
%           and the clone cells on either side (the second one is 0 if absent).
%           This can be computed from the other data.
%        generation(ei)         An integer recording the generation at which
%                               this edge was created.
%
% Other parameters:
%        colorSaturation        ...
%        splitThreshold         ...
%        colorparams            ...
%        jiggleAmount           ...
%        averagetargetarea      ...
%   uniformhue: the hue of cells having cloneindex 0.
%   uniformhuerange: the random range of hue for cells having cloneindex 0.

    if isempty(m), return; end
    [s,ok] = safemakestruct( mfilename(), varargin );
    if ~ok, return; end
    s = defaultfields( s, 'mode', 'full', ...
        'absdiam', -1, 'absarea', -1, ...
        'reldiam', -1, 'relarea', -1, ...
        'relinitarea', -1, ...
        'numcells', 20, 'sides', 12, ...
        'refinement', 1, 'add', true, 'probperFE', [], 'probpervx', [], ...
        'axisratio', 1, ...
        'allowoverlap', true, ...
        'cellcolors', [], ...
        'colors', [], ...
        'colorvariation', [], ...
        'vertexdata', [], ...
        'celldata', [], ...
        'vertexindexing', [] );
    ok = checkcommandargs( mfilename(), s, 'only', ...
        'mode', ...
        'absdiam', 'absarea', ...
        'reldiam', 'relarea', ...
        'relinitarea', ...
        'numcells', 'sides', ...
        'refinement', 'add', 'probperFE', 'probpervx', ...
        'axisratio', ...
        'allowoverlap', ...
        'cellcolors', ...
        'colors', 'colorvariation', ...
        'vertexdata', 'celldata', 'vertexindexing' );
    if ~ok, return; end
    haveexternaldata = (~isempty( s.vertexdata )) || (~isempty( s.celldata ));
    if haveexternaldata
        if isempty( s.vertexdata )
            complain( '%s: celldata was supplied without vertexdata.', ...
                mfilename() );
            ok = false;
            return;
        elseif isempty( s.celldata )
            complain( '%s: vertexdata was supplied without celldata.', ...
                mfilename() );
            ok = false;
            return;
        end
        s = rmfield( s, { ...
                'absdiam', 'absarea', 'reldiam', 'relarea', ...
                'relinitarea', 'numcells', 'sides', 'refinement', ...
                'probperFE', 'probpervx', 'axisratio', 'allowoverlap' ...
            } );
        s.mode = 'external';
        if isempty( s.vertexindexing )
            for i=1:length(s.celldata)
                if any( s.celldata{i}==0 )
                    s.vertexindexing = 0;
                    break;
                end
            end
        end
        if s.vertexindexing==0
            for i=1:length(s.celldata)
                s.celldata{i} = s.celldata{i}+1;
            end
        end
    end
    
    if ~haveexternaldata
        if s.sides < 3
            s.sides = 3;
        end

        if s.absdiam > 0
            celldiameter = s.absdiam;
        elseif s.absarea > 0
            celldiameter = sqrt(4*s.absarea/pi);
        elseif s.reldiam > 0
            celldiameter = s.reldiam * sqrt(m.globalDynamicProps.currentArea);
        elseif s.relarea > 0
            cellarea = s.relarea * m.globalDynamicProps.currentArea;
            celldiameter = sqrt(cellarea*4/pi);
      % elseif s.abssize > 0
      %     celldiameter = s.abssize * 2;
      % elseif s.relsize > 0
      %     celldiameter = s.relsize * 2 * sqrt(m.globalDynamicProps.currentArea);
        elseif s.relinitarea > 0
            cellarea = s.relinitarea * m.globalProps.initialArea;
            celldiameter = sqrt(cellarea*4/pi);
        elseif (~strcmp( s.mode, 'full' ))&& (~strcmp( s.mode, 'voronoi' ))
          % celldiameter = s.relFEsize * 2 * sqrt(m.globalDynamicProps.cellscale);
            complain( '%s: no cell size found for mode "%s": use one of absdiam, absarea, reldiam, or relarea.', ...
                mfilename(), s.mode );
            ok = false;
            return;
        end
    end

    if isempty( s.colors ) && hasSecondLayer( m )
        s.colors = m.globalProps.colors;
    end
    if isempty( s.colorvariation ) && hasSecondLayer( m )
        s.colorvariation = m.globalProps.colorvariation;
    end
    if (~isempty( s.colors )) || (~isempty( s.colorvariation ))
        if ~isempty( s.colors )
            m.globalProps.colors = s.colors;
        end
        if ~isempty( s.colorvariation )
            m.globalProps.colorvariation = s.colorvariation;
        end
        m.globalProps.colorparams = ...
            makesecondlayercolorparams( m.globalProps.colors, ...
                                        m.globalProps.colorvariation );
    end

    if (~haveexternaldata) && ~isempty(s.probpervx)
        if ischar(s.probpervx) || (numel(s.probpervx)==1)
            mgenIndex = FindMorphogenIndex( m, s.probpervx );
            if isempty(mgenIndex)
                if ischar(s.probpervx)
                    fprintf( 1, '%s: morphogen %s not found.\n', s.probpervx );
                else
                    fprintf( 1, '%s: morphogen %d not found.\n', s.probpervx );
                end
                s.probpervx = [];
                return;
            else
                s.probpervx = m.morphogens(:,mgenIndex(1));
            end
        end
    end
    
    global findFEMCellStats
    findFEMCellStats(:) = 0;

    switch s.mode
        case 'external'
            m = makeExternCells( m, s.vertexdata, s.celldata, s.add );
            m = setSecondLayerColorInfo( m, s.colors, s.colorvariation );
            if (size(s.cellcolors,1)==length(s.celldata)) && (size(s.cellcolors,2)==3)
                m.secondlayer.cellcolor = s.cellcolors;
            else
                m.secondlayer.cellcolor = [ m.secondlayer.cellcolor; ...
                    randcolor( length(s.celldata), ...
                         m.globalProps.colorparams(1,[1 2 3]), ...
                         m.globalProps.colorparams(1,[4 5 6]) ) ];
            end
            m.secondlayer.cloneindex = [ m.secondlayer.cloneindex; ...
                                         zeros( length(s.celldata), 1 ) ];
            m.secondlayer.side = [ m.secondlayer.side; ...
                true( length(s.celldata), 1 ) ];
        case 'grid'
            m = makeCellGrid( m, celldiameter );
        case 'voronoi'
            m = makeVoronoiBioA( m, s.numcells, 8, [], s.colors, s.colorvariation, s.cellcolors );
            return;
        case 'full'
            m = makeFullSecondlayer2( m, s.refinement );
        case 'single'
            if (~s.add) || (~hasSecondlayer( m ))
                m.secondlayer = newemptysecondlayer();
            end
          % m = addRandomSecondLayerClump( m, cellsize*2 );
            if s.axisratio ~= 1
                m = calcPolGrad( m );
            end
            [m,ok] = addRandomSecondLayerCell( m, celldiameter, s.axisratio, s.sides, 2, [], [], s.allowoverlap, s.cellcolors );
            if ~ok
                fprintf( 1, '%s: new cell overlapped old cells, not created.\n', mfilename() );
            end
        case 'few'
            if (~s.add) || (~hasSecondlayer( m ))
                m.secondlayer = newemptysecondlayer();
            end
            if s.axisratio ~= 1
                m = calcPolGrad( m );
            end
            if s.numcells > 1
                cells = randperm( size( m.tricellvxs, 1 ) );
                m.secondlayer.cells = allocateCells( s.numcells );
                numcreated = 0;
                maxtries = ceil(s.numcells * 1.5);
                i = 1;
                while (i <= maxtries) && (numcreated < s.numcells)
                    [m,ok] = addRandomSecondLayerCell( m, ...
                                celldiameter, s.axisratio, s.sides, 2, cells(i), i, [], s.allowoverlap, s.cellcolors(numcreated+1,:) );
                    if ok
                        numcreated = numcreated+1;
                    end
                    i = i+1;
                end
                if numcreated < s.numcells
                    fprintf( 1, '%s: %d cells requested, only %d created due to overlaps.\n', ...
                        mfilename(), s.numcells, numcreated );
                end
            else
                [m,ok] = addRandomSecondLayerCell( m, celldiameter, s.axisratio, s.sides, 2, [], [], s.allowoverlap, s.cellcolors );
                if ~ok
                    fprintf( 1, '%s: new cell overlapped old cells, not created.\n', mfilename() );
                end
            end
            checkclonesvalid( m.secondlayer );
        case 'each'
            if (~s.add) || (~hasSecondLayer( m ))
                m.secondlayer = newemptysecondlayer();
            end
            if ~isempty(s.probpervx)
                posprob = max( 0, s.probpervx );
                s.probperFE = sum( posprob(m.tricellvxs), 2 )/size(m.tricellvxs,2);
                if all(s.probperFE==0)
                    s.probperFE(:) = 1;
                end
            end
            if ~isempty(s.probperFE)
                s.probperFE = s.probperFE .* m.cellareas;
                totprob = sum(s.probperFE);
                if totprob > 0
                    cis = genrand( makedistribution(s.probperFE/totprob), s.numcells );
                else
                    cis = [];
                end
            else
                cis = genrand( makedistribution(m.cellareas/sum(m.cellareas)), s.numcells );
            end
            if s.axisratio ~= 1
                m = calcPolGrad( m );
            end
            numcreated = 0;
            for i=1:length(cis)
                if i > size(s.cellcolors,1)
                    cellcolor = [];
                else
                    cellcolor = s.cellcolors(i,:);
                end
                [m,ok] = addRandomSecondLayerCell( m, celldiameter, s.axisratio, s.sides, 1, cis(i), [], s.allowoverlap, cellcolor );
                if ok
                    numcreated = numcreated+1;
                end
            end
            if numcreated < s.numcells
                fprintf( 1, '%s: %d cells requested, only %d created due to overlaps.\n', ...
                    mfilename(), s.numcells, numcreated );
            end
            checkclonesvalid( m.secondlayer );
        otherwise
            fprintf( 1, '%s: unknown mode "%s".\n', mfilename(), s.mode );
            return;
    end
    m.secondlayer.jiggleAmount = 0; % 0.02;
    if ~isfield( m.secondlayer, 'side' ) || isempty( m.secondlayer.side )
        m.secondlayer.side = true( length(m.secondlayer.cells), 1 );
    end
    
    m = calcBioACellAreas( m );
    numcells = length(m.secondlayer.cellarea);
    m.secondlayer.areamultiple = ones(1,numcells);
    if numcells == 0
        m.secondlayer.celltargetarea = [];
    else
        m.secondlayer.celltargetarea = m.secondlayer.areamultiple * ...
            (sum(m.secondlayer.cellarea)/numcells);
    end
    m.secondlayer.averagetargetarea = ...
        sum(m.secondlayer.celltargetarea)/length(m.secondlayer.celltargetarea);
    m.secondlayer.generation = zeros( size(m.secondlayer.edges,1), 1 );
end

    
