function status = writemeshobj( filedir, filename, mesh, surface, minimal, throughfaces, bbox )
%writemesh( mesh, filedir, filename, surface, minimal )
%    Write the mesh to a file in extended OBJ format.
% If minimal is true (default is false), then only the nodes and cells will
% be written, otherwise everything in the mesh is written.
% If surface and minimal are both true, the resulting file should be a
% standard OBJ file, readable by any software that reads OBJ files.

    if nargin < 4, surface = false; end
    if nargin < 5, minimal = false; end
    fullfilename = fullfile( filedir, filename );
    if ~endsWithM(fullfilename,'.obj')
        fullfilename = strcat(fullfilename,'.obj');
    end
    fid = fopen(fullfilename,'w');
    if fid==0
        fprintf( 1, 'Cannot write to file %s.\n', fullfilename );
        status = 0;
        return;
    end
    
    mesh = makeTRIvalid( mesh );
    
    fprintf( fid, '# File written by Growth Toolbox at %s.\n', ...
        datestr(clock) );
    fprintf( fid, '# Mesh iteration %d, time stamp %.3f.\n\n', ...
        mesh.globalDynamicProps.currentIter, mesh.globalDynamicProps.currenttime );
    
    mc = meshColors(mesh);
    mcheader = '';
    mcprefix = '';
    if ~isempty(mc)
        if mesh.plotdata.pervertex
            mcheader = 'Vertex colors';
            mcprefix = 'vc';
        else
            mcheader = 'Face colors';
            mcprefix = 'fc';
        end
    end
    numMgens = size(mesh.morphogens,2);
    mgeninfo = mesh.morphogens;
    maxmgens = max(mgeninfo,[],1);
    notzeros = find(maxmgens ~= 0);
    mgeninfo(:,notzeros) = mgeninfo(:,notzeros) ./ ...
        repmat( maxmgens(notzeros), size(mesh.morphogens,1), 1 );
    if surface
        % Vertex positions.
        pts = fitToBbox( mesh.nodes, bbox );
        writearray( fid, pts, 'v', '%f %f %f' );

        % Faces.
        writemesharray( fid, mesh, 'tricellvxs', 'f', '%d %d %d' );
        
        % Vertex or face colours.
        if ~isempty(mc)
            writearray( fid, mc, mcprefix, '%f %f %f', 0, mcheader );
        end

        if ~minimal
            % Morphogen names.
            writeMgenNames( fid, mesh.mgenIndexToName );

            % Morphogen values.
            mgenfmt = [ '%f', repeatString( ' %f', numMgens-1 ) ];
            writearray( fid, mgeninfo, 'm', mgenfmt, 0, 'Morphogens' );
        end
    else
        % Vertex positions.
        pts = fitToBbox( mesh.prismnodes, bbox );
        writearray( fid, pts, 'v', '%f %f %f' );

        % Faces.
        if throughfaces
            throughedges = 1:size(mesh.edgecells,1);
        else
            throughedges = find(mesh.edgecells(:,2)==0);
        end
        throughcells = mesh.edgecells(throughedges,1);
        bsidevxs = mesh.tricellvxs*2;
        asidevxs = bsidevxs-1;
        writearray( fid, asidevxs(:,[1 3 2]), 'f', '%d %d %d', 0, 'A side' );
        writearray( fid, bsidevxs, 'f', '%d %d %d', 0, 'B side' );
        if ~isempty(throughedges)
            if throughfaces
                fprintf( fid, '# Through faces\n' );
            else
                fprintf( fid, '# Rim faces\n' );
            end
            for i=throughedges(:)' % size(mesh.edgeends,1)
                c = mesh.edgecells(i,1);
                % find edge i in cell c
                cei = find( mesh.celledges(c,:)==i, 1 );
                vs = mesh.tricellvxs( c, othersOf3( cei ) );
                % vs contains the two ends of edge i, in positive order.  These
                % are the same vertexes as mesh.edgeends(i,:), but the latter
                % are in arbitrary order.
                bvs = vs*2;
                avs = bvs-1;
                fprintf( fid, 'f %d %d %d\nf %d %d %d\n', ...
                    avs(1), avs(2), bvs(2), ... % The A (bottom) side in negative order.
                    avs(1), bvs(2), bvs(1) ); % The B (top) side in positive order.
            end
        end

        % Vertex or face colours.
        if ~isempty(mc) && ~isempty(mcprefix)
            if mesh.plotdata.pervertex
                mc = reshape( [mc'; mc'], 3, [] )';
            else
                throughcolors = mc(throughcells,:);
                throughcolors = reshape( [throughcolors';throughcolors'], 3, [] )';
                mc = [mc; mc; throughcolors];
            end
            writearray( fid, mc, mcprefix, '%f %f %f', 0, 'Vertex colors' );
        end

        if ~minimal
            % Morphogen names.
            writeMgenNames( fid, mesh.mgenIndexToName );

            % Morphogen values.
            mgeninfo = reshape( [mgeninfo'; mgeninfo'], numMgens, [] )';
            mgenfmt = [ '%f', repeatString( ' %f', numMgens-1 ) ];
            writearray( fid, mgeninfo, 'm', mgenfmt, 0, 'Morphogens' );
        end
    end

    if ~minimal
        globFields = fieldnames( mesh.globalProps );
        for i=1:length(globFields)
            if isnumeric( mesh.globalProps.(globFields{i}) )
                writeGlobalProp( fid, mesh, globFields{i} );
            end
        end
        fprintf( fid, '\n' );

        writemesharray( fid, mesh, 'celllabel', 'celllabel', '%d', 1 );
        morphogenFormat = [ '%f', repmat( ' %f', 1, size(mesh.morphogens,2)-1 ) ];;
        writemesharray( fid, mesh, 'morphogens', 'g', morphogenFormat );
        writemesharray( fid, mesh, 'morphogenclamp', 'gcl', morphogenFormat );
    end

    status = fclose(fid)==0;
    if ~status
        fprintf( 1, 'Could not close file %s after writing mesh.\n', ...
            fullfilename );
    end
end

function writeMgenNames( fid, names )
    fprintf( fid, '# Morphogen names\nmn' );
    for i=1:length(names)
        fprintf( fid, ' %s', names{i} );
    end
    fprintf( fid, '\n\n' );
end

function writeGlobalProp( fid, mesh, globField )
    if isfield( mesh.globalProps, globField )
        fprintf( fid, '%s', globField );
        fprintf( fid, ' %f', mesh.globalProps.(globField) );
        fprintf( fid, '\n' );
    end
end

function writearray( fid, array, name, fmt, indexed, message )
    if (nargin<4) || isempty(fmt)
        fmt = '%f';
    end
    if (nargin<5) || isempty(indexed)
        indexed = 0;
    end
    if (nargin<6) || isempty(message)
        message = [ '# ', name ];
    end
    fprintf( fid, '# %s\n', message );
    if indexed
        fprintf( fid, [ name, ' %d ', fmt, '\n' ], [ (1:size(array,1))', array ]' );
    else
        fprintf( fid, [ name, ' ', fmt, '\n' ], array' );
    end
    fprintf( fid, '\n' );
end

function writemesharray( fid, mesh, field, name, fmt, indexed, message )
    if isfield( mesh, field )
        if (nargin<4) || isempty(name)
            name = field;
        end
        if (nargin<5) || isempty(fmt)
            fmt = '%f';
        end
        if (nargin<6) || isempty(indexed)
            indexed = 0;
        end
        if (nargin<7) || isempty(message)
            message = name;
        end
        writearray( fid, mesh.(field), name, fmt, indexed, message )
    end
end

