function m = leaf_subdivide( m, varargin )
%m = leaf_subdivide( m, ... )
%   Subdivide every edge of m where a specified morphogen is above and/or below
%   thresholds, and the length of the current edge is at least a certain
%   value.
%
%   NB. This function is currently NOT suitable for calling from an
%   interaction function.  It will go wrong.
%
%   Note that this command will subdivide every eligible edge every time it
%   is called.  It does not remember which edges it has subdivided before
%   and refrain from subdividing them again.
%
%   Options:
%       'morphogen': The name or index of the morphogen
%       'min':       The value that the morphogen must be at least equal to.
%       'max':       The value that the morphogen must not exceed.
%       'mode':      'all' [default], 'any', or 'mid'.
%       'minabslength': A real number.  No edge shorter than this will be
%                    subdivided.
%       'minrellength': A real number.  This is a fraction of the current
%                    threshold for automatic splitting of edges.  No edge
%                    shorter than this will be subdivided.  The current
%                    threshold value is returned by currentEdgeThreshold(m).
%       'levels':    The levels of morphogen the new vertices will adopt
%                    'all','interp','none'
%       'force':     If true (the default), all of the requested
%                    subdivisions will happen.  If false, edges will not be
%                    split where this would result in excessively small
%                    cells or cell angles.
%
%   An edge will be subdivided if and only if it satisfies all of the
%   conditions that are specified.  Any combination of the arguments can be
%   given.  No arguments gives no subdivision.
%
%   'mode' is only relevant if 'min' or 'max' has been specified.
%   If mode is 'all', then each edge is split for which both ends satisfy
%   the min/max conditions.
%   If mode is 'any', each edge is split for which either edge
%   satisfies the conditions.
%   If mode if 'mid', each edge is split for which the average of the
%   morphogen values at its ends satisfies the conditions.
%
%   This command ignores the setting, that can be set through the GUI or
%   leaf_setproperty(), that enables or disables automatic splitting of
%   long edges.
%
%   Topics: Mesh editing.

    if isempty(m), return; end
    [s,ok] = safemakestruct( mfilename(), varargin );
    if ~ok, return; end
    s = defaultfields( s, ...
            'levels', 'all', ...
            'mode', 'all', ...
            'force', true );
    ok = checkcommandargs( mfilename(), s, 'only', ...
            'levels', ...
            'morphogen', ...
            'min', ...
            'max', ...
            'mode', ...
            'minabslength', ...
            'minrellength', ...
            'force' );
    if ~ok, return; end
    
    if isfield( s, 'morphogen' )
        mgenIndex = FindMorphogenIndex( m, s.morphogen, mfilename() );
        if isempty(mgenIndex)
            % FindMorphogenIndex has already complained.
            return;
        end
    else
        mgenIndex = 0;
    end
    
    if isfield(s, 'levels')
        leveltype=s.levels;
    else
        leveltype='interp';
    end

    if isfield( s, 'morphogen' )
        [kmorph_i,kmorph_p,kmorph_a,kmorph_l] = getMgenLevels( m, s.morphogen );
        current_max_morph_value = max(kmorph_p);
    end
    current_vertices = m.nodes;
    
    edgesToSplit = 1:size(m.edgeends,1);
    dosplit = false;

    if mgenIndex ~= 0
        dosplit = true;
        mgenEdgeLevels = reshape( m.morphogens( m.edgeends, mgenIndex ), size(m.edgeends) );
        switch s.mode
            case 'all'
                if isfield( s, 'min' )
                    ind_min = find( mgenEdgeLevels >= s.min );
%                     edgesToSplit = find( all( mgenEdgeLevels(edgesToSplit,:) >= s.min, 2 ) );
                end
                if isfield( s, 'max' )
                    ind_max = find( mgenEdgeLevels <= s.max );
%                     edgesToSplit = find( all( mgenEdgeLevels(edgesToSplit,:) <= s.max, 2 ) );
                end
                edgesToSplit=intersect(ind_min,ind_max);
            case 'any'
                if isfield( s, 'min' )
                    ind_min = find( mgenEdgeLevels >= s.min );
%                     edgesToSplit = find( any( mgenEdgeLevels(edgesToSplit,:) >= s.min, 2 ) );
                end
                if isfield( s, 'max' )
                    ind_max = find( mgenEdgeLevels <= s.max );
%                     edgesToSplit = find( any( mgenEdgeLevels(edgesToSplit,:) <= s.max, 2 ) );
                end
                edgesToSplit=[ind_min,ind_max];
            case 'mid'
                mgenEdgeLevels = sum(mgenEdgeLevels,2)/2;
                if isfield( s, 'min' )
                    ind_min = find( mgenEdgeLevels >= s.min );
%                     edgesToSplit = find( any( mgenEdgeLevels(edgesToSplit,:) >= s.min, 2 ) );
                end
                if isfield( s, 'max' )
                    ind_max = find( mgenEdgeLevels <= s.max );
%                     edgesToSplit = find( any( mgenEdgeLevels(edgesToSplit,:) <= s.max, 2 ) );
                end
                edgesToSplit=intersect(ind_min,ind_max);
            case 'levels'
                mgenEdgeLevels = sum(mgenEdgeLevels,2)/2;
                if isfield( s, 'min' )
                    ind_min = find( mgenEdgeLevels >= s.min );
%                     edgesToSplit = find( any( mgenEdgeLevels(edgesToSplit,:) >= s.min, 2 ) );
                end
                
                    [knor_i,knor_p,knor_a,knor_l] = getMgenLevels( m, 'KNOR' );

            otherwise
                fprintf( 1, '%s: Option ''mode'' has invalid value ''%s''.\n', ...
                    mfilename(), s.mode );
                fprintf( 1, '    Allowed values are ''min'', ''max'', or ''mid''.\n' );
                return;
        end
    end

    if isfield(s, 'minabslength') || isfield(s, 'minrellength')
        dosplit = true;
        edgevecs = m.nodes( m.edgeends( edgesToSplit, 2 ), : ) ...
                   - m.nodes( m.edgeends( edgesToSplit, 1 ), : );
        edgelensqs = sum( edgevecs.*edgevecs, 2 );
        if isfield( s, 'minabslength' )
            edgesToSplit = edgesToSplit( edgelensqs >= s.minabslength*s.minabslength );
        end
        if isfield( s, 'minrellength' )
            splitthreshold = s.minrellength * currentEdgeThreshold( m );
            edgesToSplit = edgesToSplit( edgelensqs >= splitthreshold*splitthreshold );
        end
    end
    if dosplit
        m = splitalledges( m, edgesToSplit, s.force );
    end
    
    if isfield( s, 'morphogen' )
        if strcmp(leveltype,'all')
            [temp,ind] = setdiff(m.nodes,current_vertices,'rows');
            kmorph_p(ind) = current_max_morph_value;
            m.morphogens(:,kmorph_i) = kmorph_p; %
        elseif strcmp(leveltype,'none')
            [temp,ind] = setdiff(m.nodes,current_vertices,'rows');
            kmorph_p(ind) = 0;
            m.morphogens(:,kmorph_i) = kmorph_p; %
        end % default is interpolation
    end
end

