Merge 7316f4d into cd30473

External/imtool3D_td/src/StatsGUI.m

@@ -7,7 +7,7 @@
     Nvol = size(I,5);
 end
 if ~exist('fields','var') || isempty(fields)
-    fields = cellfun(@(x) ['vol #' x],strsplit(num2str(Nvol:-1:1)),'uni',0); 
+    fields = cellfun(@(x) ['vol #' x],strsplit(num2str(1:Nvol)),'uni',0); 
 else
     fields = cellfun(@(X) X(max(1,end-30):end),fields,'uni',0);
 end
@@ -22,9 +22,9 @@
     Stats = cell(length(fields),7);
     for iii=1:length(fields)
         if iscell(I)
-            datiii = nanmean(I{length(fields)+1-iii},4); % average 4D data along time
+            datiii = nanmean(I{iii},4); % average 4D data along time
         else
-            datiii = nanmean(I(:,:,:,:,length(fields)+1-iii),4);
+            datiii = nanmean(I(:,:,:,:,iii),4);
         end
         datiii = datiii(Mask);
         Stats{iii,1} = sum(Mask(:));

External/imtool3D_td/src/nii_get_orient.m

@@ -18,6 +18,8 @@
     useForm='s';
 elseif hdr.qform_code > 0
     useForm='q';
+else
+    useForm=[];
 end
 
 affine_transform = 1;

External/imtool3D_td/src/nii_load.m

@@ -62,7 +62,11 @@
 
     % LOAD AND RESLICE
     hdr_ref = nii_reset_orient(list{1}); % Set hdr.dim back to initial orientation
-    nii = nii_xform(list{ff},hdr_ref,[],intrp);
+    if ff==1
+        nii = nii_tool('load',list{ff});
+    else
+        nii = nii_xform(list{ff},hdr_ref,[],intrp);
+    end
 
     if nargin==1 || (~isempty(untouch) && ~untouch)
         orient = nii_get_orient(nii.hdr);

qMRLab.m

@@ -212,8 +212,16 @@ function MethodMenu(hObject, eventdata, handles, Method)
 if isappdata(0,'Model') && strcmp(class(getappdata(0,'Model')),Method) % if same method, load the current class with parameters
     Model = getappdata(0,'Model');
 else % otherwise create a new object of this method
-    modelfun  = str2func(Method);
-    Model = modelfun();
+    Modeltobesaved = getappdata(0,'Model');
+    savedModel = getappdata(0,'savedModel');
+    savedModel.(class(Modeltobesaved)) = Modeltobesaved;
+    setappdata(0,'savedModel',savedModel);
+    if isfield(savedModel,Method) && ~isempty(savedModel.(Method))
+        Model = savedModel.(Method);
+    else
+        modelfun  = str2func(Method);
+        Model = modelfun();
+    end
 end
 SetAppData(Model)
 % Create empty Data
@@ -420,6 +428,9 @@ function FitGo_FitData(hObject, eventdata, handles)
 SetAppData(FileBrowserList);
 % Show results
 handles.CurrentData = FitResults;
+if exist('hdr','var')
+    handles.CurrentData.hdr = hdr;
+end
 guidata(hObject,handles);
 DrawPlot(handles);
 
@@ -506,7 +517,7 @@ function Stats_Callback(hObject, eventdata, handles)
 
 I = handles.tool.getImage(1);
 Iraw = handles.CurrentData;
-fields = setdiff(Iraw.fields,'Mask')';
+fields = setdiff(Iraw.fields,'Mask','stable')';
 Maskall = handles.tool.getMask(1);
 Color = handles.tool.getMaskColor;
 StatsGUI(I,Maskall, fields, Color);

src/Common/FitData.m

@@ -34,7 +34,7 @@
 tStart = tic;
 tsaved = 0;
 
-h=[];
+h=[]; hwarn=[];
 if moxunit_util_platform_is_octave % ismethod not working properly on Octave
     try, Model = Model.Precompute; end
     try, Model = Model.PrecomputeData(data); end
@@ -44,6 +44,19 @@
     if ismethod(Model,'PrecomputeData'), Model = Model.PrecomputeData(data); end
 end
 
+% NaN in Mask
+if isfield(data,'Mask') && (~isempty(data.Mask)) && any(isnan(data.Mask(:)))
+    data.Mask(isnan(data.Mask))=0;
+    msg = 'NaNs will be set to 0. We recommend you to check your mask.';
+    titlemsg = 'NaN values detected in the Mask';
+    if exist('wait','var') && (wait)
+        hwarn = warndlg(msg,titlemsg);
+    end
+    fprintf('\n')
+    warning(titlemsg)
+    fprintf('%s\n\n',msg)
+end
+
 if Model.voxelwise % process voxelwise
     %############################# INITIALIZE #################################
     % Get dimensions
@@ -75,9 +88,25 @@
     end
 
 
-    % Find voxels that are not empty
+
     if isfield(data,'Mask') && (~isempty(data.Mask))
-        Voxels = find(all(data.Mask & ~computed,2));
+
+        % Set NaN values to zero if there are any 
+        if any(isnan(data.Mask(:)))
+           data.Mask(isnan(data.Mask))=0;
+            msg = 'NaNs will be set to 0. We recommend you to check your mask.';
+            titlemsg = 'NaN values detected in the Mask';
+            if exist('wait','var') && (wait)
+                hwarn = warndlg(msg,titlemsg);
+            end
+            fprintf('\n')
+            warning(titlemsg)
+            fprintf('%s\n\n',msg) 
+        end
+
+        % Find voxels that are not empty
+        Voxels = find(all(data.Mask & ~computed,2));    
+
     else
         Voxels = find(~computed)';
     end
@@ -212,6 +241,7 @@
 end
 % delete waitbar
 %if (~isempty(hMSG) && not(isdeployed));  delete(hMSG); end
+if ishandle(hwarn), delete(hwarn); end
 
 Fit.Time = toc(tStart);
 Fit.Protocol = Model.Prot;

src/Common/GUI/Custom_OptionsGUI.m

@@ -316,7 +316,9 @@ function OptionsGUI_OpeningFcn(hObject, eventdata, handles, varargin)
             % Move down
             uicontrol(handles.(fields{ii}).panel,'Units','normalized','Position',[.53 0.02*N .44 .02*N],'Style','pushbutton','String','Move down','Callback',@(hObject, eventdata) PointDown_Callback(hObject, eventdata, handles,fields{ii}));
             % LOAD
-            uicontrol(handles.(fields{ii}).panel,'Units','normalized','Position',[.03 0    .94 .02*N],'Style','pushbutton','String','Load','Callback',@(hObject, eventdata) LoadProt_Callback(hObject, eventdata, handles,fields{ii}));
+            uicontrol(handles.(fields{ii}).panel,'Units','normalized','Position',[.03 0      .44 .02*N],'Style','pushbutton','String','Load','Callback',@(hObject, eventdata) LoadProt_Callback(hObject, eventdata, handles,fields{ii}));
+            % Create
+            uicontrol(handles.(fields{ii}).panel,'Units','normalized','Position',[.53 0      .44 .02*N],'Style','pushbutton','String','Create','Callback',@(hObject, eventdata) CreateProt_Callback(hObject, eventdata, handles,fields{ii}));
         end
 
     end
@@ -470,16 +472,42 @@ function DefaultProt_Callback(hObject, eventdata, handles)
 OptionsGUI_OpeningFcn(hObject, eventdata, handles, Model, handles.caller)
 
 function LoadProt_Callback(hObject, eventdata, handles, MRIinput)
-[FileName,PathName] = uigetfile({'*.mat;*.xls;*.xlsx;*.txt;*.scheme'},'Load Protocol Matrix');
+FileFormat = '*.mat;*.xls;*.xlsx;*.txt';
+if strcmp(MRIinput,'DiffusionData')
+    FileFormat = ['*.bvec;*.scheme;' FileFormat];
+end
+[FileName,PathName] = uigetfile({FileFormat},'Load Protocol Matrix');
 if PathName == 0, return; end
 fullfilepath = [PathName, FileName];
 Prot = ProtLoad(fullfilepath);
+if Prot == 0, return; end
 if ~isnumeric(Prot), errordlg('Invalid protocol file'); return; end
 set(handles.(MRIinput).table,'Data',Prot)
 Model = getappdata(0,'Model');
 Model.Prot.(MRIinput).Mat = Prot;
 UpdateProt(MRIinput,Prot,handles)
 
+function CreateProt_Callback(hObject, eventdata, handles, MRIinput)
+Model = getappdata(0,'Model');
+Fmt = Model.Prot.(MRIinput).Format; if ischar(Fmt), Fmt = {Fmt}; end
+answer = inputdlg(Fmt,'Enter values, vectors or Matlab expressions',[1 100]);
+if isempty(answer), return; end
+Prot = cellfun(@str2num,answer,'uni',0);
+Prot = cellfun(@(x) x(:),Prot,'uni',0);
+Nlines = max(cell2mat(cellfun(@length,Prot,'uni',0)));
+Model.Prot.(MRIinput).Mat = NaN(Nlines,length(Fmt));
+for ic = 1:length(Fmt)
+    if length(Prot{ic})>1
+        Lmax = length(Prot{ic}); % if vector, fill as many as possible
+    else
+        Lmax = Nlines; % if scalar, all lines get this value
+    end
+    Model.Prot.(MRIinput).Mat(1:Lmax,ic) = Prot{ic};
+end
+Prot = Model.Prot.(MRIinput).Mat;
+set(handles.(MRIinput).table,'Data',Prot)
+UpdateProt(MRIinput,Prot,handles)
+
 
 function UpdateProt(MRIinput,Prot,handles)
 if ~isnumeric(Prot)
@@ -530,10 +558,21 @@ function Default_Callback(hObject, eventdata, handles)
 oldModel = getappdata(0,'Model');
 modelfun = str2func(class(oldModel));
 Model = modelfun();
-Model.Prot = oldModel.Prot;
-setappdata(0,'Model',Model);
+
+answer = questdlg('What do you want to set to default?','Reset protocol?','Reset options','Reset options AND protocol','Reset protocol','Reset options');
+if strfind(answer,'options')
+    newModel = Model;
+    newModel.Prot = oldModel.Prot;
+else
+    newModel = oldModel;
+end
+if strfind(answer,'protocol')
+    newModel.Prot = Model.Prot;
+end
+
+setappdata(0,'Model',newModel);
 set(handles.ParametersFileName,'String','Parameters Filename');
-OptionsGUI_OpeningFcn(hObject, eventdata, handles, Model, handles.caller)
+OptionsGUI_OpeningFcn(hObject, eventdata, handles, newModel, handles.caller)
 
 % --- Executes on button press in Load.
 function Load_Callback(hObject, eventdata, handles)

src/Common/tools/GUIfun/DrawPlot.m

@@ -38,8 +38,9 @@ function DrawPlot(handles,CurrentName)
 
 % Change save as NIFTI function
 H = handles.tool.getHandles;
-set(H.Tools.maskSave,'Callback',@(hObject,evnt)saveMask(handles.tool,hObject,handles.CurrentData.hdr))
-set(H.Tools.maskLoad,'Callback',@(hObject,evnt)loadMask(handles.tool,hObject,handles.CurrentData.hdr))
+if isfield(handles.CurrentData,'hdr'), hdr = {handles.CurrentData.hdr}; else, hdr = {}; end
+set(H.Tools.maskSave,'Callback',@(hObject,evnt)saveMask(handles.tool,hObject,hdr{:}))
+set(H.Tools.maskLoad,'Callback',@(hObject,evnt)loadMask(handles.tool,hObject,hdr{:}))
 
 % Use Shortcut to Source button
 set(findobj('Name','qMRLab'),'Windowkeypressfcn', @(hobject, event) shortcutCallback(hobject, event,handles))

src/Common/tools/ProtLoad.m

@@ -1,6 +1,6 @@
 function Prot = ProtLoad(fullfilepath)
 [~,~,ext] = fileparts(fullfilepath);
-filterindex = find(cellfun(@isempty,strfind({'*.mat';'*.xls;*.xlsx';'*.txt;*.scheme'},ext))==0);
+filterindex = find(cellfun(@isempty,strfind({'*.mat';'*.xls;*.xlsx';'*.txt;*.scheme';'*.bvec'},ext))==0);
 switch filterindex
     case 1
         Prot = load(fullfilepath);
@@ -10,4 +10,21 @@
         Prot = xlsread(fullfilepath);
     case 3
         Prot = txt2mat(fullfilepath);
+    case 4 % manage bvec/bvals
+        bvecfile = fullfilepath;
+        [FileName,PathName] = uigetfile({'*.bval;*.txt'},'Load bval',strrep(fullfilepath,'.bvec','.bval'));
+        if PathName == 0, Prot=0; return; end
+        bvalfile = fullfile(PathName,FileName);
+        Model = getappdata(0,'Model');
+        if ismember(length(Model.Prot.DiffusionData.Format), [7 8])
+            fprintf(['Assuming max gradient of 40mT/m...\n'])
+            Prot = scd_schemefile_FSLconvert(bvalfile, bvecfile, 40);
+            Prot(:,4) = Prot(:,4)*1e3;
+            Prot(:,5:end) = Prot(:,5:end)*1e-3;
+        else
+            bvec = txt2mat(bvecfile); bvec = bvec';
+            bval = txt2mat(bvalfile); bval = bval(:);
+            Prot = cat(2,bvec,bval);            
+        end
+
 end

src/Models/Diffusion/dti.m

@@ -112,16 +112,16 @@
             Prot(Prot(:,4)==0,1:6) = 0;
             [~,c,ind] = consolidator(Prot(:,1:7),[],'count');
             cmax = max(c); % find images repeated more than 5 times (for relevant STD)
+            if ~strcmp(obj.options.fittingtype,'linear') && ~strcmp(obj.options.Riciannoisebias_Method,'fix sigma')
             if cmax<2
-                warndlg({'Your dataset doesn''t have 2 repeated measures (same bvec/bvals) --> you can''t estimate noise STD voxel-wise. Specify a fixed Sigma Noise in the option panel instead.'  'See Methods Noise/NoiseLevel.m to estimate the noise standard deviation.'},'Noise estimation method')
-                obj.options.Riciannoisebias_Method = 'fix sigma';
+                warndlg({'Your dataset doesn''t have 2 repeated measures (same bvec/bvals) --> you can''t estimate noise STD per voxel.' 'Using linear fitting instead (noise bias will impact your results).'},'Noise estimation method')
+                obj.options.fittingtype = 'linear';
             elseif cmax<4
-                warndlg({'Your dataset doesn''t have 4 repeated measures (same bvec/bvals) --> you can''t estimate noise STD voxel-wise accurately. Specify a fixed Sigma Noise in the option panel instead.'  'See Methods Noise/NoiseLevel.m to estimate the noise standard deviation.'},'Noise estimation method')
+                warndlg({'Your dataset doesn''t have 4 repeated measures (same bvec/bvals) --> you can''t estimate noise STD voxel-wise accurately.'},'Noise estimation method')
+            end
             end
             if strcmp(obj.options.Riciannoisebias_Method,'Compute Sigma per voxel')
                 obj.options.Riciannoisebias_value  = 'auto';
-            elseif isempty(obj.options.Riciannoisebias_value)
-                obj.options.Riciannoisebias_value=10;
             end
 
             % disable Rician noise panel if linear
@@ -179,6 +179,7 @@
                 else
                     SigmaNoise = obj.options.Riciannoisebias_value;
                 end
+                if isempty(SigmaNoise), SigmaNoise = max(data.DiffusionData/100); end
 
 
                 if ~moxunit_util_platform_is_octave && SigmaNoise

src/Models_Functions/DTIfun/scd_model_dti.m

@@ -1,6 +1,6 @@
 function [D,L,fiberdirection] = scd_model_dti(data,scheme)
 
-b0_index = ~scheme(:,4)>0;
+b0_index = scheme(:,4)<5e-6; % 5mT/m or less are considered as b=0
 scheme = scheme(~b0_index,:);
 data = data(~b0_index);
 bvec = scheme(:,[1 2 3]);

src/Models_Functions/Diffusion/scd_schemefile_FSLconvert.m

@@ -22,12 +22,12 @@
 % load bvec
 bvecs = txt2mat(bvecfile);
 
-disp(['\nCheck if bvecs is written by lines or by columns...'])
+fprintf(['\nCheck if bvecs is written by lines or by columns...\n'])
 if size(bvecs,1)==3 && size(bvecs,2)~=3
-    disp(['\nbvecs seems to be written in columns, transposing it...'])
+    fprintf(['\nbvecs seems to be written in columns, transposing it...\n'])
     bvecs=bvecs';
 else
-	disp(['.. OK! written by lines'])
+	fprintf(['\n.. OK! written by lines\n'])
 end
 
 % maximum b-value in the s/mm^2 unit