function m = splitSecondLayerCells( m, splitall )
%m = splitSecondLayerCells( m )
%   Split all second layer cells satisfying certain criteria.
%   If SPLITALL is true, all cells will be split.  Otherwise, the choice depends
%   on cell size.

    if ~hasSecondLayer( m )
        return;
    end
    setGlobals();
    global gOLD_ARREST;
    global gNEW_ARREST;
    if m.versioninfo.mgenversion==0
        ARREST_MGEN = gOLD_ARREST;
    else
        ARREST_MGEN = gNEW_ARREST;
    end

    numsecondcells = length( m.secondlayer.cells );
    if numsecondcells==0, return; end

    if m.globalProps.maxBioAcells > 0
        maxsplits = m.globalProps.maxBioAcells - numsecondcells;
        if  maxsplits <= 0
            return;
        end
    else
        maxsplits = -1;
    end

    if nargin < 2, splitall = 0; end

    numsplits = 0;
    cellsToSplit = zeros(1,numsecondcells);

    USEMINDIAM = 1; % 0;
    USEPOLGRAD = 0; % 1;
    MINEIGRATIO = 1;

    arrestThreshold = 0.99;
    haveArrest = max( m.morphogens( :, ARREST_MGEN ) ) >= arrestThreshold;
    splitcentres = zeros(0,3);
    splitdirs = zeros(0,3);

    for ci=1:numsecondcells
        secondlayer_vxs = m.secondlayer.cells(ci).vxs;
        cellcoords = m.secondlayer.cell3dcoords( secondlayer_vxs, : );
        numcellvxs = size( cellcoords,1 );
        
        % Don't split a cell if the arrest morphogen is high.
        % The level of arrest morphogen over the whole cell is estimated by
        % its average over the vertexes of the cell.
        if haveArrest
            femCells = m.secondlayer.vxFEMcell( secondlayer_vxs );
            femBCs = m.secondlayer.vxBaryCoords( secondlayer_vxs, : );
            femvxs = m.tricellvxs( femCells, : );
            arrestPerFEVx = reshape( m.morphogens( femvxs', ARREST_MGEN ), 3, [] )';
            arrestPerSecondVx = sum( arrestPerFEVx .* femBCs, 2 );
            secondCellArrest = sum( arrestPerSecondVx, 1 ) / length( secondlayer_vxs );
            if secondCellArrest >= arrestThreshold
                continue;
            end
        end
        
        dosplit = 1;
        
        if isa( m.globalProps.bioApresplitproc, 'function_handle' )
            [m,dosplit,centroid,v] = m.globalProps.bioApresplitproc( m, ci );
            splitall = 1;
        elseif USEMINDIAM
            [centroid,v,distance] = polysplitdiameter( cellcoords(:,[1 2]) );
            centroid = [ (centroid + randInCircle2( distance/4 )), 0 ];
            v = [ -v(2), v(1), 0 ];
        elseif USEPOLGRAD
            centroid = sum( cellcoords, 1 )/numcellvxs;
            pg = m.gradpolgrowth( ...
                    m.secondlayer.vxFEMcell( m.secondlayer.cells(ci).vxs(:) ), ...
                    : );
            v = sum(pg,1)'/size(pg,1);
            n = norm(v);
            covariance = cov( cellcoords );
            [eigs,e] = eig( covariance );
            de = diag(e);
            es = sort(de);
            r = es(2)/es(3);
            if (n <= 0) || (r < MINEIGRATIO)
                cellaxis = eigs(:,3);
                va = vecangle( v', cellaxis' );
                theta = abs(abs(va) - pi/2);
                if theta < pi/4
                    v = makeperp( v, cellaxis );
                	v = v/norm(v);
                elseif n > 0
                    v = v/n;
                else
                    v = cellaxis;
                end
            else
                v = v/n;
            end
        elseif 1
            centroid = sum( cellcoords, 1 )/numcellvxs;
            covariance = cov( cellcoords );
            [vv,e] = eig( covariance );
            v = vv(:,3);
        else
            maxd = -1;
            for i=1:numcellvxs-1
                for j=i+1:numcellvxs
                    vv = cellcoords(i,:) - cellcoords(j,:);
                    d = dotproc2(vv,vv);
                    if maxd < d
                        maxd = d;
                        v = vv';
                        maxi = i;
                        maxj = j;
                    end
                end
            end
            centroid = (cellcoords(maxi,:) + cellcoords(maxj,:)) * 0.5;
            v = v/norm(v);
        end
        if splitall
            % Nothing.
        elseif 1
            a = polyarea3( cellcoords );
          % splitSecondLayerCells_splitThreshold = m.secondlayer.splitThreshold^2
          % dosplit = a > m.secondlayer.splitThreshold^2;
            femcells = m.secondlayer.vxFEMcell( m.secondlayer.cells(ci).vxs );
            fembcs = m.secondlayer.vxBaryCoords( m.secondlayer.cells(ci).vxs, : );
            mgens = reshape( 1 + m.morphogens( m.tricellvxs( femcells, : ), ARREST_MGEN ), 3, [] );
            vxmgens = dot( mgens', fembcs, 2 );
            cellarrest = sum(vxmgens)/numel(vxmgens);
            if 1
                targetarea = m.secondlayer.averagetargetarea;
            elseif 1
                targetarea = m.secondlayer.areamultiple(ci) * ...
                             m.secondlayer.celltargetarea(ci);
            else
                targetarea = m.secondlayer.areamultiple(ci) * ...
                         (m.morphogens( ...
                            m.tricellvxs( ...
                                m.secondlayer.vxFEMcell( ...
                                    m.secondlayer.cells(ci).vxs(1) ), 1 ) ...
                            , ARREST_MGEN ) + 1) * ...
                         m.secondlayer.celltargetarea(ci);
            end
            dosplit = a * cellarrest > targetarea * 1.414;
        else
            p = cellcoords * v;
            celllength = max(p)-min(p);
            dosplit = celllength > m.secondlayer.splitThreshold;
        end
        if dosplit
            numsplits = numsplits+1;
            cellsToSplit(numsplits) = ci;
            v = v + 0.0*rand(size(v));
            splitdirs(numsplits,:) = v';
            splitcentres(numsplits,:) = centroid;
            if (maxsplits > 0) && (numsplits >= maxsplits)
                break;
            end
        end
    end
    cellsToSplit = cellsToSplit(1:numsplits);
    if numsplits > 0
        fprintf( 1, 'splitSecondLayerCells: %d cells to split\n', numsplits );
    end
    for i=1:numsplits
        m = splitclonecell( m, cellsToSplit(i), splitdirs(i,:), splitcentres(i,:) );
    end
%    if numsplits > 0
%        if numsplits==1, past = ' was'; else past = 's were'; end
%        fprintf( 1, '%d cell%s split in the second layer.  Total %d.\n', ...
%            numsplits, past, length( m.secondlayer.cells ) );
%    end
    ok = checkclonesvalid( m.secondlayer );
    if ~ok
        fprintf( 1, 'Invalid second layer in %s.\n', mfilename() );
        error( mfilename() );
    end
end
