WriteCellToVTKStructured(...

function WriteCellToVTKStructured(...
    filename, input_cell, coordinates, values, ...
    varargin)
% WriteCellToVTKStructured (FILENAME, CELL, COORDINATES, VALUES, ...)
% Writes the content of a cell array into an ascii vtk-file 'FILENAME'
% as a structured grid.
%
% For details about this kind of data see 
% http://www.itk.org/Wiki/ParaView/Users_Guide/VTK_Data_Model#Curvilinear_Grid_.28Structured_Grid.29
% and
% www.vtk.org/VTK/img/file-formats.pdf
%
% This function takes a 3D cell array of structs. For each struct, it then
% uses either given fields of the struct or its indices as VTK coordinates,
% and other fields as the values of the grid.
%
% Input:
% filename      destination (string)
% cell          N1xN2xN3 cell array of structs
% coordinates   Cell with three cells of strings. These define what to use
%               for the points' coordinates. (Details below)
% values        Cell with cells of strings. These define the scalar values
%               at each point. (details below)
%
% Optional input parameters of the form ..., 'name', value, ...:
% header        Name of the dataset, will be saved into the file.
% valuenames    Names of the values. The defaults are their field names.
% spacing       3x1 matrix of inter-point spacings if indices are used as
%               coordinates.
% offset        3x1 matrix of offsets for this case.
%
% Example usage:
%
% This example assumes that 'results' is a cell array of structs, where
% each struct has (at least) the following valid fields:
% 'vw', 'effects.energy' and 'effects.charge'
%
% WriteCellToVTKStructured('test.vtk', results, ...
%     {{}, {}, {'effects', 'energy'}}, ...
%     {{'vw'}, {'effects' 'charge'}}, ...
%     'offset', [1 0.4 0.1], ...
%     'spacing', [1 0.05 0.1])
%
% Then use e.g. ParaView to view the resultin file 'test.vtk'.
%

% (c) Daniel Hornung, 2012-01-23
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
% POSSIBILITY OF SUCH DAMAGE.

[N1, N2, N3] = size(input_cell);
N = numel(input_cell);

%%%%%%%%
% Get the optional input parameters
%%%
pars = inputParser;
% header        Name of the dataset, will be saved into the file.
[~, basename, ~] = fileparts(filename);
pars.addParamValue('header', basename);

% valuenames    Names of the values. The defaults are their field names.
pars.addParamValue('valuenames', {[] [] []});

% spacing       3x1 matrix of inter-point spacings if indices are used as
%               coordinates.
pars.addParamValue('spacing', [1 1 1]);

% offset        3x1 matrix of offsets.
pars.addParamValue('offset', [0 0 0]);

pars.parse(varargin{:});
header = pars.Results.header;
header = header(1:min(256, numel(header)));
valuenames= pars.Results.valuenames;
spacing = pars.Results.spacing;
offset = pars.Results.offset;

%%%%% Input params parsing finished %%%%%

%%%%%%%%
% Convenience functions for coordinates and values
%%%

f_co = cell(3,1);
for k = 1:3
    f_co{k} = @(x) x;
    for i = 1:numel(coordinates{k})
        f_co{k} = @(x) getfield(f_co{k}(x), coordinates{k}{i}); %#ok<GFLD>
    end
end

num_v = numel(values);
f_val = cell(num_v, 1);
for k = 1:num_v
    f_val{k} = @(x) x;
    for i = 1:numel(values{k})
        f_val{k} = @(x) getfield(f_val{k}(x), values{k}{i}); %#ok<GFLD>
    end
end

%%%%% Struct parsing functions finished %%%%%

fid = fopen(filename, 'w');

% Writing the header
fprintf(fid, [...
    '# vtk DataFile Version 3.0\n'...
    '%s\n'...
    'ASCII\n' ...
    'DATASET STRUCTURED_GRID\n' ...
    'DIMENSIONS %d %d %d\n' ...
    'POINTS %d double\n' ...
    ], header, N1, N2, N3, N);

% Writing the point coordinates
for i = 1:numel(input_cell)
    [sx, sy, sz] = ind2sub(size(input_cell), i);
    sub = [sx, sy, sz] - 1;
    sub = offset + spacing .* sub;
    input_struct = input_cell{i};
    x = f_co{1}(input_struct);
    if (isstruct(x))
        x = sub(1);
    end
    y = f_co{2}(input_struct);
    if (isstruct(y))
        y = sub(2);
    end
    z = f_co{3}(input_struct);
    if (isstruct(z))
        z = sub(3);
    end
    fprintf(fid, '%f %f %f\n', x, y, z);
end

fprintf(fid, 'POINT_DATA %d\n', N);

% Writing the point values
for k = 1:numel(f_val)
    if numel(valuenames{k}) == 0
        valuenames{k} = fullfile(values{k}{:});
        disp(['Using default valuename ' valuenames{k} ...
            ' for value ' num2str(k)])
    end
    fprintf(fid, [ ...
        'SCALARS %s double 1\n' ...
        'LOOKUP_TABLE default\n' ...
        ], ...
        valuenames{k} ...
    );
    for i = 1:numel(input_cell)
        if ~mod(i, size(input_cell, 1))
            fprintf(fid, '\n');
        end
        input_struct = input_cell{i};
        fprintf(fid, ...
            '%f ', ...
            f_val{k}(input_struct));
    end
end
fclose(fid);

end