include IVIM code from Dr. Amita Shukla-Dave Lab at MSKCC

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/IVIM_standard_bcin.m

@@ -0,0 +1,113 @@
+function [f_arr, D_arr, Dx_arr, s0_arr, fitted_dwi_arr, RSS, rms_val, chi, AIC, BIC, R_sq] = IVIM_standard_bcin(dwi_arr,bval_arr,sigmadwi,LB0,UB0,x0in,parallelFlag,previewMode,previewAxes)
+
+
+
+if size(bval_arr,2) < size(bval_arr,1)
+    bval_arr = bval_arr';
+end
+
+numVoxels = size(dwi_arr,2);
+numBvals = length(bval_arr);
+
+f_lb = 0;
+% f_ub = 0.5;
+f_ub = 1.0;
+% f_ub = 0.75;
+% D_lb = 0;
+D_lb = 1e-6;
+D_ub = 0.003;
+% D_ub = 4e-3;
+% Dx_lb = 0;
+Dx_lb = 1e-6;
+Dx_ub = 5e-2;
+% Dx_ub = 0.5;
+% Dx_ub = 300e-3;
+s0_lb = 0;
+
+% %Yousef bounds Jul 16 2018
+d10 = 1e-3;
+d20 = 1e-2;
+f0 = 0.2;
+
+% p0(i,:) = [    smax(i)       d10        d20    f0];
+% lb(i,:) = [0.5*smax(i)  0.25*d10   0.25*d20   0.0];
+% ub(i,:) = [2.0*smax(i)   4.0*d10    4.0*d20   1.0];
+
+if nargin < 6
+    x0in = [f0 d10 d20 1];  % entry for s (4th) is multiplicative factor, not absolute
+    if nargin < 4
+        %     LB = [f_lb D_lb Dx_lb s0_lb];
+        %     UB = [f_ub D_ub Dx_ub 0];
+        % entry for s (4th) is multiplicative factor, not absolute
+        LB0 = [0 0.25*d10 0.25*d20 0.5];
+        UB0 = [1 4*d10 4*d20 2];
+    end
+end
+
+optimization_iterations = 4;
+
+options = optimset('MaxFunEvals',400,'MaxIter',200, ...
+    'TolFun',1e-6,'TolX',1e-6,'Display','off');
+
+axisLabels = {'b-value (s/mm^2)','signal (a.u.)'};
+
+tic
+for i = 1:numVoxels
+    if rem(i,1000) == 0
+        toc
+        disp(['Voxel ' num2str(i) ' of ' num2str(numVoxels)]);
+        tic
+    end
+    s = dwi_arr(:,i);
+
+    LB = LB0;
+    UB = UB0;
+    x0 = x0in;
+
+    LB(end) = LB0(end)*max(s);
+    UB(end) = UB0(end)*max(s);
+    x0(end) = x0in(end)*max(s);
+    %     disp(['Voxel ' num2str(i) ' of ' num2str(numVoxels)]);
+    if ~parallelFlag
+        for j = 1:optimization_iterations
+            %         x0 = (LB + UB) / 2; %+ (UB - LB) * (rand - 0.5);
+            [x_fit(j,:),resnorm(j),~,~,~] = lsqcurvefit(@(x,b)ivim_modeling_noise(x,b,sigmadwi), ...
+                x0, bval_arr', s, LB, UB, options);
+        end
+    else
+        parfor j = 1:optimization_iterations
+            %         x0 = (LB + UB) / 2; %+ (UB - LB) * (rand - 0.5);
+            [x_fit(j,:),resnorm(j),~,~,~] = lsqcurvefit(@(x,b)ivim_modeling_noise(x,b,sigmadwi), ...
+                x0, bval_arr', s, LB, UB, options);
+        end
+    end
+    [~,best_fit_idx] = min(resnorm);
+    p = x_fit(best_fit_idx,:);
+
+    fitted_dwi_arr(:,i) = ivim_modeling_noise(p,bval_arr',sigmadwi);
+    f_arr(i) = p(1);
+    D_arr(i) = p(2)*(1e3);
+    if p(3) > 1e-5
+        Dx_arr(i) = p(3)*(1e3);
+    else
+        Dx_arr(i) = 1e-5;
+    end
+    s0_arr(i) = p(4);
+
+    if previewMode == 1 || previewMode == 100
+        axisLabels{3} = ['IVIM Fit (Voxel # ' num2str(i) '/' num2str(numVoxels) ')'];
+        if previewMode == 100 && i < 101
+            updatePreviewAxes(previewAxes,bval_arr',s,fitted_dwi_arr(:,i),axisLabels);
+        elseif previewMode == 1
+            updatePreviewAxes(previewAxes,bval_arr',s,fitted_dwi_arr(:,i),axisLabels);
+        end
+     end
+
+end
+
+% Quantification of Fitting Quality
+Nb = numel(bval_arr);
+Np = numel(p);
+observed_arr = dwi_arr;
+predicted_arr = fitted_dwi_arr;
+[RSS,rms_val,chi,AIC,BIC,R_sq] = qualityFit(Nb,Np,observed_arr,predicted_arr);

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/README.asv

@@ -0,0 +1,33 @@
+# MRI-QAMPER_IVIM
+
+This is the codebase for intravoxel incoherent motion (IVIM) code from the MRI-QAMPER MATLAB package, developed by Dr. Shukla-Dave's lab at Memorial Sloan Kettering Cancer Center
+
+Authors: Eve LoCastro ([email protected]), Dr. Ramesh Paudyal ([email protected]), Dr. Amita Shukla-Dave ([email protected])
+Institution: Memorial Sloan Kettering Cancer Center
+Department: Medical Physics
+Address: 321 E 61st St, New York, NY 10022
+
+The codebase and subdirectories should be added to the MATLAB path. An example usage script is provided in `demo_QAMPER_IVIM.m`.
+
+```
+% This is an example usage script for MSK Medical Physics Dave Lab QAMPER IVIM
+% Please replace the variable names below with path values for 
+%       qamper_path: path to MRI-QAMPER_IVIM folder
+%       img_nii: multi-b value DWI (4-D NIfTI)
+%       bval_file: b-value information (txt file)
+%       roi_nii: single-volume mask ROI image (NIfTI)
+
+
+
+qamper_path = 'path:\to\MRI-QAMPER_IVIM';
+addpath(genpath(qamper_path));
+
+img_nii = 'dwi.nii';
+bval_file = fullfile(qamper_path,'test_data','702-HN401D-D2019_10_08.bval');
+roi_nii = fullfile(qamper_path,'test_data','702-HN401D-D2019_10_08_BD_REDO_SV.nii.gz');
+
+save_folder = fullfile(qamper_path,'test_data');
+
+batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);
+
+```

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/README.md

@@ -0,0 +1,30 @@
+# MRI-QAMPER_IVIM
+
+This is the codebase for intravoxel incoherent motion (IVIM) code from the MRI-QAMPER MATLAB package, developed by Dr. Shukla-Dave's lab at Memorial Sloan Kettering Cancer Center.
+
+__Authors:__ Eve LoCastro ([email protected]), Dr. Ramesh Paudyal ([email protected]), Dr. Amita Shukla-Dave ([email protected]) </br>
+__Institution:__ Memorial Sloan Kettering Cancer Center </br>
+__Department:__ Medical Physics</br>
+__Address:__ 321 E 61st St, New York, NY 10022
+
+The codebase and subdirectories should be added to the MATLAB path. An example usage script is provided in `demo_QAMPER_IVIM.m`.
+
+```
+% This is an example usage script for MSK Medical Physics Dave Lab QAMPER IVIM
+% Please replace the variable names below with path values for 
+%       qamper_path: path to MRI-QAMPER_IVIM folder
+%       img_nii: multi-b value DWI (4-D NIfTI)
+%       bval_file: b-value information (txt file)
+%       roi_nii: single-volume mask ROI image (NIfTI)
+
+qamper_path = 'path:\to\MRI-QAMPER_IVIM';
+addpath(genpath(qamper_path));
+
+img_nii = 'dwi.nii';
+bval_file = 'dwi.bval';
+roi_nii = 'roi.nii';
+
+save_folder = fullfile(qamper_path,'test_data');
+
+batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);
+```

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/avgRepeatBvalVols.m

@@ -0,0 +1,16 @@
+function [uniq_bval_arr, dwi_img_new] = avgRepeatBvalVols(bval_arr,dwi_img)
+
+disp('Averaging multiple bval vols to increase signal')
+
+[uniq_bval_arr,IA,~] = unique(bval_arr);
+
+dwi_img_new = [];
+
+for i = 1:numel(uniq_bval_arr)
+    if i < numel(uniq_bval_arr)
+        dwi_avg_vol = mean(dwi_img(:,:,:,IA(i):IA(i+1)-1),4);
+    else
+        dwi_avg_vol = mean(dwi_img(:,:,:,IA(i):end),4);
+    end
+    dwi_img_new = cat(4,dwi_img_new, dwi_avg_vol); 
+end

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/buildBatchFromCSV.m

@@ -0,0 +1,107 @@
+function batchStruct = buildBatchFromCSV(img_nii,roi_nii,qamper_mode,optionS)
+
+batchStruct = [];
+batchStruct.batch_id = num2str(now);
+
+if strcmp(qamper_mode,'DWI')
+    if isfield(optionS,'DWIstruct')
+        DWIstruct = optionS.DWIstruct;
+    end
+
+    DWIstruct.runDWI = 1;
+    [droot,fbase,ext] = fileparts(img_nii);
+    batchStruct.currentNiiPath = fullfile(droot,'dwi_maps');
+
+
+    if strcmp(ext,'.gz')
+        [~,fbase,~] = fileparts(fbase);
+    end
+
+    bval_mat = fullfile(droot,[fbase '.mat']);
+
+    T = load(bval_mat);
+
+    DWIstruct.bval_arr  = T.bval_arr;
+
+    img = nii_load(img_nii,1);
+    vol_num = size(img.img,4);
+
+    DWIstruct.files{1} = img_nii; DWIstruct.vols = vol_num;
+    DWIstruct.files{2} = roi_nii;
+
+    if ~isfield(DWIstruct,'ADC')
+        DWIstruct.ADC = 0;
+    end
+    if ~isfield(DWIstruct,'IVIM')
+        DWIstruct.IVIM = 0;
+    end
+    if ~isfield(DWIstruct,'NGIVIM')
+        DWIstruct.NGIVIM = 0;
+    end
+    if ~isfield(DWIstruct,'Kurtosis')
+        DWIstruct.Kurtosis = 0;
+    end
+
+    if ~isfield(DWIstruct,'roiType')
+        DWIstruct.roiType = 'node';
+    end
+
+    if ~isfield(DWIstruct,'avgRepeatBvals')
+        DWIstruct.avgRepeatBvals = 1;
+    end
+
+    if ~isfield(DWIstruct,'channels')
+        DWIstruct.channels = 16; %
+    end
+
+    if ~isfield(DWIstruct, 'avgChannels')
+        DWIstruct.avgChannels = 2; %
+    end
+
+    if DWIstruct.IVIM
+        if ~isfield(DWIstruct,'LB_IVIM')
+            IVIMBoundsTable = [{0} {0} {0} {0}; {0.5} {0.005} {0.5} {1.0}; {0.05} {0.0005} {0.01} {0.5}];
+
+            DWIstruct.LB_IVIM = cell2mat(IVIMBoundsTable(1,:));
+            DWIstruct.UB_IVIM = cell2mat(IVIMBoundsTable(2,:));
+            DWIstruct.x0_IVIM = cell2mat(IVIMBoundsTable(3,:));
+        end
+    end
+    if DWIstruct.NGIVIM
+        if ~isfield(DWIstruct,'LB_NG')
+            NGIVIMDefaultBoundsTable = [{0} {0} {0} {0} {0}; {0.5} {0.005} {0.1} {1.0} {2.0}; {0.05} {0.0005} {0.01} {1.0} {0.5}];
+
+            DWIstruct.LB_NG = cell2mat(NGIVIMDefaultBoundsTable(1,:));
+            DWIstruct.UB_NG = cell2mat(NGIVIMDefaultBoundsTable(2,:));
+            DWIstruct.x0_NG = cell2mat(NGIVIMDefaultBoundsTable(3,:));
+        end
+    end    
+
+    if ~isfield(DWIstruct,'Parallel')
+        DWIstruct.Parallel = 1;
+    end
+
+    DWIstruct.bval_tf = ones(size(DWIstruct.bval_arr));
+%     
+%     DWIstruct.previewAxes = handles.fittingAxes;
+%     
+    DWIstruct.nanopt = 1;
+
+    if ~isfield(DWIstruct,'kernelsize')
+        DWIstruct.kernelsize = 1;
+    end
+    if ~isfield(DWIstruct,'smoothopt')
+        DWIstruct.smoothopt = 1;
+    end
+
+    DWIstruct.previewMode = 0;
+    DWIstruct.previewAxes = 0;
+
+else
+    DWIstruct.runDWI = 0;
+end
+
+batchStruct.DWIstruct = DWIstruct;
+batchStruct.T1struct.runT1 = 0;
+batchStruct.DCEstruct.runDCE = 0;
+batchStruct.T2struct.runT2 = 0;