gtVTKStructGridWriter.m

function gtVTKStructGridWriter(vol, filename, varargin)
% GTVTKSTRUCTGRIDWRITER  Write a 3D volume in a VTK file.
%
%   gtVTKStructGridWriter(vol, filename, varargin)
%   --------------------------------------------------------------------------
%
%     INPUT:
%       vol         = <3Dimage> Input 3D volume
%       filename    = <string>  Path to the output VTK file
%
%     OPTIONAL INPUT (varargin as a list of pairs, see parse_pv_pairs.m):
%       'binary'    = <bool>    Write output in binary format (true by default)
%       'spacing'   = <double>  Vector defining spacing in each space dimension
%                               ([1 1 1] by default)
%       'origin'    = <double>  Vector defining volume origin
%                               ([0 0 0] by default)
%       'header'    = <string>  VTK file header description
%       'datatype'  = <string>  Force data type to consider for input volume
%       'cmap'      = <cmap>    Color map for volume data for vizualization
%       'opacity'   = <double>  Opacity for vizualization (1 by default)
%       'verbose'   = <bool>    Verbose mode
%
%     Version 001 24-11-2012 by YGuilhem yoann.guilhem@esrf.fr

% General parameters
par.spacing = ones(1, 3);
par.origin  = zeros(1, 3);
par.header  = 'VTK structured volume created by DCT http://sourceforge.net/projects/dct/';
par.binary = true;
par.verbose = false;
par.cmap = [];
par.opacity = 1;

% Parameters for each field
par.dataname = {'Values'};
par.datatype = [];

% Parse parameters
par = parse_pv_pairs(par, varargin);

%% Check parameters

% Check spacing
if length(par.spacing) == 1
    par.spacing = par.spacing * ones(1, 3);
end

% Check origin
if length(par.origin) ~= 3
    gtError('gtVTKStructGridWriter:wrong_origin', ...
        'Origin should be a 3 coordinates vector!');
end

% Check number of fields
if ~iscell(vol)
    vol = {vol};
end
Nvol = length(vol);

% Check size of fields
if Nvol > 1
    for ivol=2:Nvol
        if any(size(vol{ivol}) ~= size(vol{1}))
            gtError('gtVTKStructGridWriter:wrong_volume_size', ...
                'All the volumes should have the same size!');
        end
    end
end
[nx, ny, nz] = size(vol{1});
ntot = nx*ny*nz;

% Check cmap
if ~isempty(par.cmap)
    for ivol=1:Nvol
        if size(par.cmap, 1) < max(max(max(vol{ivol})))
            warning('gtVTKStructGridWriter:cmap_warning', ...
                'You may have not enough values in your color map...');
        end
    end
end

%% Get file parts and check extension
[fpath, fname, fext] = fileparts(filename);
if isempty(fext)
    fext = '.vtk';
    filename = fullfile(fpath, [fname fext]);
elseif ~strcmp(fext, '.vtk')
    gtError('gtVTKMeshWriter:wrong_file_extension', ...
        'Output file extension should be ''.vtk''!');
end

%% Guess input & output datatype
if isempty(par.datatype)
    for ivol=1:Nvol
        par.datatype{ivol} = class(vol{ivol});
    end
elseif ~iscell(par.datatype)
    tmp = par.datatype;
    for ivol=1:Nvol
        if class(vol{ivol}) ~= tmp
            gtError('gtVTKStructGridWriter:wrong_data_type', ...
                'All the volumes should have the same data type!');
        end
        par.datatype{ivol} = tmp;
    end
end

outputTypeVTK = cell(Nvol);
outputType = cell(Nvol);
outputFMT = cell(Nvol);
for ivol=1:Nvol
    [outputTypeVTK{ivol}, outputType{ivol}, outputFMT{ivol}] = ...
        gtVTKConvertTypeMatlab2VTK(par.datatype{ivol});
end

if par.verbose
    disp(['Header        : ' par.header]);
    disp(['Spacing       : ' num2str(par.spacing)]);
    disp(['Origin        : ' num2str(par.origin)]);
    disp(['Color map     : ' num2str(par.cmap)]);
    for ivol=1:Nvol
        disp(['Field         : ' num2str(ivol)]);
        disp(['  Dataname      : ' par.dataname]);
        disp(['  OutputType    : ' outputType]);
        disp(['  OutputTypeVTK : ' outputTypeVTK]);
    end
end

%% Open the file
fid = fopen(filename, 'w');
if fid ==-1
    gtError('gtVTKStructGridWriter:bad_output_file', 'Can''t open the file.');
end

%% Write ASCII file header
fprintf(fid, '# vtk DataFile Version 3.0\n');
fprintf(fid, '%s\n', par.header);
if par.binary
    fprintf(fid, 'BINARY\n');
else
    fprintf(fid, 'ASCII\n');
end
fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
fprintf(fid, 'DIMENSIONS %d %d %d\n', nx, ny, nz);
fprintf(fid, 'ORIGIN  %4.3f   %4.3f  %4.3f\n', par.origin);
fprintf(fid, 'SPACING %4.3f   %4.3f  %4.3f\n', par.spacing);

for ivol=1:Nvol
    fprintf(fid, 'POINT_DATA %d\n', ntot);
    fprintf(fid, 'SCALARS %s %s 1\n', par.dataname{ivol}, outputTypeVTK{ivol});

    if ~isempty(par.cmap)
        fprintf(fid, 'LOOKUP_TABLE map\n');
    else
        fprintf(fid, 'LOOKUP_TABLE default\n');
    end

    %% Write volume data
    if par.binary
        count = fwrite(fid, reshape(vol{ivol}, 1, ntot), ...
            par.datatype{ivol}, 0, 'ieee-be');
        if count ~= numel(vol{ivol})
            warning('gtVTKStructGridWriter:structgrid_volume_warning', ...
                'Problem in writing structured grid volume.');
        end
        fprintf(fid, '\n');
    else
        fprintf(fid, [outputFMT{ivol} '\n'], reshape(vol{ivol}, 1, ntot));
    end
end

%% Write lookup table if needed
if ~isempty(par.cmap)
    mapSize = size(par.cmap);
    if mapSize(2) == 3
        output = [par.cmap'; par.opacity*ones(1, mapSize(1))];
    elseif mapSize(2) == 4
        output = par.cmap';
    else
        gtError('gtVTKStructGridWriter:wrong_cmap_size', ...
            'Color map sould be N row vector of size 3 or 4!');
    end
    fprintf(fid, 'LOOKUP_TABLE map %d\n', mapSize(1));
    if par.binary
        count = fwrite(fid, uint8(225*output), 'uchar', 0, 'ieee-be');
        if count ~= numel(output)
            warning('gtVTKMeshWriter:cmap_warning', ...
                'Problem in writing color map!');
        end
        fprintf(fid, '\n');
    else
        fprintf(fid, '%f %f %f %f\n', output);
    end
end

%% Close file and terminate
fclose(fid);

end % end of function