function m = leaf_plot( m, varargin )
%m = leaf_plot( m, ... )
%   Plot the leaf.
%   There are many options.
%
% leaf_plot stores all of the plotting options in the mesh, so that a
% subsequent call to leaf_plot2 with only the mesh as an argument will plot
% the same thing.
%
%   Equivalent GUI operation: none.  The leaf is plotted automatically.
%   Various options may be set in the "Plot options" panel, and the scroll
%   bars on the picture change the orientation.
%
%   See also: leaf_plotoptions.
%
%   Topics: Plotting.

    % Sanity check.
    if isempty(m), return; end
    
    [m,plotinfo] = leaf_plotoptions( m, varargin{:} );
    
    s = m.plotdefaults;
    
    % Eliminate invalid picture handles.
    m.pictures = m.pictures( ishandle( m.pictures ) );
    
    % If there is no picture handle, make one.
    if isempty( m.pictures )
        v = struct( 'Visible', boolchar( s.invisibleplot, 'off', 'on' ) );
        m = leaf_addpicture( m, ...
                'properties', v, ...
                'uicontrols', s.uicontrols );
        if isempty( m.plotdefaults.axisRange )
            h = guidata( m.pictures(1) );
            setaxis( h.picture, boundingbox( m, 0.1 ) );
        else
            setaxis( m.pictures(1), m.plotdefaults.axisRange );
        end
        setViewFromMesh( m );
    end
    
    % Get the axes for the first picture handle.
    pichandles = guidata( m.pictures(1) );
    theaxes = pichandles.picture;
    cla( theaxes );
    hold( theaxes, 'on' );
    if m.plotdefaults.light
        lightAxes( theaxes, true );
    end

    % We only have a colour bar when we are plotting a single quantity.
    % This is the case if all of the following hold:
    % (i) cmaptype is nonempty.
    % (ii) The colorbar object exists.
    % (iii) We are plotting exactly one scalar value -- that is, we are not
    % plotting different things on both sides or multiple morphogens, and
    % we are plotting something.
    wantColorBar = true;
    if isempty( s.cmaptype )
        % No colour map ==> no colorbar.
        wantColorBar = false;
    end
    if plotinfo.haveAplot && plotinfo.haveBplot
        % Plotting something separately on each side ==> no colorbar.
        wantColorBar = false;
    end
    havemultipleplot = iscell(s.morphogen) && (numel( s.morphogen ) > 1);
    havemultipleAplot = iscell(s.morphogenA) && (numel( s.morphogenA ) > 1);
    havemultipleBplot = iscell(s.morphogenB) && (numel( s.morphogenB ) > 1);
    if havemultipleplot || havemultipleAplot || havemultipleBplot
        % Plotting multiple morphogens ==> no colorbar.
        wantColorBar = false;
    end

    % Calculate the data to be plotted as colours.
    m.plotdata = struct();
    m = calculatePlotData( m, '' );
    m = calculatePlotData( m, 'A' );
    m = calculatePlotData( m, 'B' );
    numscalars = 0;
    for fnc = { '', 'A', 'B' }
        fn_value = ['value' fnc{1}];
        fn_mgen = ['morphogen' fnc{1}];
        fn_output = ['outputquantity' fnc{1}];
        if isfield( m.plotdata, fn_value ) && ~isempty( m.plotdata.(fn_value) ) && (size( m.plotdata.(fn_value), 2 )==1)
            numscalars = numscalars+1;
            if s.autoColorRange || isempty( s.crange ) || (s.crange(1) >= s.crange(2))
                crange = [ min(0,min(m.plotdata.(fn_value))), max(0,max(m.plotdata.(fn_value))) ];
            else
                crange = s.crange;
            end
            monocolors = []; % [1 0 0];
            if strcmp( s.cmaptype, 'monochrome' )
                if ~isempty( s.monocolors )
                    monocolors = s.monocolors;
                elseif ~isempty( s.(fn_mgen) )
                    mgens = FindMorphogenIndex( m, s.(fn_mgen) );
                    if length(mgens)==1
                        monocolors = [ m.mgenposcolors( :, mgens )'; ...
                                       m.mgennegcolors( :, mgens )' ];
                    end
                elseif ~isempty( s.(fn_output) )
                    o = regexprep( s.(fn_output), 'actual', 'resultant' );  % Hack.
                    o = regexprep( o, 'rate$', '' );
                    monocolors = m.outputcolors.(o);
                else
                    monocolors = [0 0 1; 1 0 0];
                end
                if size( monocolors, 1 )==1
                    monocolors = [ monocolors; oppositeColor( monocolors ) ];
                end
            end
            [s.cmap,s.crange] = chooseColorMap( s.cmaptype, crange, monocolors );
            m.plotdefaults.cmap = s.cmap;
            m.plotdefaults.crange = s.crange;
            if ~wantColorBar
                m.plotdata.(fn_value) = translateToColors( m.plotdata.(fn_value), crange, s.cmap );
            end
        end
    end
    if numscalars ~= 1
        wantColorBar = false;
    end
    if wantColorBar
        colormap( pichandles.picture, s.cmap );
        caxis( pichandles.picture, s.crange );
    end

    commonpatchargs = struct( 'Parent', theaxes, ...
                              'FaceLighting', s.lightmode, ...
                              'AmbientStrength', s.ambientstrength, ...
                              'FaceAlpha', s.alpha, ...
                              'EdgeColor', s.FElinecolor );
    m = findVisiblePart( m );
    vistriangles = m.tricellvxs( m.visible.cells, : );
    visnodeindexes = find(m.visible.nodes);
    renumberNodes = zeros( size(m.nodes,1), 1 );
    renumberNodes(visnodeindexes) = (1:length(visnodeindexes))';
    vistriangles = renumberNodes(vistriangles);

    % Draw the mesh
    if s.drawleaf && ~isempty(vistriangles)
        viscellindexes = find(m.visible.cells);
        if ~s.thick
            % Mid-plane
            if s.drawedges < 2
                thickness = 0;
            else
                thickness = m.plotdefaults.FEthinlinesize;
            end
            suffix = '';
            for fnc = { '', 'A', 'B' }
                fn_value = ['value' fnc{1}];
                if isfield( m.plotdata, fn_value ) && ~isempty( m.plotdata.(fn_value) )
                    suffix = fnc{1};
                    break;
                end
            end
            vxs = m.nodes( m.visible.nodes, : );
            m = makeCanvasColor( m, '' );
            plotmeshsurface( theaxes, vxs, vistriangles, ...
                    m.plotdata.(['value' suffix]), ...
                    m.plotdata.(['pervertex' suffix]), ...
                    m.visible, thickness, commonpatchargs, ...
                    struct( 'faceindexes', viscellindexes, 'face', 'M' ) );
            if s.drawedges==1
                % Draw border edges.
                plotedges( m, theaxes, '', ...
                    m.visible.borderedges, ...
                    m.plotdefaults.FEthinlinesize, ...
                    m.plotdefaults.FElinecolor );
                if s.drawseams
                    plotedges( m, theaxes, '', ...
                        m.visible.borderedges & m.seams, ...
                        m.plotdefaults.seamlinesize, ...
                        m.plotdefaults.seamlinecolor );
                end
            end
            if ~isempty( m.selection.highlightedVxs )
                plotpts( theaxes, m.nodes( m.selection.highlightedVxs, : ), ...
                    '.b', 'MarkerSize', m.plotdefaults.highlightthickness*10 );
            end
        else
            if s.drawedges < 2
                thicknessA = 0;
                thicknessB = 0;
            elseif s.decorateAside % m.globalProps.dorsaltop
                thicknessA = m.plotdefaults.FEthicklinesize;
                thicknessB = m.plotdefaults.FEthinlinesize;
            else
                thicknessA = m.plotdefaults.FEthinlinesize;
                thicknessB = m.plotdefaults.FEthicklinesize;
            end
            if isfield( m.plotdata, 'valueA' ) && ~isempty( m.plotdata.valueA )
                Asuffix = 'A';
            else
                Asuffix = '';
            end
            if isfield( m.plotdata, 'valueB' ) && ~isempty( m.plotdata.valueB )
                Bsuffix = 'B';
            else
                Bsuffix = '';
            end
            if ~isfield( m.plotdata, ['value' Asuffix] )
                m = makeCanvasColor( m, 'A' );
                Asuffix = 'A';
            end
            if ~isfield( m.plotdata, ['value' Bsuffix] )
                m = makeCanvasColor( m, 'B' );
                Bsuffix = 'B';
            end
            visprismindexes = visnodeindexes*2;

            % B side
            vxs = m.prismnodes( visprismindexes, : );
            plotmeshsurface( theaxes, vxs, vistriangles, ...
                    m.plotdata.(['value' Bsuffix]), ...
                    m.plotdata.(['pervertex' Bsuffix]), ...
                    m.visible, thicknessB, commonpatchargs, ...
                    struct( 'faceindexes', viscellindexes, 'face', 'B' ) );
            % A side
            vxs = m.prismnodes( visprismindexes-1, : );
            plotmeshsurface( theaxes, vxs, vistriangles, ...
                    m.plotdata.(['value' Asuffix]), ...
                    m.plotdata.(['pervertex' Asuffix]), ...
                    m.visible, thicknessA, commonpatchargs, ...
                    struct( 'faceindexes', viscellindexes, 'face', 'A' ) );
            % Borders
            plotborder( m, theaxes );
            
            % Seams
            if s.drawseams
                plotedges( m, theaxes, 'A', ...
                    m.visible.edges & m.seams, ...
                    m.plotdefaults.seamlinesize, ...
                    m.plotdefaults.seamlinecolor );
                plotedges( m, theaxes, 'B', ...
                    m.visible.edges & m.seams, ...
                    m.plotdefaults.seamlinesize, ...
                    m.plotdefaults.seamlinecolor );
            end
            if ~isempty( m.selection.highlightedVxs )
                ni = m.selection.highlightedVxs*2;
                ni = [(ni-1);ni];
                plotpts( theaxes, m.prismnodes( ni, : ), ...
                    '.b', 'MarkerSize', m.plotdefaults.highlightthickness*10 );
            end
        end
    end
    
    % Normal vectors.
    % Gradient vectors.
    % Tensor axes.
    sparsedistance = 0;
    if s.drawnormals || s.drawgradients || s.drawtensoraxes || s.drawdisplacements
        diam = max( max(m.nodes,[],2) - min(m.nodes,[],2) );
        sparsedistance = m.plotdefaults.sparsedistance*diam;
        if sparsedistance > 0
            decorscale = m.plotdefaults.decorscale * ...
                            max( m.globalDynamicProps.cellscale, sparsedistance );
        else
            decorscale = m.plotdefaults.decorscale * m.globalDynamicProps.cellscale;
        end
    end
    if s.drawnormals || s.drawgradients || s.drawtensoraxes
        ccentres = cellcentres( m );
        if sparsedistance > 0
            if true
                bordernodeindexes = unique( m.edgeends( m.edgecells(:,2)==0, : ) );
                bordernodes = m.nodes(bordernodeindexes,:);
                borderdistance = decorscale/2.4; %sparsedistance/2;
                if ~m.plotdefaults.staticdecor
                    m.decorFEs = [];
                    m.decorBCs = [];
                end
                if ~m.plotdefaults.thick
                    [selcc,selbc,selpts] = ...
                        randPointsOnSurface( m.nodes, ...
                            m.tricellvxs, m.cellareas, sparsedistance, ...
                            m.decorFEs(:), m.decorBCs, bordernodes, borderdistance );
                elseif m.plotdefaults.decorateAside
                    [selcc,selbc,selpts] = ...
                        randPointsOnSurface( m.prismnodes(1:2:end,:), ...
                            m.tricellvxs, m.cellareas, sparsedistance, ...
                            m.decorFEs(:), m.decorBCs, bordernodes, borderdistance );
                else
                    [selcc,selbc,selpts] = ...
                        randPointsOnSurface( m.prismnodes(2:2:end,:), ...
                            m.tricellvxs, m.cellareas, sparsedistance, ...
                            m.decorFEs(:), m.decorBCs, bordernodes, borderdistance );
                end
                if m.plotdefaults.staticdecor
                    m.decorFEs = selcc;
                    m.decorBCs = selbc;
                end
            else
                selcc = sparsifyPoints( ccentres, sparsedistance );
                selbc = repmat( [1 1 1]/3, [length(selcc), 1] );
                selpts = ccentres;
            end
        else
            numcells = size(m.tricellvxs,1);
            selcc = (1:numcells)';
            selbc = repmat( [1 1 1]/3, [numcells, 1] );
            selpts = ccentres;
        end
        selpts = selpts( m.visible.cells(selcc), : );
        selbc = selbc( m.visible.cells(selcc), : );
        selcc = selcc( m.visible.cells(selcc) );
        ccentres = ccentres(selcc,:);
    end
    
    if s.drawnormals
        if s.decorateAside
            normalsign = -1;
        else
            normalsign = 1;
        end
        m.normalhandles = myquiver3( ...
            ccentres, ...
            (normalsign * decorscale) * m.unitcellnormals(selcc,:), ...
            [], m.plotdefaults.arrowheadsize, m.plotdefaults.arrowheadratio, 1, 0, ...
            'Color', 'k', ...
            'LineWidth', m.plotdefaults.arrowthickness, ...
            'Parent', theaxes );
    end
    
    if s.drawgradients
        arrows = true;
        crosses = false;
        equalsize = true;
        m.strainhandles = plotGrowthCrosses( ...
            m, selcc, selbc, selpts, sparsedistance, ...
            arrows, crosses, equalsize, decorscale, ...
            trimStruct( m.plotdefaults, ...
                        { 'arrowthickness', ...
                          'crossthickness', ...
                          'arrowheadsize', ...
                          'arrowheadratio', ...
                          'highgradcolor', ...
                          'lowgradcolor' } ) );
    end
    
    if s.drawtensoraxes ... % && isfield( s, 'axesdrawn' ) && ~isempty(s.axesdrawn) ...
             && isfield( m.plotdata, 'frames' ) && ~isempty(m.plotdata.frames) ...
             && isfield( m.plotdata, 'selcpts' ) && ~isempty(m.plotdata.selcpts)
        cpts = m.plotdata.components( selcc, m.plotdata.selcpts );
        maxperFE = max( abs(cpts), [], 2 );
        poscpts = maxperFE > 0;
        if any( poscpts )
            cpts = abs( cpts( poscpts, : ) );
            maxperFE = maxperFE( poscpts );
            frameselcc = selcc(poscpts);
            framevecs = m.plotdata.frames(:,m.plotdata.selcpts,frameselcc);
            framevecs = reshape( permute( framevecs, [2,3,1] ), [], 3 );
            if s.unitcrosses
                cpts = cpts ./ repmat( maxperFE, 1, size(cpts,2) );
            else
                cpts = cpts / max(maxperFE);
            end
            cpts = repmat( reshape( cpts', [], 1 ), 1, 3 );
          % maxcpt = max(abs(cpts(:)));
          % cpts = cpts/maxcpt;
            framevecs = framevecs .* cpts;
            m.displacementhandles = myquiver3( ...
                reshape( repmat( selpts(poscpts,:)', length(m.plotdata.selcpts), 1 ), 3, [] )', ...
                scalevecs( framevecs, decorscale ), ...
                [], ...
                0, 0, 0.4, 0.4, ...
                'Color', m.plotdefaults.axescolor, ...
                'LineWidth', m.plotdefaults.arrowthickness, ...
                'Parent', theaxes );
        end
    end
    
    if s.drawdisplacements
        selnodes = sparsifyPoints( m.nodes, sparsedistance );
        selnodes(~m.visible.nodes(selnodes)) = [];
        selprismnodes = selnodes*2;
        selprismnodes = reshape([ selprismnodes-1, selprismnodes ]', [], 1 );
        if ~isempty(m.displacements)
            m.displacementhandles = myquiver3( ...
                m.prismnodes(selprismnodes,:), ...
                scalevecs( m.displacements(selprismnodes,:), decorscale ), ...
                [], ...
                0.5, 0.6, 1, 0, ...
                'Color', [0.4 0.2 0.1], ...
                'LineWidth', m.plotdefaults.arrowthickness, ...
                'Parent', theaxes );
        end
    end

    % Plot the second layer.
    if s.drawsecondlayer && hasSecondLayer( m )
        m.plothandles.secondlayerhandle = plotSecondLayer( m );
    else
        m.plothandles.secondlayerhandle = [];
    end
    
    drawAxisLabels( theaxes, s.bgcolor );

    hold( theaxes, 'on' );
    
    if wantColorBar
        drawColorbar( pichandles.colorbar, s.cmap, s.crange, s.cmaptype, get( pichandles.picture, 'Color' ) );
    else
        blankColorBar( pichandles.colorbar, get( pichandles.picture, 'Color' ) );
    end
    
    if ~isempty( s.userplotproc )
        try
            m = s.userplotproc( m, theaxes );
        catch
            e = lasterror();
            fprintf( 1, 'User plotting procedure %s failed.  MATLAB error message:\n%s\n', ...
                char( s.userplotproc ), e.message );
        end
    end

    for i=1:length(m.pictures)
        if i >= 2
            fig = m.pictures(i);
            pichandles = guidata( fig );
        end
        grid( pichandles.picture, 'off' );
        if s.axisVisible
            axis( pichandles.picture, 'on' );
        else
            axis( pichandles.picture, 'off' );
        end
        if s.drawlegend
            set( pichandles.legend, 'Visible', 'on' );
        else
            set( pichandles.legend, 'Visible', 'off' );
        end
        
        if i >= 2
            colormap( pichandles.picture, s.cmap );
            caxis( pichandles.picture, s.crange );
            cla(pichandles.picture);
            copyobj( get( theaxes, 'Children' ), pichandles.picture );
        end
        if isfield( pichandles, 'pictureBackground' )
            uistack( pichandles.pictureBackground, 'bottom' );
        end
    end
    
    hold( theaxes, 'off' );
end

function plotnodes( m, theaxes, nodeindexes, marksize )
    plotpts( theaxes, m.nodes( nodeindexes, : ), 'Color', 'k', 'MarkerSize', marksize );
end

function plotedges( m, theaxes, side, edges, edgethickness, edgecolor )
    [lw,ls] = basicLineStyle( edgethickness );
    switch side
        case ''
            coords = m.nodes;
        case 'A'
            coords = m.prismnodes( 1:2:end, : );
        case 'B'
            coords = m.prismnodes( 2:2:end, : );
    end
        x = [ coords( m.edgeends(edges,1), 1 ), coords( m.edgeends(edges,2), 1 ) ]';
        y = [ coords( m.edgeends(edges,1), 2 ), coords( m.edgeends(edges,2), 2 ) ]';
        z = [ coords( m.edgeends(edges,1), 3 ), coords( m.edgeends(edges,2), 3 ) ]';
    line( ...
        x, y, z, ...
        'LineWidth', lw, ...
        'LineStyle', ls, ...
        'Color', edgecolor, ...
        'Parent', theaxes, ...
        'ButtonDownFcn', @doEdgeClick );
end

function plotborder( m, theaxes )
    borderedgeends = m.edgeends(m.visible.borderedges,:)';
    bordernodelist = find(m.visible.bordernodes);
    bordernoderenumber = zeros(size(m.visible.bordernodes));
    bordernoderenumber( m.visible.bordernodes ) = (1:length(bordernodelist))';
    borderedgeendsreduced = bordernoderenumber(borderedgeends);
    bordernodepositions = m.nodes( m.visible.bordernodes, : );
    borderprismnodesB = bordernodelist*2;
    borderprismnodesA = borderprismnodesB-1;
    bordernodeApositions = m.prismnodes( borderprismnodesA, : );
    bordernodeBpositions = m.prismnodes( borderprismnodesB, : );
    edgeApositions = reshape( bordernodeApositions( borderedgeendsreduced, : )', 6, [] );
    edgeBpositions = reshape( bordernodeBpositions( borderedgeendsreduced, : )', 6, [] );
    
    if isfield( m.plotdata, 'value' ) && ~isempty( m.plotdata.value )
        bordernodeABpositions = reshape( [ bordernodeApositions'; bordernodeBpositions' ], ...
                                         3, [] )';
        paintborder( m, theaxes, ...
            bordernodeABpositions, ...
            borderedgeendsreduced, ...
            m.plotdata.pervertex, ...
            false, ...
            m.plotdata.value );
    else
        bordernodeAMpositions = reshape( [ bordernodeApositions'; bordernodepositions' ], ...
                                         3, [] )';
        bordernodeBMpositions = reshape( [ bordernodeBpositions'; bordernodepositions' ], ...
                                         3, [] )';
        paintborder( m, theaxes, ...
            bordernodeAMpositions, ...
            borderedgeendsreduced, ...
            m.plotdata.pervertexA, ...
            m.plotdefaults.taper, ...
            m.plotdata.valueA );
        paintborder( m, theaxes, ...
            bordernodeBMpositions, ...
            borderedgeendsreduced, ...
            m.plotdata.pervertexB, ...
            m.plotdefaults.taper, ...
            m.plotdata.valueB );
    end

    if m.plotdefaults.drawedges > 0
        % Draw border edges.
        [lw,ls] = basicLineStyle( m.plotdefaults.FEthinlinesize );
        borderedgehandles = line( ...
            [ bordernodeApositions( :, 1 ), bordernodeBpositions( :, 1 ) ]', ...
            [ bordernodeApositions( :, 2 ), bordernodeBpositions( :, 2 ) ]', ...
            [ bordernodeApositions( :, 3 ), bordernodeBpositions( :, 3 ) ]', ...
            'LineWidth', lw, ...
            'LineStyle', ls, ...
            'Color', m.plotdefaults.FElinecolor, ...
            'Parent', theaxes, ...
            'ButtonDownFcn', @doVxClick );

        for i=1:length(borderedgehandles)
            set( borderedgehandles(i), ...
                'UserData', struct('vertex',bordernodelist(i)) );
        end
        if m.plotdefaults.drawedges < 2
            line( edgeApositions([1 4],:), ...
                  edgeApositions([2 5],:), ...
                  edgeApositions([3 6],:), ...
                  'LineWidth', lw, ...
                  'LineStyle', ls, ...
                  'Color', m.plotdefaults.FElinecolor, ...
                  'Parent', theaxes, ...
                  'ButtonDownFcn', @doEdgeClick );
            line( edgeBpositions([1 4],:), ...
                  edgeBpositions([2 5],:), ...
                  edgeBpositions([3 6],:), ...
                  'LineWidth', lw, ...
                  'LineStyle', ls, ...
                  'Color', m.plotdefaults.FElinecolor, ...
                  'Parent', theaxes, ...
                  'ButtonDownFcn', @doEdgeClick );
        end
    end
end

function borderhandles = paintborder( m, theaxes, bordernodepositions, borderedgeends, ...
                                      pervertex, taper, data )
    quadVxIndexes = [ borderedgeends*2-1; ...
                      borderedgeends([2 1],:)*2 ];
    numedges = size( borderedgeends, 2 );
    ex1 = reshape( ...
            bordernodepositions( quadVxIndexes, 1 ), ...
            4, [] );
    ey1 = reshape( ...
            bordernodepositions( quadVxIndexes, 2 ), ...
            4, [] );
    ez1 = reshape( ...
            bordernodepositions( quadVxIndexes, 3 ), ...
            4, [] );
    if ~isempty( data )
        if pervertex
            data = data( m.visible.bordernodes, : );
        else
            data = data( m.visible.bordercells, : );
        end
    end
    if isempty( data )
        % Draw white quads.
        ec1 = 'w';
    elseif pervertex
        if size( data, 2 )==1
            if taper
                ec1 = reshape( ...
                            [ data( borderedgeends ); ...
                              zeros( 2, numedges ) ], ...
                            4, [] );
            else
                ec1 = reshape( ...
                            [ data( borderedgeends ); ...
                              data( borderedgeends( [2 1], : ) ) ], ...
                            4, [] );
            end
        elseif size( data, 2 )==3
            ec1 = reshape( data( borderedgeends, : ), 2, [], 3 );
            if taper
                ec1 = [ ec1; ones(size(ec1)) ];
            else
                ec1 = [ ec1; ec1([2 1],:,:) ];
            end
        end
    elseif size( data, 2 )==1
        ec1 = data';
        if taper
            ec1 = [ ec1; ec1; zeros(2,size(ec1,2)) ];
        end
    else
        ec1 = permute( data, [3,1,2] );
        if taper
            white = ones(size(ec1));
            ec1 = [ ec1; ec1; white; white ];
        end
    end
    borderhandles = [];
    commonPlotArgs = { 'FaceLighting', m.plotdefaults.lightmode, ...
                       'AmbientStrength', m.plotdefaults.ambientstrength, ...
                       'FaceAlpha', m.plotdefaults.alpha, ...
                       'LineStyle', 'none', ...
                       'Parent', theaxes };
    if size(ex1,2)==3
        % We can't plot all three at once, as Matlab will mistake the
        % three-element row vector of color indexes for a single
        % RGB value.
        if ~isempty(ex1)
            if ischar( ec1 )
                ec1a = ec1;
                ec1b = ec1;
            else
                ec1a = ec1([1 2]);
                ec1b = ec1(3);
            end
            borderhandles = [ ...
                patch( ex1(:,[1 2]), ey1(:,[1 2]), ez1(:,[1 2]), ec1a, commonPlotArgs{:} ); ...
                patch( ex1(:,3), ey1(:,3), ez1(:,3), ec1b, commonPlotArgs{:} ) ...
            ];
        end
    else
        borderhandles = patch( ex1, ey1, ez1, ec1, commonPlotArgs{:} );
    end
end

function m = calculatePlotData( m, side )
    fn_value = ['value' side];
    fn_pervertex = ['pervertex' side];
    fn_perelement = ['perelement' side];
    fn_tensors = ['tensor' side];
    fn_morphogen = ['morphogen' side];
    fn_outputquantity = ['outputquantity' side];
    fn_outputaxes = ['outputaxes' side];
    fn_perelementaxes = ['perelementaxes' side];
    fn_perelementcomponents = ['perelementcomponents' side];
%     fn_axesdrawn = ['axesdrawn' side];
%     fn_axesquantity = ['axesquantity' side];
    fn_components = ['components' side];
    fn_selcpts = ['selcpts' side];
    fn_frames = ['frames' side];
%     if ~isempty(m.plotdefaults.(fn_axesdrawn))
%         if ~isempty(m.plotdefaults.(fn_axesquantity))
%             m.plotdata.axes = foo;
%         elseif ~isempty(m.plotdefaults.(fn_tensors))
%         elseif ~isempty(m.plotdefaults.(fn_outputquantity))
%         end
%     end
    haveuserframes = ~isempty( m.plotdefaults.(fn_perelementaxes) );
    if haveuserframes
        m.plotdata.(fn_frames) = m.plotdefaults.(fn_perelementaxes);
        numaxes = size( m.plotdefaults.(fn_perelementaxes), 2 );
        m.plotdata.(fn_selcpts) = 1:numaxes;
        if isempty( m.plotdefaults.(fn_perelementcomponents) );
            m.plotdata.(fn_components) = ones( size( m.tricellvxs, 1 ), numaxes );
        else
            m.plotdata.(fn_components) = m.plotdefaults.(fn_perelementcomponents);
        end
    end
    if ~isempty(m.plotdefaults.(fn_pervertex))
        m.plotdata.(fn_pervertex) = true;
        m.plotdata.(fn_value) = m.plotdefaults.(fn_pervertex);
    elseif ~isempty(m.plotdefaults.(fn_perelement))
        m.plotdata.(fn_pervertex) = false;
        m.plotdata.(fn_value) = m.plotdefaults.(fn_perelement);
    elseif ~isempty( m.plotdefaults.(fn_morphogen) )
        mgenindexes = FindMorphogenIndex( m, m.plotdefaults.(fn_morphogen) );
        m.plotdata.(fn_pervertex) = true;
        if length( mgenindexes )==1
            m.plotdata.(fn_value) = m.morphogens( :, mgenindexes );
        else
            m.plotdata.(fn_value) = multicolourmgens( ...
                    m, mgenindexes, m.plotdefaults.multibrighten, [] );
        end
    elseif ~isempty( m.plotdefaults.(fn_tensors) )
        m.plotdata.(fn_pervertex) = false;
        [m.plotdata.(fn_value), ...
         m.plotdata.(fn_components), ...
         theframes, ...
         m.plotdata.(fn_selcpts)] = getTensorValues( ...
            m, ...
            m.plotdefaults.(fn_tensors), ...
            m.plotdefaults.(fn_outputaxes), ...
            m.cellFrames );
        if ~haveuserframes
            m.plotdata.(fn_frames) = theframes;
        end
    elseif ~isempty( m.plotdefaults.(fn_outputquantity) )
        m.plotdata.(fn_pervertex) = false;
        if regexp( m.plotdefaults.(fn_outputquantity), '^rotation' )
            % Set rotation data.  NOT IMPLEMENTED.
            m.plotdata.(fn_value) = zeros( size(m.tricellvxs,1), 1 );
        else
            [m.plotdata.(fn_value), ...
             m.plotdata.(fn_components), ...
             m.plotdata.(fn_frames), ...
             m.plotdata.(fn_selcpts)] = getTensorValues( ...
                m, ...
                m.plotdefaults.(fn_outputquantity), ...
                m.plotdefaults.(fn_outputaxes), ...
                m.cellFrames );
        end
    else
        % Nothing to be plotted.
    end
end

function m = makeCanvasColor( m, side )
    fn_value = ['value' side];
    fn_pervertex = ['pervertex' side];
    if ~isfield( m.plotdata, fn_value ) || isempty( m.plotdata.(fn_value) )
        if ~isfield( m.plotdata, fn_pervertex )
            m.plotdata.(fn_pervertex) = false;
        end
        if m.plotdata.(fn_pervertex)
            numdata = size( m.nodes, 1 );
        else
            numdata = size( m.tricellvxs, 1 );
        end
        m.plotdata.(fn_value) = repmat( m.plotdefaults.canvascolor, numdata, 1 );
    end
end

function [values,components,frames,selectedCpts] = getTensorValues( m, oq, oa, cf )
%[values,axisvectors] = getTensorValues( m, oq, oa, cf )
%   Calculate the tensor values specified by oq (name of an output
%   quantity), oa (set of axes), and cf (cell frames).

    values = zeros( size(cf,3), 1 );
    xx = regexp( oq, '^(?<base>.*)(?<rate>rate)$', 'names' );
    israte = ~isempty(xx);
    if israte
        oq = xx.base;
    end
    oq = regexprep( oq, 'growth', 'strain' );
    oq = regexprep( oq, 'resultant', 'actual' );
    xx = regexp( oq, '^(?<base>.*)(?<bend>bend|strain|anisotropy)$', 'names' );
    if isempty(xx)
        return;
    end
    oq = [ xx.base, 'strain' ];
    isbend = strcmp( xx.bend, 'bend' );
    isaniso = strcmp( xx.bend, 'anisotropy' );
    if ~isfield( m.outputs, oq )
        return;
    end
    polAligned = ~isaniso && (strcmp(oa,'parallel') || strcmp(oa,'perpendicular') || strcmp(oa,'normal') || strcmp(oa,'areal'));
    [components,frames] = getCpts( m.outputs.(oq), isbend, polAligned );
    if ~israte && (m.globalProps.timestep > 0)
        components = components*m.globalProps.timestep;
    end
    selectedCpts = [];
    if isaniso
        values = components(:,1) - components(:,2);
    else
        switch oa
            case {'','total'}
                selectedCpts = 1:size(components,2);
            case 'areal'
                selectedCpts = [1 2];
            case { 'parallel', 'major' }
                selectedCpts = 1;
            case { 'perpendicular', 'minor' }
                selectedCpts = 2;
            case 'normal'
                selectedCpts = 3;
        end
        values = sum( components(:,selectedCpts), 2 );
    end
    
    function [components,frames] = getCpts( tensors, bend, polAligned )
        if bend
            tensors = (tensors.A-tensors.B)/2;
        else
            tensors = (tensors.A+tensors.B)/2;
        end
        if polAligned
            [components,frames] = tensorsToComponents( tensors, m.cellFrames );
        else
            [components,frames] = tensorsToComponents( tensors );
        end
    end
end
 
function drawAxisLabels( theaxes, bgcolor )
    elementColor = contrastColor( bgcolor );
    set( theaxes, ...
        'XColor', elementColor, ...
        'YColor', elementColor, ...
        'ZColor', elementColor );
    theaxesRange = axis(theaxes);
    if length(theaxesRange) < 6
        theaxesRange = [ theaxesRange, 0, 0 ];
    end
    labeloffset = 0.075;
    xaxislabel = get(theaxes,'XLabel');
    yaxislabel = get(theaxes,'YLabel');
    zaxislabel = get(theaxes,'ZLabel');
%     xloc = get( theaxes, 'XAxisLocation' );
%     yloc = get( theaxes, 'YAxisLocation' );
    set( theaxes, ...
        'XColor', elementColor, ...
        'YColor', elementColor, ...
        'ZColor', elementColor );
    set( xaxislabel, ...
        'String','X', ...
        'Position', [ (theaxesRange(1) + theaxesRange(2))/2, ...
                      theaxesRange(3)*(1+labeloffset) - theaxesRange(4)*labeloffset, ...
                      theaxesRange(5) ], ...
        'Rotation', 0, ...
        'Color', elementColor);
    set( yaxislabel, ...
        'String','Y', ...
        'Position', [ theaxesRange(1)*(1+labeloffset) - theaxesRange(2)*labeloffset, ...
                      (theaxesRange(3) + theaxesRange(4))/2, ...
                      theaxesRange(5) ], ...
        'Rotation', 0, ...
        'Color', elementColor );
%    if length(theaxesRange) == 6
        set( zaxislabel, ...
        'String','Z', ...
        'Position', [ theaxesRange(1)*(1+labeloffset) - theaxesRange(2)*labeloffset, ...
                      theaxesRange(4)*(1+labeloffset) - theaxesRange(3)*labeloffset, ...
                      (theaxesRange(5) + theaxesRange(6))/2 ], ...
        'Rotation', 0, ...
        'Color', elementColor );
%    end
end
