% MRI3DPLOT - plot 3-D density image translucently on top of the mean MR 
%               brain image used in DIPPLOT. Plot brain slices in directions
%               'top' (axial), or 'side' (sagittal), or 'rear' (coronal).
%               Creates a new FIGURE. Smoothing uses Matlab SMOOTH3
% Usage:
%      >> [smoothed_3ddens, mriplanes] = mri3dplot(array3d, mri, 'key', 'val', ...);
%
% Input: 
%   array3d     - 3-D array to plot translucently on top of MRI image planes
%                  (e.g., as returned by DIPOLEDENSITY, unit: dipoles/cc).
%   mri         - [string or struct] base MR image structure (as returned by 
%                 dipoledensity.m or mri file (matlab format or file format read 
%                 by fcdc_read_mri. See dipplot.m help for more information.
%
% Optional inputs:
%   'mriview'   - ['top'|'side'|rear'] MR image slices to plot. 'axial',
%                    'coronal', and 'saggital' are also recognized keywords
%                    {default|[]: 'top'|'axial'}. May also be a cell array
%                    of such keyword (one per slice).
%   'mrislices' - [real vector] MR slice plotting coordinates in 'mriview' 
%                    direction {default|[]: -50 to +50 in steps of 11 mm}
%   'kernel'    - 3-D smoothing ht, width, & depth in voxels {default|0: none}
%                    > 0 -> uses gaussian kernel of std. dev. 0.65 voxels)
%   'geom'      - [rows, cols] geometry of slice array in figure output
%                    {default: smallest enclosing near-square array}
%   'rotate'    - [0|90|180|270 deg] 2-D slice image rotation {default: 90}
%   'cmap'      - [float array] colormap for plotting the 3-D array 
%                    {default: 'hot'}
%   'cmax'      - [float] color palette max value {default: array3d max}
%   'cmin'      - [float] color palette min value {default: 0}
%   'cbar'      - ['on'|'off'] plot colorbar. Default is 'on'.  
%   'subplot'   - ['on'|'off'] for single slice only, plot within a sub-plot
%                 panel. If 'on', this automatically sets 'cbar' to 'off'.
%                 Default is 'off'.  
%   'plotintersect' - ['on'|'off'] plot intersection between plotted slices.
%                 Default is 'on'.
%   'mixfact'   - [float] factor for mixing the background image with the
%                 array3d information. Default is 0.5.
%   'mixmode'   - ['add'|'overwrite'] 'add' will allow for transparency
%                 while 'overwrite' will preserve the original MRI image 
%                 and overwrite the pixel showind density changes.
%
% Outputs:
%   smoothed_3ddens - the plotted (optionally smoothed) 3-D density matrix
%   mriplanes       - cell array of the plotted MR image slices
%
% Author: Arnaud Delorme & Scott Makeig, SCCN, 10 June 2003
%
% Example:
%   dipfitdefs;
%   load('-mat', template_models(1).mrifile); % load mri variable
%   array = gauss3d(91,109,91);
%   mri3dplot(array, mri);
%
% See also: PLOTMRI

% Copyright (C) Arnaud Delorme, sccn, INC, UCSD, 2003-
% 03/29/2013 Makoto. Line 370 added to avoid negative matrix indices.
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

function [newprob3dori, mriplanes] = mri3dplot(prob3d, mri, varargin)

    % REDUCEPATCH  Reduce number of patch faces.

    if nargin < 1
        help mri3dplot;
        return;
    end

    DEFAULT_SPACING = 11;     % default MR slice interval (in mm)
    translucency    = 0.5;    % default alpha value
    mri_lim         = 85;     % +/- axis limits of MNI head image
    icadefs;
    g = finputcheck( varargin, { 'mriview'   { 'string','cell' 'real' }  { { 'sagital','axial','coronal','top','side','rear' } {} [] }   'top';
                        'mixmode'   'string'   { 'add','overwrite','min' }     'add';
                        'mrislices' 'float'    []                        [];
                        'view'      'float'    []                        [];
                        'geom'      'float'    []                        [];
                        'cmap'      'float'    []                        feval(DEFAULT_COLORMAP, 64);
                        'cmax'      'float'    []                        [];
                        'mixfact'   'float'    []                        0.5;
                        'cmin'      'float'    []                        0;
                        'plotintersect' 'string'   { 'on','off' }            'on';
                        'cbar'      'string'   { 'on','off' }            'on';
                        'subplot'   'string'   { 'on','off' }            'off';
                        'rotate'    'integer'  { 0,90,180,270 }          90;
                        'kernel'    'float'    []                        0; 
                        'addrow'    'integer'  []                        0;
                        'fighandle' ''         []                        []});
    if ischar(g), error(g); end
    if ischar(g.mriview) == 1, g.plotintersect = 'off'; end
    if strcmpi(g.mriview,'sagittal'),    g.mriview = 'side'; 
    elseif strcmpi(g.mriview,'axial'),   g.mriview = 'top'; 
    elseif strcmpi(g.mriview,'coronal'), g.mriview = 'rear';
    end
    if strcmpi(g.subplot, 'on') % plot colorbar
        g.cbar = 'off';
    end
    if ischar(mri)
         try
            mri = load('-mat', mri);
            mri = mri.mri;
         catch
            disp('Failed to read Matlab file. Attempt to read MRI file using function read_fcdc_mri');
            try
                warning off;
                mri = read_fcdc_mri(mri);
                mri.anatomy = round(gammacorrection( mri.anatomy, 0.8));
                mri.anatomy = uint8(round(mri.anatomy/max(reshape(mri.anatomy, prod(mri.dim),1))*255));
                % WARNING: if using double instead of int8, the scaling is different 
                % [-128 to 128 and 0 is not good]
                % WARNING: the transform matrix is not 1, 1, 1 on the diagonal, some slices may be 
                % misplaced
                warning on;
            catch
                error('Cannot load file using read_fcdc_mri');
            end
         end
    end
    
    % normalize prob3d for 1 to ncolors and create 3-D dim
    % ----------------------------------------------------
    if ~iscell(prob3d), prob3d = { prob3d }; end
    if length(prob3d) > 1
        if isempty(g.cmax), g.cmax = max(max(prob3d{1}(:)),max(prob3d{2}(:))); end
        [newprob3d{1}] = prepare_dens(prob3d{1}, g, 'abscolor');
        [newprob3d{2}] = prepare_dens(prob3d{2}, g, 'abscolor');
    else
        if isempty(g.cmax), g.cmax = max(prob3d{1}(:)); end
        [newprob3d{1}, maxdens1] = prepare_dens(prob3d{1}, g, 'usecmap');
    end
    fprintf('Brightest color denotes a density of: %1.6f (presumed unit: dipoles/cc)\n', g.cmax);
    
    % plot MRI slices
    % ---------------
    if isempty(g.mrislices), 
        g.mrislices = linspace(-50, 50, DEFAULT_SPACING); 
    end
    
    if strcmpi(g.cbar, 'on'), add1 = 1; else add1 = 0; end
    if isempty(g.geom)
        g.geom = ceil(sqrt(length(g.mrislices)+add1)); 
        g.geom(2) = ceil((length(g.mrislices)+add1)/g.geom)+g.addrow;
    end

    if strcmpi(g.subplot, 'off')
        if isempty(g.fighandle)
             fig = figure;
        else
             fig = g.fighandle;
             clf(fig);
        end
        
        pos = get(fig, 'position');
        set(fig, 'position', [ pos(1)+15 pos(2)+15 pos(3)/4*g.geom(1) pos(4)/3*g.geom(2) ]);
    end
    
    disp('Plotting...');
    
    % custom view for each slice
    if ~iscell( g.mriview )
        g.mriview = { g.mriview };
        g.mriview(2:length( g.mrislices )) = g.mriview(1);
    end
    
    newprob3dori = newprob3d;
    for index = 1:length( g.mrislices ) %%%%%%% for each plotted MR image slice %%%%%%%%

        % plot intersection between plotted slices
        % ----------------------------------------
        newprob3d = newprob3dori;
        if strcmpi(g.plotintersect, 'on')
            for index2 = setdiff_bc(1:length( g.mrislices ), index)
                switch g.mriview{index2}
                    case 'side', coord = [  g.mrislices(index2) 0 0 1 ]; 
                    case 'top' , coord = [  0 0 g.mrislices(index2) 1 ]; 
                    case 'rear', coord = [  0 g.mrislices(index2) 0 1 ]; 
                end
                coord = round( pinv(mri.transform)*coord' )';
                for i = 1:length(newprob3d)
                    switch g.mriview{index2}
                     case 'side', newprob3d{i}(  coord(1), :, :, :) = 0;
                     case 'top' , newprob3d{i}(  :, :, coord(3), :) = 0;
                     case 'rear', newprob3d{i}(  :, coord(2), :, :) = 0;
                    end
                end
            end
        end

        % create axis if necessary
        % ------------------------
        if strcmpi(g.subplot, 'off')
            mysubplot(g.geom(1), g.geom(2), index); % get an image slice axis
        end
        
        % find coordinate
        % ---------------
        switch g.mriview{index}
         case 'side', coord = [  g.mrislices(index) 0 0 1 ]; 
         case 'top' , coord = [  0 0 g.mrislices(index) 1 ]; 
         case 'rear', coord = [  0 g.mrislices(index) 0 1 ]; 
        end
        
        coord = round( pinv(mri.transform)*coord' )';
        
        % get MRI slice
        % -------------
        switch g.mriview{index}
         case 'side', mriplot  = squeeze( mri.anatomy(coord(1), :, :) );
         case 'top' , mriplot  = squeeze( mri.anatomy(:, :, coord(3)) );
         case 'rear', mriplot  = squeeze( mri.anatomy(:, coord(2), :) );
        end

        mriplot(:,:,2) = mriplot(:,:,1);
        mriplot(:,:,3) = mriplot(:,:,1);
        mriplot = rotatemat( mriplot, g.rotate );
        
        % get dipole density slice
        % ------------------------
        for i = 1:length(newprob3d)
            switch g.mriview{index}
             case 'side', densplot{i} = squeeze( newprob3d{i}(coord(1), :, :, :) );
             case 'top' , densplot{i} = squeeze( newprob3d{i}(:, :, coord(3), :) );
             case 'rear', densplot{i} = squeeze( newprob3d{i}(:, coord(2), :, :) );
            end
            densplot{i} = rotatemat( densplot{i}, g.rotate );
        end
  
        if isa(mriplot, 'uint8')
            % check if densplot is in uint8
            % if not - transform to uint8
            for dlen = 1:length(densplot)
                if ~isa(densplot{dlen}, 'uint8')
                    % check if multiply by 255 (when double
                    % with values ranging from 0 - 1) and
                    % then transform to uint8
                    checkval = densplot{dlen} >= 0 & densplot{dlen} <= 1;
                    checkval = sum(checkval(:)) == numel(densplot{dlen});
                    if checkval
                        densplot{dlen} = uint8(densplot{dlen} * 255); %#ok<AGROW>
                        continue
                    end
                    
                    % check if it's ok to transform
                    % straight to int8
                    checkval = densplot{dlen} >= 0 & densplot{dlen} <= 255;
                    checkval = sum(checkval(:)) == numel(densplot);
                    testint = isequal(densplot{dlen}, round(densplot{dlen}));
                    
                    if checkval && testint
                        densplot{dlen} = uint8(densplot{dlen}); %#ok<AGROW>
                    end
                    
                end
            end
            clear dlen checkval testint
        end
        
        if length(densplot) == 1
            densplot = densplot{1};
            if strcmpi(g.mixmode, 'add')
                densplot(isnan(densplot)) = 0; % do not plot colors outside the brain, as indicated by NaNs
                mriplot  = mriplot*g.mixfact + densplot*(1-g.mixfact); % Mix 1/2 MR image + 1/2 density image
            elseif strcmpi(g.mixmode, 'overwrite')
                indsnon0 = sum(densplot(:,:,:),3) > 0;
                tmpmri = mriplot(:,:,1); tmpdens = densplot(:,:,1); tmpmri(indsnon0) = tmpdens(indsnon0); mriplot(:,:,1) = tmpmri;
                tmpmri = mriplot(:,:,2); tmpdens = densplot(:,:,2); tmpmri(indsnon0) = tmpdens(indsnon0); mriplot(:,:,2) = tmpmri;
                tmpmri = mriplot(:,:,3); tmpdens = densplot(:,:,3); tmpmri(indsnon0) = tmpdens(indsnon0); mriplot(:,:,3) = tmpmri;
            elseif strcmpi(g.mixmode, 'min')
                densplot(isnan(densplot)) = 0; % do not plot colors outside the brain, as indicated by NaNs
                mriplot  = min(mriplot, densplot); % min
            end
            clear densplot;
        else
            densplot{1}(isnan(densplot{1})) = 0; % do not plot colors outside the brain, as indicated by NaNs
            densplot{2}(isnan(densplot{2})) = 0; % do not plot colors outside the brain, as indicated by NaNs
            if strcmpi(g.mixmode, 'add')
                mriplot(:,:,1) = mriplot(:,:,1)*g.mixfact + densplot{1}(:,:)*(1-g.mixfact); % min
                mriplot(:,:,3) = mriplot(:,:,3)*g.mixfact + densplot{2}(:,:)*(1-g.mixfact); % min
                mriplot(:,:,2) = mriplot(:,:,2)*g.mixfact; % min
            elseif strcmpi(g.mixmode, 'overwrite')
                indsnon01 = densplot{1}(:,:,:) > 0;
                indsnon02 = densplot{2}(:,:,:) > 0;
                tmpmri = mriplot(:,:,1); tmpdens = densplot{1}(:,:); tmpmri(indsnon01) = tmpdens(indsnon01); mriplot(:,:,1) = tmpmri;
                tmpmri = mriplot(:,:,2); tmpdens = densplot{2}(:,:); tmpmri(indsnon02) = tmpdens(indsnon02); mriplot(:,:,2) = tmpmri;
            else
                mriplot(:,:,1) = max(mriplot(:,:,1), densplot{1}); % min
                mriplot(:,:,2) = max(mriplot(:,:,2), densplot{2}); % min
            end
        end
        mriplanes{index} = mriplot;
        
        imagesc(mriplot); % plot [background MR image + density image]

        axis off; hold on;

        xl = xlim;
        yl = ylim;
        zl = zlim;
        
        % options = { 'FaceColor','texturemap', 'EdgeColor','none', 'CDataMapping', ...
        %             'scaled','facelighting', 'none', 'facealpha', translucency};
        % h = surface( [xl(1) xl(2); xl(1) xl(2)], [yl(1) yl(1); yl(2) yl(2)], ...
        %              [1 1; 1 1], densplot, options{:});

        axis equal;
        if ~isempty(g.view), view(g.view); end
        tit = title( [ int2str(g.mrislices(index)) ' mm' ]);
        if strcmpi(g.subplot, 'off'), set(tit, 'color', 'w'); end
    end
    
    % plot colorbar
    % -------------
    if strcmpi(g.cbar, 'on') % plot colorbar

        h = mysubplot(g.geom(1), g.geom(2), length(g.mrislices)+1);
        pos = get(h, 'position');
        pos(1) = pos(1)+pos(3)/3;
        pos(3) = pos(3)/6;
        pos(2) = pos(2)+pos(4)/5;
        pos(4) = pos(4)*3/5;
        axis off;
        h = axes('unit', 'normalized', 'position', pos); % position the colorbar

        if strcmpi(g.mixmode, 'add')
            tmpmap = g.cmap/2 + ones(size(g.cmap))/4;  % restrict range to [1/4, 3/4] of cmap
        else
            tmpmap = g.cmap; %g.cmap/2 + ones(size(g.cmap))/4;  % restrict range to [1/4, 3/4] of cmap
        end
        colormap(tmpmap);
        cbar(h, [1:length(g.cmap)], [g.cmin g.cmax]);
        box off;
        set(h, 'ycolor', [0.7 0.7 0.7]);
    end
    
    fprintf('\n');
    if exist('fig') == 1 
        if length(prob3d) == 1
            set(fig,'color', g.cmap(1,:)/2);
        else
            set(fig,'color', [0.0471 0.0471 0.0471]/1.3);
        end
    end
return;

function mat = rotatemat(mat, angle);

    if angle == 0, return; end
    if ndims(mat) == 2
        if angle >= 90,  mat = rot90(mat); end
        if angle >= 180, mat = rot90(mat); end
        if angle >= 270, mat = rot90(mat); end
    else
        if angle >= 90,
            newmat(:,:,1) = rot90(mat(:,:,1));
            newmat(:,:,2) = rot90(mat(:,:,2));
            newmat(:,:,3) = rot90(mat(:,:,3));
            mat = newmat;
        end
        if angle >= 180,
            newmat(:,:,1) = rot90(mat(:,:,1));
            newmat(:,:,2) = rot90(mat(:,:,2));
            newmat(:,:,3) = rot90(mat(:,:,3));
            mat = newmat;
        end
        if angle >= 270,
            newmat(:,:,1) = rot90(mat(:,:,1));
            newmat(:,:,2) = rot90(mat(:,:,2));
            newmat(:,:,3) = rot90(mat(:,:,3));
            mat = newmat;
        end
    end
    
function [newprob3d maxdens] = prepare_dens(prob3d, g, col);

    if g.kernel ~= 0
        disp('Smoothing...');
        smoothprob3d    = smooth3(prob3d, 'gaussian', g.kernel);
        prob3d          = smoothprob3d;
    end

    maxdens = max(prob3d(:));
    ncolors = size(g.cmap,1);
    
    prob3d    = round((prob3d-g.cmin)/(g.cmax - g.cmin)*(ncolors-1))+1; % project density image into the color space: [1:ncolors]
    prob3d( find(prob3d > ncolors) ) = ncolors;
    prob3d( find(prob3d < 1))        = 1; % added by Makoto
    newprob3d = zeros(size(prob3d,1), size(prob3d,2), size(prob3d,3), 3);

    outOfBrainMask = find(isnan(prob3d)); % place NaNs in a mask, NaNs are assumed for points outside the brain
    prob3d(outOfBrainMask) = 1;
    
    if strcmpi(col, 'abscolor')
        newprob3d = prob3d/ncolors;
    else
        tmp = g.cmap(prob3d,1); newprob3d(:,:,:,1) = reshape(tmp, size(prob3d));
        tmp = g.cmap(prob3d,2); newprob3d(:,:,:,2) = reshape(tmp, size(prob3d));
        tmp = g.cmap(prob3d,3); newprob3d(:,:,:,3) = reshape(tmp, size(prob3d));
    end
    
function h = mysubplot(geom1, geom2, coord);
    
    coord = coord-1;
    horiz_border = 0;
    vert_border  = 0.1;
    
    coordy = floor(coord/geom1);
    coordx = coord - coordy*geom1;
    
    posx   = coordx/geom1+horiz_border*1/geom1/2;
    posy   = 1-(coordy/geom2+vert_border*1/geom2/2)-1/geom2;
    width  = 1/geom1*(1-horiz_border);
    height = 1/geom2*(1- vert_border);
    
    h = axes('unit', 'normalized', 'position', [ posx posy width height ]);
    %h = axes('unit', 'normalized', 'position', [ coordx/geom1 1-coordy/geom2-1/geom2 1/geom1 1/geom2 ]);
    
