function varargout = bst_mriviewer(varargin)
%BST_MRIVIEWER - Application M-file for bst_mriviewer.fig
% function varargout = bst_mriviewer(varargin)
%    FIG = BST_MRIVIEWER launch bst_mriviewer GUI.
%    BST_MRIVIEWER('callback_name', ...) invoke the named callback.

%<autobegin> ---------------------- 23-Sep-2004 10:29:59 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
%   publictoolbox\mri_toolbox_v2p0\avw_img_read.m
%   toolbox\avw2brainstorm.m
%   toolbox\bst_color_scheme.m
%   toolbox\bst_layout.m
%   toolbox\bst_message_window.m
%   toolbox\bst_mriviewer.m  NOTE: Routine calls itself explicitly
%   toolbox\bst_wavedata_display.m
%   toolbox\file_selection_win.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%   toolbox\load_bstmri.m
%   toolbox\makeuswait.m
%   toolbox\modulo.m
%   toolbox\mri2scs.m
%   toolbox\mri_parameters.m
%   toolbox\pcs_params.m
%   toolbox\read_ctf_mri.m
%   toolbox\save_fieldnames.m
%   toolbox\scs2mri.m
%
% Subfunctions in this file, in order of occurrence in file:
%   [] = clearAll(fig);
%   varargout = slider_Callback(h, eventdata, handles, varargin)
%   XYZ = wherearewe(handles)
%   varargout = LoadBrainStorm_Callback(h, eventdata, handles, varargin)
%   varargout = ctfmri_Callback(h, eventdata, handles, varargin)
%   varargout = fifmri_Callback(h, eventdata, handles, varargin)
%   varargout = importANALYZE_Callback(h, eventdata, handles, varargin)
%   varargout = gismri_Callback(h, eventdata, handles, varargin)
%   varargout = dicommri_Callback(h, eventdata, handles, varargin)
%   varargout = rawmri_Callback(h, eventdata, handles, varargin)
%   varargout = rawmri_done
%   varargout = VolumeParameters_Callback(h, eventdata, handles, varargin)
%   varargout = SaveAsBrainStorm_Callback(h, eventdata, handles, varargin)
%   varargout = DefineSCS_Callback(h, eventdata, handles, varargin)
%   varargout = DefineSCSdone_Callback(h, eventdata, handles, varargin)
%   varargout = ShowFiducials_Callback(h, eventdata, handles, varargin)
%   [] = isvisible(MRI,hObjects)
%   varargout = LoadSourceMap_Callback(h, eventdata, handles, varargin)
%   varargout = SaveScouts_Callback(h, eventdata, handles, varargin)
%   varargout = LoadScouts_Callback(h, eventdata, handles, varargin)
%   varargout = cortic_remap_Callback(h, eventdata, handles, varargin)
%   varargout = ImportSPM_Callback(h, eventdata, handles, varargin)
%   varargout = Talairach_Callback(h, eventdata, handles, varargin)
%   varargout = TalairachTransform_Callback(h, eventdata, handles, varargin)
%   varargout = DefineScouts_Callback(h, eventdata, handles, varargin)
%   varargout = DefineCorticalScouts_Callback(h, eventdata, handles, varargin)
%   varargout = Quit_Callback(h, eventdata, handles, varargin)
%   [ObjectHandles] = plot_MRI(MRI,pointsCoord,varargin);
%   [MRI] = show_MRI(MRI, varargin);
%   [] = reset_axis(H)
%   MRI = FlipUpDown
%   MRI = FlipLeftRight
%   MRI = Transpose
%   mriCoord = vox2mm(MRI,sliceIndices);
%   sliceIndices = mm2vox(MRI,mriCoord);
%   varargout = LockCursor_Callback(h, eventdata, handles, varargin)
%   [] = stopmoving
%   [] = move2cursor
%   varargout = AddCrosshair_Callback(h, eventdata, handles, varargin)
%   varargout = Zoom_Callback(h, eventdata, handles, varargin)
%   varargout = PointList_Callback(h, eventdata, handles, varargin)
%   varargout = SelectPoint_Callback(h, eventdata, handles, varargin)
%   varargout = DeletePoint_Callback(h, eventdata, handles, varargin)
%   varargout = RenamePoint_Callback(h, eventdata, handles, varargin)
%   varargout = SeeAllPoints_Callback(h, eventdata, handles, varargin)
%   varargout = HideAllPoints_Callback(h, eventdata, handles, varargin)
%   varargout = Brighten_Callback(h, eventdata, handles, varargin)
%   varargout = Darken_Callback(h, eventdata, handles, varargin)
%   varargout = LittleEndian_Callback(h, eventdata, handles, varargin)
%
% Application data and their calls in this file:
%   'MRI'
%   'TileType'
%   'mriResults'
%   'tag'
%   
%   setappdata(MRIViewer,'MRI',MRI);
%   setappdata(MRIfig,'MRI',MRI);
%   setappdata(fig,'MRI',MRI);
%   setappdata(fig,'TileType','T')
%   setappdata(findobj(0,'tag','MRIViewer','type','figure'),'MRI',MRI)
%   setappdata(gcbf,'MRI',MRI)
%   setappdata(gcbf,'MRI',MRI);
%   setappdata(handles.MRIViewer,'MRI',MRI)
%   setappdata(handles.MRIViewer,'MRI',MRI);
%   setappdata(mriFig,'MRI',MRI)
%   setappdata(mriFig,'mriResults',mriResults)
%   setappdata(timeseriesFig,'TileType','T');
%   setappdata(topoFig,'TileType','T');
%   
%   MRI = getappdata(MRIfig,'MRI');
%   MRI = getappdata(fig,'MRI');
%   MRI = getappdata(gcbf,'MRI');
%   MRI = getappdata(handles.MRIViewer,'MRI');
%   MRI = getappdata(mriFig,'MRI');
%   MRI.SCS2Landmarks = getappdata(findobj(0,'Tag','pcsparams'),'LookUpFiducials');
%   allData = getappdata(fig);
%   mriResults = getappdata(mriFig,'mriResults');
%   tmpMRI = getappdata(handles.MRIViewer,'tmpMRI');
%
% Figure Files opened by this function:
%   mfilename
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure bst_mriviewer.fig
%   figure::MRIViewer "" uses DeleteFcn for bst_mriviewer('Quit_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:listbox:PointList "" uses Callback for <automatic>
%   uicontrol:pushbutton:Brighten "Brighten" uses Callback for <automatic>
%   uicontrol:pushbutton:Darken "Darken" uses Callback for <automatic>
%   uicontrol:pushbutton:DeletePoint "Delete" uses Callback for <automatic>
%   uicontrol:pushbutton:Quit "Quit" uses Callback for <automatic>
%   uicontrol:pushbutton:RenamePoint "Rename" uses Callback for <automatic>
%   uicontrol:pushbutton:SelectPoint "Select" uses Callback for <automatic>
%   uicontrol:slider:slider_x "" uses Callback for
%     bst_mriviewer('slider_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:slider:slider_y "" uses Callback for
%     bst_mriviewer('slider_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:slider:slider_z "" uses Callback for
%     bst_mriviewer('slider_Callback',gcbo,[],guidata(gcbo))
%   uimenu::AddCrosshair "" uses Callback for <automatic>
%   uimenu::cortic_remap "" uses Callback for <automatic>
%   uimenu::ctfmri "" uses Callback for <automatic>
%   uimenu::DefineCorticalScouts "" uses Callback for <automatic>
%   uimenu::DefineSCS "" uses Callback for <automatic>
%   uimenu::dicommri "" uses Callback for <automatic>
%   uimenu::fifmri "" uses Callback for <automatic>
%   uimenu::gismri "" uses Callback for <automatic>
%   uimenu::HideAllPoints "" uses Callback for <automatic>
%   uimenu::importANALYZE "" uses Callback for <automatic>
%   uimenu::ImportSPM "" uses Callback for <automatic>
%   uimenu::LittleEndian "" uses Callback for
%     bst_mriviewer('LIttleEndian_Callback',gcbo,[],guidata(gcbo))
%   uimenu::LoadBrainStorm "" uses Callback for <automatic>
%   uimenu::LoadScouts "" uses Callback for <automatic>
%   uimenu::LoadSourceMap "" uses Callback for <automatic>
%   uimenu::LockCursor "" uses Callback for <automatic>
%   uimenu::rawmri "" uses Callback for <automatic>
%   uimenu::SaveAsBrainStorm "" uses Callback for <automatic>
%   uimenu::SaveScouts "" uses Callback for <automatic>
%   uimenu::SeeAllPoints "" uses Callback for <automatic>
%   uimenu::ShowFiducials "" uses Callback for <automatic>
%   uimenu::TalairachTransform "" uses Callback for <automatic>
%   uimenu::Untitled_1 "" uses Callback for colormapeditor
%   uimenu::VolumeParameters "" uses Callback for <automatic>
%   uimenu::Zoom "" uses Callback for <automatic>
%
% At Check-in: $Author: Mosher $  $Revision: 61 $  $Date: 9/23/04 8:58a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 23-Sep-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 23-Sep-2004 10:29:59 -----------------------


% Last Modified by GUIDE v2.0 24-May-2004 17:32:33


% /---Script Authors-------------------------------------\
% |  *** Sylvain Baillet, Ph.D.                          |
% |  Cognitive Neuroscience & Brain Imaging Laboratory   |
% |  CNRS UPR640 - LENA                                  | 
% |  Hopital de la Salpetriere, Paris, France            |
% |  sylvain.baillet@chups.jussieu.fr                    |
% |                                                      |
% | John C. Mosher, Ph.D. 
% \------------------------------------------------------/
%  
% Script History ---------------------------------------------------------------
% SB  10-Mar-2004  Creation (BrainStorm v.1 update)
% JCM 26-Apr-2004  Comments cleanup, making SCS more automated for outside
%                  calling.
% JCM 25-May-2004  Fixed listing of the solution dipole numbers in the listbox
% JCM 16-Jun-2004  in launch GUI, always call bst_layout, since openfig('reuse')
%                  appears to shift location
% JCM 17-Jun-2004  Added menubar 'figure' to the Timeseries and Topography
%                  figures
% Script History ---------------------------------------------------------------

if nargin == 0  % LAUNCH GUI
   
   fig = openfig(mfilename,'reuse');
   
   % Use BST color scheme for figure:
   setappdata(fig,'TileType','T')
   bst_color_scheme(fig)
   bst_layout('align',fig,1,2,1);
   set(fig,'visible','on')
   
   
    % DoubleBuffer avoids flickering
    set(fig,'DoubleBuffer','on')
    
    % Generate a structure of handles to pass to callbacks, and store it. 
    handles = guihandles(fig);
    guidata(fig, handles);
    
    
    
    % Basic diplay features
    % Assign proper look to viewing axis
    reset_axis([handles.axc, handles.axs, handles.axa]);
    
    if nargout > 0
        varargout{1} = fig;
    end
    
elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK
    
    try
        if (nargout)
            [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
        else
            feval(varargin{:}); % FEVAL switchyard
        end
    catch
        disp(lasterr);
    end
    
end


% --------------------------------------------------------------------
function [] = clearAll(fig);
% Delete all appdata linked to fingure fig
allData = getappdata(fig);
allDataFields = fieldnames(allData);
for k = 1: length(allDataFields)
    if ~strcmp(lower(allDataFields{k}),'tiletype')
        rmappdata(fig,allDataFields{k});
    end
end
h = guihandles(fig);
set(h.LockCursor,'Checked','off')
set(h.Zoom,'Checked','off')

% --------------------------------------------------------------------
function varargout = slider_Callback(h, eventdata, handles, varargin)
% Browse through slices

MRI = getappdata(gcbf,'MRI');

% where are we now ?
newIndices = mm2vox(MRI,wherearewe(handles)); % coordinates in Voxels
%newIndices(newIndices == 0) = 1;

MRI = show_MRI(MRI,'update',newIndices);

% % Update crosshair location if necessary
% if strcmp(get(findobj(gcbf,'tag','AddCrosshair'),'checked'),'on')
%     bst_mriviewer('AddCrosshair_Callback',findobj(gcbf,'tag','AddCrosshair'),'',handles,1); 
% end

% --------------------------------------------------------------------
function XYZ = wherearewe(handles)
%Return MRI coordinates (in mm) of current point in volume

XYZ(1) = get(handles.slider_x,'Value');
XYZ(2) = get(handles.slider_y,'Value');
XYZ(3) = get(handles.slider_z,'Value');

% --------------------------------------------------------------------
function varargout = LoadBrainStorm_Callback(h, eventdata, handles, varargin)
% Load a MRI volume in BrainStorm format

gcbf = findobj(0,'Tag','MRIViewer','Type','figure');
clearAll(gcbf) % Clear all appdata currently stored in figure

if nargin == 3 % No file specified - let user choose
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    % Select file
    [filebs, pathbs] = uigetfile('*subjectimage*.mat','Select a .mat file containing the subject''s MRI.');
    if filebs==0
        return
    end
    cd(pathbs)
    filebs = fullfile(pathbs,filebs);
else
    filebs = varargin{1};
end

% Load it
makeuswait('start')
MRI = load_BsTMRI(filebs);

% Show it
MRI = show_MRI(MRI);

% Store it 
setappdata(gcbf,'MRI',MRI) 

if isfield(MRI,'PCS') % Older MRI brainstorm file
    bst_message_window('wrap',{...
            'Older BrainStorm file format detected',...
            'This is an older BrainStorm MRI file. You may need to re-define the subject''s coordinate system.',...
            'Transforming image cube to new BrainStorm coordinates. . .'})
    set(gcbf,'CurrentAxes',handles.axa)
    MRI = FlipUpDown;
    set(gcbf,'CurrentAxes',handles.axc)
    MRI = FlipLeftRight;
    
    % Fill out new SCS structure
    MRI.SCS.FiducialName = MRI.PCS.FiducialName;
    
    MRI.SCS.CubeFiducial = MRI.PCS.CubeFiducial([2 1 3],:);
    MRI.SCS.System = MRI.PCS.Comment; % Type of SCS. 
    
    MRI.SCS.CubeFiducial(2,:) = size(MRI.Cube,2) + 1 - MRI.SCS.CubeFiducial(2,:);
    MRI.SCS.CubeFiducial(3,:) = size(MRI.Cube,3) + 1 - MRI.SCS.CubeFiducial(3,:);
    
    % Move to mm coordinates
    for k= 1:size(MRI.SCS.CubeFiducial,2)
        MRI.SCS.mmCubeFiducial(:,k) = vox2mm(MRI,MRI.SCS.CubeFiducial(:,k)');
    end
    
    MRI = rmfield(MRI,'PCS');
    
    % Now compute R,T transform 
    [transf] = mri2scs(MRI);
    MRI.SCS.R = transf.R;
    MRI.SCS.T = transf.T;
    MRI.SCS.Origin = transf.Origin;
    
    setappdata(gcbf,'MRI',MRI) 
    
    set(findobj(0,'tag','ShowFiducials'),'enable','on')        
    bst_mriviewer('SaveAsBrainStorm_Callback','','',handles)
    
    bst_message_window('wrap',{...
            '->Done'})
end

% Enable to Show fiducials if they exist
if isfield(MRI,'SCS')
    % Store to Landmark strucure
    if isfield( MRI,'Landmarks')
        MRI.Landmarks.Names = {MRI.SCS.FiducialName{:},MRI.Landmarks.Names{:}};
        MRI.Landmarks.MRImmXYZ = [MRI.SCS.mmCubeFiducial,MRI.Landmarks.MRImmXYZ];
    else
        MRI.Landmarks.Names = MRI.SCS.FiducialName;
        MRI.Landmarks.MRImmXYZ = [MRI.SCS.mmCubeFiducial];
    end
    
    MRI.SCS2Landmarks = [1 2 3 4]; % SCS fiducials are first 3 landmark points 
    % Draw graphic object for each fiducial
    %setappdata(gcbf,'MRI',MRI) 
    
    MRI.Landmarks.Handles = plot_MRI(MRI,MRI.Landmarks.MRImmXYZ);
    
    set(findobj(0,'tag','ShowFiducials'),'enable','on','checked','off')        
    
    setappdata(gcbf,'MRI',MRI) 
    bst_mriviewer('ShowFiducials_Callback','','','')
else
    set(findobj(0,'tag','ShowFiducials'),'enable','off')        
end


% That's it
makeuswait('stop')


% --------------------------------------------------------------------
function varargout = ctfmri_Callback(h, eventdata, handles, varargin)
% Proper import of CTF .mri file format 

% Machine format can be a problem here 
% (especially when importing from other plateforms)

clearAll(gcbf) % Clear all appdata currently stored in figure

% Is the MRIParameters window open ?
MRIParamFig = findobj(0,'tag','MRIParameters','type','figure');
if isempty(MRIParamFig) % No - open it
    mri_parameters('Create')
    % Now here is the figure
    MRIParamFig = findobj(0,'tag','MRIParameters','type','figure');
end

MRIParamHandles = guihandles(MRIParamFig);
set(MRIParamHandles.DataFormat,'Value',2); % CTF scans are sagittal
switch get(MRIParamHandles.MachineFormat,'value')
    case 1 % Big Endian
        MACHINEFORMAT = 'b';
    case 2 % Little Endian
        MACHINEFORMAT = 'l';
end

Current = get_user_directory;
cd(Current.SUBJECTS)
[filebs, pathbs] = uigetfile('*.mri','Please select a CTF .mri file');
if filebs==0,    return,end
cd(pathbs); 

makeuswait('start')
bst_message_window('Reading CTF MRI file. . .')
[MRI.Cube, Header] = read_ctf_mri(filebs,MACHINEFORMAT);

bst_message_window({...
        char(Header.identifierString),...
        'Reading CTF MRI file -> DONE'})

MRI.Voxsize = [Header.mmPerPixel_sagittal Header.mmPerPixel_coronal Header.mmPerPixel_axial];

%permVect: Permutation vector to apply to the original cube
% and associated parameters to get the MR volume in the proper 
% orientation of the MRViwer and BrainStorm format.
permVect = [3,1,2]; 

% reArrange voxel location to move to BrainStorm MRI coordinate system
MRI.Cube = permute(MRI.Cube,permVect);
MRI.Voxsize = MRI.Voxsize(permVect);

MRI.Cube = MRI.Cube(:,end:-1:1,end:-1:1);

MRI.header = 0; % Now for BrainStorm format
tmpString = {'uchar','ushort'};
MRI.DataClass = tmpString{Header.dataSize};
MRI.scantype = 'sagittal';
% Show it
MRI = show_MRI(MRI);

setappdata(findobj(0,'tag','MRIViewer','type','figure'),'MRI',MRI)
mri_parameters('Create','','','', findobj(0,'tag','MRIViewer','type','figure'))


% If defined, fill-out SCS information
try %... to read the headmodel information
    if isfield(Header,'HeadModel_Info') % store in BrainStorm format
        bst_message_window('Importing fiducial locations. . .')
        MRI.SCS.FiducialName = {'Nasion'  'LeftEar'  'RightEar'};
        MRI.SCS.CubeFiducial = [Header.HeadModel_Info.Nasion_Sag,Header.HeadModel_Info.LeftEar_Sag,Header.HeadModel_Info.RightEar_Sag;...
                Header.HeadModel_Info.Nasion_Cor,Header.HeadModel_Info.LeftEar_Cor,Header.HeadModel_Info.RightEar_Cor;...
                Header.HeadModel_Info.Nasion_Axi,Header.HeadModel_Info.LeftEar_Axi,Header.HeadModel_Info.RightEar_Axi...
            ];
        MRI.SCS.System = 'CTF'; % Type of SCS. 
        
        MRI.SCS.CubeFiducial(2,:) = size(MRI.Cube,2) + 1 - MRI.SCS.CubeFiducial(2,:);
        MRI.SCS.CubeFiducial(3,:) = size(MRI.Cube,3) + 1 - MRI.SCS.CubeFiducial(3,:);
        
        % Move to mm coordinates
        for k= 1:size(MRI.SCS.CubeFiducial,2)
            MRI.SCS.mmCubeFiducial(:,k) = vox2mm(MRI,MRI.SCS.CubeFiducial(:,k)');
        end
        
        % Store to Landmark strucure
        if isfield( MRI,'Landmarks')
            MRI.Landmarks.Names = {MRI.SCS.FiducialName{:},MRI.Landmarks.Names{:}};
            MRI.Landmarks.MRImmXYZ = [MRI.SCS.mmCubeFiducial,MRI.Landmarks.MRImmXYZ];
        else
            MRI.Landmarks.Names = MRI.SCS.FiducialName;
            MRI.Landmarks.MRImmXYZ = [MRI.SCS.mmCubeFiducial];
        end
        
        MRI.SCS2Landmarks = 1:length(MRI.SCS.FiducialName); % SCS fiducials are first 3 landmark points 
        % Draw graphic object for each fiducial
        %setappdata(gcbf,'MRI',MRI) 
        
        MRI.Landmarks.Handles = plot_MRI(MRI,MRI.Landmarks.MRImmXYZ);
        
        
        set(findobj(0,'tag','ShowFiducials'),'enable','on')        
    end
    
    setappdata(findobj(0,'tag','MRIViewer','type','figure'),'MRI',MRI)
    
    set(findobj(0,'Tag','ShowFiducials'),'checked','off','enable','on')
    % Show/Hide fiducials in the MR volume 
    ShowFiducials_Callback([],[], handles)
    DefineSCS_Callback([], [], handles)
    
catch
    % do nothing    
    set(findobj(0,'tag','ShowFiducials'),'enable','off')        
end

setappdata(findobj(0,'tag','MRIViewer','type','figure'),'MRI',MRI)

bst_message_window({...
        'That''s it.'})


makeuswait('stop')


% --------------------------------------------------------------------
function varargout = fifmri_Callback(h, eventdata, handles, varargin)



% --------------------------------------------------------------------
function varargout = importANALYZE_Callback(h, eventdata, handles, varargin)

% Which machine format is this ?
handles = guihandles(gcbf);
switch (get(handles.LittleEndian,'checked'))
    case 'on' % Little endian (default)
        MACHINEFORMAT = 'ieee-le';
    otherwise
        MACHINEFORMAT = 'ieee-be'; % Big_endian
end

set(0,'CurrentFigure',handles.MRIViewer)

Current = get_user_directory;
cd(Current.SUBJECTS)   

[anatfile, anatpath] = uigetfile('*.img','Please select ANALYZE MRI file (.img)' );    
if (anatfile == 0), return, end
cd(anatpath);

% Call D. Weber's tool
makeuswait('start')
[ avw, machine ] = avw_img_read(strtok(anatfile,'.'),'',MACHINEFORMAT);

% Convert to BrainStorm format
MRI = avw2brainstorm(avw, machine);

% Show it
MRI = show_MRI(MRI);

setappdata(handles.MRIViewer,'MRI',MRI)

makeuswait('stop')

%-----------------------------------------------------------------------------------------------

function varargout = gismri_Callback(h, eventdata, handles, varargin)
% Import GIS format (an MRI file format for better interchange with BrainVisa software [http://brainvisa.info])

Users = get_user_directory;

% Look for .ima extension of GIS file pair (one .dim and one .ima)
[Files, Filenames] = file_selection_win({'gis image'},Users.SUBJECTS);

if isempty(Files{1}), return, end
if length(Files{1} > 1), 
    Files = Files{1}; % Consider only first file
end

makeuswait('start')
% Read .dim header file
try 
    [nx, ny, nslices,datatype,dx,dy,dz,BO,OM] = ...
    textread(strrep(Files,'.ima','.dim'),'%d%d%d\n-type%s\n -dx%f -dy%f -dz%f\n-bo%s\n-om%s');
catch
    [nx, ny, nslices,nt, datatype,dx,dy,dz,dt,BO,OM] = ...
    textread(strrep(Files,'.ima','.dim'),'%d%d%d%d\n-type%s\n -dx%f -dy%f -dz%f -dt%f\n-bo%s\n-om%s');
end


% Store required information in MRI structure

MRI.header = 0; 
MRI.Voxsize = [dx dy dz]; clear dx dy dz
MRI.scantype = 'axial'; % Default (not specified in .dim header)
MRI.DataClass = strrep(datatype{1},'S','int');


switch(get(handles.LittleEndian,'checked'))
    case 'on'
        file = fopen(strrep(Files,'.dim','.ima'),'rb','l'); % 
    case 'off'
        file = fopen(strrep(Files,'.dim','.ima'),'rb','b'); %  big-endian byte ordering
end

MRI.Cube = int16(fread(file, prod([nx ny nslices]),MRI.DataClass)); % Read the whole volume / store in a huge vector   
fclose(file);
MRI.Cube = reshape(MRI.Cube,nx,ny,nslices); % reshape it into a 3-D array.

permVect = [3,1,2]; 

% reArrange voxel location to move to BrainStorm MRI coordinate system
MRI.Cube = permute(MRI.Cube,permVect);
MRI.Voxsize = MRI.Voxsize(permVect);

MRI.Cube = MRI.Cube(:,end:-1:1,end:-1:1);


% Show it
MRI = show_MRI(MRI);

setappdata(findobj(0,'tag','MRIViewer','type','figure'),'MRI',MRI)

makeuswait('stop')
% --------------------------------------------------------------------
function varargout = dicommri_Callback(h, eventdata, handles, varargin)


% --------------------------------------------------------------------
function varargout = rawmri_Callback(h, eventdata, handles, varargin)
% Import MRI raw scans

% GUI for input of MRI parameters
mri_parameters('Create',h,'','',gcbf)

paramFig = findobj(0,'Tag','MRIParameters','type','figure');
paramHandles = guihandles(paramFig);
% Change Button String
% set(paramHandles.Quit,'String','Done') % JCM: did it in the GUI
bst_message_window('wrap',...
    {,' ',' -> Fill-out the basic parameter values for the MRI images to be imported into BrainStorm, then click on the ''Done'' button.'});

% --------------------------------------------------------------------

function varargout = rawmri_done
% Go on once scan parameters are entered (called by mri_parameters.m)
MRIViewer = findobj(0,'Tag','MRIViewer','type','figure');
if isempty(MRIViewer), return, end

handles = guihandles(MRIViewer);

if ~isappdata(handles.MRIViewer,'tmpMRI') % User did killed the mri_parameters window before filling out all fields
    return
else
    tmpMRI = getappdata(handles.MRIViewer,'tmpMRI'); % Fetch image and file parameters
end

% Ask for file extension
[fileType] = inputdlg('Type the extension of the MRI file(s) to import (e.g. *.mri)','',1,{'*.*'});
if isempty(fileType{:})
    return
end

Users = get_user_directory;
[Files, Filenames] = file_selection_win(fileType,Users.SUBJECTS);

% Note: if several files were selected, scans are concatenated

islice = 0;
% Now read all files
makeuswait('start')
MRI.scantype = tmpMRI.scantype;

for file = Files
    try 
        file_mri = fopen(file{1},'r',tmpMRI.MachineFormat);
        if file_mri < 1
            bst_message_window('wrap',...
                {' ',sprintf('-> Invalid File Name %s',file{1}),' '})
            return;
        end
        
        fseek(file_mri,tmpMRI.HeaderSize,'bof'); % Skip header
        
        switch MRI.scantype % Axial/Sagittal/Coronal
            case 'axial'
                MRI.Cube = (zeros(tmpMRI.ImageSize(1),tmpMRI.ImageSize(2),tmpMRI.NumberOfSlices)); % Initialze
            case 'sagittal'
                MRI.Cube = (zeros(tmpMRI.ImageSize(1),tmpMRI.NumberOfSlices,tmpMRI.ImageSize(2))); % Initialze
            case 'coronal'
                MRI.Cube = (zeros(tmpMRI.ImageSize(2),tmpMRI.ImageSize(1),tmpMRI.NumberOfSlices)); % Initialze
        end
        
        hwait = waitbar(0,'Please wait while BrainStorm is importing the scans...');
        for i = 1:tmpMRI.NumberOfSlices
            islice = islice+1; % Global slice index across all files 
            if feof(file_mri)
                errordlg('Please Check Image size / Number of Slices : do not match file size')
                return
            end

            switch MRI.scantype % Axial/Sagittal/Coronal
                case 'axial'
                    MRI.Cube(:,:,islice) = (fread(file_mri,tmpMRI.ImageSize,tmpMRI.DataFormat));
                case 'sagittal'
                    MRI.Cube(:,islice,:) = (fread(file_mri,tmpMRI.ImageSize,tmpMRI.DataFormat));
                case 'coronal'
                    MRI.Cube(islice,:,:) = (fread(file_mri,tmpMRI.ImageSize,tmpMRI.DataFormat));        
            end
                                 
            if ~rem(islice,10)
                waitbar(i/tmpMRI.NumberOfSlices)
            end
            
        end   
        fclose(file_mri);
        close(hwait);
        
    catch
        % Error occured
        if ~isempty(Files{1})
            errordlg(sprintf('Wrong parameter specification for scan(s) in file %s', file{1}))
        end
        makeuswait('stop')
        return
    end
end

% Fill-out remaining parameters
MRI.DataClass = 'ushort';
MRI.Voxsize = [tmpMRI.FOV(1)/size(MRI.Cube,1), tmpMRI.FOV(2)/size(MRI.Cube,2), tmpMRI.SliceThickness];



% and store
setappdata(handles.MRIViewer,'MRI',MRI);

% House cleaning
rmappdata(handles.MRIViewer,'tmpMRI');

% Display
show_MRI(MRI);

makeuswait('stop')
% --------------------------------------------------------------------
function varargout = VolumeParameters_Callback(h, eventdata, handles, varargin)

% Call specific function
mri_parameters('Create',h,'','',gcbf)

% --------------------------------------------------------------------
function varargout = SaveAsBrainStorm_Callback(h, eventdata, handles, varargin)

% Save MRI in BrainStorm -mat format
MRI = getappdata(handles.MRIViewer,'MRI');

if ~isfield(MRI,'FileName')
    Current = get_user_directory;
    PATHSTR = Current.SUBJECTS;
    NAME = 'subjectimage.mat';
else
    [PATHSTR,NAME,EXT,VERSN] = fileparts(MRI.FileName);
end

cd(PATHSTR)
[file,path] = uiputfile(NAME,'Save Anatomical Scans in...');
if file == 0,     return,end
if isempty(findstr(file,'subjectimage'))
    file = strtok(file,'.');
    file = [file,'_subjectimage.mat'];
end

cd(path)
makeuswait('start')

% Clean-up MRI structure (useless or old fieldnames)
Fields2BDeleted  = {'FOV','Origin','sag','ax','cor','hFiducials','header','Comment','filename'};
SCSFields2BDeleted  = {'Origin','Comment'};

for k = 1:length(Fields2BDeleted)
    if isfield(MRI,Fields2BDeleted{k})
        MRI = rmfield(MRI,Fields2BDeleted{k});
    end
end

if isfield(MRI,'SCS')
    for k = 1:length(SCSFields2BDeleted)
        if isfield(MRI.SCS,SCSFields2BDeleted{k})
            MRI.SCS = rmfield(MRI.SCS,SCSFields2BDeleted{k});
        end
    end
end

if isfield(MRI,'SCS2Landmarks') % don't need to be stored in file
    tmpLandmarks = MRI.Landmarks;
    MRI.Landmarks.MRImmXYZ = MRI.Landmarks.MRImmXYZ(:,setdiff([1:end],MRI.SCS2Landmarks));
    MRI.Landmarks.Names = MRI.Landmarks.Names(setdiff([1:end],MRI.SCS2Landmarks));
    tmpHandles = MRI.Landmarks.Handles;
    MRI.Landmarks = rmfield(MRI.Landmarks,'Handles');
    tmpSCS2Landmarks = MRI.SCS2Landmarks; 
    MRI = rmfield(MRI,'SCS2Landmarks');
end

save_fieldnames(MRI,file)

if exist('tmpLandmarks','var') % don't need to be stored in file
    MRI.Landmarks = tmpLandmarks;
    MRI.SCS2Landmarks = tmpSCS2Landmarks;
    MRI.Landmarks.Handles = tmpHandles;
end

MRI.FileName = fullfile(path,file);
setappdata(handles.MRIViewer,'MRI',MRI);

bst_message_window('wrap',...
    sprintf('MRI saved in %s ',fullfile(path,file)))


makeuswait('stop')


% --------------------------------------------------------------------
function varargout = DefineSCS_Callback(h, eventdata, handles, varargin)

% Call external function and GUI
pcs_params('create','','','',get(handles.PointList,'string'));

%waitfor(findobj(0,'type','figure','tag','pcsparams')) % Get info about fiducials


% --------------------------------------------------------------------
function varargout = DefineSCSdone_Callback(h, eventdata, handles, varargin)

% Launched by pressing the "done" button in pcs_params, requires that
% pcs_params still be present with the information

% Get fiducials coordinates
% Since this routine may have been launched from outside, don't assume on
% input.
fig = findobj(0,'type','figure','tag','MRIViewer'); % find the handle to the bst_mriviewer
handles = guihandles(fig); % reload the handles

MRI = getappdata(fig,'MRI'); % get the current MRI structure

if isfield(MRI,'SCS')
    if isfield(MRI.SCS,'CubeFiducial') % Deprecated field; should not be saved in file
        MRI.SCS = rmfield(MRI.SCS,'CubeFiducial');
    end
end

% Get fiducial information
MRI.SCS2Landmarks = getappdata(findobj(0,'Tag','pcsparams'),'LookUpFiducials');
paramHandles = guihandles(findobj(0,'Tag','pcsparams'));
FiduSystem = get(paramHandles.checkboxnmag,'Value');
switch (FiduSystem)
    case 1
        FiduSystem = 'NEUROMAG';
    otherwise
        FiduSystem = 'CTF';
end
MRI.SCS.System = FiduSystem;
MRI.SCS.mmCubeFiducial = MRI.Landmarks.MRImmXYZ(:,MRI.SCS2Landmarks);
% Remove fiducial points from landmarks
MRI.SCS.FiducialName = MRI.Landmarks.Names(MRI.SCS2Landmarks);

[transf] = mri2scs(MRI);

MRI.SCS.R = transf.R;
MRI.SCS.T = transf.T;
MRI.SCS.Origin = transf.Origin;

MRI.SCS.mmCubeFiducial(:,end+1) = MRI.SCS.Origin;
MRI.SCS.FiducialName{end+1} = sprintf('%s SCS origin',MRI.SCS.System);

setappdata(fig,'MRI',MRI);

% Display 
plot_MRI(MRI,MRI.SCS.Origin,[],{sprintf('%s SCS origin',MRI.SCS.System)});
show_MRI(MRI,'update',mm2vox(MRI,MRI.SCS.Origin));

% Save
bst_mriviewer('SaveAsBrainStorm_Callback','','',handles);



% --------------------------------------------------------------------
function varargout = ShowFiducials_Callback(h, eventdata, handles, varargin)
% Show/Hide fiducials in the MR volume 

MRIfig = findobj(0,'Tag','MRIViewer','Type','figure');
MRI = getappdata(MRIfig,'MRI');
if ~isfield(MRI,'SCS'),
    errordlg('No fiducials were defined.','Unvalid operation')
    return
end

switch(get(findobj(0,'Tag','ShowFiducials'),'checked'))
    case 'off' % Show fiducials
        % Plot fiducials in volume
        %MRI.hFiducials = plot_MRI(MRI,MRI.SCS.mmCubeFiducial,[],MRI.SCS.FiducialName); % in meters
        %setappdata(MRIfig,'MRI',MRI)
        set(MRI.Landmarks.Handles(MRI.SCS2Landmarks),'visible','on')
        
        % Show fiducial list in lisbox
        set(findobj(MRIfig,'tag','PointList'),'String',MRI.Landmarks.Names,'Value',1);
        % Plot fiducial points
        bst_mriviewer('PointList_Callback','','', guihandles(MRIfig), '');
        
        
        % Move to MR slices that contain the Nasion point, i.e. point #1
        show_MRI(MRI,'update',mm2vox(MRI,MRI.SCS.mmCubeFiducial(:,1)'));
        
        
        set(findobj(0,'Tag','ShowFiducials'),'checked','on')
        
    case 'on' % Hide fiducials
        % Delete associated patch objects
        % 		if ishandle(MRI.hFiducials(1))
        % 			delete(MRI.hFiducials)
        % 		end
        set(MRI.Landmarks.Handles(MRI.SCS2Landmarks,:),'visible','off')		
        %MRI = rmfield(MRI,'hFiducials');
        %setappdata(MRIfig,'MRI',MRI)
        
        % Show all landmarks but those defined as SCS fiducials
        set(findobj(MRIfig,'tag','PointList'),'String',MRI.Landmarks.Names(setdiff([1:end],MRI.SCS2Landmarks)),'Value',1);
        
        set(findobj(0,'Tag','ShowFiducials'),'checked','off')
end


% --------------------------------------------------------------------

function [] = isvisible(MRI,hObjects)
%ISVISIBLE Check which object is visible in current slices
% isvisible(MRI,MRIplot)

MRIViewer = findobj(0,'tag','MRIViewer','type','figure');
handles = guihandles(MRIViewer);
currCoord = mm2vox(MRI,wherearewe(handles)); % In voxel units

hObjects = hObjects(:);
hObjects = hObjects(hObjects>0);
hObjects = hObjects(ishandle(hObjects));

ntol = 2; % Apply a n-voxel ntol tolerance on selected slice
for k = 1:length(hObjects)
    
    XYZObject = mm2vox(MRI,get(hObjects(k),'Userdata')');
    % Full 3D coordinates need to be stored in handles' Userdata
    if isempty(XYZObject) % Object is outside the slices
        %errordlg('Selected object is not within MRI volume. Please check its coordinates or redefine fiducial points.')
        return
    end
    
    switch(get(hObjects(k),'Parent')) % Now, depending on which slice the object belongs to
        case(handles.axa)
            if ismember(XYZObject(3),[currCoord(3)-ntol:currCoord(3)+ntol]) 
                set(hObjects(k),'Visible','on')
            else
                set(hObjects(k),'Visible','off')
            end
        case(handles.axc)
            if ismember(XYZObject(2),[currCoord(2)-ntol:currCoord(2)+ntol]) % Apply a one-voxel tolerance on selected slice
                set(hObjects(k),'Visible','on')
            else
                set(hObjects(k),'Visible','off')
            end
        case(handles.axs)
            if ismember(XYZObject(1),[currCoord(1)-ntol:currCoord(1)+ntol]) % Apply a one-voxel tolerance on selected slice
                set(hObjects(k),'Visible','on')
            else
                set(hObjects(k),'Visible','off')
            end
    end
    
end


% --------------------------------------------------------------------
function varargout = LoadSourceMap_Callback(h, eventdata, handles, varargin)
% Load a results file into the viewer
% Optional: filename is stored in varargin{1}

makeuswait('start')
VERBOSE = 1;

if nargin == 3
    Current = get_user_directory;
    cd(Current.STUDIES)
    % Select file
    [resFile, pathbs] = uigetfile('*result*.mat','Select a BrainStorm result file.');
    if resFile==0
        return
    end
    cd(pathbs)
    resFile = fullfile(pathbs,resFile);
else
    
    resFile = varargin{1}; % Result file
    % Load it
    Current = get_user_directory;
    resFile = fullfile(Current.STUDIES,resFile);
    cd(fileparts(resFile));
end

Results = load(resFile);

if VERBOSE,
    try
        bst_message_window('wrap',...
            {'Results file from'...
                sprintf('%s',fullfile(Current.STUDIES,Results.GUI.DataName))...
                'loaded into results visualization tool'...
                '------------------'});
    end
    
    %    bst_message_window({'Click on the list of sources to place a source on the MRI scan'...
    %          'To append list of sources choose LoadSources\Load Parametric '...
    %          'To delete a source from the list right click '...
    %          'while the cursor is over the source list to get a drop-down menu'...
    %          'and choose an appropriate action'...
    %          '------------------'});
end


% Is MRIViewer opened ?
mriFig = findobj(0,'Tag','MRIViewer','type','figure');
if isempty(mriFig)
    mriFig = bst_mriviewer;
end

% File location of subject's anatomy
try
    % if list is empty open MRI tool window 
    if(isfield(Results.Subject,'SubjectImage'))
        SubjectImageFile = fullfile(Current.SUBJECTS,Results.Subject.SubjectImage);
        %         elseif(isfield(Results.Subject,'Anatomy'))
        %             SubjectImageFile = fullfile(Current.SUBJECTS,Results.Subject.Anatomy);
    elseif isfield(Results.StudySubject,'Subject')
        Current = get_user_directory;
        load(fullfile(Current.SUBJECTS,Results.StudySubject.Subject),'Anatomy')
        SubjectImageFile = fullfile(Current.SUBJECTS,Anatomy);
    else
        warndlg('<Anatomy> or <SubjectImage> fields are not specified in the Results.Subject object');
        
        bst_message_window({'Warning!: <Anatomy> or <SubjectImage> fields are'...
                'not specified in the Results.Subject object'...
                '------------------'});
        return;
    end;
catch
    
    warndlg('Results structure does not contain all necessary fields');
    
end;

% Are there results already in memory of MRIViewer ?
if ~isappdata(mriFig,'mriResults') % No : load MRI  
    
    % Load subjectimage MRI
    bst_mriviewer('LoadBrainStorm_Callback','','','',SubjectImageFile);
    
    % Store for future usage
    mriResults.SubjectImageFile = SubjectImageFile;
    
else
    mriResults = getappdata(mriFig,'mriResults');
end

% Are current results from current subject in MRIViewer ?
if ~strcmp(SubjectImageFile,mriResults.SubjectImageFile) % No 
    mriResults.SubjectImageFile = SubjectImageFile;
    % Load subjectimage MRI
    bst_mriviewer('LoadBrainStorm_Callback','','','',mriResults.SubjectImageFile);
end

MRI = getappdata(mriFig,'MRI'); % Get MRI Information

% Do we append current sources to some former session ?
if(isfield(mriResults,'Source'))
    n_old = length(mriResults.Source);
    mriResults.nloads = mriResults.nloads + 1;
else
    n_old = 0;
    mriResults.nloads = 1;
end;
mriResults.ChannelFlag = Results.ChannelFlag;

nsource = length(Results.SourceLoc); % How many sources

% Start storing source information in mriResults structure
for n=1:nsource
    
    mriResults.Source(n_old+n).ChannelFlag = Results.ChannelFlag;
    mriResults.Source(n_old+n).Name= num2str(n);
    mriResults.Source(n_old+n).Loc = Results.SourceLoc{n};
    mriResults.Source(n_old+n).TimeSeries = Results.TimeSeries(:,n);
    
    
    if isfield(Results,'DataFlag') % Original data: either MEG, EEG or Fusion 
        mriResults.Source(n_old+n).DataFlag = Results.DataFlag;
    else
        mriResults.Source(n_old+n).DataFlag = 'MEG'; % Assume default
    end
    
    Data = load(Results.GUI.DataName,'Time'); % Load original data time vector
    if(~isempty(Data))
        mriResults.Source(n_old+n).Time = Data.Time(Results.Time);
    else
        mriResults.Source(n_old+n).Time =  1:length(mriResults.Source(n_old+n).TimeSeries);
    end;
    
    %    % JCM 15-Nov-2003 for some reason, AOS removed the mean. Reversed that
    %    %  and multiplied the time series  by 1e9 to yield nA-m
    %    if 0,
    %       mriResults.Source(n_old+n).TCPlotHandle = plot(mriResults.Source(n_old+n).Time*1000,...
    %          mriResults.Source(n_old+n).TimeSeries-...
    %          mean(mriResults.Source(n_old+n).TimeSeries),...
    %          'LineWidth',1,'Color','b','Visible','off',...
    %          'Parent',handles.sourcetimeseries);
    %    else
    %       % time in milliseconds, amplitude in nA-m
    %       mriResults.Source(n_old+n).TCPlotHandle = plot(mriResults.Source(n_old+n).Time * 1000,...
    %          mriResults.Source(n_old+n).TimeSeries * 1e9,...
    %          'LineWidth',1,'Color','b','Visible','off',...
    %          'Parent',handles.sourcetimeseries);
    %    end
    
    
    mriResults.Source(n_old+n).SourceCorrelation = Results.SourceCorrelation(n);      
    mriResults.Source(n_old+n).IndepTopo = Results.IndepTopo(:,n);
    
    if(isfield(Results.StudySubject,'Channel')) % Channel is a field from StudySubject
        mriResults.Source(n_old+n).Channel=Results.StudySubject.Channel;      
    elseif(isfield(Results.Study,'Channel'))
        mriResults.Source(n_old+n).Channel=Results.Study.Channel;
    else
        mriResults.Source(n_old+n).Channel=[];
        bst_message_window({sprintf('Default channel file %s not found',ChannelTest),...
                'Source topography will not be available'});
    end;
    
    
    mriResults.Source(n_old+n).SourceOrder = Results.SourceOrder(n);
    mriResults.Source(n_old+n).Comment = Results.Comment;
    % mriResults.Source(n_old+n).LocMRI = transform('rap2mri',Results.SourceLoc{n},PCS);
    
    % Move to mri coordinates
    [mriResults.Source(n_old+n).LocMRI] = scs2mri(MRI,Results.SourceLoc{n}*1000);
    
    if(mriResults.Source(n_old+n).SourceOrder==-1)
        tmp = scs2mri(MRI,1000*Results.SourceOrientation{n});                                         
        mriResults.Source(n_old+n).OrientMRI = tmp/norm(tmp);
    elseif(mriResults.Source(n_old+n).SourceOrder==2)
        tmp = scs2mri(MRI,1000*Results.SourceOrientation{n}(1:3)); 
        mriResults.Source(n_old+n).OrientMRI(:,1) = tmp/norm(tmp);
        tmp = scs2mri(MRI,1000*Results.SourceOrientation{n}(4:6));
        mriResults.Source(n_old+n).OrientMRI(:,2) = tmp/norm(tmp);
    else
        mriResults.Source(n_old+n).OrientMRI = [0 0 0]';
    end;
    
    
    % See source indep topos
    if n == 1;
        topoFig = figure;
        %set(topoFig,'visible','off'); % hide while we build
        % give it a unique name for now
        set(topoFig,'name',sprintf('Topography Fig %.0f',topoFig));
        setappdata(topoFig,'TileType','T');
        bst_color_scheme(topoFig)
        set(topoFig,'menubar','figure')
        bst_layout('align',topoFig,2,2,2);
        %set(topoFig,'visible','on'); 
        
        nrows = ceil(sqrt(nsource));
        ncols = ceil(sqrt(nsource));
        
        
    end
    figure(topoFig), hold on
    OPTIONS.Axes = subplot(nrows,ncols,n);
    OPTIONS.DisplayTypes = [0 0 5]; % 2D disc
    OPTIONS.Colorbar = 0;
    Channel = load(mriResults.Source(n_old+n).Channel);
    Channel = Channel.Channel;
    OPTIONS.SelectedChannels{1} = good_channel(Channel,...
        mriResults.Source(n_old+n).ChannelFlag,mriResults.Source(n_old+n).DataFlag); % CHEAT Hardwired 'MEG' results
    OPTIONS.Channel =  Channel(OPTIONS.SelectedChannels{1}); % CHEAT Hardwired 'MEG' results
    OPTIONS.SelectedChannels{1} = 1:length(OPTIONS.Channel);
    OPTIONS.Modality = 1; % CHEAT Hardwired 'MEG' results
    OPTIONS.ModalityLabel = {sprintf('Source %d',n)};
    OPTIONS.SensorMarkers = 0;
    bst_wavedata_display({mriResults.Source(n_old+n).IndepTopo},OPTIONS)
    
end;

timeseriesFig = figure; 
setappdata(timeseriesFig,'TileType','T');
%set(timeseriesFig,'visible','off'); % hide while we build
% give it a unique name for now
set(timeseriesFig,'name',sprintf('Timeseries Fig %.0f',timeseriesFig));
bst_color_scheme(timeseriesFig)
set(timeseriesFig,'menubar','figure');
bst_layout('align',timeseriesFig,2,2,4);
%set(timeseriesFig,'visible','on');
OPTIONS.DisplayTypes = [1 0 0]; % Overlay
OPTIONS.DataType = 'source';
OPTIONS.Modality = 1; % CHEAT Hardwired 'MEG' results
OPTIONS.ModalityLabel = {'Source(s)'};
OPTIONS.Time = mriResults.Source(n_old+1).Time;
OPTIONS.WaveLabel = {mriResults.Source(n_old+1:end).Name};
bst_wavedata_display({mriResults.Source(n_old+1:end).TimeSeries},OPTIONS)


setappdata(mriFig,'mriResults',mriResults)

% See sources (with orientations)
handles = guihandles(mriFig);
nptlist = length(get(handles.PointList,'String'));

ObjectHandles = plot_MRI(MRI,[mriResults.Source([n_old+1:end]).LocMRI],...
    [mriResults.Source([n_old+1:end]).OrientMRI],...
    {mriResults.Source(n_old+1:end).Name});

% Update MRI.Landmarks structure
if ~isfield(MRI,'Landmarks')
    MRI.Landmarks.MRImmXYZ = [];
    MRI.Landmarks.Names = {};
    MRI.Landmarks.Handles = [];
end

MRI.Landmarks.MRImmXYZ = [MRI.Landmarks.MRImmXYZ,[mriResults.Source([n_old+1:end]).LocMRI]];
MRI.Landmarks.Handles = [MRI.Landmarks.Handles;ObjectHandles];

if isempty(MRI.Landmarks.Names)
    MRI.Landmarks.Names = cellstr(num2str([1:size(ObjectHandles,1)]'));
else
    MRI.Landmarks.Names(end+1:end+size(ObjectHandles,1)) = cellstr(num2str([1:size(ObjectHandles,1)]'));
end

setappdata(mriFig,'MRI',MRI)

set(handles.PointList,'String',MRI.Landmarks.Names,'Value',nptlist+1)
bst_mriviewer('AddCrosshair_Callback','','','',1);
bst_mriviewer('PointList_Callback','','','')

makeuswait('stop')

% --------------------------------------------------------------------
function varargout = SaveScouts_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = LoadScouts_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = cortic_remap_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = ImportSPM_Callback(h, eventdata, handles, varargin)
% Import SPM functional maps and overlay with anatomical
% Go on once scan parameters are entered (called by mri_parameters.m)

makeuswait('start')

MRIViewer = findobj(0,'Tag','MRIViewer','type','figure');
if isempty(MRIViewer), return, end

handles = guihandles(MRIViewer);

if ~isappdata(handles.MRIViewer,'tmpMRI') % User did killed the mri_parameters window before filling out all fields
    tmpMRI.MachineFormat = 'ieee-be'; % Default machine format is big-endian
else
    tmpMRI = getappdata(handles.MRIViewer,'tmpMRI'); % Fetch image and file parameters
end


%% --- Load ANALYZE files
Users = get_user_directory;
cd(Users.SUBJECTS)   

[anatfile, anatpath] = uigetfile('*.img','Please select the aMRI file (.img)' );    
if (anatfile == 0), return, end
cd(anatpath);

% Is SPM installed on this machine ?
tmp = which('spm_vol');
if isempty(tmp)
    errordlg('SPM toolbox was not found in Matlab search paths. Action is aborted.','SPM toolbox not found');
    return
end

eval('V = spm_vol(anatfile);','ok=0;');
Vorig = V.mat;
clear V

[blobfile, blobpath] = uigetfile('*.img','Please select the fMRI map');    
if (blobfile == 0), return, end
cd(blobpath);

eval('V = spm_vol(blobfile);','ok=0;');


%---------------------------------------
% Read img volumes
% Anatomy

% Call D. Weber's tool
[ avw, machine ] = avw_img_read(strtok(anatfile,'.'),'', tmpMRI.MachineFormat );

% Convert to BrainStorm format
MRI = avw2brainstorm(avw, machine);

clear avw machine

% fMRI
bst_message_window('Reading functional files...');
file = fopen(V.fname,'rb',tmpMRI.MachineFormat );
%fseek(file,0,'bof');
fMRI.Cube = zeros(V.dim(1),V.dim(2),V.dim(3));
    
for i= 1:V.dim(3)%:-1:1 % Lower Slice first
    if feof(file)
        msgbox('Please check that data is of ''uint8'' type', 'BrainStorm Message')
        errordlg('Please Check Image size / Number of Slices : do not match file size')
        return
    end
    fMRI.Cube(:,:,i) = (fread(file,[V.dim(1),V.dim(2)],'uint8'));
    %waitbar((mri_scan.blob.nslice-i)/mri_scan.blob.nslice);
end
fclose(file);
fMRI.Cube(isnan(fMRI.Cube)) = 0;
bst_message_window('Reading functional files... -> Done');


%---------------------------------------

% 
% disp('Reading functional files...');
% file = fopen(V.fname,'rb',MACHINEFORMAT);
% fseek(file,header,'bof');
% for i= 1:mri_scan.blob.nslice%:-1:1 % Lower Slice first
%     if feof(file)
%         msgbox('Please check that data is of ''uint8'' type', 'BrainStorm Message')
%         errordlg('Please Check Image size / Number of Slices : do not match file size')
%         return
%     end
%     blob(:,:,i) = (fread(file,[siz(1),siz(2)],mri_scan.code));
%     %waitbar((mri_scan.blob.nslice-i)/mri_scan.blob.nslice);
% end
% % fclose(file);
% % disp('Reading functional files... -> Done');
% 
% mri_scan.blob.scan = blob;
% mri_scan.blob.V = V;  
% 
% eval('mri_scan.blob')
% 
% %% --- reading of the functional file completed
% 
% 
% %% -- Interpolation to the resolution of the anatomical MR scan
% 
% pix_siz = mri_scan.pix_siz ;  % Voxel size of the anatomical MR scan 

%% -- Reslicing
bst_message_window('wrap',...
    'Interpolating functional scans into anatomical volume...')

% [XI,YI] = meshgrid(linspace(1,size(fMRI.Cube,1),size(MRI.Cube,1)),...
%     linspace(1,size(fMRI.Cube,2),size(MRI.Cube,2)));
% iblob = zeros(size(MRI.Cube,2),size(MRI.Cube,1),size(fMRI.Cube,3));
% for i=1:size(fMRI.Cube,3)
%     iblob(:,:,i) = interp2(fMRI.Cube(:,:,i),XI,YI,'*linear');
% end

% Interpolate along the z-axis
[qo,ro] = qr(Vorig(1:3,1:3));
[q,r] = qr(V.mat(1:3,1:3));

scale = abs(r(9)/ro(9)); % Scaling along the z-axis
n_slice = floor(scale * size(fMRI.Cube,3)); 

[XI,YI,ZI] = meshgrid(linspace(1,size(fMRI.Cube,2),size(MRI.Cube,2)),...
    linspace(1,size(fMRI.Cube,1),size(MRI.Cube,1)),linspace(1,size(fMRI.Cube,3),n_slice));
iiblob = interp3(fMRI.Cube,XI,YI,ZI,'*linear');
iiblob(isnan(iiblob)) = 0;

scale_func = size(fMRI.Cube); clear fMRI

% Position of the functional blobs in mm
[iorig] = find(iiblob);
t = iiblob(iorig); % Effects 
[i,j,k] = ind2sub(size(iiblob),iorig);
xyz = [i';j';k'];
clear i j k 

% Rescale to original voxel units
scale_anat= (size(iiblob));
scale = diag( scale_func./scale_anat);
xyz = scale*xyz ;
xyz =  [V.mat] * [xyz;ones(1,size(xyz,2))] ; % affine transformation mapping from voxel coordinates to real world coordinates (fMRI space)
xyz = xyz(1:3,:);
xyz =  inv(Vorig) * [xyz;ones(1,size(xyz,2))] ; % affine transformation mapping from real world coordinates to voxel coordinates (MRI space)
xyz = round(xyz(1:3,:));
%xyz = round(ROTmm2mri * xyz + repmat(TRANSmm2mri,1, size(xyz,2)));


mri_scan.CubeOrigin = Vorig(1:3,4);%round([ORI./scale]'); % Indices of the slices containing the origin

mri_scan.blob.Cube = xyz;


%iiblob = zeros(size(mri_scan.mri,1),size(mri_scan.mri,2),size(mri_scan.mri,3));
sizz = size(MRI.Cube);
tt = sub2ind([sizz],mri_scan.blob.Cube(1,:),mri_scan.blob.Cube(2,:),mri_scan.blob.Cube(3,:));

% % For Kevin ----- flip U/D
% tmp = zeros(size(iiblob(:,:,1)));
% for k = 1:size(iiblob,3)
%     tmp = flipud(iiblob(:,:,k));
%     iiblob(:,:,k) = tmp; 
% end
% clear tmp
[iorig] = find(iiblob);
t = iiblob(iorig); % Effects 
clear iiblob 
ngray = floor(max(MRI.Cube(:))/2);
MRI.Cube(tt(t>0*mean(t))) =  ngray+(t(t>0*mean(t))+eps)*ngray/max(t);

% Show result
setappdata(MRIViewer,'MRI',MRI);
show_MRI(MRI);



% ax = findobj(MAINWINDOW,'Tag','colorbar');
% axes(ax)
% colormapp = split(65:end,:);
% n = size(colormapp,1);
% image([0, 1], [0 10],(n+1:2*n)');
% set(ax,'Ydir','normal','Tag','colorbar')
% set(ax,'xtick',[],'ycolor',frontcolor,'xcolor',frontcolor)


% Colormap
split = [bone(ngray);hot(ngray+10)];
colormap(split)



% --------------------------------------------------------------------
function varargout = Talairach_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = TalairachTransform_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = DefineScouts_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = DefineCorticalScouts_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = Quit_Callback(h, eventdata, handles, varargin)

delete(findobj(0,'Tag','MRIViewer','Type','figure'))
delete(findobj(0,'Tag','pcsparams','Type','figure'))

% --------------------------------------------------------------------
function [ObjectHandles] = plot_MRI(MRI,pointsCoord,varargin);
%plot_MRI - Plot points in MR volume using BST_MRIVIEWER
% hFiducials = plot_MRI(MRI,pointsCoord);
% hFiducials = plot_MRI(MRI,pointsCoord,pointOrient);
% hFiducials = plot_MRI(MRI,pointsCoord,pointOrient, pointsLabel);
% MRI is a valid BrainStorm MRI structure (see subjectimage file format)
% pointsCoord is a 3x1 vector or a 3xN matrix which contains the coordinates of the points (in mm) to be plotted 
% pointsOrient is a 3x1 vector or a 3xN matrix which contains the orientation of the sources 
% pointsLabel is a cell array of strings containing a label for each landmark
%
% ObjectHandles a Nx3 (or Nx6 when orientations are specified) matrix of patch (and line) object handles one for each point.
% ObjectHandles(k,i) is the handle of the kth point in the coronal (i=1), sagittal (i=2), axial (i=3) plane 
% (i=4:6 are the same for line objects associated to orientations in all 3 planes, if specified in pointsOrient)

if size(pointsCoord,1) ~= 3
    pointsCoord = pointsCoord';
end

if nargin > 2
    pointsOrient = varargin{1};
    if size(pointsOrient,1) ~= 3
        pointsOrient = pointsOrient';
    end
    % 	pointsLabel = cellstr(num2str(1:size(pointsCoord,2))); % Default labelling
    % 	if nargin == 4
    % 		pointsLabel = varargin{2};
    % 	end
else
    pointsOrient = [];
end

ObjectHandles = zeros(size(pointsCoord,2),6);

fig = findobj(0,'tag','MRIViewer','type','figure');
handles = guihandles(fig);

Colors = get(0,'DefaultAxesColorOrder'); % Default line colors

for k = 1:size(pointsCoord,2) % For each point
    % Axial plane
    axes(handles.axa), hold on
    % Handle to patch object in axial plane
    ObjectHandles(k,3) = scatter(pointsCoord(1,k),pointsCoord(2,k),'o','filled');
    set(ObjectHandles(k,3),'Parent',handles.axa,'visible','off');
    if ~isempty(pointsOrient) % Draw line
        pointsOrient = 3*pointsOrient; % plot scaling factor
        ObjectHandles(k,6) = line([pointsCoord(1,k),pointsCoord(1,k)+pointsOrient(1,k)],[pointsCoord(2,k),pointsCoord(2,k)+pointsOrient(2,k)]);
    end
    
    % same story for next 2 planes
    axes(handles.axc), hold on
    ObjectHandles(k,1) = scatter(pointsCoord(1,k),pointsCoord(3,k),'o','filled');
    set(ObjectHandles(k,1),'Parent',handles.axc,'visible','off');
    if ~isempty(pointsOrient) % Draw line
        ObjectHandles(k,4) = line([pointsCoord(1,k),pointsCoord(1,k)+pointsOrient(1,k)],[pointsCoord(3,k),pointsCoord(3,k)+pointsOrient(3,k)]);
    end
    
    axes(handles.axs), hold on
    ObjectHandles(k,2) = scatter(pointsCoord(2,k),pointsCoord(3,k),'o','filled');
    set(ObjectHandles(k,2),'Parent',handles.axs,'visible','off');
    if ~isempty(pointsOrient) % Draw line
        ObjectHandles(k,5) = line([pointsCoord(2,k),pointsCoord(2,k)+pointsOrient(2,k)],[pointsCoord(3,k),pointsCoord(3,k)+pointsOrient(3,k)]);
    end
    
    set([ObjectHandles(k,1:3)],...
        'Userdata',pointsCoord(:,k),'markerfacecolor',Colors(modulo(k,size(Colors,1)),:),...
        'MarkerEdgeColor',Colors(modulo(k,size(Colors,1)),:)) % Store full 3D coordinates in handles' Userdata for isvisible check
    
    if ~isempty(pointsOrient) % Draw line
        set([ObjectHandles(k,4:6)],...
            'Userdata',pointsCoord(:,k),'color',Colors(modulo(k,size(Colors,1)),:),'linewidth',2) % Store full 3D coordinates in handles' Userdata for isvisible check
    end    
end




% --------------------------------------------------------------------

function [MRI] = show_MRI(MRI, varargin);
% show_MRI - Show a MRI volume in BST_MRIVIEWER
% function [] = show_MRI(MRI, 'update', [iX,iY,iZ]);
% iX (res. iY, iZ) is the index of the coronal (res. sagittal, axial) slice to be displayed

fig = findobj(0,'type','figure','tag','MRIViewer');

if nargin > 1
    
    handles = guihandles(fig);
    action = varargin{1};
    switch(action)
        case 'update' % Update display with new slices
            newIndices = varargin{2};
            set(findobj(handles.axa,'type','image'),'CData',MRI.Cube(:,:,newIndices(3))','CdataMapping','scaled')
            set(findobj(handles.axc,'type','image'),'CData',squeeze(MRI.Cube(:,newIndices(2),:))','CdataMapping','scaled')
            set(findobj(handles.axs,'type','image'),'CData',squeeze(MRI.Cube(newIndices(1),:,:))','CdataMapping','scaled')
            
            mriCoord = vox2mm(MRI,newIndices);
            
    end
    
    
else  % Create new display
    
    
    fig = bst_mriviewer;
    handles = guihandles(fig);
    
    MRI.FOV = size(MRI.Cube).*MRI.Voxsize;
    
    ngray = 256; % Default grayscale 
    
    % Delete previous graphic objects
    delete(findobj([get(handles.axa,'children'),get(handles.axc,'children'),get(handles.axs,'children')],...
        'type','patch'))
    delete(findobj([get(handles.axa,'children'),get(handles.axc,'children'),get(handles.axs,'children')],...
        'type','line'))
    delete(findobj([get(handles.axa,'children'),get(handles.axc,'children'),get(handles.axs,'children')],...
        'type','text'))
    delete(findobj([get(handles.axa,'children'),get(handles.axc,'children'),get(handles.axs,'children')],...
        'type','image'))
    
    %     if isappdata(fig,'MRIPoint')
    %         rmappdata(fig,'MRIPoint');
    %     end
    set(handles.PointList,'String','')
    
    %     delete(findobj(get(handles.axc,'children'),'type','patch'))
    %     delete(findobj(get(handles.axs,'children'),'type','patch'))
    %     delete(findobj(get(handles.axa,'children'),'type','line'))
    %     delete(findobj(get(handles.axc,'children'),'type','line'))
    %     delete(findobj(get(handles.axs,'children'),'type','line'))
    %     delete(findobj(get(handles.axa,'children'),'type','image'))
    %     delete(findobj(get(handles.axc,'children'),'type','image'))
    %     delete(findobj(get(handles.axs,'children'),'type','image'))
    
    if 1 % deprecated code
        % Scale gray levels to 255 to save the MRI in uint8
        % Memory requirememts necessitate a slice-by-slice loop on most machines
        M = double(max(MRI.Cube(:))); % May be rather time consuming if limited RAM is availlable
        if M > (ngray-1)    
            hw = waitbar(0,'Scaling the MRI to 256 gray levels...');
            for k = 1:size(MRI.Cube,3)
                if ~rem(k,50)
                    waitbar(k/size(MRI.Cube,3),hw);
                end
                
                tmp = double(MRI.Cube(:,:,k)) ;
                MRI.Cube(:,:,k) = int16(tmp * ((ngray - 1 )/M)) ;
            end
            
            MRI.Cube = uint8(MRI.Cube);
            delete(hw)          
            drawnow
            setappdata(fig,'MRI',MRI);
        end
        
        if M < ngray - 1
            tmp = double(MRI.Cube) ;
            MRI.Cube = int16(tmp * ((ngray - 1 )/M)) ;
        end
        
        clear tmp
    end
    
    % Default slice index values
    iX  = round(size(MRI.Cube,1)/2);
    iY  = round(size(MRI.Cube,2)/2);
    iZ  = round(size(MRI.Cube,3)/2);
    
    % Display 
    
    %% -- Axial View
    axes(handles.axa)
    MRI.ax = image([0, MRI.FOV(1)],[0,MRI.FOV(2)],MRI.Cube(:,:,iZ)','Parent',handles.axa);% DO NOT REMOVE transpose IF YOU KEEP ON USING the image COMMAND
    V = axis;
    set(handles.axa,'Tag','axa','Xtick',[],'Ytick',[])
    text(V(1),V(4),'A','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'P','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'L','horizontalalignment','left',...
        'verticalalignment','cap','color','w');
    text(V(2),V(3),'R','horizontalalignment','right',...
        'verticalalignment','cap','color','w');
    
    caxis([0 ngray-1])
    %axis([0, max(MRI.FOV) 0, max(MRI.FOV)])
    %% -- Axial view Done
    
    %% -- Sagittal view
    axes(handles.axs) 
    MRI.sag = image([0 MRI.FOV(2)],[0 MRI.FOV(3)],squeeze(MRI.Cube(iX,:,:))','Parent',handles.axs);% DO NOT REMOVE transpose IF YOU KEEP ON USING the image COMMAND
    V = axis;
    set(handles.axs,'Tag','axs','Xtick',[],'Ytick',[])
    text(V(1),V(4),'T','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'B','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'P','horizontalalignment','left',...
        'verticalalignment','cap','color','w');
    text(V(2),V(3),'A','horizontalalignment','right',...
        'verticalalignment','cap','color','w');
    
    caxis([0 ngray-1])
    %axis([0, max(MRI.FOV) 0, max(MRI.FOV)])
    %% -- Sagittal view Done
    
    %% -- Coronal view
    axes(handles.axc)
    MRI.cor = image([0 MRI.FOV(1)],[0 MRI.FOV(3)],squeeze(MRI.Cube(:,iY,:))','Parent',handles.axc);% DO NOT REMOVE transpose IF YOU KEEP ON USING the image COMMAND
    V = axis;
    set(handles.axc,'Tag','axc','Xtick',[],'Ytick',[]);
    text(V(1),V(4),'T','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'B','horizontalalignment','right',...
        'verticalalignment','baseline','color','w');
    text(V(1),V(3),'L','horizontalalignment','left',...
        'verticalalignment','cap','color','w');
    text(V(2),V(3),'R','horizontalalignment','right',...
        'verticalalignment','cap','color','w');
    caxis([0 ngray-1])
    %axis([0, max(MRI.FOV) 0, max(MRI.FOV)])
    %% -- Coronal view Done
    
    % Reset axes look 
    reset_axis([handles.axc, handles.axs, handles.axa]);
    set([MRI.ax,MRI.cor,MRI.sag],'uicontextmenu',get(handles.axc,'uicontextmenu'))
    
    % Update navigation handles
    mriCoord = vox2mm(MRI,[iX,iY,iZ]);
    set(handles.slider_x,'Visible','on','Min', 0,'Max',MRI.FOV(1))
    set(handles.slider_y,'Visible','on','Min', 0,'Max',MRI.FOV(2))
    set(handles.slider_z,'Visible','on','Min', 0,'Max',MRI.FOV(3))
    
    % Apply colormap
    colormap(bone(ngray))
end

set(handles.X,'fontunits','points','fontsize',7,'string',sprintf('X = %3.1f',mriCoord(1)))
set(handles.Y,'fontunits','points','fontsize',7,'string',sprintf('Y = %3.1f',mriCoord(2)))
set(handles.Z,'fontunits','points','fontsize',7,'string',sprintf('Z = %3.1f',mriCoord(3)))
if isfield(MRI,'SCS')
    [transf,scsCoord] = mri2scs(MRI,mriCoord');
    set(handles.scsX,'fontunits','points','fontsize',7,'string',sprintf('X = %3.1f',scsCoord(1)))
    set(handles.scsY,'fontunits','points','fontsize',7,'string',sprintf('Y = %3.1f',scsCoord(2)))
    set(handles.scsZ,'fontunits','points','fontsize',7,'string',sprintf('Z = %3.1f',scsCoord(3)))
else
    set([handles.scsX,handles.scsY,handles.scsZ],'fontunits','points','fontsize',7,'string',' - ')
end

set(handles.slider_x,'Value',mriCoord(1))
set(handles.slider_y,'Value',mriCoord(2))
set(handles.slider_z,'Value',mriCoord(3))

% Update crosshair location if necessary
if strcmp(get(findobj(fig,'tag','AddCrosshair'),'checked'),'on')
    bst_mriviewer('AddCrosshair_Callback',findobj(fig,'tag','AddCrosshair'),'',handles,1); 
end

% Check which object is visible in current slices
if isfield(MRI,'Landmarks')
    if isfield(MRI.Landmarks,'Handles')
        isvisible(MRI,MRI.Landmarks.Handles)
    end
end



% --------------------------------------------------------------------

function   [] = reset_axis(H)
% Properly reset the appearance of the MRIViewer axes

line_color = .8*[1 1 1]; % Axis Borders

set(H,'visible','on','Xcolor',line_color ,'Ycolor',line_color,'Color','k',...
    'box','on','layer','top','xgrid','off','ygrid','off','Ydir','normal','gridlinestyle','-',...
    'XtickLabel','','YTickLabel','', 'Xtick',[], 'Ytick',[]);   

axis(H, 'equal')

set(0,'CurrentFigure',findobj(0,'Tag','MRIViewer','Type','figure'))
axisMenu = uicontextmenu;
uimenu(axisMenu,'Label','Flip Up/Down','Callback',...
    'bst_mriviewer(''FlipUpDown'')');
uimenu(axisMenu,'Label','Flip Left/Right','Callback',...
    'bst_mriviewer(''FlipLeftRight'')');
uimenu(axisMenu,'Separator','on','Label','Transpose','Callback',...
    'bst_mriviewer(''Transpose'')');


set(H,'uicontextmenu',axisMenu)




% --------------------------------------------------------------------

function MRI = FlipUpDown
% Flip MRI images up/down, based on on currentaxes

% Get MRI volume
MRI = getappdata(gcbf,'MRI');
handles = guihandles(gcbf);

switch(gca) % Depending on current axes
    case (handles.axc)
        MRI.Cube = MRI.Cube(:,:,[end:-1:1]);
    case (handles.axa)
        MRI.Cube = MRI.Cube(:,[end:-1:1],:);
    case (handles.axs)
        MRI.Cube = MRI.Cube(:,:,[end:-1:1]);
end

XYZ = wherearewe(handles);
setappdata(gcbf,'MRI',MRI)
MRI = show_MRI(MRI,'update',mm2vox(MRI,XYZ));

if isfield(MRI,'Landmarks') % permute also landmarks locations
    tmp = size(MRI.Cube) .* MRI.Voxsize;
    switch(gca) % Depending on current axes
        case {handles.axc,handles.axs}
            MRI.Landmarks.MRImmXYZ(3,:) = tmp(1)-MRI.Landmarks.MRImmXYZ(3,:);
        case (handles.axa)
            MRI.Landmarks.MRImmXYZ(2,:) = tmp(2)-MRI.Landmarks.MRImmXYZ(2,:);
    end
    
    delete(MRI.Landmarks.Handles(MRI.Landmarks.Handles~=0))
    MRI.Landmarks = rmfield(MRI.Landmarks,'Handles');
    MRI.Landmarks.Handles = plot_MRI(MRI,MRI.Landmarks.MRImmXYZ,[],MRI.Landmarks.Names);
    isvisible(MRI,MRI.Landmarks.Handles)
end


setappdata(gcbf,'MRI',MRI)

% --------------------------------------------------------------------

function MRI = FlipLeftRight
% Flip MRI images Left/Right, based on on currentaxes

% Get MRI volume
MRI = getappdata(gcbf,'MRI');
handles = guihandles(gcbf);

switch(gca) % Depending on current axes
    case (handles.axc)
        MRI.Cube = MRI.Cube([end:-1:1],:,:);
    case (handles.axa)
        MRI.Cube = MRI.Cube([end:-1:1],:,:);
    case (handles.axs)
        MRI.Cube = MRI.Cube(:,[end:-1:1],:);
end

XYZ = wherearewe(handles);
setappdata(gcbf,'MRI',MRI)
MRI = show_MRI(MRI,'update',mm2vox(MRI,XYZ));

if isfield(MRI,'Landmarks') % permute also landmarks locations
    tmp = size(MRI.Cube) .* MRI.Voxsize;
    switch(gca) % Depending on current axes
        case {handles.axc,handles.axa}
            MRI.Landmarks.MRImmXYZ(1,:) = tmp(1)-MRI.Landmarks.MRImmXYZ(1,:);
        case (handles.axs)
            MRI.Landmarks.MRImmXYZ(2,:) = tmp(2)-MRI.Landmarks.MRImmXYZ(2,:);
    end
    
    delete(MRI.Landmarks.Handles(MRI.Landmarks.Handles~=0))
    MRI.Landmarks = rmfield(MRI.Landmarks,'Handles');
    MRI.Landmarks.Handles = plot_MRI(MRI,MRI.Landmarks.MRImmXYZ,[],MRI.Landmarks.Names);
    isvisible(MRI,MRI.Landmarks.Handles)
end


setappdata(gcbf,'MRI',MRI)

%bst_mriviewer('slider_Callback','','',guihandles(gcbf))

% --------------------------------------------------------------------

function MRI = Transpose
% Transpose MRI images, based on on currentaxes

makeuswait('start')
% Get MRI volume
MRI = getappdata(gcbf,'MRI');
handles = guihandles(gcbf);

switch(gca) % Depending on current axes
    case (handles.axc)
        perm = [3 2 1]; % permutation of indices
    case (handles.axa)
        perm = [2 1 3]; % permutation of indices
    case (handles.axs)
        perm = [1 3 2]; % permutation of indices
end

set(gcbf,'pointer','watch')
MRI.Cube = permute(MRI.Cube, perm);

% update voxel dimensions & FOV
MRI.Voxsize = MRI.Voxsize(perm);
MRI.FOV = size(MRI.Cube) .* MRI.Voxsize ;

XYZ = wherearewe(handles);
setappdata(gcbf,'MRI',MRI)
MRI = show_MRI(MRI,'update',mm2vox(MRI,XYZ));

if isfield(MRI,'Landmarks') % permute also landmarks locations
    MRI.Landmarks.MRImmXYZ = MRI.Landmarks.MRImmXYZ(perm,:);
    delete(MRI.Landmarks.Handles(MRI.Landmarks.Handles~=0))
    MRI.Landmarks = rmfield(MRI.Landmarks,'Handles');
    MRI.Landmarks.Handles = plot_MRI(MRI,MRI.Landmarks.MRImmXYZ,[],MRI.Landmarks.Names);
    isvisible(MRI,MRI.Landmarks.Handles)
end

setappdata(gcbf,'MRI',MRI)

% if isfield(MRI,'SCS')
%     % Finally, show fiducials if requested
%     bst_mriviewer('ShowFiducials_Callback','','','')
% end

makeuswait('stop')

% --------------------------------------------------------------------
function mriCoord = vox2mm(MRI,sliceIndices);
% Coordinate transform from slice indices to mm (MR coordinates)

if size(sliceIndices,1) == 3
    sliceIndices = sliceIndices';
end

mriCoord = ((size(MRI.Cube).*MRI.Voxsize)./(size(MRI.Cube)-1)) .* (sliceIndices -1);
if ~isempty(find(mriCoord > (size(MRI.Cube).*MRI.Voxsize))) % Point effectively outside image
    mriCoord = [];    
end

% --------------------------------------------------------------------
function sliceIndices = mm2vox(MRI,mriCoord);
% Coordinate transform from mm (MR coordinates) indices to indices

if size(mriCoord,1) == 3
    mriCoord = mriCoord';
end

sliceIndices = round((((size(MRI.Cube)-1)./(size(MRI.Cube).*MRI.Voxsize)) .* mriCoord) +1);
if ~isempty(find(sliceIndices > size(MRI.Cube))) | ~isempty(find(sliceIndices <= 0))
    sliceIndices = []; % Point effectively outside image
end


% --------------------------------------------------------------------
function varargout = LockCursor_Callback(h, eventdata, handles, varargin)

switch(get(h,'checked'))
    case 'on'
        set(h,'checked','off')
        set([findobj(handles.axa,'type','image'),findobj(handles.axc,'type','image'),findobj(handles.axs,'type','image')],...
            'Buttondownfcn','')
        set(gcbf,'WindowButtonMotionFcn','','WindowButtonDownFcn','','WindowButtonUpFcn','');
    case 'off'
        set(h,'checked','on')
        set(handles.Zoom,'Checked','off'); zoom off
        set([findobj(handles.axa,'type','image'),findobj(handles.axc,'type','image'),findobj(handles.axs,'type','image')],...
            'Buttondownfcn','bst_mriviewer(''move2cursor'')');
        set(gcbf,'WindowButtonMotionFcn','','WindowButtonDownFcn','bst_mriviewer(''move2cursor'')','WindowButtonUpFcn','bst_mriviewer(''stopmoving'')');
end

% --------------------------------------------------------------------
function [] = stopmoving
set(gcbf,'WindowButtonMotionFcn','','WindowButtonDownFcn','bst_mriviewer(''move2cursor'')')

% --------------------------------------------------------------------
function [] = move2cursor

handles = guihandles(gcbf);
MRI = getappdata(gcbf,'MRI');


switch get(gcbf,'CurrentAxes')
    case handles.axa
        
        tmp = get(handles.axa,'Currentpoint');
        tmp = tmp(1,1:2);
        Z = get(handles.slider_z,'Value');
        mriCoord = [tmp,Z];
        
    case handles.axs
        
        tmp = get(handles.axs,'Currentpoint');
        tmp = tmp(1,1:2);
        X = get(handles.slider_x,'Value');
        mriCoord =[X,tmp];
        
    case handles.axc
        
        tmp = get(handles.axc,'Currentpoint');
        tmp = tmp(1,1:2);
        Y = get(handles.slider_y,'Value');
        mriCoord = [tmp(1),Y,tmp(2)];
        
    otherwise
        return
end

newIndices = mm2vox(MRI,mriCoord); 

if isempty(newIndices), return, end  % Point is outisde image

show_MRI(MRI,'update',newIndices);

% % Update crosshair location if necessary
% if strcmp(get(findobj(gcbf,'tag','AddCrosshair'),'checked'),'on')
%     bst_mriviewer('AddCrosshair_Callback',findobj(gcbf,'tag','AddCrosshair'),'',handles,1); 
% end

set(gcbf,'WindowButtonMotionFcn','bst_mriviewer(''move2cursor'')',...
    'WindowButtonDownFcn','bst_mriviewer(''move2cursor'')')%,'DoubleBuffer','on')

% --------------------------------------------------------------------
function varargout = AddCrosshair_Callback(h, eventdata, handles, varargin)
% Add crosshair cursor to orthogonal views 
%if varargin == 1: show crosshair cursor anyhow


fig = findobj('Tag','MRIViewer','type','figure');
handles = guihandles(fig);
MRI = getappdata(fig,'MRI');
h = findobj(fig,'tag','AddCrosshair');


if isempty(varargin)
    switch(get(h,'checked'))
        case 'on'
            set(h,'checked','off') % Remove cursor
            
            axaXtick = [];
            axaYtick = [];
            
            axcXtick = [];
            axcYtick = [];
            
            axsXtick = [];
            axsYtick = [];
            varargin{1} = 0;  % we want off
        case 'off'
            set(h,'checked','on')
            varargin{1} = 1; % we want on
    end
end

if varargin{1} == 1
    set(h,'checked','on')
end

axaXtickLabel = '';
axaYtickLabel = '';
axcXtickLabel = '';
axcYtickLabel = '';
axsXtickLabel = '';
axsYtickLabel = '';
        
if strcmp(get(h,'checked'),'on')% Update crosshair location
    
    X = get(handles.slider_x,'Value');
    Y = get(handles.slider_y,'Value');
    Z = get(handles.slider_z,'Value');
    mriCoord = [X,Y,Z];
    
    axaXtick = [mriCoord(1)];
    axaXtickLabel = '';
    axaYtick = [mriCoord(2)];
    axaYtickLabel = '';
    
    axcXtick = [mriCoord(1)];
    axcXtickLabel = '';
    axcYtick = mriCoord(3);
    axcYtickLabel = '';
    
    axsXtick = mriCoord(2);
    axsXtickLabel = '';
    axsYtick = mriCoord(3);
    axsYtickLabel = '';
    
end

set(handles.axa,'xgrid','on','ygrid','on','Xtick',axaXtick,'XtickLabel',axaXtickLabel,...
    'Ytick',axaYtick,'YtickLabel',axaYtickLabel)
set(handles.axc,'xgrid','on','ygrid','on','Xtick',axcXtick,'XtickLabel',axcXtickLabel,...
    'Ytick',axcYtick,'YtickLabel',axcYtickLabel)
set(handles.axs,'xgrid','on','ygrid','on','Xtick',axsXtick,'XtickLabel',axsXtickLabel,...
    'Ytick',axsYtick,'YtickLabel',axsYtickLabel)


% --------------------------------------------------------------------
function varargout = Zoom_Callback(h, eventdata, handles, varargin)
% Toggle zoom state and adjust views (axes limits) in orthogonal views

% Check / uncheck zoom menu item
switch(get(handles.Zoom,'Checked'))
    case 'on'
        zoom off
        set(handles.Zoom,'Checked','off')
        set(handles.MRIViewer,'WindowButtonMotionFcn','','WindowButtonUpFcn','')
        %reset_axis([handles.axa,handles.axc,handles.axs])
    case 'off'
        set(handles.Zoom,'Checked','on')
        %set([handles.axa,handles.axc,handles.axs],'UIContextMenu',[])
        set(handles.LockCursor,'Checked','off')         
        zoom on
        set(handles.MRIViewer,'WindowButtonMotionFcn',' nowzooming','WindowButtonUpFcn',' nowzooming')
        
end


% --------------------------------------------------------------------

function varargout = PointList_Callback(h, eventdata, handles, varargin)

fig = findobj(0,'Tag','MRIViewer','Type','figure');

MRI = getappdata(fig,'MRI');
handles = guihandles(fig);
cPoint = get(handles.PointList,'Value');

% Show appropriate slices
try
    show_MRI(MRI,'update',mm2vox(MRI,MRI.Landmarks.MRImmXYZ(:,cPoint)'));
catch
    errordlg('Selected point is outside the MRI volume. Please check its coordinates or redefine fiducial points.')
    return
end


% Where are we ?
mmPoint = wherearewe(guihandles(fig)); % Current coordinates in mm
voxPoint = mm2vox(MRI,mmPoint);
bst_message_window({...
        'Current point:',...
        sprintf('%3.1f %3.1f %3.1f (mm)',mmPoint),...
        sprintf('%3.1f %3.1f %3.1f (voxels)',voxPoint)...
    })




% --------------------------------------------------------------------
function varargout = SelectPoint_Callback(h, eventdata, handles, varargin)
% Select current point in MR views and store as a landmark

MRIfig = findobj(0,'Tag','MRIViewer','Type','figure');

MRI = getappdata(MRIfig,'MRI');
mmPoint = wherearewe(handles); % Current coordinates in mm
voxPoint = mm2vox(MRI,mmPoint);
bst_message_window({...
        'Current point:',...
        sprintf('%3.1f %3.1f %3.1f (mm)',mmPoint),...
        sprintf('%3.1f %3.1f %3.1f (voxels)',voxPoint)...
    })

prompt  = {'Enter a label for current point:'}; 
title   = 'Please assign a label to the selected point';
lines = 1;
def     = {'Freak out !','Wazoo','Pumpkin',...
        'Zoot Allure','Call Any Vegetable','Creamcheese','Ruben','King Kong','Uncle Bernie'};
idef = randperm(length(def)); idef = idef(1);
answer  = inputdlg(prompt,title,lines,def(idef));

if isempty(answer)
    return
end     

% Update display
if size(mmPoint,1) ~= 3
    mmPoint = mmPoint';
end
ObjectHandles = plot_MRI(MRI,mmPoint);

% Update MRI.Landmarks structure
if ~isfield(MRI,'Landmarks')
    MRI.Landmarks.MRImmXYZ = [];
    MRI.Landmarks.Names = {};
    MRI.Landmarks.Handles = [];
end

for k = 1:size(mmPoint,2) % For each point
    tmp_XYZ = mm2vox(MRI,mmPoint(:,k))';
    if isempty(tmp_XYZ)
        warndlg(sprintf('Source # %d is outside the MRI cube.',k),'Source is discarded')
    else
        MRI.Landmarks.MRImmXYZ = [MRI.Landmarks.MRImmXYZ ,mmPoint(:,k)]; % 3 x N matrix    
    end
end

MRI.Landmarks.Handles = [MRI.Landmarks.Handles;ObjectHandles];

if isempty(MRI.Landmarks.Names)
    MRI.Landmarks.Names = answer;
else
    MRI.Landmarks.Names(end+1:end+length(answer)) = answer;
end

%switch(get(findobj(0,'Tag','ShowFiducials'),'checked'))
set(handles.PointList,'String',MRI.Landmarks.Names,'Value', length(MRI.Landmarks.Names))

setappdata(MRIfig,'MRI',MRI); 

bst_mriviewer('PointList_Callback','','', guihandles(MRIfig), '');

% --------------------------------------------------------------------
function varargout = DeletePoint_Callback(h, eventdata, handles, varargin)

% Delete current point in list

MRI = getappdata(gcbf,'MRI');
handles = guihandles(gcbf);
cPoint = get(handles.PointList,'Value');


% Update structure
%MRIPoint = getappdata(gcbf,'MRIPoint'); 
MRI.Landmarks.MRImmXYZ = MRI.Landmarks.MRImmXYZ(:,[1:end]~=cPoint);
%MRI.Landmarks.XYZvox = MRI.Landmarks.XYZvox(:,[1:end]~=cPoint);
MRI.Landmarks.Names = MRI.Landmarks.Names([1:end]~=cPoint);

% delete patch corresponding patch object in all 3 views
ndx = find(MRI.Landmarks.Handles(cPoint,:)); % non-zero handles
delete(MRI.Landmarks.Handles(cPoint,ndx))
MRI.Landmarks.Handles = MRI.Landmarks.Handles([1:end]~=cPoint,:);

setappdata(gcbf,'MRI',MRI);

set(handles.PointList,'String',MRI.Landmarks.Names,'value',length(MRI.Landmarks.Names))
bst_mriviewer('PointList_Callback','','','')



% --------------------------------------------------------------------
function varargout = RenamePoint_Callback(h, eventdata, handles, varargin)

% Rename current point in list
MRI = getappdata(gcbf,'MRI');
handles = guihandles(gcbf);
cPoint = get(handles.PointList,'Value');

prompt  = {'Enter a label for current point:'}; 
title   = 'Please assign a label to the selected point';
lines = 1;
answer  = inputdlg(prompt,title,lines,MRI.Landmarks.Names(cPoint));

if isempty(answer)
    return
end     

MRI.Landmarks.Names{cPoint} = answer{1};
set(handles.PointList,'String',MRI.Landmarks.Names)

setappdata(gcbf,'MRI',MRI)




% --------------------------------------------------------------------
function varargout = SeeAllPoints_Callback(h, eventdata, handles, varargin)

% Sell all patch points in figure
set(findobj(handles.MRIViewer,'type','patch'),'Visible','on')
set(findobj(handles.MRIViewer,'type','line'),'Visible','on')
set(handles.SeeAllPoints,'checked','on')
set(handles.HideAllPoints,'checked','off')


% --------------------------------------------------------------------
function varargout = HideAllPoints_Callback(h, eventdata, handles, varargin)

% Hide all patch points in figure
set(findobj(handles.MRIViewer,'type','patch'),'Visible','off')
set(findobj(handles.MRIViewer,'type','line'),'Visible','off')
set(handles.SeeAllPoints,'checked','off')
set(handles.HideAllPoints,'checked','on')


% --------------------------------------------------------------------
function varargout = Brighten_Callback(h, eventdata, handles, varargin)

clim = get(findobj(handles.MRIViewer,'type','axes'),'CLim');
set(findobj(handles.MRIViewer,'type','axes'),'CLim',.9*clim{1})


% --------------------------------------------------------------------
function varargout = Darken_Callback(h, eventdata, handles, varargin)

clim = get(findobj(handles.MRIViewer,'type','axes'),'CLim');
set(findobj(handles.MRIViewer,'type','axes'),'CLim',1.1*clim{1})



% --------------------------------------------------------------------
function varargout = LittleEndian_Callback(h, eventdata, handles, varargin)
% Flag for specification of machine format - default is little endian

switch get(h,'checked')
    case 'on'
        set(h,'checked','off')
    otherwise
        set(h,'checked','on')
end

