use fit instead of lsqcurvefit in AIFbiexpcon + weighting support, mi…

…sc fixes

.gitignore

@@ -179,4 +179,5 @@ email_preferences.txt
 dce_last_state.mat
 *info.log
 dce/*info.log
+test_data/BIDS_test/derivatives
 

dce/AIFbiexpcon.m

@@ -3,28 +3,27 @@
 %% This models the AIF as a biexponential function convolved with a
 %% rectangular function to model the injection
 
-function Cp = AIFbiexpcon(x, xdata)%A, B, c, d, T1, step)
+% function Cp = AIFbiexpcon(x, xdata)
+function Cp = AIFbiexpcon(A, B, c, d, T1, step, fittingAU, baseline)
 
-A = x(1);
-B = x(2);
-c = x(3);
-d = x(4);
-T1 = xdata{1}.timer;
-step=xdata{1}.step;
-maxer=xdata{1}.maxer;
+% A = x(1);
+% B = x(2);
+% c = x(3);
+% d = x(4);
+% T1 = xdata{1}.timer;
+% step=xdata{1}.step;
+% maxer=xdata{1}.maxer;
 
 T1 = T1(:);
 OUT = find(step > 0);
 
-
-
 %Set as model A, MacGrath, MRM 2009
-
 %B = maxer-A;
 
-if isfield(xdata{1}, 'raw') && xdata{1}.raw == true && isfield(xdata{1}, 'baseline')
-    baseline = xdata{1}.baseline;
-else
+% if exist('xdata', 'var') && isfield(xdata{1}, 'raw') && xdata{1}.raw == true && isfield(xdata{1}, 'baseline')
+if ~fittingAU
+%     baseline = baseline;
+% else
     baseline = 0;
 end
 
@@ -40,3 +39,4 @@
         Cp(j) = A.*(exp(-c.*(T1(j) - T1(OUT(end))))) + B.*(exp(-d.*(T1(j) - T1(OUT(end)))));
     end
 end
+Cp = Cp';

dce/AIFbiexpfithelp.m

@@ -58,14 +58,25 @@
 options = optimset('lsqcurvefit');
 
 %increase the number of function evaluations for more accuracy
-options.MaxFunEvals = MaxFunEvals;
-options.MaxIter     = MaxIter;
-options.TolFun      = TolFun;
-options.TolX        = TolX;
-options.Diagnostics = 'off';
-options.Display     = 'off';
-options.Algorithm   = 'levenberg-marquardt';
-options.Robust      = Robust;
+% options.MaxFunEvals = MaxFunEvals;
+% options.MaxIter     = MaxIter;
+% options.TolFun      = TolFun;
+% options.TolX        = TolX;
+% options.Diagnostics = 'off';
+% options.Display     = 'off';
+% options.Algorithm   = 'levenberg-marquardt';
+% options.Robust      = Robust;
+options = fitoptions('Method', 'NonlinearLeastSquares',...
+    'Algorithm', 'Levenberg-Marquardt',...
+    'MaxIter', MaxIter,...
+    'MaxFunEvals', MaxFunEvals,...
+    'TolFun', TolFun,...
+    'TolX', TolX,...
+    'Display', 'off',...
+    'Lower',lower_limits,...
+    'Upper', upper_limits,...
+    'StartPoint', initial_values,...
+    'Robust', Robust);
 
 % Choose upper and lower bounds only for trust-region methods.
 % lb = [0 0 0 0];
@@ -97,15 +108,19 @@
 step(start_index:end_index) = 1;
 xdata{1}.step = step;
 
-% W is the weighting matrix, should you want to emphasise certain
-% datapoints
-W = ones(size(Cp));
+
 [~, max_index] = max(Cp.*step);
 % WW= sort(Cp.*step, 'descend');
 % ind(1) = find(Cp == WW(1));
 % ind(2) = find(Cp == WW(2));
 % ind(3) = find(Cp == WW(3));
 
+% W is the weighting matrix, should you want to emphasise certain
+% datapoints
+W = ones(size(Cp)) .* 10;
+W(1:max_index) = 0;
+options.Weights = W;
+
 step(max_index+1:end) = 0;
 xdata{1}.step = step;
 
@@ -133,7 +148,7 @@
 
 maxer = Cp(max_index);
 xdata{1}.maxer = maxer;
-Cp = Cp.*W;
+% Cp = Cp.*W;
 
 end_baseline = find(xdata{1}.step > 0);
 baseline = mean(Cp(1:end_baseline(1)));
@@ -142,24 +157,31 @@
 upper_limits(2) = maxer*2;
 initial_values(1) = maxer*0.5;
 initial_values(2) = maxer*0.5;
+options.Upper = upper_limits;
+options.StartPoint = initial_values;
 if verbose>0
     disp('Fitting AIF values, limits and initial values adjusted');
     fprintf('lower_limits = %s\n',num2str(lower_limits));
     fprintf('upper_limits = %s\n',num2str(upper_limits));
     fprintf('initial_values = %s\n',num2str(initial_values));
 end
 % Currently, we use AIF
-[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(@AIFbiexpcon, ...
-    initial_values, xdata, ...
-    Cp',lower_limits,upper_limits,options);
+% [x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(@AIFbiexpcon, ...
+%     initial_values, xdata, ...
+%     Cp',lower_limits,upper_limits,options);
+
+ft = fittype('AIFbiexpcon(A, B, c, d, T1, step, fittingAU, baseline)', ...
+    'independent', {'T1', 'step'}, 'coefficients', {'A', 'B', 'c', 'd'}, 'problem', {'fittingAU', 'baseline'});
+[f, gof, output] = fit([timer, step], Cp, ft, options, 'problem', {xdata{1}.fittingAU, xdata{1}.baseline});
 
 xdata{1}.timer = oldt;
-rsquare = 1 - resnorm / norm(Cp-mean(Cp))^2;
+% rsquare = 1 - resnorm / norm(Cp-mean(Cp))^2;
 if verbose>0
-    disp(['R^2 of AIF fit = ' num2str(rsquare)]);
+    disp(['Adjusted R^2 of AIF fit = ' num2str(gof.adjrsquare)]);
+    % disp(['Adjusted R^2 of AIF fit = ' num2str(rsquare)]);
 end
 
-out = AIFbiexpcon(x, xdata);
-
-
-
+% out = AIFbiexpcon(x, xdata);
+out = AIFbiexpcon(f.A, f.B, f.c, f.d, timer, step, xdata{1}.fittingAU, xdata{1}.baseline)';
+x = [f.A, f.B, f.c, f.d];
+rsquare = gof.adjrsquare;

dce/B_AIF_fitting_func.m

@@ -230,6 +230,7 @@
 aif_name = '';
 if isempty(import_aif_path)
     if(fit_aif==1)
+        xdata{1}.fittingAU = false;
         [Cp_fitted xAIF xdataAIF] = AIFbiexpfithelp(xdata, 1);
         Cp_use = Cp_fitted;
 
@@ -239,6 +240,7 @@
         %Fit raw data curve
         Cptemp = xdata{1}.Cp;
         xdata{1}.Cp = xdata{1}.Stlv;
+        xdata{1}.fittingAU = true;
         [Stlv_fitted, ~, ~] = AIFbiexpfithelp(xdata, 1);
         xdata{1}.Cp = Cptemp;
         Stlv_use = Stlv_fitted;

dce/dce_auto_aif.m

@@ -1,7 +1,4 @@
 function [end_ss, aif_index, dynamic_output] = dce_auto_aif(dynamic_input, mask_index,dimx,dimy,dimz,injection_duration)
-% load('C:\Users\sbarnes\Documents\data\6 DCE Stroke\aging\Sam Analysis\Raw AIF\1 young - Copy\A-preAIFcalc.mat')
-% load('C:\Users\sbarnes\Documents\data\6 DCE Stroke\aging\Sam Analysis\Fitted Auto AIF\1 young\A-preAIFcalcPostDrift.mat')
-% load('C:\Users\sbarnes\Documents\data\6 DCE Stroke\aging\Sam Analysis\Fitted Auto AIF\2 old\A-preAutoFitWithDrift.mat')
 
 % Only calculate the end of steady state, stop here
 if nargout == 1

dce/dce_preferences.txt

@@ -26,12 +26,12 @@ autoaif_sobel_threshold = 0.005
 % Initial Upper limits and initial values for A and B are auto-set
 % values order = A B c d
 aif_lower_limits = 0 0 0 0
-aif_upper_limits = 5 5 5 5
-aif_initial_values = 1 1 1 1
+aif_upper_limits = 5 5 50 50
+aif_initial_values = 1 1 1 0.01
 
 % Advanced Options
 aif_TolFun = 10^(-20)
-aif_TolX = 10^(-20)
+aif_TolX = 10^(-23)
 aif_MaxIter = 1000
 aif_MaxFunEvals = 1000
 %robust options no quotes: "off" "LAR" "Bisquare"