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

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/.gitignore

@@ -0,0 +1 @@
+*.asv

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.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 (locastre@mskcc.org), Dr. Ramesh Paudyal (paudyalr@mskcc.org), Dr. Amita Shukla-Dave (davea@mskcc.org) </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/bvalOrderFix.m

@@ -0,0 +1,15 @@
+function [bval_arr_new, dwi_img_new] = bvalOrderFix(bval_arr,dwi_img)
+
+disp(['Reordering B0 to start of imaging sequence'])
+
+if size(bval_arr,1) > size(bval_arr,2)
+    bval_arr = bval_arr';
+end
+
+% bval_arr_new = [bval_arr(end) bval_arr(1:end-1)];
+[bval_arr_new,idx] = sort(bval_arr);
+
+% dwi_img_new = cat(4,dwi_img(:,:,:,end),dwi_img(:,:,:,1:end-1));
+dwi_img_new = dwi_img(:,:,:,idx);
+
+disp(bval_arr_new);

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

@@ -0,0 +1,30 @@
+function [dwi_arr, bval_arr, sigmadwi, SNR, noise] = dataPrepDWI(dwi_img,bval_arr,ROI,channels,N_avgs)
+
+bgMask = zeros(size(dwi_img,1),size(dwi_img,2));
+bgMask(10:30,10:30) = 1;
+% bgidx = sub2ind(size(bgMask),10:30,10:30);
+bgidx = find(bgMask);
+
+maskidx = find(ROI);
+
+dwi_arr = map2arr(dwi_img,ROI);
+dwi_arr(:,end+1) = mean(dwi_arr,2);
+
+noise_square = bgMask.*dwi_img(:,:,1,1);
+noise_arr = noise_square(bgidx);noise_arr;
+
+%n = 8; % need to update with new protocol
+%navg = 4;
+n = channels;
+navg = N_avgs;
+
+% if data were converted to nifti before conv_panel started to extract
+% channels and N_avgs, then those numbers default to this
+if isempty(n)
+    n = 16;
+    navg = 2;
+end
+
+sigmadwi = sqrt(navg*norm(noise_arr(:))^2/length(noise_arr(:)))/sqrt(2*n);
+noise = mean(noise_arr);
+SNR = mean(dwi_arr(1,:))/sigmadwi;

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

@@ -0,0 +1,30 @@
+% Authors: Eve LoCastro (locastre@mskcc.org), Dr. Ramesh Paudyal (paudyalr@mskcc.org), Dr. Amita Shukla-Dave (davea@mskcc.org)
+% Institution: Memorial Sloan Kettering Cancer Center
+% Address: 321 E 61st St, New York, NY 10022
+
+% 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'; %path where the MRI-QAMPER folder is located
+addpath(genpath(qamper_path));
+
+% Uncomment below and substitute the names of your files
+% img_nii = '702-D2019_10_08.nii';
+% bval_file = fullfile(qamper_path,'test_data','702-D2019_10_08.bval');
+% roi_nii = fullfile(qamper_path,'test_data','702-D2019_10_08_ROI.nii.gz');
+
+% Optional function below: download the OSIPI test data from Zenodo
+[img_nii, bval_file, roi_nii] = getTestData;
+
+save_folder = fullfile(qamper_path,'output_files'); %change this location to where you want the output files to be saved
+if ~exist(save_folder,'dir')
+    mkdir(save_folder);
+end
+
+% Run the analysis
+batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);

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

@@ -0,0 +1,16 @@
+function [img_nii,bval_file,roi_nii] = getTestData
+
+[d,~,~] = fileparts(which('getTestData'));
+test_data_folder = fullfile(d,'test_data');
+if ~exist(test_data_folder,'dir')
+    mkdir(test_data_folder);
+end
+
+zenodo_url = 'https://zenodo.org/records/14605039/files/OSIPI_TF24_data_phantoms.zip';
+save_zip = fullfile(test_data_folder,'OSIPI_TF24_data_phantoms.zip');
+websave(save_zip,zenodo_url);
+unzip(save_zip, test_data_folder);
+
+img_nii = fulllfile(test_data_folder,'Data','brain.nii.gz');
+bval_file = fulllfile(test_data_folder,'Data','brain.bval');
+roi_nii = fulllfile(test_data_folder,'Data','brain_mask_gray_matter.nii.gz');

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

@@ -0,0 +1,17 @@
+function s = ivim_modeling_noise(x,b,sigma)
+%IVIM analysis
+%s = s0*((1-f)exp(-b.D) + f.exp(- b(D + Dstar)))
+%x = [f,D,Dstar];
+
+
+f = x(1);
+D = x(2);
+Dstar = x(3);
+s0 = x(4);
+
+s = s0*((1-f)*exp(-b*D) + f*exp(-b*Dstar));
+% n = 8;
+% navg = 4;
+% s = sqrt(s.^2 + 2*n*sigma.^2/navg);
+
+s = sqrt(s.^2 + sigma.^2);

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/qa_fitting/outputHNMPAResults.m

@@ -0,0 +1,57 @@
+function current_results_out = outputHNMPAResults(current_results,csv_outdir,modality,roiType,resultsCell)
+
+
+current_results_out = current_results;
+
+struct_string = [modality '_' roiType];
+
+%         current_results = outputHNMPAResults(current_results,'IVIM',roiType,{bval_arr,dwi_arr,[f_arr;D_arr;Dx_arr;s0_arr],fitted_dwi_arr,RSS,rms_val,chi,AIC,BIC});
+%                 current_results = outputHNMPAResults(current_results,'IVIM',roiType,{bval_arr,dwi_arr,[f_arr;D_arr;Dx_arr;s0_arr],fitted_dwi_arr,RSS,rms_val,chi,AIC,BIC,SNR,sigma,LB_IVIM,UB_IVIM,ROI});
+eval(['current_results_out.' struct_string '_b_value_fit = resultsCell{1};']); bval_arr = resultsCell{1};
+eval(['current_results_out.' struct_string '_dwi_data = resultsCell{2};']);
+
+params_out = resultsCell{3};
+
+eval(['current_results_out.' struct_string '_params_f = params_out(1,:);']); f = params_out(1,:);
+eval(['current_results_out.' struct_string '_features_f = returnFeatures(params_out(1,:));']);
+
+eval(['current_results_out.' struct_string '_params_D = params_out(2,:);']); D = params_out(2,:);
+eval(['current_results_out.' struct_string '_features_D = returnFeatures(params_out(2,:));']);
+
+eval(['current_results_out.' struct_string '_params_Dx = params_out(3,:);']); Dx = params_out(3,:);
+eval(['current_results_out.' struct_string '_features_Dx = returnFeatures(params_out(3,:));']);
+
+eval(['current_results_out.' struct_string '_params_s0 = params_out(4,:);']);
+eval(['current_results_out.' struct_string '_features_s0 = returnFeatures(params_out(4,:));']);
+
+eval(['current_results_out.' struct_string '_fitted_dwi_data1 = resultsCell{4};']);
+%Fit quality
+eval(['current_results_out.' struct_string '_quality_RSS = resultsCell{5};']); RSS = resultsCell{5};
+eval(['current_results_out.' struct_string '_quality_rms_val = resultsCell{6};']);
+eval(['current_results_out.' struct_string '_quality_chi = resultsCell{7};']);
+eval(['current_results_out.' struct_string '_quality_AIC = resultsCell{8};']);
+eval(['current_results_out.' struct_string '_quality_BIC = resultsCell{9};']);
+eval(['current_results_out.' struct_string '_quality_SNR = resultsCell{10};']); SNR = resultsCell{10};
+eval(['current_results_out.' struct_string '_quality_sigma = resultsCell{11};']); sigmadwi = resultsCell{11};
+eval(['current_results_out.' struct_string '_LB_IVIM = resultsCell{12};']); LB_IVIM = resultsCell{12};
+eval(['current_results_out.' struct_string '_UB_IVIM = resultsCell{13};']); UB_IVIM = resultsCell{13};
+eval(['current_results_out.' struct_string '_quality_Rsq = resultsCell{end-1};']); R_sq = resultsCell{end-1};
+eval(['current_results_out.' struct_string '_ROI_map = resultsCell{end};']); ROI = resultsCell{end};
+
+outstringheader='pat_ID,analysis_batch,tissue,D_Mean,D_Median,D_StDev,D_Skewness,D_Avg_Signal,D*_Mean,D*_Median,D*_StDev,D*_Skewness,D*_Avg_Signal,f_Mean,f_Median,f_StDev,f_Skewness,f_Avg_Signal,volume,voxels,RSS,b-values,SNR,sigma,R_sq,D_Kurtosis,D*_Kurtosis,f_Kurtosis';
+
+outcsv = [csv_outdir filesep 'dwi_IVIM.csv']; voxdim = current_results.dwinii.hdr.dime.pixdim(2) * current_results.dwinii.hdr.dime.pixdim(3) * current_results.dwinii.hdr.dime.pixdim(4);
+voxnum = numel(find(ROI));
+outstring = [current_results.pat_id ',' num2str(current_results.batch_id) ',' roiType ',' num2str(mean(D(1:end-1))) ',' num2str(median(D(1:end-1))) ',' num2str(std(D(1:end-1))) ',' num2str(skewness(D(1:end-1))) ',' num2str(D(end)) ','  ...
+    num2str(mean(Dx(1:end-1))) ',' num2str(median(Dx(1:end-1))) ',' num2str(std(Dx(1:end-1))) ',' num2str(skewness(Dx(1:end-1))) ',' num2str(Dx(end)) ','  ...
+    num2str(mean(f(1:end-1))) ',' num2str(median(f(1:end-1))) ',' num2str(std(f(1:end-1))) ',' num2str(skewness(f(1:end-1))) ',' num2str(f(end)) ','  ...
+    num2str(voxdim * voxnum) ',' num2str(voxnum) ',' num2str(mean(RSS)) ',' num2str(numel(bval_arr)) ',' num2str(SNR) ',' num2str(sigmadwi) ',' num2str(R_sq) ',' num2str(kurtosis(D(1:end-1))) ',' num2str(kurtosis(Dx(1:end-1))) ',' num2str(kurtosis(f(1:end-1)))];
+if ~exist(outcsv,'file')
+    fid = fopen(outcsv,'w');
+    fprintf(fid,'%s\n',outstringheader);
+else
+    fid = fopen(outcsv,'a');
+end
+fprintf(fid,'%s\n',outstring);
+fclose(fid);
+eval(['current_results_out.' struct_string '_CSV = outstring ;']);

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/qa_fitting/qualityFit.m

@@ -0,0 +1,31 @@
+function [RSS,rms_val,chi,AIC,BIC,R_sq_mean] = qualityFit(Nb,Np,observed_arr,predicted_arr)
+
+numVoxels = size(observed_arr,2);
+RSS = zeros(numVoxels,1);
+rms_val = zeros(numVoxels,1);
+chi = zeros(numVoxels,1);
+AIC = zeros(numVoxels,1);
+BIC = zeros(numVoxels,1);
+
+for i = 1:numVoxels
+    observed_signal = observed_arr(:,i);
+    predicted_signal = predicted_arr(:,i);
+
+    signal_diff = observed_signal - predicted_signal;
+
+    obs_mean = mean(observed_signal);
+    SST(i) = sum((observed_signal - obs_mean).^2);
+
+    RSS(i) = sum(signal_diff.^2);  %SSE
+
+    R_sq(i) = 1 - (RSS(i)/SST(i));
+
+    rms_val(i) = sqrt(RSS(i) / Nb);
+    chi(i) = (sum(signal_diff.^2 ./ predicted_signal)) / Nb;
+    AIC(i) = 2 * Np + Nb * log(RSS(i) / Nb) + (2 * Np * (Np + 1)) / (Nb - Np - 1); % AICc
+    BIC(i) = Nb * log(RSS(i)/Nb) + Np * log(Nb);
+end
+
+disp(['mean Rsq ' num2str(mean(R_sq(find(R_sq > 0))))]);
+
+R_sq_mean = mean(R_sq(find(R_sq > 0)));

src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/qa_fitting/returnFeatures.m

@@ -0,0 +1,10 @@
+function arrFeatures = returnFeatures(param_arr)
+
+measure_set = {'Mean','Median','StDev','Skewness','Avg_Signal'};
+
+arrFeatures = zeros(1,5);
+
+arrFeatures(1) = mean(param_arr);
+arrFeatures(2) = median(param_arr);
+arrFeatures(3) = std(param_arr);
+arrFeatures(4) = skewness(param_arr);