function mritool_cb(action) 
%MRITOOL_CB - OBSOLETE Switchyard corresponding to the OLDER MRITOOL GUI (mribox.fig)
% function mritool_cb(action) 
% Switchyard corresponding to the MRITOOL GUI (mribox.m)

%<autobegin> ---------------------- 14-Jun-2004 17:10:40 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_color_scheme.m
%   toolbox\bst_layout.m
%   toolbox\colnorm.m
%   toolbox\get_user_directory.m
%   toolbox\ginput2.m
%   toolbox\good_channel.m
%   toolbox\leaflet.m
%   toolbox\mritool_cb.m  NOTE: Routine calls itself explicitly
%   toolbox\norlig.m
%   toolbox\tri_interp.m
%   toolbox\vertices_connectivity.m
%
% Subfunctions in this file, in order of occurrence in file:
%   visu_pt(pthandle)
%   A = get_gain(fid,ndx);
%
% Application data and their calls in this file:
%   'TileType'
%   
%   setappdata(MAINWINDOW,'TileType','T');
%   
%
% Figure Files opened by this function:
%   'mribox.fig'
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure mribox.fig
%   text::Axes1Text1 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text1 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text2 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text2 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text3 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text3 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text4 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   text::Axes1Text4 "" uses ButtonDownFcn for ctlpanel SelectMoveResize
%   uicontrol:checkbox:big-endian "Big-endian" uses Callback for
%     mutincomp([findobj(gcbf,'Tag','big-endian'), findobj(gcbf,'Tag','little-endian')]);
%   uicontrol:checkbox:flip_lr_a "" uses Callback for mritool_cb flip_lr_a
%   uicontrol:checkbox:flip_ud_c "" uses Callback for mritool_cb flip_ud_c
%   uicontrol:checkbox:flip_ud_s "" uses Callback for mritool_cb flip_ud_s
%   uicontrol:checkbox:gotopoint "'Go to Point' mode" uses Callback for mritool_cb gotopoint
%   uicontrol:checkbox:little-endian "Little-endian" uses Callback for
%     mutincomp([findobj(gcbf,'Tag','big-endian'), findobj(gcbf,'Tag','little-endian')]);
%   uicontrol:checkbox:mouse_slct "Select Slice with Mouse" uses Callback for mritool_cb mouse_slct
%   uicontrol:checkbox:show_coordinates "Show Coordinates" uses Callback for
%     mritool_cb show_coordinates
%   uicontrol:checkbox:zoom "Activate Zoom" uses Callback for mritool_cb zoom
%   uicontrol:listbox:listbox "" uses Callback for mritool_cb current_point
%   uicontrol:pushbutton:BLOBS "Load SPM Clusters" uses Callback for mritool_cb addblobs
%   uicontrol:pushbutton:deletepoint "Delete" uses Callback for mritool_cb deletepoint
%   uicontrol:pushbutton:Hscalppoints "Load MR Volume" uses Callback for mritool_cb LoadBrainStorm
%   uicontrol:pushbutton:LEAFLET "Leaflet" uses Callback for mritool_cb leaflet
%   uicontrol:pushbutton:newlabel "New Label" uses Callback for mritool_cb newlabel
%   uicontrol:pushbutton:next ">" uses Callback for mritool_cb current_point
%   uicontrol:pushbutton:PCS "Define SCS" uses Callback for mritool_cb pcs
%   uicontrol:pushbutton:previous "<" uses Callback for mritool_cb current_point
%   uicontrol:pushbutton:pushbutton1 "Cortical Scout" uses Callback for mritool_cb CorticalScout
%   uicontrol:pushbutton:Quit "Quit" uses Callback for mritool_cb quit
%   uicontrol:pushbutton:scalppoints "Rendering" uses Callback for rendering_cb create
%   uicontrol:pushbutton:seelabels "Show Label" uses Callback for mritool_cb seelabels
%   uicontrol:pushbutton:seepoints "See/Hidel Points" uses Callback for mritool_cb seepoints
%   uicontrol:pushbutton:showfiducials "Show Fiducials" uses Callback for mritool_cb showfiducials
%   uicontrol:pushbutton:slctpoint "Brighten" uses Callback for mritool_cb brighten
%   uicontrol:pushbutton:slctpoint2 "Darken" uses Callback for mritool_cb darken
%   uicontrol:pushbutton:slctpoint20 "Select" uses Callback for mritool_cb slctpoint
%   uicontrol:pushbutton:slctpoint3 "Process fMRI Files" uses Callback for mritool_cb concatenate
%   uicontrol:pushbutton:slctpoint5 "Load Generators" uses Callback for mritool_cb readresults
%   uicontrol:radiobutton:ZoomButton1 "" uses Callback for mritool_cb ZoomButton
%   uicontrol:radiobutton:ZoomButton2 "" uses Callback for mritool_cb ZoomButton
%   uicontrol:radiobutton:ZoomButton3 "" uses Callback for mritool_cb ZoomButton
%   uicontrol:slider:1 "" uses ButtonDownFcn for mritool_cb sliders
%   uicontrol:slider:1 "" uses Callback for mritool_cb sliders
%   uicontrol:slider:2 "" uses ButtonDownFcn for mritool_cb sliders
%   uicontrol:slider:2 "" uses Callback for mritool_cb sliders
%   uicontrol:slider:3 "" uses ButtonDownFcn for mritool_cb sliders
%   uicontrol:slider:3 "" uses Callback for mritool_cb sliders
%   uimenu:: "" uses Callback for mritool_cb fifmri
%   uimenu::FileImport MR Scans.mri formatuimenu1 "" uses Callback for mritool_cb mrislice
%   uimenu::FileImport MRI Data.mri format1 "" uses Callback for mritool_cb mriformat
%   uimenu::FileRead MR Scansuimenu1 "" uses Callback for mritool_cb saveaxial
%   uimenu::Fileuimenu3 "" uses Callback for mritool_cb SaveAsBrainstorm
%   uimenu::Fileuimenu4 "" uses Callback for mritool_cb LoadBrainStorm
%   uimenu::img format "" uses Callback for mritool_cb imgformat
%   uimenu::Import/Export Cortical Remappinguimenu1 "" uses Callback for mritool_cb saveas_points
%   uimenu::Import/Export Import MEG Cortical Mapsuimenu1 "" uses Callback for
%     mritool_cb cortic_remap
%   uimenu::Import/Export Read MR ScansLoad .mri formatuimenu1 "" uses Callback for mritool_cb ctfmri
%   uimenu::Import/Export uimenu1 "" uses Callback for mritool_cb MEGBlobs
%   uimenu::IMPORT_mridata "" uses Callback for mritool_cb IMPORT_mridata 
%
% At Check-in: $Author: Mosher $  $Revision: 22 $  $Date: 6/14/04 3:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-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> ------------------------ 14-Jun-2004 17:10:40 -----------------------

%-- Sylvain Baillet, Ph.D. 
%-- Signal and Image Processing Institute
%-- University of Southern Califonia
%-- Los Angeles - USA

% History ----------------------------------------------------------------------
% 10/25/00 JCM    added alphabetization to mri slice routine
% 04/20/01 SB -   Modified routine to import CTF format  
% JCM 08-Sep-2003 Commenting
% SB  15-Dec-2003 Updated color scheme
% JCM 19-May-2004 Change use of variable "dir" to "rootdir" line 4053
% ------------------------------------------------------------------------------

frontcolor = [.8 .8 0];
TASKBAR = findobj(0,'Tag','TASKBAR');

DEFAULTslice_thickness = 2;
DEFAULTsiz = 256;
line_color = [.7 .7 .0]; % Line color for coordinates display
FORMATtick = '%2.1f';
NGRAY_MAX = 255; % Max gray level value for the images

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

switch action
    
case 'create'
    MAINWINDOW = [];%findobj(0,'Tag','mainwindow');
    if isempty(MAINWINDOW)
        MAINWINDOW = openfig('mribox.fig','reuse');
        %MAINWINDOW = findobj(0,'Tag','mainwindow');
        % bst_color_scheme(MAINWINDOW); % doesn't work well
        setappdata(MAINWINDOW,'TileType','T');
        bst_color_scheme(MAINWINDOW)
        bst_layout('align',MAINWINDOW,1,1,1);
        set(MAINWINDOW,'visible','on')
    end
    
    %-----------------------------------------------------------------------------------------------
    
case 'read_file_scan'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    listboxx = findobj(0,'Tag','listbox');
    set(listboxx,'visible','off')
    listboxx_title = findobj(MAINWINDOW,'Tag','listboxx_title');
    set(listboxx_title,'visible','off')
    legend_sources = findobj(MAINWINDOW,'Tag','legend_sources');
    set(legend_sources,'visible','off')
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    ax = findobj(MAINWINDOW,'Tag','colorbar');
    if ~isempty(ax)
        set(ax,'Visible','off')
    end
    
    set(0,'CurrentFigure',   MAINWINDOW )
    mri_scan = get(gcf,'Userdata');
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    mri_scan.ptname = {};
    mri_scan.otherfig = [];
    mri_scan.ptcoord = [];
    mri_scan.pthandle  = [];
    mri_scan.ptlabels  = [];
    mri_scan.gotopoint = 0;
    mri_scan.seepoints = 0;
    mri_scan.seelabels = 1;
    %  mri_scan.pcs = 0;
    
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    filename = mri_scan.filemri;
    set(file_mri,'Visible','on','String',filename);
    
    header =  mri_scan.header ;
    set(Header,'String',int2str(header))
    
    siz = mri_scan.siz;
    
    set(image_size,'String',num2str(siz,'%1d %1d'));
    if isempty(siz)
        siz = [DEFAULTsiz,DEFAULTsiz];
        h = warndlg('Image size field is empty - default value is [256 256]');
        uiwait(h);
    end
    if length(siz)==1
        siz = [siz,siz];
    end
    
    scan_type = mri_scan.scantype;
    
    fov = mri_scan.fov ;
    set(FOV,'String',num2str([fov(1),fov(2)],'%1.1f %1.1f'));
    if isempty(fov)
        fov = siz;
    end
    
    slice_thickness = mri_scan.slice_thick;
    set(SLICE_thickness ,'String',num2str(slice_thickness));
    if isempty(slice_thickness)
        slice_thickness  = DEFAULTslice_thickness;
        warndlg(['Slice thickness field empty - using default',num2str(DEFAULTslice_thickeness)])   
    end
    
    
    nslice = mri_scan.nslice;
    set(NSLICE,'String',int2str(nslice));
    
    
    if isempty(nslice)
        errordlg('Please enter the number of slices')
        return
    end
    set(NSLICE,'String',num2str(nslice));
    
    siz_tmp = siz;
    
    clear siz_tmp scan_type_tmp
    
    fov(3) = nslice*slice_thickness;
    
    code = mri_scan.code;
    
    
    if (~mri_scan.readslice)& ~isfield(mri_scan,'mri')
        
        file_mri = fopen(filename,'r',MACHINEFORMAT);
        if file_mri < 1
            disp('Invalid File Name')
            return;
        end
        
        fseek(file_mri,header,'bof');
        
        mri_scan.mri = uint8(zeros(siz(1),siz(2),nslice));
        
        h= waitbar(0,'Please wait...');
        for i = 1:nslice
            if feof(file_mri)
                errordlg('Please Check Image size / Number of Slices : do not match file size')
                return
            end
            mri_scan.mri(:,:,i) = uint8(fread(file_mri,[siz],code));
            waitbar(i/nslice)   ;
            
        end   
        fclose(file_mri);
        delete(h);
    else
        %   mri = mri_scan.mri;
    end
    
    pix_siz = [fov(1)/(siz(1)-1) fov(2)/(siz(2)-1) slice_thickness];
    
    if strcmp(scan_type,'Coronal')
        mri_scan.mri = permute(mri_scan.mri,[1,3,2]);   
    elseif strcmp(scan_type,'Sagittal')   
        mri_scan.mri = permute(mri_scan.mri,[3,1,2]);   
    end
    
    tmp = size(mri_scan.mri);
    siz = tmp(1:2);
    nslice = tmp(3);
    
    
    if strcmp(scan_type,'Coronal')
        pix_siz  = pix_siz([3,2,1]);   
        fov = fov([3,2,1]);
    elseif strcmp(scan_type,'Sagittal')
        pix_siz  = pix_siz([1,3,2]);   
        fov = fov([1,3,2]);
    end
    
    mri_scan.pix_siz = pix_siz;
    mri_scan.fov = fov;
    
    ic = round(siz(2)/2); % Image indices for coronal,sagittal and axial slices respectively
    is = round(siz(1)/2);
    ia = round(nslice/2);
    
    
    mri_scan.ic = ic;
    mri_scan.is = is;
    mri_scan.ia = ia;
    
    
    % Delete all previous points
    delete(findobj(get(axa,'children'),'type','patch'))
    delete(findobj(get(axc,'children'),'type','patch'))
    delete(findobj(get(axs,'children'),'type','patch'))
    delete(findobj(get(axa,'children'),'type','line'))
    delete(findobj(get(axc,'children'),'type','line'))
    delete(findobj(get(axs,'children'),'type','line'))
    
    
    ngray = 256;
    %M = (max(max(mri_scan.mri(:,:,round(size(mri_scan.mri,3)/2)))));
    
    M = double(max(mri_scan.mri(:))); % May be rather time consuming if limited RAM is availlable
    
    if M > (ngray-1)    
        % Scale gray levels to 255 to save the MRI in uint8
        % Memory requirememts necessitate a slice-by-slice loop on most machines
        hw = waitbar(0,'Scaling the MRI to 256 gray levels...');
        for k = 1:size(mri_scan.mri,3)
            waitbar(k/size(mri_scan.mri,3),hw);
            tmp = double(mri_scan.mri(:,:,k)) ;
            mri_scan.mri(:,:,k) = int16(tmp * ((ngray - 1 )/M)) ;
        end
        
        mri_scan.mri = uint8(mri_scan.mri);
        delete(hw)          
        drawnow
    end
    
    
    %% -- Axial View
    axes(axa)
    ax = imagesc([0, fov(1)],[0,fov(2)],(mri_scan.mri(:,:,ia)),'Parent',axa);    
    set(axa,'Tag','axa','Units','Normal');
    set(axa,'visible','on','xcolor',line_color ,'Ycolor',line_color,'Color','k');   
    mri_scan.ax = ax;
    mri_scan.axa = axa;
    set(axa,'box','on','layer','top','xgrid','on','xtick',[],'ygrid','on','ytick',[],'gridlinestyle','-')
    axis([0 fov(1), 0 fov(2)])  
    set(ax,'erasemode','none');
    caxis([0 ngray-1])
    axis equal
    %% -- Axial view Done
    
    
    
    %% -- Sagittal view
    axes(axs) 
    sag = imagesc([0 fov(1)],[0 fov(3)],((fliplr(squeeze(mri_scan.mri(is,:,:))))'),'Parent',axs);
    set(axs,'Tag','axs','Units','Normal');
    mri_scan.sag = sag;
    mri_scan.axs = axs;
    set(axs,'color',[0 0 0],'box','on','layer','top','xgrid','on','xtick',[],'ygrid','on','ytick',[],'gridlinestyle','-')
    set(axs,'visible','on','xcolor',line_color ,'Ycolor',line_color,'DataAspectRatioMode','auto' );   
    axis([0 fov(1), 0 fov(3)])  
    set(sag,'erasemode','none');
    set(axs,'Xlim', get(axa,'Xlim'))
    caxis([0 ngray-1])
    axis equal
    %% -- Sagittal view Done
    
    %% -- Coronal view
    axes(axc)
    cor = imagesc([0 fov(2)],[0 fov(3)],(fliplr(flipud(squeeze(mri_scan.mri(:,ic,:))'))),'Parent',axc);
    set(axc,'Tag','axc','Units','Normal');
    mri_scan.cor = cor;
    mri_scan.axc = axc;
    set(axc,'color',[0 0 0],'box','on','layer','top','xgrid','on','xtick',[],'ygrid','on','ytick',[],'gridlinestyle','-')
    set(axc,'visible','on','xcolor',line_color ,'Ycolor',line_color,'DataAspectRatioMode','auto');   
    axis([0 fov(2), 0 fov(3)])  
    set(cor,'erasemode','none');
    set(axc,'Xlim', get(axa,'Ylim'),'Ylim', get(axs,'Ylim'))
    caxis([0 ngray-1])
    axis equal
    
    axes(axa)
    set(axa,'Xlim', get(axs,'Xlim'),'Ylim', get(axc,'Xlim'))
    axis equal
    axes(axs)
    set(axs,'Xlim', get(axa,'Xlim'),'Ylim', get(axc,'Ylim'))
    axis equal
    axes(axc)
    set(axc,'Xlim', get(axa,'Ylim'),'Ylim', get(axs,'Ylim'))
    axis equal
    %% -- Coronal view Done
    
    eval(['colormap bone(',int2str(ngray),');']);
    
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    set(s1,'Visible','on','Min', 0,'Max',fov(2),'Value',is*pix_siz(2))
    set(s2,'Visible','on','Min', 0,'Max',fov(3),'Value',ia*pix_siz(3))
    set(s3,'Visible','on','Min', 0,'Max',fov(1),'Value',ic*pix_siz(1))
    
    
    ind = 1;
    
    set(MAINWINDOW,'Userdata',mri_scan)   
    set([sag,ax,cor],'ButtonDownFcn','mritool_cb gotopoint')   
    
    %-----------------------------------------------------------------------------------------------
    
case 'brighten'
    MAINWINDOW = gcbf;
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    
    clim = get(axs,'CLim');
    set(axs,'CLim',.9*clim)
    set(axa,'CLim',.9*clim)
    set(axc,'CLim',.9*clim)
    
    %-----------------------------------------------------------------------------------------------
    
case 'darken'
    
    MAINWINDOW = gcbf;
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    
    clim = get(axs,'CLim');
    set(axs,'CLim',1.1*clim)
    set(axa,'CLim',1.1*clim)
    set(axc,'CLim',1.1*clim)
    
    
    %-----------------------------------------------------------------------------------------------
    
case 'importANALYZE' % Import anatomical image in ANALYZE format  - still needs to be done (I've just cut and paste the importBlobs section below)
    
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    MAINWINDOW = gcbf;
    set(0,'CurrentFigure',MAINWINDOW )
    mri_scan = get(gcbf,'Userdata');
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    
    % -----
    
    fig = MAINWINDOW;
    
    %set(fig,'colormap',bone(64))
    
    %% ---
    Current = get_user_directory;
    cd(Current.SUBJECTS)   
    
    [anatfile, anatpath] = uigetfile('*.img','Please select the high-resolution MRI file (.img)' );    
    if (anatfile == 0), return, end
    cd(anatpath);
    eval('V = spm_vol(anatfile);','ok=0;');
    Vorig = V.mat;
    clear V
    
    [blobfile, blobpath] = uigetfile('*.img','Please select the file containing functional blobs');    
    if (blobfile == 0), return, end
    cd(blobpath);
    
    eval('V = spm_vol(blobfile);','ok=0;');
    
    
    %---------------------------------------
    % Read img volume
    
    header = 0;
    mri_scan.blob.header = header;
    
    siz = [V.dim(1), V.dim(2)];
    mri_scan.blob.siz = siz;
    
    mri_scan.blob.scantype= 'Axial';    
    
    mri_scan.blob.nslice = V.dim(3);
    
    mri_scan.code = 'uint8'; % That's all folks !
    
    blob = zeros(siz(1),siz(2),mri_scan.blob.nslice);
    
    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
    disp('Reslicing functional scans...')
    [XI,YI] = meshgrid(linspace(1,size(blob,1),mri_scan.siz(1)),linspace(1,size(blob,2),mri_scan.siz(2)));
    iblob = zeros(mri_scan.siz(2),mri_scan.siz(1),size(blob,3));
    for i=1:size(blob,3)
        iblob(:,:,i) = interp2(blob(:,:,i),XI,YI,'*linear');
    end
    scale_func = (size(blob));
    clear blob
    % 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(iblob,3)); 
    
    iiblob = zeros(size(iblob,1),size(iblob,2),n_slice);
    iiblob(:,:,1:floor(n_slice/2)) = repmat(iblob(:,:,1),[1,1,length(1:floor(n_slice/2))]);
    iiblob(:,:,floor(n_slice/2)+1:(n_slice)) = repmat(iblob(:,:,2),[1,1,length(floor(n_slice/2)+1:(n_slice))]);
    
    clear iblob
    
    % 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 
    
    ROTmri2mm = (V.mat(1:3,1:3));
    TRANSmri2mm = V.mat(1:3,4);
    
    Vorig = inv(Vorig);
    ROTmm2mri = Vorig(1:3,1:3);
    TRANSmm2mri = Vorig(1:3,4);
    
    
    % Rescale in original voxel units
    %scale_func = (size(blob));
    scale_anat= (size(iiblob));
    scale = diag( scale_func./scale_anat);
    xyz = scale*xyz ;
    
    xyz =  ROTmri2mm* xyz + repmat((TRANSmri2mm),1, size(xyz,2));
    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;
    
    
    %iiiblob = zeros(size(mri_scan.mri,1),size(mri_scan.mri,2),size(mri_scan.mri,3));
    sizz = size(mri_scan.mri);
    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 
    %iiiblob(tt(t>0*mean(t))) =  64+(t(t>0*mean(t))+eps)*64/max(t);
    mri_scan.mri(tt(t>0*mean(t))) =  64+(t(t>0*mean(t))+eps)*64/max(t);
    
    %%%----- Superimpose to regular MRI
    %  mri_scan.mri = mri_scan.mri + iiiblob;
    
    
    [Mt,ind] = max(t);
    tmp = mri_scan.blob.Cube(:,ind);
    
    clear blob iblob iiblob iiiblob
    
    ia = tmp(3);
    ic = tmp(2);
    is = tmp(1);
    
    ax = mri_scan.ax;
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    set(sag ,'CData',(flipud(squeeze(mri_scan.mri(is,:,:))')),'erasemode','normal');
    set(get(sag,'Parent'),'CLim',[0 128])
    set(cor ,'CData',(fliplr(flipud(squeeze(mri_scan.mri(:,ic,:))'))),'erasemode','normal');
    set(get(cor,'Parent'),'CLim',[0 128])
    set(ax ,'CData',(mri_scan.mri(:,:,ia)),'erasemode','normal');    
    set(get(ax,'Parent'),'CLim',[0 128])
    
    
    set(MAINWINDOW,'Userdata',mri_scan);
    
    split(1:64,:) = bone(64);
    split(65:128,:) = hot(64);
    
    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)
    
    %-----------------------------------------------------------------------------------------------
    
    
    
case 'addblobs' % Add functional blobs as written in a .img (SPM compliant) format
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    MAINWINDOW = gcbf;
    set(0,'CurrentFigure',MAINWINDOW )
    mri_scan = get(gcbf,'Userdata');
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    
    % -----
    
    fig = MAINWINDOW;
    
    %set(fig,'colormap',bone(64))
    
    %% ---
    Current = get_user_directory;
    cd(Current.SUBJECTS)   
    
    [anatfile, anatpath] = uigetfile('*.img','Please select the high-resolution MRI file (.img)' );    
    if (anatfile == 0), return, end
    cd(anatpath);
    eval('V = spm_vol(anatfile);','ok=0;');
    Vorig = V.mat;
    clear V
    
    [blobfile, blobpath] = uigetfile('*.img','Please select the file containing functional blobs');    
    if (blobfile == 0), return, end
    cd(blobpath);
    
    eval('V = spm_vol(blobfile);','ok=0;');
    
    
    %---------------------------------------
    % Read img volume
    
    header = 0;
    mri_scan.blob.header = header;
    
    siz = [V.dim(1), V.dim(2)];
    mri_scan.blob.siz = siz;
    
    mri_scan.blob.scantype= 'Axial';    
    
    mri_scan.blob.nslice = V.dim(3);
    
    mri_scan.code = 'uint8'; % That's all folks !
    
    blob = zeros(siz(1),siz(2),mri_scan.blob.nslice);
    
    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
    disp('Reslicing functional scans...')
    [XI,YI] = meshgrid(linspace(1,size(blob,1),mri_scan.siz(1)),linspace(1,size(blob,2),mri_scan.siz(2)));
    iblob = zeros(mri_scan.siz(2),mri_scan.siz(1),size(blob,3));
    for i=1:size(blob,3)
        iblob(:,:,i) = interp2(blob(:,:,i),XI,YI,'*linear');
    end
    scale_func = (size(blob));
    clear blob
    % 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(iblob,3)); 
    
    iiblob = zeros(size(iblob,1),size(iblob,2),n_slice);
    iiblob(:,:,1:floor(n_slice/2)) = repmat(iblob(:,:,1),[1,1,length(1:floor(n_slice/2))]);
    iiblob(:,:,floor(n_slice/2)+1:(n_slice)) = repmat(iblob(:,:,2),[1,1,length(floor(n_slice/2)+1:(n_slice))]);
    
    clear iblob
    
    % 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 
    
    ROTmri2mm = (V.mat(1:3,1:3));
    TRANSmri2mm = V.mat(1:3,4);
    
    Vorig = inv(Vorig);
    ROTmm2mri = Vorig(1:3,1:3);
    TRANSmm2mri = Vorig(1:3,4);
    
    
    % Rescale in original voxel units
    %scale_func = (size(blob));
    scale_anat= (size(iiblob));
    scale = diag( scale_func./scale_anat);
    xyz = scale*xyz ;
    
    xyz =  ROTmri2mm* xyz + repmat((TRANSmri2mm),1, size(xyz,2));
    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;
    
    
    %iiiblob = zeros(size(mri_scan.mri,1),size(mri_scan.mri,2),size(mri_scan.mri,3));
    sizz = size(mri_scan.mri);
    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 
    %iiiblob(tt(t>0*mean(t))) =  64+(t(t>0*mean(t))+eps)*64/max(t);
    mri_scan.mri(tt(t>0*mean(t))) =  64+(t(t>0*mean(t))+eps)*64/max(t);
    
    %%%----- Superimpose to regular MRI
    %  mri_scan.mri = mri_scan.mri + iiiblob;
    
    
    [Mt,ind] = max(t);
    tmp = mri_scan.blob.Cube(:,ind);
    
    clear blob iblob iiblob iiiblob
    
    ia = tmp(3);
    ic = tmp(2);
    is = tmp(1);
    
    ax = mri_scan.ax;
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    set(sag ,'CData',(flipud(squeeze(mri_scan.mri(is,:,:))')),'erasemode','normal');
    set(get(sag,'Parent'),'CLim',[0 128])
    set(cor ,'CData',(fliplr(flipud(squeeze(mri_scan.mri(:,ic,:))'))),'erasemode','normal');
    set(get(cor,'Parent'),'CLim',[0 128])
    set(ax ,'CData',(mri_scan.mri(:,:,ia)),'erasemode','normal');    
    set(get(ax,'Parent'),'CLim',[0 128])
    
    
    set(MAINWINDOW,'Userdata',mri_scan);
    
    split(1:64,:) = bone(64);
    split(65:128,:) = hot(64);
    
    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)
    
    %-----------------------------------------------------------------------------------------------
    
case 'cortic_remap'
    
    %%% ---- Cotical Remapping of parametric sources on the cortex 
    
    MAINWINDOW = gcbf;
    mri_scan = get(MAINWINDOW ,'Userdata');
    if isempty(mri_scan), return, end
    
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    
    
    User = get_user_directory;
    cd(User.STUDIES)
    
    % Load result file
    
    [result_file,result_dir] = uigetfile('*result*.mat','Load Result File...'); 
    if result_file == 0, return, end
    
    cd(result_dir)
    Results = load(result_file);
    
    nsource = length(Results.SourceLoc);
    
    % Look for synchronous sources
    isync = find(Results.SourceOrder==2 | Results.SourceOrder==3 |Results.SourceOrder==4);
    insync = setdiff(1:nsource,isync);
    SourceOrder_old = Results.SourceOrder;
    
    trans = zeros(size(Results.SourceOrder,2),nsource);
    
    i = 1;
    for src = 1:length(SourceOrder_old)
        if ~isempty(isync) & ~isempty(find(isync==src))
            Results.SourceOrder(i) = SourceOrder_old(src);
            Results.SourceOrder(i+1) = SourceOrder_old(src);
            ptname{i} = int2str(src);
            ptname{i+1} = int2str(src);
            SourceLoc(:,i:i+1) = Results.SourceLoc{src};
            COLOR(i,:) = [min(.3+abs(rand),1) min(.3+abs(rand),1) min(.3+abs(rand),1) ];
            COLOR(i+1,:) = COLOR(i,:);
            i = i+2;
        else
            Results.SourceOrder(i) = SourceOrder_old(src);
            ptname{i} = int2str(src);
            SourceLoc(:,i) = Results.SourceLoc{src};
            COLOR(i,:) = [min(.3+abs(rand),1) min(.3+abs(rand),1) min(.3+abs(rand),1) ];
            i = i+1;
        end
    end
    
    nsource = size(SourceLoc,2);
    
    j = 1;
    
    for i = 1:nsource
        ptcoord(i,:) = 1000*SourceLoc(:,i)';
        icoord(i,:) = [0 0 0];
        
        if mri_scan.pcs == 0
            orig = findstr(ptname{i},'originPCS');
            if ~isempty(orig) % This set of points can be registered owing to a MEG-EEG PCS.
                mri_scan.pcs = 1; 
                mri_scan.pcs_type = ptname{i}(findstr(ptname{i},'_')+1:end);
                mri_scan.center = icoord;
                center = icoord(i,:);
                ptcenter = ptcoord(i,:)
            else 
                mri_scan.pcs = 0; 
                mri_scan.pcs_type = 'NONE';
                center = [0 0 0];
            end
        end
        
    end
    fclose all;
    disp('')
    
    % Load result file : END
    
    % Load tessellation information
    cd(User.SUBJECTS)
    load(Results.StudySubject.Subject)
    
    % Load Head Model: is a cortical gain matrix available ?
    cd(User.STUDIES)
    load(Results.StudySubject.HeadModel,'ImageGain')
    if isempty([ImageGain{:}]), errodlg('No cortical gain matrix was defined - remapping aborted'), return, end
    flag = zeros(1,length(ImageGain));
    for k = 1:length(ImageGain)
        flag(k) = ~isempty(ImageGain{k});
    end
    igain = find(flag);
    ImageGain = ImageGain{k};
    
    if isempty(Tesselation), errordlg('Subject has no cortical tessellation - remapping aborted'), return, end
    cd(User.SUBJECTS)
    load(Tesselation)
    
    disp(' ')
    fprintf('%s\t',Comment{igain})
    fprintf('%\n')
    disp(' ')
    %itess = input(['Which cortical tessellation would you like to use here ? (',deblank(num2str(igain)),') : ']);
    itess = 5; % CHEAT for Kevin  -for silvin : 5
    
    tess.vertices = Vertices{itess}'; clear Vertices
    tess.faces = Faces{itess}; clear Faces
    Comment = Comment{itess};
    
    % Load ImageGain matrix
    [PATH,NAME,EXT,VER] = fileparts( Results.StudySubject.Study); 
    file = fopen(fullfile(User.STUDIES,PATH,ImageGain{:} ),'r',MACHINEFORMAT);
    if file <0, errordlg('Unable to open the head model file'), return, end
    G = get_gain(file, 1:size(tess.vertices,1));
    
    % Channels
    load(fullfile(User.STUDIES,Results.StudySubject.Channel))
    ndx = good_channel(Channel,Results.ChannelFlag,'meg'); % CHEAT meg only here
    
    G = G(ndx,:);
    verbose = 1; 
    % Vertex Connectivity
    if isempty(VertConn)
        disp('Computing Vertex Connectivity on for this tessellation')
        drawnow
        vertconn = vertices_connectivity(tess,verbose);
        disp('Saving Connectivity Information...')
        drawnow
        [PATH,NAME,EXT,VER] = fileparts(Results.StudySubject.Subject); 
        VertConn = fullfile(PATH,'VertConn.mat');
        save(fullfile(User.SUBJECTS,VertConn),'vertconn') ;
        disp('Updating Subject Information...')
        drawnow
        save(Results.StudySubject.Subject,'VertConn','-append')
        disp('Updating Subject Information...DONE')
        drawnow
    else
        load(VertConn)
    end
    
    
    % Cortical Remapping
    disp(' ')
    fprintf('Remapping Sources on Tessellation %s \n',Comment)
    nseeds = 1;
    for k = 1:size(Results.IndepTopo,2)
        [amps{k},csources{k}] = corti_remap(tess,vertconn,G,Results.IndepTopo(:,k),nseeds,verbose);
    end
    amps = [amps{:}];
    
    mri_scan.M = max(abs(Results.TimeSeries(:)));
    
    mri_scan.climm = ceil(max(get(axa,'Clim')));
    Results.TimeSeries = Results.TimeSeries*mri_scan.climm/mri_scan.M;
    
    for k = 1:length(csources)
        mri_scan.TimeSeries{k} = Results.TimeSeries(:,k);
    end
    
    icsources = [csources{:}];
    
    
    %%%% Interpolation of the Cortical Current Density Image
    if ~exist('interp_mat.mat','file')
        % Transform PCS coordinates in the MRI volume reference system
        switch mri_scan.pcs_type
        case {'NEUROMAG','CTF'}
            if ~isfield(mri_scan,'rot')
                errordlg('Please Define a PCS First')
                return
            end
            center = mri_scan.center;
            rot = mri_scan.rot;
            mri_scan.mript = (inv(rot)*tess.vertices')';      
        end
        
        mri_scan.mript (:,1) = mri_scan.mript (:,1) + center(1);
        mri_scan.mript (:,2) = mri_scan.mript (:,2) + center(2);
        mri_scan.mript (:,3) = mri_scan.mript (:,3) + center(3);
        clear icoord
        icoord(:,1) = round((get(s1,'Max')-mri_scan.mript (:,2))/mri_scan.pix_siz(2)+1);
        icoord(:,2) = round(mri_scan.mript (:,1)/mri_scan.pix_siz(1)+1);
        icoord(:,3) = round((get(s2,'Max')-mri_scan.mript (:,3))/mri_scan.pix_siz(3)+1);
        [H] = tri_interp(size(mri_scan.mri),icoord, tess.faces)   ;
        disp('Saving Interpolation Matrix')
        User = get_user_directory;
        save interp_mat H
        clear tess vertconn G Results
    else
        clear tess vertconn G Results
        load interp_mat
    end
    
    H = H(:,icsources);
    
    mri_scan.GridAmp = ones(size(icsources))';%abs(amps(csources)); clear amps 
    %mri_scan.GridAmp  = norlig(mri_scan.GridAmp ); % See integrated activity at each vertex along time
    mri_scan.csources = csources;
    
    mri_scan.M = max(mri_scan.GridAmp);
    
    mri_scan.mri(mri_scan.mri >= mri_scan.climm) = mri_scan.climm;
    mri_scan.GridAmp = mri_scan.GridAmp*mri_scan.climm/mri_scan.M;
    mri_scan.thres = 0*mri_scan.climm; % Threshold for visualization of current density
    
    %CInterp = sparse(size(H,1),1,1000*length(csources));
    CInterp = sparse(H*mri_scan.GridAmp);% OK If enough memory is available
    
    if 0 
        tmp = sum(H,2);
        ind = find(tmp);
        blocks = unique(round(linspace(1,length(ind),100)));
        hwait = waitbar(0,'Interpolating the volumic curreny density...' );
        for k = 1:length(blocks)
            waitbar(blocks(k)/blocks(end),hwait)
            if k > 1
                vec = ind(blocks(k-1)+1:blocks(k));
            else
                vec = ind(1:blocks(1));
            end
            CInterp(vec) = sparse(H(vec , :)*mri_scan.GridAmp);
        end
        delete(hwait)
        ind2 =find(CInterp(ind)>=mri_scan.thres);
        mri_scan.mri(ind(ind2)) = 1.01*mri_scan.climm + nonzeros(CInterp(ind(ind2)));
    end
    
    ind =find(CInterp>mri_scan.thres);
    mri_scan.indmodif = ind; % indices of the voxels with modified grey levels
    mri_scan.gl_orig = mri_scan.mri(ind);
    
    %mri_scan.mri(ind) = round(1.01*mri_scan.climm + nonzeros(CInterp(ind)));
    mri_scan.mri(ind) = double(mri_scan.mri(ind)) + .9*nonzeros(CInterp(ind));
    minn = double(min(mri_scan.mri(ind)));
    mri_scan.mri(ind) =  double(mri_scan.mri(ind)) - minn + mri_scan.climm+1;
    maxx = double(max(mri_scan.mri(ind)));
    mri_scan.climm = round(maxx/2);
    %mri_scan.mri(ind) = log(double(mri_scan.mri(ind))+1); % Make the cortical remaps look more yellow
    
    % Assign time series to vertices
    mri_scan.GridAmp  = zeros(size(H,2),size(mri_scan.TimeSeries{1},1));
    
    for k = 1:length(mri_scan.csources)
        mri_scan.GridAmp(mri_scan.csources{k},:) = repmat(mri_scan.TimeSeries{k}',length(mri_scan.csources{k}),1);
    end
    mri_scan.GridAmp =  mri_scan.GridAmp([mri_scan.csources{:}],:);
    
    mags = mri_scan.GridAmp;
    
    save kevin_remap mags csources
    
    clear mags
    
    clear H mri_scan.GridAmp
    
    [M tmp] = max(mri_scan.mri(ind));
    [is ic ia] = ind2sub(size(mri_scan.mri),ind((tmp)));
    
    ax = mri_scan.ax;
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    set(sag ,'CData',(flipud(squeeze(mri_scan.mri(is,:,:))')),'erasemode','none');
    set(get(sag,'Parent'),'CLim',[0 2*mri_scan.climm])
    set(cor ,'CData',(fliplr(flipud(squeeze(mri_scan.mri(:,ic,:))'))),'erasemode','none');
    set(get(cor,'Parent'),'CLim',[0 2*mri_scan.climm])
    set(ax ,'CData',(mri_scan.mri(:,:,ia)),'erasemode','none');    
    set(get(ax,'Parent'),'CLim',[0  2*mri_scan.climm])
    
    split(1:mri_scan.climm,:) = bone(mri_scan.climm);
    split(mri_scan.climm+1:2*mri_scan.climm,:) = hot(mri_scan.climm);
    
    mri_scan.split = split;
    
    set(MAINWINDOW,'Userdata',mri_scan);
    
    split(1:mri_scan.climm,:) = bone(mri_scan.climm);
    split(mri_scan.climm+1:2*mri_scan.climm,:) = hot(mri_scan.climm);
    
    ax = findobj(MAINWINDOW,'Tag','colorbar');
    axes(ax)
    colormapp = split( mri_scan.climm+1: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)
    
    %-----------------------------------------------------------------------------------------------
    
case 'cortic_remap_update'
    
    %%% ---- Update of the Cortical Remapping of parametric sources on the cortex 
    
    MAINWINDOW = gcf;
    mri_scan = get(MAINWINDOW ,'Userdata');
    if isempty(mri_scan), return, end
    
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    mri_scan.is = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2));
    mri_scan.ia= round(get(s2,'Value')/pix_siz(3));
    mri_scan.ic = round(get(s3,'Value')/pix_siz(1));
    
    % Load interpolation matrix  
    User = get_user_directory;
    cd(User.SUBJECTS)
    load interp_mat
    H = H(:,[mri_scan.csources{:}]);
    
    mri_scan.thres = 0*mri_scan.climm; % Threshold for visualization of current density
    
    %CInterp = sparse(size(H,1),size(mri_scan.GridAmp,2));
    % clear H  
    
    ax = mri_scan.ax;
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    %figure
    %plot([mri_scan.TimeSeries{:}])
    M = norlig(mri_scan.GridAmp)';
    CInterp = sparse(H*M);
    indd = find(CInterp>.2*max(M));
    CInterp = sparse(H(indd,:)*abs(mri_scan.GridAmp));
    
    fig = figure;
    set(gcf,'Color','k','Units','normal','position',[.05 .5 .80 .3])
    
    axa2 = subplot(133);
    axc2 = subplot(132);
    axs2 = subplot(131);
    
    
    pax = 1*get(axa,'position');
    pax2 = get(axa2,'position');
    paxc = 1*get(axc,'position');
    paxc2 = get(axc2,'position');
    paxs = 1*get(axs,'position');
    paxs2 = get(axs2,'position');
    
    cop = copyobj([axs axc axa],gcf);
    set(cop,'units','normal')
    set(cop(3),'position',[pax2(1) pax2(2) pax(3) pax(4)])
    set(cop(2),'position',[paxc2(1) paxc2(2) paxc(3) paxc(4)])
    set(cop(1),'position',[paxs2(1) paxs2(2) paxs(3) paxs(4)])
    
    delete([axa2 axc2 axs2])
    
    colormap(mri_scan.split)
    
    pause
    
    
    for time = 1:1:size(mri_scan.GridAmp,2)
        time
        %CInterp = sparse(H*abs(mri_scan.GridAmp(:,time)));
        % Reset MRI to original 
        mri_scan.mri( mri_scan.indmodif) =  mri_scan.gl_orig;
        ind = find(CInterp(:,time)>mri_scan.thres);
        mri_scan.indmodif= ind; % indices of the voxels with modified grey levels
        mri_scan.gl_orig = mri_scan.mri(ind);
        mri_scan.mri(indd(ind)) = round(1.01*mri_scan.climm + nonzeros(CInterp(ind,time)));
        set(findobj(cop(1),'Type','image') ,'CData',(flipud(squeeze(mri_scan.mri(mri_scan.is,:,:))')),'erasemode','none');
        set(findobj(cop(2),'Type','image') ,'CData',(fliplr(flipud(squeeze(mri_scan.mri(:,mri_scan.ic,:))'))),'erasemode','none');
        set(findobj(cop(3),'Type','image') ,'CData',(mri_scan.mri(:,:,mri_scan.ia)),'erasemode','none');    
        drawnow
        pause(.1)
    end
    
    set(MAINWINDOW,'Userdata',mri_scan)
    
    %-----------------------------------------------------------------------------------------------
    
case 'CorticalScout'
    
    %%% ---- Initialize visulization
    MAINWINDOW = gcbf;
    mri_scan = get(MAINWINDOW ,'Userdata');
    
    if isempty(mri_scan) %~exist('mri_scan','var')
        errordlg('Please load an MR file first')
        return
    end
    
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
        
    % Delete all previous points
    delete(findobj(get(axa,'children'),'type','patch'))
    delete(findobj(get(axc,'children'),'type','patch'))
    delete(findobj(get(axs,'children'),'type','patch'))
    
    %%% ---- END: Initialize visulization
    
%     TASKBAR = findobj(0,'Tag','TASKBAR');
%     if isempty(TASKBAR), return, end
    Current = get_user_directory;
    
    cd(Current.STUDIES)
    
    [result_file,result_dir] = uigetfile('CorticalScout.mat','File containing the Cortical Scouts...');
    if result_file == 0, return, end
    cd(result_dir)
    Results = load(result_file);
    Results.SourceLoc = Results.CorticalScouts.CorticalMarkers';
    
    nsource = size(Results.SourceLoc,2);
    
    for src = 1:nsource
        Results.SourceOrder(src) = -1;
        ptname{src} = Results.CorticalScouts.CorticalMarkersLabels{src};
        SourceLoc(:,src) = Results.CorticalScouts.CorticalMarkers(src,:)';
        COLOR(src,:) = [min(.3+abs(rand),1) min(.3+abs(rand),1) min(.3+abs(rand),1) ];
    end
    
    % Is it EEG or MEG
    if ~isfield(Results,'DataFlag')
        Results.DataFlag = 'MEG'; %Assume older data types are MEG
    end
    
    global EEG MEG
    if isempty(strmatch('EEG', Results.DataFlag));
        MEG = 1;
        EEG = 0;
    else
        MEG = 0;
        EEG = 1;
    end
    
 
    nsource = size(SourceLoc,2);
    
    j = 1;
    
    for i = 1:nsource
        
        ptcoord(i,:) = 1000*SourceLoc(:,i)';
        icoord(i,:) = [0 0 0];
        
        if mri_scan.pcs == 0
            orig = findstr(ptname{i},'originPCS');
            if ~isempty(orig) % This set of points can be registered owing to a MEG-EEG PCS.
                mri_scan.pcs = 1; 
                mri_scan.pcs_type = ptname{i}(findstr(ptname{i},'_')+1:end);
                mri_scan.center = icoord;
                center = icoord(i,:);
                ptcenter = ptcoord(i,:)
            else 
                mri_scan.pcs = 0; 
                mri_scan.pcs_type = 'NONE';
                center = [0 0 0];
            end
        end
        
    end
    fclose all;
    
    % Transform PCS coordinates in the MRI volume reference system
    switch mri_scan.pcs_type
    case {'NEUROMAG','CTF'}
        if ~isfield(mri_scan,'rot')
            errordlg('Please Define a PCS First')
            return
        end
        center = mri_scan.center;
        rot = mri_scan.rot;
        ptcoord = (inv(rot)*ptcoord')';      
    end
    
    ptcoord(:,1) = ptcoord(:,1) + center(1);
    ptcoord(:,2) = ptcoord(:,2) + center(2);
    ptcoord(:,3) = ptcoord(:,3) + center(3);
    
    icoord(:,1) = round(ptcoord(:,1)/pix_siz(2));
    icoord(:,2) = round(ptcoord(:,2)/pix_siz(1));
    icoord(:,3) = round(ptcoord(:,3)/pix_siz(3));
    mri_scan.pthandle = [];
    
    for i = 1:size(ptcoord,1)
        x = ptcoord(i,1);
        y = ptcoord(i,2);
        z = ptcoord(i,3);
        
        axes(axs)
        hold on
        switch Results.SourceOrder(i)
        case -1 % Current Dipole
            marker = 'o';
        case 0 % Magnetic Dipole
            marker = 's';
        case 1 % Multipole
            marker = 'd';
        case 2 % Synchronous Current Dipoles
            marker = 'v';
        case 3 % Synchronous Current Dipoles
            marker = '^';
        case 4 % Synchronous Multipoles
            marker = '*'; 
        end
        
        pt_s = scatter(x,z,marker,'filled');
        set(pt_s,'Parent',axs,'userdata',y)
        pt_slabel = text(x,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_slabel,'Parent',axs,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',y)
        
        axes(axa)
        hold on
        pt_a = scatter(x,get(s1,'Max')-y,marker,'filled');
        set(pt_a,'Parent',axa,'userdata',z)
        pt_alabel = text(x,get(s1,'Max')-y,num2str(size(mri_scan.pthandle,1)));
        set(pt_alabel,'Parent',axa,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',z)
        
        axes(axc)
        hold on
        pt_c = scatter(y,z,marker,'filled');
        set(pt_c,'Parent',axc,'userdata',x)
        pt_clabel = text(y,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_clabel,'Parent',axc,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',x)
        
        set([pt_s,pt_a,pt_c],'markerfacecolor','b')
        mri_scan.pthandle = [mri_scan.pthandle;[pt_s,pt_a,pt_c]];
        mri_scan.ptlabels= [mri_scan.ptlabels;[pt_slabel,pt_alabel,pt_clabel]];
        
        set([pt_s,pt_a,pt_c],'markerfacecolor',COLOR(i,:),'MarkerEDGEColor',COLOR(i,:))
    end
    
    seepoints = findobj(MAINWINDOW,'Tag','seepoints');
    seelabels = findobj(MAINWINDOW,'Tag','seelabels');
    LEAFLET = findobj(MAINWINDOW,'Tag','LEAFLET');
    set(seelabels,'enable','on')
    set(seepoints,'enable','on')
    set(LEAFLET,'enable','on')
    
    set(findobj(MAINWINDOW,'Tag','deletepoint'),'enable','on')
    
    listboxx = findobj(0,'Tag','listbox');
    set(listboxx,'visible','on')
    set(listboxx,'Value',1,'String',ptname(:))
    listboxx_title = findobj(MAINWINDOW,'Tag','listboxx_title');
    set(listboxx_title,'visible','on')
    legend_sources = findobj(MAINWINDOW,'Tag','legend_sources');
    set(legend_sources,'visible','on')
    
    
    RENDERING = findobj(0,'Tag','rendering');
    if ~isempty(RENDERING)
        listboxx_rend = findobj(RENDERING,'Tag','listobx_rend');
        set(listboxx_rend,'Value',1,'String',ptname)
    end
    
    mri_scan.ptcurrent = size( mri_scan.ptcoord,1);
    mri_scan.ptname = ptname ;
    mri_scan.ptcoord = ptcoord ;
    
    visu_pt(mri_scan.pthandle)
    
    set(MAINWINDOW,'Userdata',mri_scan)
    
    %-----------------------------------------------------------------------------------------------------------
    
case 'MEGBlobs' % Visualize local cortical current density estimates within the MR volume
    % As those obtained from MNE or non-linear imaging methods
    MAINWINDOW = gcbf;
    set(0,'CurrentFigure',MAINWINDOW )
    mri_scan = get(gcbf,'Userdata');
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    axs = findobj(MAINWINDOW,'Tag','axs');
    axc = findobj(MAINWINDOW,'Tag','axc');
    axa = findobj(MAINWINDOW,'Tag','axa');
    gca = axs;
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    
    %Get Results file
    Current = get_user_directory;
    cd(Current.STUDIES)   
    [resfile, path] = uigetfile('*.mat','Pick-up a Result file');
    cd(path)
    load(resfile)
    
    if isempty(ImageGridAmp)
        errordlg('No Cortical Current Density Images Available for this Study')
        return
    end
    
    
    % Regular way
    % Tessellation file
    cd(Current.SUBJECTS)
    load(StudySubject.Subject)
    if ~exist(Tesselation,'file')
        errordlg('No Cortical Tessellation is Available for that Subject ')
        return
    end
    load(Tesselation,'Vertices','Faces');
    
    % CHEAT
    %load('g:\brainstorm\subjects\huangmx\SUBJECTESS_SMOOTH_SCS_red.MAT','Vertices')
    %load g:\brainstorm\subjects\subjects_old\whitekevin\subjectess_smooth_scs_rap.mat   
    %load G:\Brainstorm\subjects\BailletS\SUBJECtTESS_SMOOTH_SCS_red_CTF.MAT.old -mat
    
    mri_scan.mript = Vertices{GUI.iGrid}';
    clear Vertices;
    Faces = Faces{GUI.iGrid};
    
    
    % See the whole tessellation
    %GridAmp = 10*ones(size(mri_scan.mript)); 
    %ImageGridAmp = GridAmp ;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    mri_scan.GridAmp = norlig(ImageGridAmp)'; % See integrated activity at each vertex along time
    tmax = 1;% See integrated activity at each vertex along time
    mri_scan.M = max(mri_scan.GridAmp );
    
    
    %mri_scan.GridAmp = abs(ImageGridAmp); % see regular mapping
    %[mri_scan.M , tmax] = max(max(abs(ImageGridAmp)));% see regular mapping
    
    
    mri_scan.climm = ceil(max(get(gca,'Clim')));
    
    mri_scan.mri(mri_scan.mri >= mri_scan.climm) = mri_scan.climm;
    mri_scan.GridAmp = mri_scan.GridAmp*mri_scan.climm/mri_scan.M;
    
    mri_scan.thres = .2*mri_scan.climm; % Threshold for visualization of current density
    
    %mri_scan.ivisu = find(mri_scan.GridAmp >= mri_scan.thres); % Regular
    %mri_scan.ivisu = 1:size(mri_scan.mript,1); % CHEAT - see all points
    
    rooot = findstr(StudySubject.Study,'brainstormstudy.mat');
    
    if isempty(rooot)
        errodlg('Study file name should be of the form ''*brainstormstudy.mat''')
        return
    end
    
    rooot = StudySubject.Study(1:rooot-2);
    BrainStorm.rooot = rooot;
    
    mri_scan.mript = 1000*mri_scan.mript; % CHEAT
    %length(mri_scan.ivisu)
    
    
    %--------
    
    % Transform PCS coordinates in the MRI volume reference system
    switch mri_scan.pcs_type
    case {'NEUROMAG','CTF'}
        if ~isfield(mri_scan,'rot')
            errordlg('Please Define a PCS First')
            return
        end
        center = mri_scan.center;
        rot = mri_scan.rot;
        mri_scan.mript = (inv(rot)*mri_scan.mript')';      
    end
    
    mri_scan.mript (:,1) = mri_scan.mript (:,1) + center(1);
    mri_scan.mript (:,2) = mri_scan.mript (:,2) + center(2);
    mri_scan.mript (:,3) = mri_scan.mript (:,3) + center(3);
    
    icoord(:,1) = round((get(s1,'Max')-mri_scan.mript (:,2))/mri_scan.pix_siz(2)+1);
    icoord(:,2) = round(mri_scan.mript (:,1)/mri_scan.pix_siz(1)+1);
    icoord(:,3) = round((get(s2,'Max')-mri_scan.mript (:,3))/mri_scan.pix_siz(3)+1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %ind = sub2ind(size(mri_scan.mri),icoord(:,1),icoord(:,2),icoord(:,3));
    %[uind,ii] = unique(ind);
    
    %%%% Interpolation of the Cortical Current Density Image
    if ~exist('interp_mat.mat','file')
        [H] = tri_interp(size(mri_scan.mri),icoord, Faces)   ;
        disp('Saving Interpolation Matrix')
        User = get_user_directory;
        save interp_mat H
    else
        load interp_mat
    end
    
    CInterp = sparse(H*mri_scan.GridAmp(:,tmax));
    mri_scan.csources = {[1:size(H,2)]};
    clear H mri_scan.GridAmp
    
    ind =find(CInterp>=mri_scan.thres);
    
    mri_scan.indmodif= ind; % indices of the voxels with modified grey levels
    mri_scan.gl_orig = mri_scan.mri(ind);
    
    %  mri_scan.mri(ind) = 1.01*mri_scan.climm + nonzeros(CInterp(ind));
    mri_scan.mri(ind) = double(mri_scan.mri(ind)) + nonzeros(CInterp(ind));
    minn = double(min(mri_scan.mri(ind)));
    mri_scan.mri(ind) =  double(mri_scan.mri(ind)) - minn + mri_scan.climm;
    
    %%%----- Superimpose to regular MRI
    %  mri_scan.mri = mri_scan.mri + iiiblob;
    
    [M tmp] = max(mri_scan.mri(ind));
    [is ic ia] = ind2sub(size(mri_scan.mri),ind(tmp));
    
    ax = mri_scan.ax;
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    set(sag ,'CData',(flipud(squeeze(mri_scan.mri(is,:,:))')),'erasemode','none');
    set(get(sag,'Parent'),'CLim',[0 2*mri_scan.climm])
    set(cor ,'CData',(fliplr(flipud(squeeze(mri_scan.mri(:,ic,:))'))),'erasemode','none');
    set(get(cor,'Parent'),'CLim',[0 2*mri_scan.climm])
    set(ax ,'CData',(mri_scan.mri(:,:,ia)),'erasemode','none');    
    set(get(ax,'Parent'),'CLim',[0  2*mri_scan.climm])
    
    split(1:mri_scan.climm,:) = bone(mri_scan.climm);
    split(mri_scan.climm+1:2*mri_scan.climm,:) = hot(mri_scan.climm);
    mri_scan.split = split;
    
    set(MAINWINDOW,'Userdata',mri_scan);
    
    
    ax = findobj(MAINWINDOW,'Tag','colorbar');
    axes(ax)
    colormapp = split( mri_scan.climm+1: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)
    
    
    %-----------------------------------------------------------------------------------------------------------
    
    
case 'MEGBlob_update' % Update the display of the cortical current density
    
    disp('   ')
    
case 'SaveAsBrainstorm' % Saves the Image cube according to the Brainstorm data structure 
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(gcf,'Userdata');
    if isempty(mri_scan)
        errordlg('Please Load a Scan First','The BrainStorm Toolbox')
        return
    end
    
    Cube = uint8(mri_scan.mri); % Convert to int8 data type 
    
    Voxsize = mri_scan.fov./(size(Cube));
    Voxsize = mri_scan.fov([2 1 3])./(size(Cube)); % To be compliant with the reference system of the MRI volume here.
    %Voxsize(3) = mri_scan.fov([3])./size(Cube,3); % Different for slic thickness
    
    if isfield(mri_scan,'Segment')
        Segment = mri_scan.Segment;
    else
        Segment = [];
    end
    
    if isfield(mri_scan,'pcs') % A PCS has been defined 
        PCS.R = inv(mri_scan.rot); % According to BrainStorm convention
        PCS.t = mri_scan.center; % According to BrainStorm convention
        PCS.Comment = mri_scan.pcs_type;
        PCS.PCSFiducial = mri_scan.PCSFiducial;
        PCS.CubeFiducial = mri_scan.CubeFiducial;
        PCS.FiducialName = mri_scan.FiducialName;
        Comment = mri_scan.pcs_type;
    else 
        PCS = [];
        Comment = '';
    end
    
    % Current = get_user_directory;
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    [file,path] = uiputfile('subjectimage.mat','Save Anatomical Scans in...');
    if file == 0
        return
    end
    cd(path)
    eval(['save ',file,' Cube Voxsize Comment Segment PCS -mat'])
    mri_scan.filename = file;
    set(MAINWINDOW,'Userdata',mri_scan)
    
    
    %-----------------------------------------------------------------------------------------------------------
    
case 'LoadBrainStorm'
    % Pretty much the same as 'READ_HEADER', even more simple though.
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    set(0,'CurrentFigure',MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    axes(findobj(MAINWINDOW,'Tag','axs')); 
    cla 
    axes(findobj(MAINWINDOW,'Tag','axc'));
    cla
    axes(findobj(MAINWINDOW,'Tag','axa'));
    cla
    
    clear mri_scan
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    [filebs, pathbs] = uigetfile('*subjectimage.mat','Select a .mat file containing the subject image');
    if filebs==0
        return
    end
    cd(pathbs); 
    eval(['load ',filebs,' -mat'])
    
    mri_scan.filename = fullfile(pathbs,filebs);
    mri_scan.mri = Cube;
    clear Cube
    mri_scan.pix_size = Voxsize([2 1 3]);
    siz = size(mri_scan.mri);
    siz = siz([2 1 3]);
    
    %mri_scan.pix_size = mri_scan.pix_size./(1-1./siz);
    
    mri_scan.fov = Voxsize([2 1 3]).*(siz); 
    mri_scan.scantype = 'Axial'; 
    mri_scan.filemri = Comment; 
    mri_scan.header = 0; 
    mri_scan.siz = siz(1:2); 
    mri_scan.nslice= siz(3); 
    mri_scan.slice_thick = Voxsize(3); 
    mri_scan.code= 'char'; 
    mri_scan.readslice = 0;
    
    if isfield(PCS,'R')% A PCS has been defined 
        mri_scan.pcs = 1;
        mri_scan.rot = inv(PCS.R);
        mri_scan.center = PCS.t ;
        mri_scan.pcs_type = PCS.Comment ;
        mri_scan.PCSFiducial = PCS.PCSFiducial ;
        mri_scan.CubeFiducial = PCS.CubeFiducial ;
        mri_scan.FiducialName = PCS.FiducialName;
        set(findobj(MAINWINDOW,'Tag','showfiducials'),'Enable','on')
        
        if exist('scale_fiducials','var') % User uses the dummy head
            mri_scan.scale_fiducials = scale_fiducials;
            mri_scan.rot = mri_scan.scale_fiducials*mri_scan.rot ;
        end
        
    else
        mri_scan.pcs = 0;
        set(findobj(MAINWINDOW,'Tag','showfiducials'),'Enable','off')
    end
    
    set(MAINWINDOW,'Userdata',mri_scan) 
    mritool_cb read_file_scan
    
case 'mrislice'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    set(0,'CurrentFigure',MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    delete(findobj(MAINWINDOW,'Type','image'))
    clear mri_scan
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    C = get(findobj(MAINWINDOW,'Tag','CurrentDrive'),'String');
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    switch lower(computer)
    case 'pcswin'
        mri_scan.cdrive = C;
        cd(C);
    end
    
    d = dir;   
    [ignore,sorti] = sort(lower({d.name})); % sort
    str = {d.name}; % unsorted list
    str = str(sorti); % sorted without respect to capitalization
    [s,v] = listdlg('PromptString','Select a file:',...
        'SelectionMode','multiple',...
        'ListString',str);
    
    if v == 0
        return
    end
    
    while length(s)==1
        cd(str{s});
        d = dir;
        [ignore,sorti] = sort(lower({d.name})); % sort
        str = {d.name}; % unsorted list
        str = str(sorti); % sorted without respect to capitalization
        [s,v] = listdlg('PromptString','Select a file:',...
            'SelectionMode','multiple',...
            'ListString',str);
        if v==0
            return
        end
        
    end
    
    mri_scan.s = s;
    mri_scan.str = str;
    nslice = length(s);
    
    header = str2num(get(Header,'String'));
    mri_scan.header = header;
    
    siz = str2num(get(image_size,'String'));
    if length(siz)==1
        siz = [siz,siz];
    end
    mri_scan.siz = siz;
    
    scan_type_tmp = get(scan_type,'String');
    tmp = get(scan_type,'Value');
    scan_type = strtok(char(scan_type_tmp(tmp)));
    mri_scan.scantype= scan_type;    
    
    fov = str2num(get(FOV,'String'));
    if isempty(fov)
        fov = siz;
    end
    mri_scan.fov = fov;
    
    slice_thickness = str2num(get(SLICE_thickness ,'String'));
    if isempty(slice_thickness)
        slice_thickness  = DEFAULTslice_thickness;
        warndlg(['Slice thickness field empty - using default',num2str(DEFAULTslice_thickeness)])   
    end
    mri_scan.slice_thick = slice_thickness;
    
    set(NSLICE,'String',nslice);
    if isempty(nslice)
        errordlg('Please enter the number of slices')
        return
    end
    
    siz_tmp = siz;
    
    clear siz_tmp scan_type_tmp
    
    mri_scan.nslice = nslice;
    
    pix_siz = [fov(1)/(siz(1)-1) fov(2)/(siz(2)-1) slice_thickness];
    mri_scan.pix_siz = pix_siz;
    
    code = (get(data_format,'String'));   
    tmp = get(data_format,'Value');
    code = strtok(char(code(tmp)));
    mri_scan.code = code;
    
    
    mri = zeros(siz(1),siz(2),nslice);
    hh = waitbar(0,'Reading files...');
    %   str = str(3:end);
    str = str(s);
    for i=1:nslice
        file = fopen(str{i},'rb',MACHINEFORMAT);
        fseek(file,header,'bof');
        if feof(file)
            errordlg('Please Check Image size / Number of Slices : do not match file size')
            return
        end
        mri(:,:,i) = fread(file,[siz(1),siz(2)],code);
        fclose(file);
        waitbar(i/nslice);
    end
    close(hh);
    
    mri_scan.mri = mri;
    clear mri
    
    ngray = 256;
    M = max(mri_scan.mri(:)); % May be rather time consuming if limited RAM is available
    
    if M > (ngray-1)    
        mri_scan.mri = double(mri_scan.mri) ;
        M = double(M);
        mri_scan.mri = (mri_scan.mri) * ((ngray - 1 )/M) ;
        mri_scan.mri = uint8(mri_scan.mri);
    end
    
    mri_scan.readslice = 1;
    mri_scan.filemri = str{i};  
    
    set(gcf,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    
    %-----------------------------------------------------------------------------------------------------------
case 'ctfmri' % Test for the proper reading of CTF .mri file formats 
    
    % CTF MRI File Version 2.1
    % Pretty much the same as 'READ_HEADER', even simpler though.
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    set(0,'CurrentFigure',MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    axes(findobj(MAINWINDOW,'Tag','axs')); 
    cla 
    axes(findobj(MAINWINDOW,'Tag','axc'));
    cla
    axes(findobj(MAINWINDOW,'Tag','axa'));
    cla
    
    clear mri_scan
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    [filebs, pathbs] = uigetfile('*.mri','Please select a CTF .mri file');
    if filebs==0
        return
    end
    cd(pathbs); 
    
    file = fopen(filebs,'rb',MACHINEFORMAT); 
    if file < 1
        errordlg(['Error when opening : ',fullfile(pathbs, filebs)]);
        return
    end
    
    Header.identifierString = char(fread(file,32,'char')');
    Header.imageSize = fread(file,1,'short');
    Header.dataSize = fread(file,1,'short');
    Header.clippingRange = fread(file,1,'short');
    Header.imageOrientation= fread(file,1,'short');
    Header.mmPerPixel_sagittal = fread(file,1,'float32');
    Header.mmPerPixel_coronal = fread(file,1,'float32');
    Header.mmPerPixel_axial = fread(file,1,'float32');
    
    % Head Model Info
    Header.HeadModel_Info.Nasion_Sag = fread(file,1,'short');
    Header.HeadModel_Info.Nasion_Cor = fread(file,1,'short');
    Header.HeadModel_Info.Nasion_Axi = fread(file,1,'short');
    Header.HeadModel_Info.LeftEar_Sag = fread(file,1,'short');
    Header.HeadModel_Info.LeftEar_Cor = fread(file,1,'short');
    Header.HeadModel_Info.LeftEar_Axi = fread(file,1,'short');
    Header.HeadModel_Info.RightEar_Sag = fread(file,1,'short');
    Header.HeadModel_Info.RightEar_Cor = fread(file,1,'short');
    Header.HeadModel_Info.RightEar_Axi = fread(file,1,'short');
    Header.HeadModel_Info.defaultSphereX= fread(file,1,'float32');
    Header.HeadModel_Info.defaultSphereY= fread(file,1,'float32');
    Header.HeadModel_Info.defaultSphereZ= fread(file,1,'float32');
    Header.HeadModel_Info.defaultSphereRadius= fread(file,1,'float32');
    % Head Model Info ********************************************
    
    % Image Info
    Header.Image_Info.modality = fread(file,1,'short');
    Header.Image_Info.manufacturerName = fread(file,64,'char');
    Header.Image_Info.instituteName = fread(file,64,'char');
    Header.Image_Info.patientID = fread(file,32,'char');
    Header.Image_Info.dateAndTime = fread(file,32,'char');
    Header.Image_Info.scanType = fread(file,32,'char');
    Header.Image_Info.contrastAgent = fread(file,32,'char');
    Header.Image_Info.imagedNucleus = fread(file,32,'char');
    Header.Image_Info.Frequency = fread(file,1,'float32');
    Header.Image_Info.FieldStrengh = fread(file,1,'float32');
    Header.Image_Info.EchoTime = fread(file,1,'float32');
    Header.Image_Info.RepetitionTime = fread(file,1,'float32');
    Header.Image_Info.InversionTime = fread(file,1,'float32');
    Header.Image_Info.FlipAngle = fread(file,1,'float32');
    Header.Image_Info.NoExcitations = fread(file,1,'short');
    Header.Image_Info.NoAcquisitions = fread(file,1,'short');
    Header.Image_Info.commentString = fread(file,256,'char');
    Header.Image_Info.forFutureUse = fread(file,64,'char');
    % Image Info ********************************************
    
    
    Header.headOrigin_sagittal = fread(file,1,'float32');
    Header.headOrigin_coronal = fread(file,1,'float32');
    Header.headOrigin_axial  = fread(file,1,'float32');
    
    Header.rotate_coronal = fread(file,1,'float32');
    Header.headOrigin_sagittal = fread(file,1,'float32');
    Header.headOrigin_axial  = fread(file,1,'float32');
    Header.orthogonalFlag = fread(file,1,'short');
    Header.unused = fread(file,272,'uchar');
    
    fseek(file,1028,-1); % Skip the whole header by padding to 1028
    
    % Header Info  ********************************************
    
    %Read Image
    mri_scan.nslice = Header.imageSize;
    mri_scan.siz =  [256 256 mri_scan.nslice];
    
    %mri_scan.mri = zeros(mri_scan.siz(1),mri_scan.siz(2),mri_scan.nslice);
    
    
    switch(Header.dataSize)
    case 1
        mri_scan.code = 'uchar';
    case 2
        mri_scan.code = 'ushort';
    otherwise
        errordlg('Unrecognized data format - Please check the file format');
        return
    end
    
    %    hh = waitbar(0,['Reading ', char(Header.identifierString),' ...']);
    %    for i=1:mri_scan.nslice
    %       mri_scan.mri(:,:,i) = uint8(fread(file,[mri_scan.siz(1),mri_scan.siz(2)],mri_scan.code));
    %       waitbar(i/mri_scan.nslice);
    %    end
    
    hmsg = msgbox('Reading file...');
    set(MAINWINDOW,'Pointer','watch')
    drawnow
    
    mri_scan.mri = int16(fread(file, prod(mri_scan.siz),mri_scan.code)); % Read the whole volume / store in a huge vector   
    fclose(file);
    mri_scan.mri = reshape(mri_scan.mri,mri_scan.siz(1),mri_scan.siz(2),mri_scan.siz(3)); % reshape it into a 3-D array.
    mri_scan.clippingRange = Header.clippingRange ;
    
    delete(hmsg)
    set(MAINWINDOW,'Pointer','arrow')
    drawnow
    
    
    %    ngray = 256;
    %    M = max(mri_scan.mri(:)); % May be rather time consuming if limited RAM is available
    %    
    %    if M > (ngray-1)    
    %       mri_scan.mri = double(mri_scan.mri) ;
    %       M = double(M);
    %       mri_scan.mri = (mri_scan.mri) * ((ngray - 1 )/M) ;
    % 	   mri_scan.mri = uint8(mri_scan.mri);
    %    end
    %    
    mri_scan.scantype = 'Sagittal';
    mri_scan.fov = [mri_scan.siz].*[Header.mmPerPixel_sagittal Header.mmPerPixel_coronal Header.mmPerPixel_axial];
    mri_scan.slice_thick = Header.mmPerPixel_sagittal;
    
    mri_scan.readslice = 1;
    mri_scan.filemri = fullfile(pathbs, filebs);  
    mri_scan.header = 0;
    
    set(MAINWINDOW,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    mritool_cb flip_ud_s
    mritool_cb flip_lr_a
    
    %-----------------------------------------------------------------------------------------------------------
case 'fifmri' % Load FIF MRI format (4D NeuroImaging) 
    
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    set(0,'CurrentFigure',MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    axes(findobj(MAINWINDOW,'Tag','axs')); 
    cla 
    axes(findobj(MAINWINDOW,'Tag','axc'));
    cla
    axes(findobj(MAINWINDOW,'Tag','axa'));
    cla
    
    clear mri_scan
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    
    if ~exist('loadmri.m','file') % meg_1.2 routine not available on this machine
        errordlg('meg_1.2 routines not available on this machine. Please check the following link: http://boojum.hut.fi/~kuutela/meg-pd/');
        return
    end
    [filebs, pathbs] = uigetfile('*.fif','Please select a FIF MRI file');
    if filebs==0
        return
    end
    cd(pathbs); 
    disp('Loading MRI volume...' )
    drawnow
    [M,coords]=loadmri(filebs);
    disp('Loading MRI volume...-> DONE' )
    drawnow
    
    %Check if fiducials are defined
    try 
        T = loadtrans(filebs,'HEAD','MRI');
    catch
        msgdlg('Fiducials were note defined for this MRI - Please proceed manually within BrainStorm');
    end
    
    mri_scan.PCS.T = T;
    
    if isempty(M) 
        errordlg(['Error when opening : ',fullfile(pathbs, filebs)]);
        return
    end
    
    % Header Info  ********************************************
    
    %Read Image
    mri_scan.nslice = size(M,3);
    mri_scan.siz =  [size(M,1) size(M,2) mri_scan.nslice];
    
    kk = whos('M');
    mri_scan.code = kk.class;
    
    %mri_scan.mri = reshape(mri_scan.mri,mri_scan.siz(1),mri_scan.siz(2),mri_scan.siz(3)); % reshape it into a 3-D array.
    mri_scan.mri = M; clear M
    
    %mri_scan.clippingRange = Header.clippingRange ;
    
    
    mri_scan.scantype = 'Sagittal';
    mri_scan.fov = 1000*abs([mri_scan.siz].*[coords(find(coords(:,1)),1), coords(find(coords(:,2)),2),...
            coords(find(coords(:,3)),3)]);
    
    mri_scan.slice_thick = 1000*abs(coords(find(coords(:,3)),3));
    
    mri_scan.readslice = 1;
    mri_scan.filemri = fullfile(pathbs, filebs);  
    mri_scan.header = 0;
    mri_scan.PCS.t = [];;
    
    set(MAINWINDOW,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    %-----------------------------------------------------------------------------------------------------------
case 'imgformat' % Test for the proper reading of CTF (?) .img file formats - headers of these files seem to be of variable size from on file to the other
    
    
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    set(0,'CurrentFigure',MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    delete(findobj(MAINWINDOW,'Type','image'))
    clear mri_scan
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    C = get(findobj(MAINWINDOW,'Tag','CurrentDrive'),'String');
    
    mri_scan.cdrive = C;
    cd(C);
    
    d = dir;
    str = {d.name};
    [s,v] = listdlg('PromptString','Select a file:',...
        'SelectionMode','multiple',...
        'ListString',str);
    
    if v == 0
        return
    end
    
    while length(s)==1
        cd(str{s});
        d = dir;
        str = {d.name};
        [s,v] = listdlg('PromptString','Select a file:',...
            'SelectionMode','multiple',...
            'ListString',str);
        if v==0
            return
        end
        
    end
    
    mri_scan.s = s;
    mri_scan.str = str;
    nslice = length(s);
    
    header = str2num(get(Header,'String'));
    mri_scan.header = header;
    
    siz = str2num(get(image_size,'String'));
    if length(siz)==1
        siz = [siz,siz];
    end
    mri_scan.siz = siz;
    
    scan_type_tmp = get(scan_type,'String');
    tmp = get(scan_type,'Value');
    scan_type = strtok(char(scan_type_tmp(tmp)));
    mri_scan.scantype= scan_type;    
    
    fov = str2num(get(FOV,'String'));
    if isempty(fov)
        fov = siz;
    end
    mri_scan.fov = fov;
    
    slice_thickness = str2num(get(SLICE_thickness ,'String'));
    if isempty(slice_thickness)
        slice_thickness  = DEFAULTslice_thickness;
        warndlg(['Slice thickness field empty - using default',num2str(DEFAULTslice_thickeness)])   
    end
    mri_scan.slice_thick = slice_thickness;
    
    set(NSLICE,'String',nslice);
    if isempty(nslice)
        errordlg('Please enter the number of slices')
        return
    end
    
    siz_tmp = siz;
    
    clear siz_tmp scan_type_tmp
    
    mri_scan.nslice = nslice;
    
    pix_siz = [fov(1)/(siz(1)-1) fov(2)/(siz(2)-1) slice_thickness];
    mri_scan.pix_siz = pix_siz;
    
    code = (get(data_format,'String'));   
    tmp = get(data_format,'Value');
    code = strtok(char(code(tmp)));
    mri_scan.code = code;
    
    mri = zeros(siz(1),siz(2),nslice);
    hh = waitbar(0,'Reading files...');
    %   str = str(3:end);
    str = str(s);
    
    for i=1:nslice
        file = fopen(str{i},'rb',MACHINEFORMAT);
        % Assume here that code is 'int16'
        status = fseek(file,-siz(1)*siz(2)*2,'eof');
        if status == -1 
            errordlg('Please Check Image size / Number of Slices : do not match file size')
            return
        end
        mri(:,:,i) = fread(file,[siz(1),siz(2)],code);
        fclose(file);
        waitbar(i/nslice);
    end
    
    close(hh);
    
    mri_scan.mri = mri;
    clear mri
    
    ngray = 256;
    M = max(mri_scan.mri(:)); % May be rather time consuming if limited RAM is available
    
    if M > (ngray-1)    
        mri_scan.mri = double(mri_scan.mri) ;
        M = double(M);
        mri_scan.mri = (mri_scan.mri) * ((ngray - 1 )/M) ;
        mri_scan.mri = uint8(mri_scan.mri);
    end
    
    mri_scan.readslice = 1;
    mri_scan.filemri = str{i};  
    
    set(gcf,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    
    
    
    
case 'IMPORT_mridata'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    set(0,'CurrentFigure',   MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    delete(findobj(MAINWINDOW,'Type','image'))
    clear mri_scan
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    [filemri, pathmri] = uigetfile('*.*','Select a MRI image file for visualization');
    if filemri==0
        return
    end
    
    filename = findobj(gcbf,'Tag','mr_filename');
    
    set(filename,'String',filemri)
    
    
    filemri = [pathmri,filemri];
    cd(pathmri); 
    
    mri_scan.filemri = filemri;
    mri_scan.pathmri = pathmri;
    
    header = str2num(get(Header,'String'));
    mri_scan.header = header;
    
    siz = str2num(get(image_size,'String'));
    if length(siz)==1
        siz = [siz,siz];
    end
    mri_scan.siz = siz;
    
    
    scan_type_tmp = get(scan_type,'String');
    tmp = get(scan_type,'Value');
    scan_type = strtok(char(scan_type_tmp(tmp)));
    mri_scan.scantype= scan_type;    
    
    fov = str2num(get(FOV,'String'));
    if isempty(fov)
        fov = siz;
    end
    mri_scan.fov = fov;
    
    slice_thickness = str2num(get(SLICE_thickness ,'String'));
    if isempty(slice_thickness)
        slice_thickness  = DEFAULTslice_thickness;
        warndlg(['Slice thickness field empty - using default',num2str(DEFAULTslice_thickeness)])   
    end
    mri_scan.slice_thick = slice_thickness;
    
    
    nslice = str2num(get(NSLICE,'String'));
    if isempty(nslice)
        errordlg('Please enter the numer of slices')
        return
    end
    
    siz_tmp = siz;
    
    mri_scan.nslice = nslice;
    code = (get(data_format,'String'));   
    tmp = get(data_format,'Value');
    code = strtok(char(code(tmp)));
    mri_scan.code = code;
    
    mri_scan.readslice = 0;   
    set(gcf,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    
case 'mriformat'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    set(0,'CurrentFigure',   MAINWINDOW )
    set(MAINWINDOW,'Userdata',[]);
    set(findobj(MAINWINDOW,'Type','image'),'erasemode','normal')
    delete(findobj(MAINWINDOW,'Type','image'))
    clear mri_scan
    
    Header = findobj(gcbf,'Tag','header');
    image_size = findobj(gcbf,'Tag','image_size');
    SLICE_thickness = findobj(gcbf,'Tag','slice_thickness');
    FOV = findobj(gcbf,'Tag','fov');
    NSLICE = findobj(gcbf,'Tag','nslice');
    data_format = findobj(gcbf,'Tag','data_format');  
    file_mri= findobj(gcbf,'Tag','file_mri');  
    scan_type= findobj(gcbf,'Tag','scan_type');  
    
    be = findobj(MAINWINDOW,'Tag','big-endian');
    le = findobj(MAINWINDOW,'Tag','little-endian');
    
    if get(le,'value');
        MACHINEFORMAT = 'l';
    else
        MACHINEFORMAT = 'b';
    end
    
    Current = get_user_directory;
    cd(Current.SUBJECTS)
    
    [filemri, pathmri] = uigetfile('*.mri','Select a .mri image file for visualization');
    if filemri==0
        return
    end
    
    filemri = [pathmri,filemri];
    cd(pathmri); 
    
    mri_scan.filemri = filemri;
    mri_scan.pathmri = pathmri;
    
    
    file_mri = fopen(filemri,'r',MACHINEFORMAT);
    if file_mri < 1
        disp('Invalid File Name')
        return;
    end
    
    mri_scan.header = 512;
    
    for i=1:22
        fgetl(file_mri);
    end
    
    siz = fgetl(file_mri);
    [tmp,siz] = strtok(siz,'=');
    siz = str2num(siz(2:end));
    
    nslice = siz(3);
    siz = siz(1:2);
    if length(siz)==1
        siz = [siz,siz];
    end
    
    scan_type_tmp = get(scan_type,'String');
    tmp = get(scan_type,'Value');
    scan_type = strtok(char(scan_type_tmp(tmp)));
    mri_scan.scantype= scan_type;    
    
    
    cmpix = fgetl(file_mri);
    [tmp,cmpix] = strtok(cmpix,'=');
    cmpix = str2num(cmpix(2:end));
    
    fov = siz.*cmpix(1:2);
    if isempty(fov)
        fov = siz;
    end
    mri_scan.fov = fov;
    
    slice_thickness = cmpix(3);
    if isempty(slice_thickness)
        slice_thickness  = DEFAULTslice_thickness;
        warndlg(['Slice thickness field empty - using default',num2str(DEFAULTslice_thickeness)])   
    end
    mri_scan.slice_thick = slice_thickness;
    
    
    if isempty(nslice)
        errordlg('Please enter the number of slices')
        return
    end
    
    siz_tmp = siz;
    if strcmp(scan_type,'Coronal')
        siz(2) = nslice; 
        nslice = siz_tmp(2);
        
    elseif strcmp(scan_type,'Sagittal')
        siz(1) = nslice;
        siz(2) = siz_tmp(1);
        nslice = siz_tmp(2);  
    elseif strcmp(scan_type,'Axial')
        
    end
    
    clear siz_tmp scan_type_tmp
    
    
    mri_scan.nslice = nslice;
    mri_scan.siz = siz;
    
    pix_siz = [fov(1)/(siz(1)-1) fov(2)/(siz(2)-1) slice_thickness];
    mri_scan.pix_siz = pix_siz;
    
    code = 'uint8';
    set(data_format,'value',8)   
    
    mri_scan.code = code;
    
    mri_scan.readslice = 0;
    
    set(gcf,'Userdata',mri_scan)   
    
    mritool_cb read_file_scan
    
    
case 'mouse_slct'
    
    mri_scan = get(gcf,'Userdata');
    
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    
    MOUSE_SLCT = findobj(gco,'tag','mouse_slct');
    val_MOUSE_SLCT = get(MOUSE_SLCT,'value');
    
    ZOOM = findobj(gcbf,'Tag','zoom');
    set(ZOOM,'Value',0);
    set(gcbf,'Pointer','Arrow')
    zoom off
    
    gotopoint = findobj(gcbf,'Tag','gotopoint');
    set(gotopoint,'Value',0)
    
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    nslice = mri_scan.nslice;
    siz = mri_scan.siz;
    slice_thick = mri_scan.slice_thick ;
    fov = mri_scan.fov;
    
    pix_siz = mri_scan.pix_siz ;
    
    while val_MOUSE_SLCT
        [x,y,button] = ginput2(1);
        clear x y;
        current_ax = gco;
        switch current_ax
        case sag
            imax = size(mri_scan.mri,1);
            ii = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2));
            is = ii;
            axi = sag;
            str = 'ii,:,:';
            sldr = s1;   
        case ax
            imax = size(mri_scan.mri,3);
            ii = round(get(s2,'Value')/pix_siz(3));
            ia = ii;
            axi = ax;
            str = ':,:,ii';
            sldr = s2;
        case cor
            imax = size(mri_scan.mri,2);
            ii = round(get(s3,'Value')/pix_siz(1));
            ic = ii;
            axi = cor;
            str = ':,ii,:';
            sldr = s3;
            
        case MOUSE_SLCT
            set(MOUSE_SLCT,'value',0)   
            return
            
        otherwise
            set(MOUSE_SLCT,'value',0)   
            val_MOUSE_SLCT = get(MOUSE_SLCT,'value');   
        end
        
        if button== 1
            ii =ii+1;
        else
            ii = ii-1;
        end
        
        show_coordinates = findobj(gcbf,'Tag','show_coordinates');
        if get(show_coordinates,'Value')
            SHOW_COORDINATES = 1;
        else
            SHOW_COORDINATES = 0;
        end
        
        switch current_ax
        case {sag,cor} 
            if current_ax == sag
                imax = size(mri_scan.mri,1);
                is = ii;
                set(s1,'Value',get(s1,'Max')-is*pix_siz(2))   
                eval(['mri_tmp = (flipud(squeeze(mri_scan.mri(',str,'))''));']);
            elseif current_ax == cor
                imax = size(mri_scan.mri,2);
                eval(['mri_tmp = (fliplr(flipud(squeeze(mri_scan.mri(',str,'))'')));']);
                ic = ii;
                set(s3,'Value',ic*pix_siz(1))   
                
            end
            
        case ax
            %         eval(['mri_tmp = flipud(squeeze(mri(',str,')));']);
            eval(['mri_tmp = (squeeze(mri_scan.mri(',str,')));']);
            imax = size(mri_scan.mri,3);
            ia = ii;
            set(s2,'Value',ia*pix_siz(3))   
            
        end
        
        if (ii>0)&(ii<=imax) 
            set(axi,'Cdata',mri_tmp);
            drawnow
            %  set(sldr,'Value',ii)   
            switch current_ax
            case sag
                if SHOW_COORDINATES 
                    set(axc,'Xtick',get(s1,'Value'),'Xticklabel',num2str(get(s1,'Value'),FORMATtick))
                    set(axa,'Ytick',get(s1,'Max')-get(s1,'Value'));%,'Yticklabel',num2str(get(s1,'Value'),FORMATtick))
                    %set(axs,'Visible','on','Xtick',ic,'Ytick',ia)
                    
                end
            case cor
                if SHOW_COORDINATES 
                    set(axs,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
                    set(axa,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
                    %set(axa,'Visible','on','Xtick',is,'Ytick',ia)
                end
                
            case ax
                if SHOW_COORDINATES 
                    set(axc,'Ytick',get(s2,'Max')-get(s2,'Value'),'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
                    set(axs,'Ytick',get(s2,'Max')-get(s2,'Value'));%'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
                    % set(ax,'Visible','on','Xtick',ic,'Ytick',is)
                end
            end
        elseif is == 0
            is = 1;
            
        end
        
        val_MOUSE_SLCT = get(MOUSE_SLCT,'value');
        
    end
    
    %-----------------------------------------------------------------------------------------------
case 'sliders'
    
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    set(0,'CurrentFigure',   MAINWINDOW )
    mri_scan = get(gcf,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    nslice = mri_scan.nslice;
    siz = mri_scan.siz;
    ic = mri_scan.ic ;
    is = mri_scan.is ;
    ia = mri_scan.ia ;
    pix_siz = mri_scan.pix_siz ;
    
    
    gtpt = mri_scan.gotopoint;
    
    if gtpt 
        current_sldr = 'all';
    else
        current_sldr = gco;
    end
    
    
    s1 = findobj(   MAINWINDOW ,'Tag','1');
    s2 = findobj(   MAINWINDOW ,'Tag','2');
    s3 = findobj(   MAINWINDOW ,'Tag','3');
    
    visu_pt(mri_scan.pthandle)
    
    show_coordinates = findobj(   MAINWINDOW ,'Tag','show_coordinates');
    if get(show_coordinates,'Value')
        SHOW_COORDINATES = 1;
    else
        SHOW_COORDINATES = 0;
    end
    
    fov = mri_scan.fov;
    
    switch current_sldr 
        
    case s1
        imax = size(mri_scan.mri,1);
        ii = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2)+1);
        is = ii;
        is = max([1,min([size(mri_scan.mri,1),is])]);
        axi = sag;
        
        str = 'is,:,:';
        gca = axs;
        if SHOW_COORDINATES 
            set(axc,'Xtick',get(s1,'Value'),'Xticklabel',num2str(get(s1,'Value'),FORMATtick))
            set(axa,'Ytick',get(s1,'Max')-get(s1,'Value'));%,'Yticklabel',num2str(get(s1,'Value'),FORMATtick))
        end
        
        
    case s2
        imax = size(mri_scan.mri,3);
        ii = round(get(s2,'Value')/pix_siz(3)+1);
        ia = ii;
        ia = max([1,min([size(mri_scan.mri,3),ia])]);
        axi = ax;
        imax = size(mri_scan.mri,3);
        str = ':,:,ia';
        gca = axa;
        if SHOW_COORDINATES 
            set(axc,'Ytick',get(s2,'Max')-get(s2,'Value'),'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
            set(axs,'Ytick',get(s2,'Max')-get(s2,'Value'));%'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
        end
        
    case s3
        imax = size(mri_scan.mri,2);
        ii = round(get(s3,'Value')/pix_siz(1)+1);
        ic = ii;
        ic = max([1,min([size(mri_scan.mri,2),ic])]);
        axi = cor;
        imax = size(mri_scan.mri,2);
        str = ':,ic,:';
        gca = axc;
        if SHOW_COORDINATES 
            set(axs,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
            set(axa,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
        end
        
    case 'all'
        is = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2)+1);
        if SHOW_COORDINATES 
            set(axc,'Xtick',get(s1,'Value'),'Xticklabel',num2str(get(s1,'Value'),FORMATtick))
            set(axa,'Ytick',get(s1,'Max')-get(s1,'Value'));%,'Yticklabel',num2str(get(s1,'Value'),FORMATtick))
        end
        
        
        ia = round(get(s2,'Value')/pix_siz(3)+1);
        ia = max([1,min([size(mri_scan.mri,3),ia])]);
        if SHOW_COORDINATES 
            set(axc,'Ytick',get(s2,'Max')-get(s2,'Value'),'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
            set(axs,'Ytick',get(s2,'Max')-get(s2,'Value'));%'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
        end
        
        ic = round(get(s3,'Value')/pix_siz(1)+1);
        ic = max([1,min([size(mri_scan.mri,2),ic])]);
        if SHOW_COORDINATES 
            set(axs,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
            set(axa,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
        end
        ii = 1;
        imax = 2;
    end
    
    
    
    mri_scan.ic = ic;
    mri_scan.is = is;
    mri_scan.ia = ia;
    
    
    switch current_sldr
    case s1 
        eval(['mri_tmp = (flipud(squeeze(mri_scan.mri(',str,'))''));']);
    case s3 
        eval(['mri_tmp = (fliplr(flipud(squeeze(mri_scan.mri(',str,'))'')));']);
        
    case s2
        eval(['mri_tmp = (squeeze(mri_scan.mri(',str,')));']);
        %     eval(['mri_tmp = flipud(squeeze(mri(',str,')));']);
        
    end
    
    
    if (ii>0)&(ii<=imax) 
        switch current_sldr
            
        case 'all'
            strs = 'is,:,:';
            stra = ':,:,ia';
            strc = ':,ic,:';
            eval(['mri_tmp = (flipud(squeeze(mri_scan.mri(',strs,'))''));']);
            set(sag,'Cdata',(mri_tmp));
            eval(['mri_tmp = fliplr(flipud(squeeze(mri_scan.mri(',strc,'))''));']);
            set(cor,'Cdata',(mri_tmp));
            eval(['mri_tmp = squeeze(mri_scan.mri(',stra,'));']);
            %    eval(['mri_tmp = flipud(squeeze(mri(',str,')));']);
            set(ax,'Cdata',(mri_tmp));
            
            RENDERING = findobj(0,'Tag','rendering');
            if ~isempty(RENDERING)
                listboxx_rend = findobj(RENDERING,'Tag','listobx_rend');
                set(listboxx_rend,'Userdata',[is,ia,ic])
                
            end
            
            
        otherwise
            
            set(axi,'Cdata',(mri_tmp));
            drawnow
            gca = get(axi,'Parent');
            
        end
        
    elseif is == 0    
        is = 1;
        
    end
    
    %-----------------------------------------------------------------------------------------------
    
case 'show_coordinates'
    show_coordinates = findobj(gcbf,'Tag','show_coordinates');
    mri_scan = get(gcf,'Userdata');
    axa = mri_scan.axa;
    axs = mri_scan.axs;
    axc = mri_scan.axc;
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    ic = (get(s3,'Value'));
    is = (get(s1,'Value'));
    ia = (get(s2,'Value'));
    pix_siz = mri_scan.pix_siz ;
    
    if get(show_coordinates,'Value')
        set(axa,'Visible','on','Layer','Top','Xtick',ic,'Xticklabel',num2str(ic,FORMATtick),'Ytick',get(s1,'Max')-get(s1,'Value'),'Yticklabel','')
        set(axc,'Visible','on','Layer','Top','Xtick',is,'Xticklabel',num2str(is,FORMATtick),'Ytick',get(s2,'Max')-ia,'Yticklabel',num2str(ia,FORMATtick))
        set(axs,'Visible','on','Layer','Top','Xtick',ic,'Ytick',get(s2,'Max')-ia,'Xticklabel',num2str(ic,FORMATtick),'Yticklabel','')
    else
        set(axa,'Visible','on','Xtick',[],'Ytick',[])
        set(axc,'Visible','on','Xtick',[],'Ytick',[])
        set(axs,'Visible','on','Xtick',[],'Ytick',[])
    end
    
    
    
case 'zoom'
    mri_scan = get(gcf,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    ZOOM = findobj(gcbf,'Tag','zoom');
    valZOOM = get(ZOOM,'Value');
    MOUSE_SLCT = findobj(gcbf,'tag','mouse_slct');
    gotopoint = findobj(gcbf,'tag','gotopoint');
    
    val_MOUSE_SLCT = get(MOUSE_SLCT,'Value');
    val_gotopoint = get(gotopoint,'Value');
    
    %mutincomp([ZOOM MOUSE_SLCT])
    
    if valZOOM
        set(MOUSE_SLCT,'Value',0)
        set(gotopoint ,'Value',0)
        zoom on
        set(gcbf,'Pointer','Circle')
    else
        zoom off
        set(gcbf,'Pointer','arrow')
    end
    
    %-----------------------------------------------------------------------------------------------
case 'slctpoint'  
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    show_coordinates = findobj(MAINWINDOW ,'Tag','show_coordinates');
    set(show_coordinates ,'Value',1)
    mritool_cb show_coordinates
    drawnow
    
    listboxx = findobj(0,'Tag','listbox');
    set(listboxx,'visible','on')
    listboxx_title = findobj(0,'Tag','listboxx_title');
    set(listboxx_title,'visible','on')
    
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    nslice = mri_scan.nslice;
    siz = mri_scan.siz;
    ic = mri_scan.ic ;
    is = mri_scan.is ;
    ia = mri_scan.ia ;
    pix_siz = mri_scan.pix_siz ;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    
    x = get(s3,'Value');
    y = get(s1,'Value');
    z = get(s2,'Max')-get(s2,'Value');
    
    
    set(axs,'Selected','on')
    hold on
    pt_s = scatter(x,z,'o','filled');
    set(pt_s,'Parent',axs,'Userdata',y)
    pt_slabel = text(x,z,num2str(size(mri_scan.pthandle,1)));
    set(pt_slabel,'Parent',axs,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
        'Userdata',y)
    
    
    pt_a = scatter(x,get(s1,'Max')-y,'o','filled');
    set(pt_a,'Parent',axa,'Userdata',z)
    pt_alabel = text(x,get(s1,'Max')-y,num2str(size(mri_scan.pthandle,1)));
    set(pt_alabel,'Parent',axa,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
        'Userdata',z)
    
    pt_c = scatter(y,z,'o','filled');
    set(pt_c,'Parent',axc,'Userdata',x)
    pt_clabel = text(y,z,num2str(size(mri_scan.pthandle,1)));
    set(pt_clabel,'Parent',axc,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
        'Userdata',x)
    
    set([pt_s,pt_a,pt_c],'markerfacecolor','c')
    
    prompt  = {'Enter point name:'}; 
    title   = 'Please assign a label to the selected point';
    lines= 1;
    def     = {'LeftEar'};
    answer  = inputdlg(prompt,title,lines,def);
    
    if isempty(answer)
        delete([pt_s pt_c pt_a pt_slabel pt_clabel pt_alabel   ])
        clear pt_s pt_c pt_a pt_slabel pt_clabel pt_alabel
        return
    end
    hold on
    
    if ~isempty(mri_scan.ptname)
        mri_scan.ptname(end+1) = answer;
    else
        mri_scan.ptname(1) = answer;
    end
    
    set(pt_s,'tag',char(answer));   set(pt_a,'tag',char(answer));
    set(pt_c,'tag',char(answer));
    
    mri_scan.seepoints = 0;
    seepoints = findobj(gcbf,'Tag','seepoints');
    set(seepoints,'String', 'Hide Points');
    
    mri_scan.ptcoord = [mri_scan.ptcoord;[x,y,z]];
    mri_scan.pthandle = [mri_scan.pthandle;[pt_s,pt_a,pt_c]];
    mri_scan.ptcurrent = size( mri_scan.ptcoord,1);
    mri_scan.ptlabels = [mri_scan.ptlabels;[pt_slabel,pt_alabel,pt_clabel]];
    
    set(gcbf,'Userdata',mri_scan)
    
    listboxx = findobj(listboxx,'Tag','listbox');
    set(listboxx,'Value',1,'String',mri_scan.ptname)
    
    seepoints = findobj(MAINWINDOW ,'Tag','seepoints');
    seelabels = findobj(MAINWINDOW ,'Tag','seelabels');
    newlabel = findobj(gcbf,'Tag','newlabel');
    LEAFLET = findobj(MAINWINDOW ,'Tag','LEAFLET');
    set(seelabels,'enable','on')
    set(newlabel,'enable','on')
    set(seepoints,'enable','on')
    set(LEAFLET,'enable','on')
    
    set(findobj(MAINWINDOW ,'Tag','deletepoint'),'enable','on')
    
    if size(mri_scan.ptlabels,1) > 2
        set(findobj(MAINWINDOW,'Tag','PCS'),'enable','on')
    else
        set(findobj(MAINWINDOW,'Tag','PCS'),'enable','off')
    end
    
    %-----------------------------------------------------------------------------------------------
    
case 'seepoints'
    mri_scan = get(gcbf,'Userdata');
    pthandle = mri_scan.pthandle ;
    seepoints = findobj(gcbf,'Tag','seepoints');
    
    if mri_scan.seepoints == 1
        set(seepoints,'String', 'Hide Points');
        set(pthandle,'Visible','on');
        mri_scan.seepoints = 0;
    else
        set(seepoints,'String', 'See All Points');
        set(pthandle,'Visible','off');
        mri_scan.seepoints = 1;
    end
    
    set(gcf,'Userdata',mri_scan);
    
    %-----------------------------------------------------------------------------------------------
    
case 'seelabels'
    mri_scan = get(gcbf,'Userdata');
    ptlabels = mri_scan.ptlabels;
    seelabels = findobj(gcbf,'Tag','seelabels');
    
    if mri_scan.seelabels == 1
        set(seelabels,'String', 'Hide Labels');
        set(ptlabels,'Visible','on');
        mri_scan.seelabels = 0;
    else
        set(seelabels,'String','Show Labels');
        set(ptlabels,'Visible','off');
        mri_scan.seelabels = 1;
    end
    
    set(gcf,'Userdata',mri_scan);
    
    %-----------------------------------------------------------------------------------------------
    
case 'newlabel' % Assign a new label to the current point in the listbox
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    
    mri_scan = get(gcbf,'Userdata');
    ptlabels = mri_scan.ptlabels;
    ptname = mri_scan.ptname;
    seelabels = findobj(gcbf,'Tag','seelabels');
    
    prompt  = {'Enter a new label for the current point:'}; 
    title   = 'Please assign a label to the selected point';
    lines = 1;
    def     = {'I''d rather go to the beach'};
    answer  = inputdlg(prompt,title,lines,def);
    
    if isempty(answer)
        return
    else
        listboxx = findobj(MAINWINDOW,'Tag','listbox');
        current_pt = get(listboxx,'Value');
        ptname(current_pt) = answer;
        set(listboxx,'String',ptname)
        mri_scan.ptname = ptname;
        set(gcf,'Userdata',mri_scan);
    end
    
    %-----------------------------------------------------------------------------------------------   
case 'current_point'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    mri = mri_scan.mri;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    x = get(s3,'Value');
    y = get(s1,'Value');
    z = get(s2,'Max')-get(s2,'Value');
    ptname =  mri_scan.ptname ;
    ptcoord =  mri_scan.ptcoord;
    pthandle =  mri_scan.pthandle ;
    ptcurrent = mri_scan.ptcurrent;
    
    switch get(gcbo,'Tag')
        
    case 'previous'
        if ptcurrent > 1
            ptcurrent = ptcurrent - 1;
        end
    case 'next'
        if ptcurrent < size(ptcoord,1)
            ptcurrent = ptcurrent+1;
        end
    case 'listbox'
        ptcurrent = get(gcbo,'Value');
    end
    
    listboxx = findobj(MAINWINDOW,'Tag','listbox');
    set(listboxx,'Value',ptcurrent)
    
    
    mri_scan.ptcurrent = ptcurrent;
    x = ptcoord(ptcurrent,1);   
    y = ptcoord(ptcurrent,2);   
    z = ptcoord(ptcurrent,3);   
    
    set(s1,'Value',y);
    set(s2,'Value',get(s2,'Max') - z);
    set(s3,'Value',x);
    
    mri_scan.gotopoint = 1;
    set(MAINWINDOW ,'Userdata',mri_scan);
    
    mritool_cb sliders   
    
    mri_scan.gotopoint = 0;
    set(MAINWINDOW ,'Userdata',mri_scan);
    %-----------------------------------------------------------------------------------------------
case 'deletepoint'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    mri = mri_scan.mri;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    x = get(s3,'Value');
    y = get(s1,'Value');
    z = get(s2,'Max')-get(s2,'Value');
    ptname =  mri_scan.ptname ;
    ptcoord =  mri_scan.ptcoord;
    pthandle =  mri_scan.pthandle ;
    ptcurrent = mri_scan.ptcurrent;
    
    
    
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    x = get(s3,'Value');
    y = get(s1,'Value');
    z = get(s2,'Max')-get(s2,'Value');
    ptname =  mri_scan.ptname ;
    ptcoord =  mri_scan.ptcoord;
    pthandle =  mri_scan.pthandle ;
    ptcurrent = mri_scan.ptcurrent;
    ptlabels = mri_scan.ptlabels;
    
    delete(ptlabels(ptcurrent,:))
    delete(pthandle(ptcurrent,:))
    
    ptname(ptcurrent) = [];
    ptcoord(ptcurrent,:) = [];
    pthandle(ptcurrent,:) = [];
    ptlabels(ptcurrent,:)= [];
    
    
    
    ptcurrent = 1;
    
    if isempty(ptlabels)
        set(findobj(MAINWINDOW,'Tag','deletepoint'),'enable','off')
    end
    
    listboxx = findobj(MAINWINDOW,'Tag','listbox');
    set(listboxx,'Value',1,'String',ptname)
    set(listboxx,'Value',ptcurrent)
    
    mri_scan.ptname = ptname ;
    mri_scan.ptcoord   = ptcoord;
    mri_scan.pthandle = pthandle ;
    mri_scan.ptcurrent = ptcurrent;
    mri_scan.ptlabels = ptlabels;
    
    set(MAINWINDOW ,'Userdata',mri_scan);
    
    if size(mri_scan.ptlabels,1) > 2
        set(findobj(MAINWINDOW,'Tag','PCS'),'enable','on')
    else
        set(findobj(MAINWINDOW,'Tag','PCS'),'enable','off')
    end
    
    %-----------------------------------------------------------------------------------------------
case 'flip_lr_a'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    cor = mri_scan.cor;
    flip = get(findobj(gcbf,'Tag','flip_lr_a'),'Value');
    
    axa  = mri_scan.axa;
    
    
    mri_scan.mri = mri_scan.mri([size(mri_scan.mri,1):-1:1],:,:);
    
    s2 = findobj(   MAINWINDOW ,'Tag','2');
    s3 = findobj(   MAINWINDOW ,'Tag','3');
    pix_siz = mri_scan.pix_siz;
    ii = round(get(s2,'Value')/pix_siz(3)+1);
    str = ':,:,ii';
    eval(['mri_tmp = squeeze(mri_scan.mri(',str,'));']);
    set(ax,'Cdata',(mri_tmp));
    
    
    ic = round(get(s3,'Value')/pix_siz(1)+1);
    strc = ':,ic,:';
    eval(['mri_tmp = fliplr(flipud(squeeze(mri_scan.mri(',strc,'))''));']);
    set(cor,'Cdata',(mri_tmp));
    
    set(MAINWINDOW,'Userdata',mri_scan);
    
    %-----------------------------------------------------------------------------------------------
    
case 'flip_ud_s'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    axs = mri_scan.axs;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    sag = mri_scan.sag;
    cor = mri_scan.cor;
    
    %flip = get(findobj(gcbf,'Tag','flip_ud_s'),'Value');
    mri_scan.mri = mri_scan.mri(:,:,[size(mri_scan.mri,3):-1:1]);
    
    
    set(MAINWINDOW ,'Userdata',mri_scan);
    
    ii = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2)+1);
    str = 'ii,:,:';
    eval(['mri_tmp = flipud(squeeze(mri_scan.mri(',str,'))'');']);
    set(sag,'Cdata',(mri_tmp));
    
    ic = round(get(s3,'Value')/pix_siz(1)+1);
    
    strc = ':,ic,:';
    eval(['mri_tmp = fliplr(flipud(squeeze(mri_scan.mri(',strc,'))''));']);
    set(cor,'Cdata',(mri_tmp));
    
    %-----------------------------------------------------------------------------------------------
    
case 'flip_ud_c'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    axc = mri_scan.axc;
    flip = get(findobj(gcbf,'Tag','flip_ud_c'),'Value');
    if flip 
        set(axc,'Ydir','normal')
    else
        set(axc,'Ydir','reverse')
    end
    
    %-----------------------------------------------------------------------------------------------
    
case 'saveas_points'
    
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    if isempty(mri_scan)
        return
    end
    
    ptname =  mri_scan.ptname ;
    ptcoord =  mri_scan.ptcoord;
    
    if isempty(ptcoord)
        return
    end
    
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    x = ptcoord(:,1);
    y = ptcoord(:,2);
    z = ptcoord(:,3);
    
    
    pix_size = mri_scan.pix_size;
    CubeSize = size(mri_scan.mri);
    
    %VoxSize = mri_scan.pix_siz([2 1 3]);
    
    VoxSize = mri_scan.fov./CubeSize([2 1 3]);
    
    pix = VoxSize([2 1 3]);
    fov  = pix.* CubeSize;
    
    ix = ((fov(2)-y)/pix_size(2)+1);
    iy = (x/pix_size(1)+1);
    iz = ((fov(3)-z)/pix_size(3)+1);
    
    if mri_scan.pcs == 1    % If a PCS exists - get its paramters
        center = mri_scan.center;
        pcs = mri_scan.pcs_type;
        
        switch pcs
        case {'NEUROMAG','CTF'}
            if 0 
                rot = mri_scan.rot;
                tmp = rot*[x';y';z'];
                x = tmp(1,:);
                y = tmp(2,:);
                z = tmp(3,:);
            end
            PCS.R = inv(mri_scan.rot);
            PCS.t = mri_scan.center;
            
            [XYZpcs] = mri2pcs([ix,iy,iz], PCS, VoxSize([2 1 3]),CubeSize);
            
            if 0
                x = x - PCS.t(1);
                y = y - PCS.t(2);
                z = z - PCS.t(3);
                XYZpcs = PCS.R*([x,y,z]');
                XYZpcs = XYZpcs';
            end
        end
    end  
    
    [fpt_file,fpt_dir] = uiputfile('*.fpt','Save point locations in...');
    
    file = fopen([fpt_dir,fpt_file],'w+');
    if ~file
        errordlg('I/O Error',['Error trying to create ',fpt_file])
        return
    end
    
    for i=1:size(ptcoord,1)
        tmp = char(ptname{i});
        fprintf(file,'%d\t%s\t\t\t%6.2f\t%6.2f\t%6.2f\t%d\t%d\t%d\n',i,tmp,[XYZpcs(i,1),XYZpcs(i,2),XYZpcs(i,3)],round([ix(i),iy(i),iz(i)]));
    end
    fclose(file);
    
    mri_scan.fpt_file = [fpt_dir,fpt_file];
    
    %-----------------------------------------------------------------------------------------------
    
case 'leaflet' % display every fiducial or focus in a window with corresponding MR slices and coordinates
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    color = get(MAINWINDOW,'colormap');
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    ptname =  mri_scan.ptname ;
    ptcoord =  mri_scan.ptcoord;
    pthandle =  mri_scan.pthandle ;
    ptcurrent = mri_scan.ptcurrent;
    mri_scan.ptcurrent = ptcurrent;
    
    mri_scan.leaflet = 1;   
    toto = 0;
    
    
    % Determine the color of the source time series according to the MRI markers' colors
    h = findobj(flipud(get(axa,'Children')),'Type','patch');
    colors = get(h,'MarkerfaceColor');
    if ~iscell(colors) % Basically when there's just one source
        colors = {colors};
    end
    colors =  reshape([colors{:}],3,length([colors{:}])/3)';
    
    for ptcurrent = 1:size(pthandle,1)
        x = ptcoord(ptcurrent,1);   
        y = ptcoord(ptcurrent,2);   
        z = ptcoord(ptcurrent,3);   
        set(s1,'Value',y);
        set(s2,'Value',get(s2,'Max') - z);
        set(s3,'Value',x);
        mri_scan.gotopoint = 1;
        set(MAINWINDOW ,'Userdata',mri_scan);
        mritool_cb sliders   
        
        leaflet
        leaflets = findobj(0,'Tag','FIG_LEAFLET');
        set(0,'CurrentFigure',leaflets(1));
        set(gcf,'Name',['Leaflet - ',upper(char(ptname{ptcurrent}))],'colormap',color)
        l_axs = findobj(gcf,'Tag','leaf_axs');
        pos = get(l_axs,'Position');
        delete(l_axs)
        l_axs = copyobj(axs,gcf);
        set(l_axs,'Position',pos)
        
        l_axc = findobj(gcf,'Tag','leaf_axc');
        pos = get(l_axc,'Position');
        delete(l_axc)
        l_axc = copyobj(axc,gcf);
        set(l_axc,'Position',pos)
        
        l_axa = findobj(gcf,'Tag','leaf_axa');
        pos = get(l_axa,'Position');
        delete(l_axa)
        l_axa = copyobj(axa,gcf);
        set(l_axa,'Position',pos)
        
        h = copyobj(flipud(get(axs,'Children')),l_axs); 
        % flipud to plot the MRI scan first then the points above
        % otherwise points are not visible (sent to the back of the picture)
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axs)
        %axis equal
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        set(h(i_hh),'Markersize',2)   
        
        
        h = copyobj(flipud(get(axa,'Children')),l_axa); 
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axa)
        %axis equal
        axis ij
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        set(h(i_hh),'Markersize',2)   
        
        h = copyobj(flipud(get(axc,'Children')),l_axc); 
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axc)
        %axis equal
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        
        
        ax_time = findobj(gcf,'Tag','timeseries');
        if isempty(findobj(0,'Tag','mri_src_mapping'))
            set(ax_time,'Visible','off')
            couleur{ptcurrent} = 'c';
        else
            Results  = get(findobj(0,'Tag','mri_src_mapping'),'Userdata');
            S = Results.TimeSeries;
            
            % Get the Time Values
            User = get_user_directory;
            cd(User.STUDIES)
            load(Results.StudySubject.Data{1},'Time'); % in msec
            Time = 1000*Time(Results.Time);% in msec
            
            t = Time;
            
            Sn = colnorm(S);
            Sn = S./Sn(ones(1,size(S,1)),:);
            % each waveform is now unity in norm. Weight
            %  by its total contribution to the field
            axes(ax_time)
            %tmp = set(ax_time,'colororder');
            set(ax_time,'colororder',colors,'color',get(gcf,'color'))
            grid on
            hold on
            hh = plot(t,S);
            if ptcurrent > 1
                if strcmp(ptname{ptcurrent},ptname{ptcurrent-1}) % 2 synchronous sources
                    toto = toto - 1;
                else
                    toto = toto + 1;
                end
            else
                toto = 1;
            end
            
            
            set(hh(toto),'linewidth',3)
            set(ax_time,'xcolor',frontcolor,'ycolor',frontcolor)
            grid on
            %col = get(hh,'color');
            %col = col{toto};
            %set(findobj(gcf,'Type','patch'),'MarkerFaceColor',col,'MarkerEdgeColor',col)   
            %         couleur{ptcurrent} = col;
            tmp =  get(findobj(gcf,'Type','patch'),'MarkerFaceColor');       
            couleur{ptcurrent}  = tmp{1};
        end
        
        
        
    end
    
    if size(pthandle,1) > 1
        seepoints = findobj(gcbf,'Tag','seepoints');
        set(seepoints,'String', 'Hide Points');
        set(pthandle,'Visible','on');
        
        if 0
            for i = 1:size(pthandle,1)
                col = couleur{i};
                %set(pthandle(i,:),'MarkerFaceColor',col,'MarkerEdgeColor',col)
            end
        end
        
        
        mri_scan.seepoints = 0;
        
        leaflet
        leaflets = findobj(0,'Tag','FIG_LEAFLET');
        set(0,'CurrentFigure',leaflets(1));
        set(gcf,'Name','Leaflet - All Sources','colormap',color)
        l_axs = findobj(gcf,'Tag','leaf_axs');
        pos = get(l_axs,'Position');
        delete(l_axs)
        l_axs = copyobj(axs,gcf);
        set(l_axs,'Position',pos)
        
        l_axc = findobj(gcf,'Tag','leaf_axc');
        pos = get(l_axc,'Position');
        delete(l_axc)
        l_axc = copyobj(axc,gcf);
        set(l_axc,'Position',pos)
        
        l_axa = findobj(gcf,'Tag','leaf_axa');
        pos = get(l_axa,'Position');
        delete(l_axa)
        l_axa = copyobj(axa,gcf);
        set(l_axa,'Position',pos)
        
        h = copyobj(flipud(get(axs,'Children')),l_axs); 
        % flipud to plot the MRI scan first then the points above
        % otherwise points are not visible (sent to the back of the picture)
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axs)
        %axis equal
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        set(h(i_hh),'Markersize',2)   
        
        
        h = copyobj(flipud(get(axa,'Children')),l_axa); 
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axa)
        %axis equal
        axis ij
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        set(h(i_hh),'Markersize',2)   
        
        h = copyobj(flipud(get(axc,'Children')),l_axc); 
        h = findobj(h,'Type','patch');
        if isempty(h), return, end
        set(gcf,'CurrentAxes',l_axc)
        %axis equal
        a = get(h,'tag');
        i_hh = strmatch(char(ptname{ptcurrent}),[a(:)],'exact');
        set(h(i_hh),'visible','off')
        set(h(i_hh),'visible','on')   
        set(h(i_hh),'Markersize',2)   
        
        ax_time = findobj(gcf,'Tag','timeseries');
        if isempty(findobj(0,'Tag','mri_src_mapping'))
            set(ax_time,'Visible','off')
        else
            Results  = get(findobj(0,'Tag','mri_src_mapping'),'Userdata');
            S = Results.TimeSeries;
            
            t = Time;
            
            Sn = colnorm(S);
            Sn = S./Sn(ones(1,size(S,1)),:);
            % each waveform is now unity in norm. Weight
            %  by its total contribution to the field
            axes(ax_time)
            grid on
            set(ax_time,'colororder',colors,'color',get(gcf,'color'))
            hold on
            hh = plot(t,S);
            set(hh,'linewidth',2)
            set(ax_time,'xcolor',frontcolor,'ycolor',frontcolor)
            grid on
            
        end
        
    end
    
    
    mri_scan.leaflet = 0;
    mri_scan.gotopoint = 0;
    set(MAINWINDOW ,'Userdata',mri_scan);
    
    %-----------------------------------------------------------------------------------------------
    
case 'showfiducials'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    
    ptname = mri_scan.FiducialName;      
    ptcoord = [mri_scan.ptcoord;mri_scan.PCSFiducial'];
    
    % Transform PCS coordinates in the MRI volume reference system
    switch mri_scan.pcs_type
    case {'NEUROMAG','CTF'}
        if ~isfield(mri_scan,'rot')
            errordlg('Please Define a PCS First')
            return
        end
        center = mri_scan.center;
        rot = mri_scan.rot;
        ptcoord = (inv(rot)*ptcoord')';      
    end
    
    ptcoord(:,1) = ptcoord(:,1) + center(1);
    ptcoord(:,2) = ptcoord(:,2) + center(2);
    ptcoord(:,3) = ptcoord(:,3) + center(3);
    
    switch mri_scan.pcs_type
    case {'NEUROMAG','CTF'}
        icoord(:,1) = round(ptcoord(:,1)/pix_siz(3)+1);
        icoord(:,2) = round(ptcoord(:,2)/pix_siz(1)+1);
        icoord(:,3) = round(ptcoord(:,3)/pix_siz(2)+1);
    end
    
    for i = 1:size(ptcoord,1)
        x = ptcoord(i,1);
        y = ptcoord(i,2);
        z = ptcoord(i,3);
        
        axes(axs)
        hold on
        pt_s = scatter(x,z,'o','filled');
        set(pt_s,'Parent',axs,'userdata',y)
        pt_slabel = text(x,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_slabel,'Parent',axs,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',y)
        
        axes(axa)
        hold on
        pt_a = scatter(x,get(s1,'Max')-y,'o','filled');
        set(pt_a,'Parent',axa,'userdata',z)
        pt_alabel = text(x,get(s1,'Max')-y,num2str(size(mri_scan.pthandle,1)));
        set(pt_alabel,'Parent',axa,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',z)
        
        axes(axc)
        hold on
        pt_c = scatter(y,z,'o','filled');
        set(pt_c,'Parent',axc,'userdata',x)
        pt_clabel = text(y,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_clabel,'Parent',axc,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',x)
        
        set([pt_s,pt_a,pt_c],'markerfacecolor','c')
        mri_scan.pthandle = [mri_scan.pthandle;[pt_s,pt_a,pt_c]];
        mri_scan.ptlabels= [mri_scan.ptlabels;[pt_slabel,pt_alabel,pt_clabel]];
        
    end
    
    seepoints = findobj(gcbf,'Tag','seepoints');
    seelabels = findobj(gcbf,'Tag','seelabels');
    newlabel = findobj(gcbf,'Tag','newlabel');
    
    LEAFLET = findobj(gcbf,'Tag','LEAFLET');
    set(seelabels,'enable','on')
    set(newlabel,'enable','on')
    
    set(seepoints,'enable','on')
    set(LEAFLET,'enable','on')
    
    set(findobj(gcbf,'Tag','deletepoint'),'enable','on')
    
    listboxx = findobj(MAINWINDOW,'Tag','listbox');
    set(listboxx,'Visible','on')
    listboxx_title = findobj(MAINWINDOW,'Tag','listboxx_title');
    set(listboxx_title,'visible','on')
    listboxx = findobj(listboxx,'Tag','listbox');
    set(listboxx,'Value',1,'String',ptname)
    
    mri_scan.ptcurrent = size( mri_scan.ptcoord,1);
    
    mri_scan.ptname = ptname ;
    mri_scan.ptcoord = ptcoord ;
    
    visu_pt(mri_scan.pthandle)
    
    set(MAINWINDOW,'Userdata',mri_scan)
    
    
    
    points'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    pix_siz = mri_scan.pix_siz;
    s1 = findobj(MAINWINDOW ,'Tag','1');
    s2 = findobj(MAINWINDOW ,'Tag','2');
    s3 = findobj(MAINWINDOW ,'Tag','3');
    
    % Delete all previous points
    delete(findobj(get(axa,'children'),'type','patch'))
    delete(findobj(get(axc,'children'),'type','patch'))
    delete(findobj(get(axs,'children'),'type','patch'))
    
    [fpt_file,fpt_dir] = uigetfile('*.fpt','Load fiducial coordinates from...');
    if isempty(fpt_file), return, end
    
    
    file = fopen([fpt_dir,fpt_file],'r');
    if ~file
        errordlg('I/O Error',['Error while loading from ',fpt_file])
        return
    end
    i = 0;
    mri_scan.pthandle = [];
    mri_scan.ptlabels = [];
    %mri_scan.pcs = 0;
    
    while ~isempty(i)
        i = fscanf(file,'%d',1);      
        if isempty(i), break, end
        ptname{i} = fscanf(file,'%s',1);      
        ptcoord(i,:) = fscanf(file,'%f',3)';
        icoord(i,:) = fscanf(file,'%d',3)';
        
        if mri_scan.pcs == 0
            orig = findstr(ptname{i},'originPCS');
            if ~isempty(orig) % This set of points can be registered owing to a MEG-EEG PCS.
                mri_scan.pcs = 1; 
                mri_scan.pcs_type = ptname{i}(findstr(ptname{i},'_')+1:end);
                mri_scan.center = icoord;
                center = icoord(i,:);
                ptcenter = ptcoord(i,:)
            else 
                mri_scan.pcs = 0; 
                mri_scan.pcs_type = 'NONE';
                center = [0 0 0];
            end
        end
        
    end
    fclose all;
    
    % Transform PCS coordinates in the MRI volume reference system
    switch mri_scan.pcs_type
    case 'NEUROMAG'
        %xtemp = ptcoord(:,1);
        %%%%%%ptcoord(:,1) = -ptcoord(:,2);
        %ptcoord(:,1) = ptcoord(:,2);
        %ptcoord(:,2) = xtemp;      
        if ~isfield(mri_scan,'rot')
            errordlg('Please Define a PCS First')
            return
        end
        center = mri_scan.center;
        rot = mri_scan.rot;
        ptcoord = (inv(rot)*ptcoord')';      
    end
    
    ptcoord(:,1) = ptcoord(:,1) + center(1);
    ptcoord(:,2) = ptcoord(:,2) + center(2);
    ptcoord(:,3) = ptcoord(:,3) + center(3);
    
    switch mri_scan.pcs_type
    case 'NEUROMAG'
        %      ptcoord(end,:) = ptcenter;
    end
    
    
    icoord(:,1) = round(ptcoord(:,1)/pix_siz(3)+1);
    icoord(:,2) = round(ptcoord(:,2)/pix_siz(1)+1);
    icoord(:,3) = round(ptcoord(:,3)/pix_siz(2)+1);
    
    
    for i = 1:size(ptcoord,1)
        x = ptcoord(i,1);
        y = ptcoord(i,2);
        z = ptcoord(i,3);
        
        %    set(axs,'Selected','on')
        
        axes(axs)
        hold on
        pt_s = scatter(x,z,'o','filled');
        set(pt_s,'Parent',axs,'userdata',y)
        pt_slabel = text(x,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_slabel,'Parent',axs,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',y)
        
        axes(axa)
        hold on
        pt_a = scatter(x,get(s1,'Max')-y,'o','filled');
        set(pt_a,'Parent',axa,'userdata',z)
        pt_alabel = text(x,get(s1,'Max')-y,num2str(size(mri_scan.pthandle,1)));
        set(pt_alabel,'Parent',axa,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',z)
        
        axes(axc)
        hold on
        pt_c = scatter(y,z,'o','filled');
        set(pt_c,'Parent',axc,'userdata',x)
        pt_clabel = text(y,z,num2str(size(mri_scan.pthandle,1)));
        set(pt_clabel,'Parent',axc,'VerticalAlignment','top','Visible','off','Color','y','FontW','Bold',...
            'Userdata',x)
        
        set([pt_s,pt_a,pt_c],'markerfacecolor','c')
        mri_scan.pthandle = [mri_scan.pthandle;[pt_s,pt_a,pt_c]];
        mri_scan.ptlabels= [mri_scan.ptlabels;[pt_slabel,pt_alabel,pt_clabel]];
        
    end
    
    seepoints = findobj(gcbf,'Tag','seepoints');
    seelabels = findobj(gcbf,'Tag','seelabels');
    newlabel = findobj(gcbf,'Tag','newlabel');
    
    LEAFLET = findobj(gcbf,'Tag','LEAFLET');
    set(seelabels,'enable','on')
    set(newlabel,'enable','on')
    set(seepoints,'enable','on')
    set(LEAFLET,'enable','on')
    
    set(findobj(gcbf,'Tag','deletepoint'),'enable','on')
    
    listboxx = findobj(MAINWINDOW,'Tag','listbox');
    set(listboxx,'Visible','on')
    listboxx_title = findobj(MAINWINDOW,'Tag','listboxx_title');
    set(listboxx_title,'visible','on')
    legend_sources = findobj(MAINWINDOW,'Tag','legend_sources');
    set(legend_sources,'visible','on')
    listboxx = findobj(listboxx,'Tag','listbox');
    set(listboxx,'Value',1,'String',ptname)
    
    mri_scan.ptcurrent = size( mri_scan.ptcoord,1);
    
    mri_scan.ptname = ptname ;
    mri_scan.ptcoord = ptcoord ;
    
    visu_pt(mri_scan.pthandle)
    
    set(MAINWINDOW,'Userdata',mri_scan)
    
    %-----------------------------------------------------------------------------------------------
    
case 'readresults'
    
    mri_src_mapping 
    
    %-----------------------------------------------------------------------------------------------
    
case 'concatenate'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    Header = findobj(gcbf,'Tag','header');
    header = str2num(get(Header,'String'));
    
    C = get(findobj(MAINWINDOW,'Tag','CurrentDrive'),'String');
    cd(C);
    
    d = dir;
    str1 = {d.name};
    [s1,v] = listdlg('PromptString','First Set of Files :',...
        'SelectionMode','multiple',...
        'ListString',str1);
    
    
    if v == 0
        return
    end
    
    while length(s1)==1
        cd(str1{s1});
        d = dir;
        str1 = {d.name};
        
        [s1,v] = listdlg('ListString',str1);
        
        if v==0
            return
        end
        
    end
    
    dir1 = pwd;
    nfile1 = length(s1);
    
    %-------------------------
    
    d = dir;
    str2 = {d.name};
    [s2,v] = listdlg('PromptString','Second Set of Files :',...
        'SelectionMode','multiple',...
        'ListString',str2);
    
    if v == 0
        return
    end
    
    while length(s2)==1
        cd(str2{s2});
        d = dir;
        str2 = {d.name};
        [s2,v] = listdlg('PromptString','Select a file :',...
            'SelectionMode','multiple',...
            'ListString',str2);
        if v==0
            return
        end
        
    end
    dir2 = pwd;
    nfile2 = length(s2);
    
    if nfile1~=nfile2
        errodlg('Please Select the Same Number of Files for Both Sets')
        return
    end
    
    data_format = findobj(gcbf,'Tag','data_format');  
    code = (get(data_format,'String'));   
    tmp = get(data_format,'Value');
    code = strtok(char(code(tmp)));
    
    [rootfile,dirroot] = uiputfile('toto.img','Enter a Root File Name...');
    hh = waitbar(0,'Concatenating files...');   
    for i=1:nfile1
        root = strtok(rootfile,'.');
        if isempty(root)
            root= rootfile;
        end
        
        if i<10
            str = '00';
        elseif i<100
            str = '0';
        else
            str = '';
        end
        
        f1 = fopen([dir1,'\',str1{s1(i)}],'r');
        fseek(f1,header,'bof');
        f2 = fopen([dir2,'\',str2{s2(i)}],'r');
        fseek(f2,header,'bof');
        file = fopen([dirroot,root,str,num2str(i),'.img'],'w+');
        
        if ~file
            errordlg('I/O Error',['Error trying to create ',root])
            delete(hh)
            return
        end   
        
        mri1 = fliplr(fread(f1,[128,128],code)); % Convention spm
        mri2 = fliplr(fread(f2,[128 128],code));
        
        fwrite(file,[mri1,mri2],code);
        
        fclose all;
        waitbar(i/nfile1);
        
    end
    delete(hh)
    
    %-----------------------------------------------------------------------------------------------
    
    
case 'saveaxial'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata'); 
    nslice = mri_scan.nslice;
    
    data_format = findobj(gcbf,'Tag','data_format');  
    code = (get(data_format,'String'));   
    tmp = get(data_format,'Value');
    code = strtok(char(code(tmp)));
    
    
    [rootfile,rootdir] = uiputfile('*.img','Enter a Root File Name...');
    
    
    hh = waitbar(0,'Writing files...');   
    for i=1:size(mri_scan.mri,3)
        root = strtok(rootfile,'.');
        if isempty(root)
            root= rootfile;
        end
        
        if i<10
            str = '00';
        elseif i<100
            str = '0';
        else
            str = '';
        end
        
        file = fopen([rootdir,root,str,num2str(i),'.img'],'w+');
        if ~file
            errordlg('I/O Error',['Error trying to create ',root])
            delete(hh)
            return
        end   
        
        fwrite(file,(mri_scan.mri(:,:,i)),code);
        
        fclose(file);
        waitbar(i/nslice);
        
    end
    delete(hh)
    
    %-----------------------------------------------------------------------------------------------
    
case 'gotopoint'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata'); 
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    nslice = mri_scan.nslice;
    fov = mri_scan.fov;
    siz = mri_scan.siz;
    ic = mri_scan.ic ;
    is = mri_scan.is ;
    ia = mri_scan.ia ;
    pix_siz = mri_scan.pix_siz ;
    s1 = findobj(gcbf,'Tag','1');
    s2 = findobj(gcbf,'Tag','2');
    s3 = findobj(gcbf,'Tag','3');
    show_coordinates = findobj(   MAINWINDOW ,'Tag','show_coordinates');
    if get(show_coordinates,'Value')
        SHOW_COORDINATES = 1;
    else
        SHOW_COORDINATES = 0;
    end
    pix_siz = mri_scan.pix_siz ;
    
    
    MOUSE_SLCT = findobj(MAINWINDOW ,'tag','mouse_slct');
    set(MOUSE_SLCT,'value',0);
    gotopoint = findobj(MAINWINDOW,'Tag','gotopoint');
    
    if ~get(gotopoint,'Value')
        return
    end
    
    
    switch get(MAINWINDOW,'CurrentObject')
    case ax
        tmp = get(axa,'Currentpoint');
        pt = tmp(1,1:2);
        x = pt(1);
        y = pt(2);
        set(s1,'Value',get(s1,'Max')-y);
        set(s3,'Value',x);
        ic = round(get(s3,'Value')/pix_siz(1)+1);
        
        if SHOW_COORDINATES 
            set(axs,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
            set(axa,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
        end
        
        is = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2)+1);
        if SHOW_COORDINATES 
            set(axc,'Xtick',get(s1,'Value'),'Xticklabel',num2str(get(s1,'Value'),FORMATtick))
            set(axa,'Ytick',get(s1,'Max')-get(s1,'Value'));
        end
        
        strs = 'is,:,:';
        strc = ':,ic,:';
        eval(['mri_tmp = flipud(squeeze(mri_scan.mri(',strs,'))'');']);
        set(sag,'Cdata',(mri_tmp));
        eval(['mri_tmp = fliplr(flipud(squeeze(mri_scan.mri(',strc,'))''));']);
        set(cor,'Cdata',(mri_tmp));
        
    case cor
        tmp = get(axc,'Currentpoint');
        pt = tmp(1,1:2);
        y = pt(1);
        z = pt(2);
        set(s1,'Value',y);
        set(s2,'Value',get(s2,'Max')-z);
        is = round((get(s1,'Max')-get(s1,'Value'))/pix_siz(2)+1);
        if SHOW_COORDINATES 
            set(axc,'Xtick',get(s1,'Value'),'Xticklabel',num2str(get(s1,'Value'),FORMATtick))
            set(axa,'Ytick',get(s1,'Max')-get(s1,'Value'));
        end
        
        ia = round(get(s2,'Value')/pix_siz(3)+1);
        if SHOW_COORDINATES 
            set(axc,'Ytick',get(s2,'Max') - get(s2,'Value'),'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
            set(axs,'Ytick',get(s2,'Max') - get(s2,'Value'));
        end
        
        strs = 'is,:,:';
        stra = ':,:,ia';
        eval(['mri_tmp = flipud(squeeze(mri_scan.mri(',strs,'))'');']);
        set(sag,'Cdata',(mri_tmp));
        eval(['mri_tmp = squeeze(mri_scan.mri(',stra,'));']);
        set(ax,'Cdata',(mri_tmp));
        
    case sag
        tmp = get(axs,'Currentpoint');
        pt = tmp(1,1:2);
        x = pt(1);
        z = pt(2);
        set(s3,'Value',x);
        set(s2,'Value',get(s2,'Max')-z);
        ic = round(get(s3,'Value')/pix_siz(1)+1);
        if SHOW_COORDINATES 
            set(axs,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
            set(axa,'Xtick',get(s3,'Value'),'Xticklabel',num2str(get(s3,'Value'),FORMATtick))
        end
        
        ia = round(get(s2,'Value')/pix_siz(3)+1);
        if SHOW_COORDINATES 
            set(axc,'Ytick',get(s2,'Max') - get(s2,'Value'),'Yticklabel',num2str(get(s2,'Value'),FORMATtick ))
            set(axs,'Ytick',get(s2,'Max') - get(s2,'Value'));
        end
        
        stra = ':,:,ia';
        strc = ':,ic,:';
        eval(['mri_tmp = fliplr(flipud(squeeze(mri_scan.mri(',strc,'))''));']);
        set(cor,'Cdata',(mri_tmp));
        eval(['mri_tmp = squeeze(mri_scan.mri(',stra,'));']);
        set(ax,'Cdata',(mri_tmp));
        
    end
    
    %-----------------------------------------------------------------------------------------------
    
case 'ZoomButton'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata'); 
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    
    fig = figure;
    if isempty(mri_scan.otherfig )
        mri_scan.otherfig = fig;
    else
        mri_scan.otherfig = [  mri_scan.otherfig,fig];
    end
    
    gcob = get(MAINWINDOW,'CurrentObject');
    switch get(gcob,'Tag')
    case 'ZoomButton1'
        handle = axs;
    case 'ZoomButton2'
        handle = axa;
    case 'ZoomButton3'
        handle = axc;
    end
    
    color = get(MAINWINDOW,'colormap');
    handle = copyobj([handle],fig);
    
    set(handle,'Units','Normal','Position',[.01 .01 .98 .98])
    set(MAINWINDOW ,'Userdata' , mri_scan ); 
    set(0,'CurrentFigure',fig)
    set(fig,'colormap',color,'color',color(1,:),'Numbertitle','off','Name',['Zooooooom ',int2str(length( mri_scan.otherfig))]);
    set(fig,'Deletefcn','mritool_cb DeleteOtherFig')   
    colormenu
    set(gcob,'Value',0)   
    axes(handle)
    set(handle,'color',color(1,:));
    
    axis equal
    %-----------------------------------------------------------------------------------------------
    
case 'quit'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');    
    if ~isempty(mri_scan)
        if ~isempty(mri_scan.otherfig)
            delete(mri_scan.otherfig)
        end
    end
    
    
    close(MAINWINDOW)
    close(findobj(0,'Tag','fig_listofpoints'))
    
    %-----------------------------------------------------------------------------------------------
    
case 'DeleteOtherFig'
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata');    
    fig = (gcbf);
    mri_scan.otherfig = mri_scan.otherfig(mri_scan.otherfig ~= fig);
    set(MAINWINDOW,'Userdata',mri_scan);
    %-----------------------------------------------------------------------------------------------
    
    
case 'pcs'
    % Change of coordinate system -  to Neuromag, BTi, Ctf  -
    MAINWINDOW = findobj(0,'Tag','mainwindow');
    mri_scan = get(MAINWINDOW ,'Userdata'); 
    ax = mri_scan.ax;
    axa = mri_scan.axa;
    sag = mri_scan.sag;
    axs = mri_scan.axs;
    cor = mri_scan.cor;
    axc = mri_scan.axc;
    nslice = mri_scan.nslice;
    fov = mri_scan.fov;
    siz = mri_scan.siz;
    ic = mri_scan.ic ;
    is = mri_scan.is ;
    ia = mri_scan.ia ;
    pix_siz = mri_scan.pix_siz ;
    s1 = findobj(gcf,'Tag','1');
    s2 = findobj(gcf,'Tag','2');
    s3 = findobj(gcf,'Tag','3');
    pthandle = mri_scan.pthandle;
    ptlabels = mri_scan.ptlabels;
    ptcoord =  mri_scan.ptcoord;
    ptname = mri_scan.ptname;
    pcs_choice
    PCSWin = gcf;
    quit = findobj(PCSWin,'Tag','QUIT');
    set(gcf,'CurrentObject',findobj(gcf,'Tag','NEUROMAG'))
    while get(gcf,'CurrentObject') ~= quit
        waitforbuttonpress
        set(quit,'String',get(gco,'String'))
        if get(gcf,'CurrentObject') ~= quit
            pcs = get(get(gcf,'CurrentObject'),'Tag');
        else 
            delete(gcf)
            break   
        end
    end
    
    listboxx = findobj(0,'Tag','listbox');
    
    dummy = axes;
    set(dummy,'Visible','off')
    set(MAINWINDOW,'CurrentObject',dummy)
    
    h = msgbox('Please Select the Nasion Point in the List');
    waitfor(h)
    waitfor(MAINWINDOW,'CurrentObject',listboxx)
    NAS = get(listboxx,'Value');
    
    h = msgbox('Please Select the RIGHT Preauricular Point in the List');
    waitfor(h)
    waitfor(listboxx,'Value')
    RPA = get(listboxx,'Value');
    
    h = msgbox('Please Select the LEFT Preauricular Point in the List');
    waitfor(h)
    waitfor(listboxx,'Value')
    LPA = get(listboxx,'Value');
    delete(dummy)
    
    for i=1:3
        set(pthandle(NAS,i),'Markerfacecolor','y')
        set(pthandle(LPA,i),'Markerfacecolor','y')
        set(pthandle(RPA,i),'Markerfacecolor','y')
    end
    
    % Get the coordinates of the fiducials
    
    str_listboxx = get(listboxx,'String');
    mri_scan.FiducialName = {str_listboxx{NAS},str_listboxx{LPA},str_listboxx{RPA},'Origin'};
    
    NAS = ptcoord(NAS,:);
    LPA = ptcoord(LPA,:);
    RPA = ptcoord(RPA,:);
    
    scale = [pix_siz(2),pix_siz(1),pix_siz(3)];
    
    mri_scan.CubeFiducial = round([NAS./scale;LPA./scale;RPA./scale]'); % Indices of the slices containing the fiducials
    
    
    % New Coordinate System - PCS
    switch pcs
    case 'NEUROMAG'
        LR = LPA - RPA; % Left\Right Vector
        LR = LR/norm(LR);  
        LRp = [0 -1 0;1 0 0;0 0 0]* LR'; % A vector orthogonal to LR that will define a plane perpendicular to the NLR plane
        LRp = LRp'/norm(LRp);
        % note: The LR line is defined according to the intersection between 2 orthogonal planes.
        
        nLR= cross(NAS-RPA,LPA-RPA);  % Vector normal to the NLR plane
        nLR = nLR/norm(nLR);
        
        A = [nLR;LRp;LR];
        B = sum(A.*[LPA;LPA;NAS],2);
        
        center = A\B;
        
        % Rotation associee 
        vy = (RPA-center')'/norm((RPA-center'));
        vx = (NAS-center')'/norm((NAS-center'));
        rot = [0 1 0;1 0 0;0 0 1]*inv([vx,vy,cross(vx,vy)]);
        mri_scan.pcs_type = 'NEUROMAG';
        
    case 'CTF'
        old = 0;
        
        if old
            LR = LPA - RPA; % Left\Right Vector
            LR = LR/norm(LR);  
            LRp = [0 -1 0;1 0 0;0 0 0]* LR'; % A vector orthogonal to LR that will define a plane perpendicular to the NLR plane
            LRp = LRp'/norm(LRp);
            % note: The LR line is defined according to the intersection between 2 orthogonal planes.
            
            nLR= cross(NAS-RPA,LPA-RPA);  % Vector normal to the NLR plane
            nLR = nLR/norm(nLR);
            
            A = [nLR;LRp;LR];
            B = sum(A.*[LPA;LPA;NAS],2);
            
            center = A\B;
            
            % Rotation associee 
            vy = (RPA-center')'/norm((RPA-center'));
            vx = (NAS-center')'/norm((NAS-center'));
        else
            center = .5*(LPA+RPA)';
            nLR= cross(NAS-RPA,LPA-RPA);  % Vector normal to the NLR plane
            nLR = nLR/norm(nLR);
            
            vx = (NAS-center')'/norm((NAS-center'));
            vy = inv([vx';nLR;cross(vx',nLR)])*[0 0 1]'; % y direction
            vy = vy/norm(vy);
            
        end
        
        rot = [0 1 0;1 0 0;0 0 1]*inv([vx,vy,cross(vx,vy)]);
        rot = [0 1 0;-1 0 0;0 0 1]* rot;
        mri_scan.pcs_type = 'CTF';
        
    end
    
    mri_scan.CubeFiducial = [mri_scan.CubeFiducial,round(center./scale')];
    
    x = center(1);
    y = center(2);
    z = center(3);
    
    center_s = scatter(x,z,'o','filled');
    set(center_s,'Parent',axs,'Userdata',y,'Visible','on','MarkerfaceColor','r')
    
    
    center_a = scatter(x,get(s1,'Max')-y,'o','filled');
    set(center_a,'Parent',axa,'Userdata',z,'Visible','on','MarkerfaceColor','r')
    
    center_c = scatter(y,z,'o','filled');
    set(center_c,'Parent',axc,'Userdata',x,'Visible','on','MarkerfaceColor','r')
    
    mri_scan.pcs = 1;
    mri_scan.center = center;
    mri_scan.RPA = RPA;
    mri_scan.LPA = LPA;
    mri_scan.NAS = NAS;
    mri_scan.rot = rot; 
    
    
    % Coordinates of the fiducials in the PCS
    fiduc(:,1) = rot*(NAS'-center);
    fiduc(:,2) = rot*(LPA'-center);
    fiduc(:,3) = rot*(RPA'-center);
    mri_scan.PCSFiducial = [fiduc,[0 0 0]'];
    
    set(MAINWINDOW,'Userdata',mri_scan)
    
    mritool_cb SaveAsBrainstorm
    
end

%%----------------------------------------------------------------------------------%%
%%----------------------------------------------------------------------------------%%
%%----------------------------------------------------------------------------------%%
%%                                Subroutines
%%----------------------------------------------------------------------------------%%
%%----------------------------------------------------------------------------------%%
%%------------------------------------------------------------------------------

function visu_pt(pthandle)
% Check whether there is any selected point in the current slices
% If no : set them visible | off
% If yes : set them visible | on
% Author : Sylvain Baillet, Ph.D.
% Signal and Image Processing Institute - University of Southern California
% Los Angeles - USA
% The BrainStorm Toolbox
% October 7, 1998.

if isempty(pthandle)
    return
end
MAINWINDOW = findobj(0,'Tag','mainwindow');
mri_scan = get(MAINWINDOW ,'Userdata');
s1 = findobj(MAINWINDOW ,'Tag','1');
s2 = findobj(MAINWINDOW ,'Tag','2');
s3 = findobj(MAINWINDOW ,'Tag','3');

axa = mri_scan.axa;
axs = mri_scan.axs;
axc = mri_scan.axc;

x = get(s3,'Value');
y = get(s1,'Value');
z = get(s2,'Max')-get(s2,'Value');


obj_c = findobj(MAINWINDOW ,'Parent',axc,'Userdata',x);
if ~isempty(obj_c)
    if mri_scan.seelabels == 0
        set(findobj(obj_c,'type','text'),'Visible','on');
    end
    set(findobj(MAINWINDOW ,'Parent',axc,'type','patch'),'visible','off');
    set(findobj(obj_c,'type','patch'),'Visible','on');
    
else
    set(findobj(MAINWINDOW ,'Parent',axc,'Type','patch'),'Visible','off');
    set(findobj(MAINWINDOW ,'Parent',axc,'Type','text'),'Visible','off');
end

obj_s = findobj(MAINWINDOW ,'Parent',axs,'Userdata',y);
if ~isempty(obj_s)
    if mri_scan.seelabels == 0
        set(findobj(obj_s,'type','text'),'Visible','on');
    end
    set(findobj(MAINWINDOW ,'Parent',axs,'type','patch'),'visible','off');
    set(findobj(obj_s,'type','patch'),'Visible','on');
    
    
else
    tmp = findobj(MAINWINDOW ,'Parent',axs,'Type','patch');
    set(tmp,'Visible','off');
    set(findobj(MAINWINDOW ,'Parent',axs,'Type','text'),'Visible','off');
    
end

obj_a = findobj(MAINWINDOW ,'Parent',axa);
val = get(obj_a,'Userdata');
val = [val{:}];

i = find(abs(val-z)/z < 1e-10);
obj_a = obj_a(i);

if ~isempty(obj_a)
    
    if mri_scan.seelabels == 0
        
        set(findobj(obj_a,'type','text'),'Visible','on');
        
    end
    
    set(findobj(MAINWINDOW ,'Parent',axa,'type','patch'),'visible','off');
    
    set(findobj(obj_a,'type','patch'),'Visible','on');
    
    
    
else
    
    set(findobj(MAINWINDOW ,'Parent',axa,'Type','patch'),'Visible','off');
    
    set(findobj(MAINWINDOW ,'Parent',axa,'Type','text'),'Visible','off');
    
end


% --------------------- SUBFUNCTION TO READ BINARY DATA ------------------------
function A = get_gain(fid,ndx);

%GET_GAIN - Get a set of source foward fields from a .bin gain matrix file
% function A = get_gain(fid,ndx);
% Get a set of source foward fields A from a .bin gain matrix file
% INPUTS
%       fid : a file identifier following a fopen call on a valid binary gain matrix file
%       ndx : vector array of source indices for which the forward fields are requested
% OUTPUTS
%       A   : array of requested forward fields
% See also READ_GAIN

% \------------------------------------------------------/
% | *** Sylvain Baillet, Ph.D.                           |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   | 
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% \------------------------------------------------------/
%  
% Date of creation: October 2002 
% Script History -----------------------------------------------------------------------------------
%---------------------------------------------------------------------------------------------------


frewind(fid);

rows = fread(fid,1,'uint32');

if nargin == 2 
   cols = length(ndx);
   
   A = zeros(rows,cols);
   
   for i = 1:cols,
      % 4 bytes per element, find starting point
      offset = 4 + (ndx(i)-1)*rows*4;
      status = fseek(fid,offset,'bof');
      if(status == -1),
         error(sprintf('Error reading file at column %.0f',i));
      end
      
      A(:,i) = fread(fid,[rows,1],'float32'); 
      
   end
   
   return
   
else % ndx is not specified: read the whole matrix
   
   offset = 4;
   status = fseek(fid,offset,'bof');
   fbegin = ftell(fid);
   status = fseek(fid,0,'eof');
   fend = ftell(fid);
   
   frewind(fid);
   offset = 4;
   status = fseek(fid,offset,'bof');
   
   cols = (fend - fbegin + 1) / (4* rows); % Number of sources / columns
   
   A = fread(fid,[rows, cols],'float32');
   
end
