function hs = mcxpreview(cfg, varargin)
%
% Format:
%    hs=mcxpreview(cfg);
%
% Preview the simulation configuration for both MCXLAB and MMCLAB
%
% Author: Qianqian Fang <q.fang at neu.edu>
%
% Input:
%    cfg: a struct, or struct array. Each element of cfg defines
%         the parameters associated with a simulation. Please run
%         'help mcxlab' or 'help mmclab' to see the details.
%         mcxpreview supports the cfg input for both mcxlab and mmclab.
%
% Output:
%    hs: a struct to the list of handles of the plotted domain elements.
%
% Dependency:
%    this function depends on the Iso2Mesh toolbox (http://iso2mesh.sf.net)
%
% This function is part of Monte Carlo eXtreme (MCX) URL: http://mcx.space
%
% License: GNU General Public License version 3, please read LICENSE.txt for details
%

if (nargin == 0)
    error('input field cfg must be defined');
end

if (~isstruct(cfg))
    error('cfg must be a struct or struct array');
end

if (~exist('latticegrid', 'file'))
    error('you must install the iso2mesh toolbox from https://github.com/fangq/iso2mesh first');
end

len = length(cfg);
if (nargout > 0)
    hs = cell(len, 1);
end

rngstate = rand ('state');

randseed = hex2dec('623F9A9E');
rand('state', randseed);
surfcolors = rand(1024, 3);
isholdplot = ishold;

for i = 1:len
    if (~isfield(cfg(i), 'vol') && ~isfield(cfg(i), 'node') && ~isfield(cfg(i), 'shapes'))
        error('cfg.vol or cfg.node or cfg.shapes is missing');
    end

    if (i > 1)
        figure;
    end
    hold on;

    voxelsize = 1;
    if (isfield(cfg(i), 'unitinmm'))
        voxelsize = cfg(i).unitinmm;
    end

    offset = 1;
    if (isfield(cfg(i), 'issrcfrom0') && cfg(i).issrcfrom0 == 1 || isfield(cfg(i), 'node'))
        offset = 0;
    end

    hseg = [];
    if (isfield(cfg(i), 'vol') && ~isfield(cfg(i), 'node'))
        % render mcxlab voxelated domain
        dim = size(cfg(i).vol);
        if (ndims(cfg(i).vol) == 4) % for spatially varying medium
            dim = dim(2:end);
            cfg(i).vol = ones(dim);
        end
        [bbxno, bbxfc] = latticegrid(0:dim(1):dim(1), 0:dim(2):dim(2), 0:dim(3):dim(3));
        hbbx = plotmesh((bbxno + offset) * voxelsize, bbxfc, 'facecolor', 'none');

        val = unique(cfg(i).vol(:));
        val(val == 0) = [];
        padvol = zeros(dim + 2);
        padvol(2:end - 1, 2:end - 1, 2:end - 1) = cfg(i).vol;

        if (length(val) > 1)
            hseg = zeros(length(val), 1);
            for id = 1:length(val)
                [no, fc] = binsurface(padvol == val(id));
                hseg(id) = plotmesh((no - 1) * voxelsize, fc, 'facealpha', 0.3, 'linestyle', 'none', 'facecolor', surfcolors(val(id) + 1, :), varargin{:});
            end
        end
    elseif (isfield(cfg(i), 'node') && isfield(cfg(i), 'elem'))
        % render mmclab mesh domain

        elemtype = ones(size(cfg(i).elem, 1), 1);
        if (isfield(cfg(i), 'elemprop'))
            elemtype = cfg(i).elemprop;
        else
            if (size(cfg(i).elem, 2) > 4)
                elemtype = cfg(i).elem(:, 5);
            end
        end
        etypes = unique(elemtype);
        no = cfg(i).node * voxelsize;
        hseg = zeros(length(etypes), 1);
        for id = 1:size(etypes, 1)
            hseg(id) = plotmesh(no, [], cfg(i).elem(elemtype == etypes(id), :), 'facealpha', 0.3, 'linestyle', 'none', 'facecolor', surfcolors(etypes(id) + 1, :), varargin{:});
        end
    end

    if (isfield(cfg(i), 'shapes'))
        if (isfield(cfg(i), 'vol'))
            hseg = mcxplotshapes(cfg(i).shapes, size(cfg(i).vol), offset, hseg, voxelsize, varargin{:});
        else
            hseg = mcxplotshapes(cfg(i).shapes, [60, 60, 60], offset, hseg, voxelsize, varargin{:});
        end
    end

    % rendering source position and direction
    if (~isfield(cfg(i), 'srcpos') || ~isfield(cfg(i), 'srcdir'))
        error('cfg.srcpos or cfg.srcdir is missing');
    end

    cfg(i).srcpos = cfg(i).srcpos;
    srcpos = cfg(i).srcpos * voxelsize;
    hsrc = plotmesh(srcpos, 'r*');
    srcvec = cfg(i).srcdir * 10 * voxelsize;
    headsize = 1e2;
    if (isoctavemesh)
        headsize = 0.5;
    end
    hdir = quiver3(srcpos(1), srcpos(2), srcpos(3), srcvec(1), srcvec(2), srcvec(3), 'linewidth', 3, 'color', 'r', 'MaxHeadSize', headsize, 'AutoScaleFactor', 1, varargin{:});

    % rendering area-source aperature
    hsrcarea = [];
    if (isfield(cfg(i), 'srctype'))
        if (strcmp(cfg(i).srctype, 'disk') || strcmp(cfg(i).srctype, 'gaussian') || strcmp(cfg(i).srctype, 'zgaussian') || strcmp(cfg(i).srctype, 'ring'))
            if (~isfield(cfg(i), 'srcparam1'))
                error('cfg.srcparam1 is missing');
            end
            [ncyl, fcyl] = meshacylinder(srcpos, srcpos + cfg(i).srcdir(1:3) * 1e-5, cfg(i).srcparam1(1) * voxelsize, 0, 0);
            hsrcarea = plotmesh(ncyl, fcyl{end - 1}, 'facecolor', 'r', 'linestyle', 'none');
            if (cfg(1).srcparam1(2) > 0)
                [ncyl, fcyl] = meshacylinder(srcpos, srcpos + cfg(i).srcdir(1:3) * 1e-5, cfg(i).srcparam1(2) * voxelsize, 0, 0);
                hsrcarea = plotmesh(ncyl, fcyl{end - 1}, 'facecolor', 'k', 'linestyle', 'none');
            end
        elseif (strcmp(cfg(i).srctype, 'planar') || strcmp(cfg(i).srctype, 'pattern') || strcmp(cfg(i).srctype, 'fourier') || ...
                strcmp(cfg(i).srctype, 'fourierx') || strcmp(cfg(i).srctype, 'fourierx2d') || strcmp(cfg(i).srctype, 'pencilarray'))
            if (~isfield(cfg(i), 'srcparam1') || ~isfield(cfg(i), 'srcparam2'))
                error('cfg.srcparam2 or cfg.srcparam2 is missing');
            end
            if (strcmp(cfg(i).srctype, 'fourierx') || strcmp(cfg(i).srctype, 'fourierx2d'))
                vec2 = cross(cfg(i).srcdir, cfg(i).srcparam1(1:3) * voxelsize);
            else
                vec2 = cfg(i).srcparam2(1:3) * voxelsize;
            end
            nrec = [0 0 0; cfg(i).srcparam1(1:3) * voxelsize; cfg(i).srcparam1(1:3) * voxelsize + vec2; vec2];
            hsrcarea = plotmesh(nrec + repmat(srcpos(:)', [4, 1]), {[1 2 3 4 1]});
        elseif (strcmp(cfg(i).srctype, 'pattern3d'))
            dim = cfg(i).srcparam1(1:3);
            [bbxno, bbxfc] = latticegrid(0:dim(1):dim(1), 0:dim(2):dim(2), 0:dim(3):dim(3));
            hbbx = plotmesh(((bbxno + repmat(cfg(i).srcpos(1:3), size(bbxno, 1), 1)) + offset) * voxelsize, bbxfc, 'facecolor', 'y', 'facealpha', 0.3);
        elseif (strcmp(cfg(i).srctype, 'slit') || strcmp(cfg(i).srctype, 'line'))
            if (~isfield(cfg(i), 'srcparam1'))
                error('cfg.srcparam1 is missing');
            end
            hsrcarea = plotmesh([srcpos(1:3); cfg(i).srcparam1(1:3) * voxelsize], [1 2], 'linewidth', 3, 'color', 'r');
        end
    end

    % rendering detectors
    hdet = [];
    if (isfield(cfg(i), 'detpos'))
        hdet = zeros(size(cfg(i).detpos, 1), 1);
        detpos = cfg(i).detpos(:, 1:3) * voxelsize;
        if (size(cfg(i).detpos, 2) == 4)
            detpos(:, 4) = cfg(i).detpos(:, 4) * voxelsize;
        else
            detpos(:, 4) = 1;
        end
        for id = 1:size(detpos, 1)
            [sx, sy, sz] = sphere;
            hdet(id) = surf(sx * detpos(id, 4) + (detpos(id, 1)), ...
                            sy * detpos(id, 4) + (detpos(id, 2)), ...
                            sz * detpos(id, 4) + (detpos(id, 3)), ...
                            'facealpha', 0.3, 'facecolor', 'g', 'linestyle', 'none');
        end
    end

    % combining all handles
    if (nargout > 0)
        hs{i} = struct('bbx', hbbx, 'seg', hseg, 'src', hsrc, 'srcarrow', hdir, 'srcarea', hsrcarea, 'det', hdet);
    end
end

if (~isholdplot)
    hold off;
end

rand ('state', rngstate);