Skip to content

Commit c150706

Browse files
oliverchampionoliverchampion
authored
Merge pull request #52 from locastre/main
include IVIM code from Dr. Amita Shukla-Dave Lab at MSKCC.
2 parents dd869b9 + bd3dbf6 commit c150706

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

72 files changed

+15545
-0
lines changed

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
*.asv

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

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
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)
2+
3+
4+
5+
if size(bval_arr,2) < size(bval_arr,1)
6+
bval_arr = bval_arr';
7+
end
8+
9+
numVoxels = size(dwi_arr,2);
10+
numBvals = length(bval_arr);
11+
12+
f_lb = 0;
13+
% f_ub = 0.5;
14+
f_ub = 1.0;
15+
% f_ub = 0.75;
16+
% D_lb = 0;
17+
D_lb = 1e-6;
18+
D_ub = 0.003;
19+
% D_ub = 4e-3;
20+
% Dx_lb = 0;
21+
Dx_lb = 1e-6;
22+
Dx_ub = 5e-2;
23+
% Dx_ub = 0.5;
24+
% Dx_ub = 300e-3;
25+
s0_lb = 0;
26+
27+
% %Yousef bounds Jul 16 2018
28+
d10 = 1e-3;
29+
d20 = 1e-2;
30+
f0 = 0.2;
31+
32+
% p0(i,:) = [ smax(i) d10 d20 f0];
33+
% lb(i,:) = [0.5*smax(i) 0.25*d10 0.25*d20 0.0];
34+
% ub(i,:) = [2.0*smax(i) 4.0*d10 4.0*d20 1.0];
35+
36+
if nargin < 6
37+
x0in = [f0 d10 d20 1]; % entry for s (4th) is multiplicative factor, not absolute
38+
if nargin < 4
39+
% LB = [f_lb D_lb Dx_lb s0_lb];
40+
% UB = [f_ub D_ub Dx_ub 0];
41+
% entry for s (4th) is multiplicative factor, not absolute
42+
LB0 = [0 0.25*d10 0.25*d20 0.5];
43+
UB0 = [1 4*d10 4*d20 2];
44+
end
45+
end
46+
47+
optimization_iterations = 4;
48+
49+
options = optimset('MaxFunEvals',400,'MaxIter',200, ...
50+
'TolFun',1e-6,'TolX',1e-6,'Display','off');
51+
52+
axisLabels = {'b-value (s/mm^2)','signal (a.u.)'};
53+
54+
tic
55+
for i = 1:numVoxels
56+
if rem(i,1000) == 0
57+
toc
58+
disp(['Voxel ' num2str(i) ' of ' num2str(numVoxels)]);
59+
tic
60+
end
61+
s = dwi_arr(:,i);
62+
63+
LB = LB0;
64+
UB = UB0;
65+
x0 = x0in;
66+
67+
LB(end) = LB0(end)*max(s);
68+
UB(end) = UB0(end)*max(s);
69+
x0(end) = x0in(end)*max(s);
70+
% disp(['Voxel ' num2str(i) ' of ' num2str(numVoxels)]);
71+
if ~parallelFlag
72+
for j = 1:optimization_iterations
73+
% x0 = (LB + UB) / 2; %+ (UB - LB) * (rand - 0.5);
74+
[x_fit(j,:),resnorm(j),~,~,~] = lsqcurvefit(@(x,b)ivim_modeling_noise(x,b,sigmadwi), ...
75+
x0, bval_arr', s, LB, UB, options);
76+
end
77+
else
78+
parfor j = 1:optimization_iterations
79+
% x0 = (LB + UB) / 2; %+ (UB - LB) * (rand - 0.5);
80+
[x_fit(j,:),resnorm(j),~,~,~] = lsqcurvefit(@(x,b)ivim_modeling_noise(x,b,sigmadwi), ...
81+
x0, bval_arr', s, LB, UB, options);
82+
end
83+
end
84+
[~,best_fit_idx] = min(resnorm);
85+
p = x_fit(best_fit_idx,:);
86+
87+
fitted_dwi_arr(:,i) = ivim_modeling_noise(p,bval_arr',sigmadwi);
88+
f_arr(i) = p(1);
89+
D_arr(i) = p(2)*(1e3);
90+
if p(3) > 1e-5
91+
Dx_arr(i) = p(3)*(1e3);
92+
else
93+
Dx_arr(i) = 1e-5;
94+
end
95+
s0_arr(i) = p(4);
96+
97+
if previewMode == 1 || previewMode == 100
98+
axisLabels{3} = ['IVIM Fit (Voxel # ' num2str(i) '/' num2str(numVoxels) ')'];
99+
if previewMode == 100 && i < 101
100+
updatePreviewAxes(previewAxes,bval_arr',s,fitted_dwi_arr(:,i),axisLabels);
101+
elseif previewMode == 1
102+
updatePreviewAxes(previewAxes,bval_arr',s,fitted_dwi_arr(:,i),axisLabels);
103+
end
104+
end
105+
106+
end
107+
108+
% Quantification of Fitting Quality
109+
Nb = numel(bval_arr);
110+
Np = numel(p);
111+
observed_arr = dwi_arr;
112+
predicted_arr = fitted_dwi_arr;
113+
[RSS,rms_val,chi,AIC,BIC,R_sq] = qualityFit(Nb,Np,observed_arr,predicted_arr);

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

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
# MRI-QAMPER_IVIM
2+
3+
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.
4+
5+
__Authors:__ Eve LoCastro (locastre@mskcc.org), Dr. Ramesh Paudyal (paudyalr@mskcc.org), Dr. Amita Shukla-Dave (davea@mskcc.org) </br>
6+
__Institution:__ Memorial Sloan Kettering Cancer Center </br>
7+
__Department:__ Medical Physics</br>
8+
__Address:__ 321 E 61st St, New York, NY 10022
9+
10+
The codebase and subdirectories should be added to the MATLAB path. An example usage script is provided in `demo_QAMPER_IVIM.m`.
11+
12+
```
13+
% This is an example usage script for MSK Medical Physics Dave Lab QAMPER IVIM
14+
% Please replace the variable names below with path values for
15+
% qamper_path: path to MRI-QAMPER_IVIM folder
16+
% img_nii: multi-b value DWI (4-D NIfTI)
17+
% bval_file: b-value information (txt file)
18+
% roi_nii: single-volume mask ROI image (NIfTI)
19+
20+
qamper_path = 'path:\to\MRI-QAMPER_IVIM';
21+
addpath(genpath(qamper_path));
22+
23+
img_nii = 'dwi.nii';
24+
bval_file = 'dwi.bval';
25+
roi_nii = 'roi.nii';
26+
27+
save_folder = fullfile(qamper_path,'test_data');
28+
29+
batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);
30+
```

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

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function [uniq_bval_arr, dwi_img_new] = avgRepeatBvalVols(bval_arr,dwi_img)
2+
3+
disp('Averaging multiple bval vols to increase signal')
4+
5+
[uniq_bval_arr,IA,~] = unique(bval_arr);
6+
7+
dwi_img_new = [];
8+
9+
for i = 1:numel(uniq_bval_arr)
10+
if i < numel(uniq_bval_arr)
11+
dwi_avg_vol = mean(dwi_img(:,:,:,IA(i):IA(i+1)-1),4);
12+
else
13+
dwi_avg_vol = mean(dwi_img(:,:,:,IA(i):end),4);
14+
end
15+
dwi_img_new = cat(4,dwi_img_new, dwi_avg_vol);
16+
end

‎src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/bvalOrderFix.m

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
function [bval_arr_new, dwi_img_new] = bvalOrderFix(bval_arr,dwi_img)
2+
3+
disp(['Reordering B0 to start of imaging sequence'])
4+
5+
if size(bval_arr,1) > size(bval_arr,2)
6+
bval_arr = bval_arr';
7+
end
8+
9+
% bval_arr_new = [bval_arr(end) bval_arr(1:end-1)];
10+
[bval_arr_new,idx] = sort(bval_arr);
11+
12+
% dwi_img_new = cat(4,dwi_img(:,:,:,end),dwi_img(:,:,:,1:end-1));
13+
dwi_img_new = dwi_img(:,:,:,idx);
14+
15+
disp(bval_arr_new);

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

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
function [dwi_arr, bval_arr, sigmadwi, SNR, noise] = dataPrepDWI(dwi_img,bval_arr,ROI,channels,N_avgs)
2+
3+
bgMask = zeros(size(dwi_img,1),size(dwi_img,2));
4+
bgMask(10:30,10:30) = 1;
5+
% bgidx = sub2ind(size(bgMask),10:30,10:30);
6+
bgidx = find(bgMask);
7+
8+
maskidx = find(ROI);
9+
10+
dwi_arr = map2arr(dwi_img,ROI);
11+
dwi_arr(:,end+1) = mean(dwi_arr,2);
12+
13+
noise_square = bgMask.*dwi_img(:,:,1,1);
14+
noise_arr = noise_square(bgidx);noise_arr;
15+
16+
%n = 8; % need to update with new protocol
17+
%navg = 4;
18+
n = channels;
19+
navg = N_avgs;
20+
21+
% if data were converted to nifti before conv_panel started to extract
22+
% channels and N_avgs, then those numbers default to this
23+
if isempty(n)
24+
n = 16;
25+
navg = 2;
26+
end
27+
28+
sigmadwi = sqrt(navg*norm(noise_arr(:))^2/length(noise_arr(:)))/sqrt(2*n);
29+
noise = mean(noise_arr);
30+
SNR = mean(dwi_arr(1,:))/sigmadwi;

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

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
% Authors: Eve LoCastro (locastre@mskcc.org), Dr. Ramesh Paudyal (paudyalr@mskcc.org), Dr. Amita Shukla-Dave (davea@mskcc.org)
2+
% Institution: Memorial Sloan Kettering Cancer Center
3+
% Address: 321 E 61st St, New York, NY 10022
4+
5+
% This is an example usage script for MSK Medical Physics Dave Lab QAMPER IVIM
6+
% Please replace the variable names below with path values for
7+
% qamper_path: path to MRI-QAMPER_IVIM folder
8+
% img_nii: multi-b value DWI (4-D NIfTI)
9+
% bval_file: b-value information (txt file)
10+
% roi_nii: single-volume mask ROI image (NIfTI)
11+
12+
13+
qamper_path = 'path:\to\MRI-QAMPER_IVIM'; %path where the MRI-QAMPER folder is located
14+
addpath(genpath(qamper_path));
15+
16+
% Uncomment below and substitute the names of your files
17+
% img_nii = '702-D2019_10_08.nii';
18+
% bval_file = fullfile(qamper_path,'test_data','702-D2019_10_08.bval');
19+
% roi_nii = fullfile(qamper_path,'test_data','702-D2019_10_08_ROI.nii.gz');
20+
21+
% Optional function below: download the OSIPI test data from Zenodo
22+
[img_nii, bval_file, roi_nii] = getTestData;
23+
24+
save_folder = fullfile(qamper_path,'output_files'); %change this location to where you want the output files to be saved
25+
if ~exist(save_folder,'dir')
26+
mkdir(save_folder);
27+
end
28+
29+
% Run the analysis
30+
batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);

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

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function [img_nii,bval_file,roi_nii] = getTestData
2+
3+
[d,~,~] = fileparts(which('getTestData'));
4+
test_data_folder = fullfile(d,'test_data');
5+
if ~exist(test_data_folder,'dir')
6+
mkdir(test_data_folder);
7+
end
8+
9+
zenodo_url = 'https://zenodo.org/records/14605039/files/OSIPI_TF24_data_phantoms.zip';
10+
save_zip = fullfile(test_data_folder,'OSIPI_TF24_data_phantoms.zip');
11+
websave(save_zip,zenodo_url);
12+
unzip(save_zip, test_data_folder);
13+
14+
img_nii = fulllfile(test_data_folder,'Data','brain.nii.gz');
15+
bval_file = fulllfile(test_data_folder,'Data','brain.bval');
16+
roi_nii = fulllfile(test_data_folder,'Data','brain_mask_gray_matter.nii.gz');

‎src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/inpaint_nans.m

Lines changed: 1 addition & 0 deletions
Large diffs are not rendered by default.

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

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
function s = ivim_modeling_noise(x,b,sigma)
2+
%IVIM analysis
3+
%s = s0*((1-f)exp(-b.D) + f.exp(- b(D + Dstar)))
4+
%x = [f,D,Dstar];
5+
6+
7+
f = x(1);
8+
D = x(2);
9+
Dstar = x(3);
10+
s0 = x(4);
11+
12+
s = s0*((1-f)*exp(-b*D) + f*exp(-b*Dstar));
13+
% n = 8;
14+
% navg = 4;
15+
% s = sqrt(s.^2 + 2*n*sigma.^2/navg);
16+
17+
s = sqrt(s.^2 + sigma.^2);

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

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
function current_results_out = outputHNMPAResults(current_results,csv_outdir,modality,roiType,resultsCell)
2+
3+
4+
current_results_out = current_results;
5+
6+
struct_string = [modality '_' roiType];
7+
8+
% 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});
9+
% 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});
10+
eval(['current_results_out.' struct_string '_b_value_fit = resultsCell{1};']); bval_arr = resultsCell{1};
11+
eval(['current_results_out.' struct_string '_dwi_data = resultsCell{2};']);
12+
13+
params_out = resultsCell{3};
14+
15+
eval(['current_results_out.' struct_string '_params_f = params_out(1,:);']); f = params_out(1,:);
16+
eval(['current_results_out.' struct_string '_features_f = returnFeatures(params_out(1,:));']);
17+
18+
eval(['current_results_out.' struct_string '_params_D = params_out(2,:);']); D = params_out(2,:);
19+
eval(['current_results_out.' struct_string '_features_D = returnFeatures(params_out(2,:));']);
20+
21+
eval(['current_results_out.' struct_string '_params_Dx = params_out(3,:);']); Dx = params_out(3,:);
22+
eval(['current_results_out.' struct_string '_features_Dx = returnFeatures(params_out(3,:));']);
23+
24+
eval(['current_results_out.' struct_string '_params_s0 = params_out(4,:);']);
25+
eval(['current_results_out.' struct_string '_features_s0 = returnFeatures(params_out(4,:));']);
26+
27+
eval(['current_results_out.' struct_string '_fitted_dwi_data1 = resultsCell{4};']);
28+
%Fit quality
29+
eval(['current_results_out.' struct_string '_quality_RSS = resultsCell{5};']); RSS = resultsCell{5};
30+
eval(['current_results_out.' struct_string '_quality_rms_val = resultsCell{6};']);
31+
eval(['current_results_out.' struct_string '_quality_chi = resultsCell{7};']);
32+
eval(['current_results_out.' struct_string '_quality_AIC = resultsCell{8};']);
33+
eval(['current_results_out.' struct_string '_quality_BIC = resultsCell{9};']);
34+
eval(['current_results_out.' struct_string '_quality_SNR = resultsCell{10};']); SNR = resultsCell{10};
35+
eval(['current_results_out.' struct_string '_quality_sigma = resultsCell{11};']); sigmadwi = resultsCell{11};
36+
eval(['current_results_out.' struct_string '_LB_IVIM = resultsCell{12};']); LB_IVIM = resultsCell{12};
37+
eval(['current_results_out.' struct_string '_UB_IVIM = resultsCell{13};']); UB_IVIM = resultsCell{13};
38+
eval(['current_results_out.' struct_string '_quality_Rsq = resultsCell{end-1};']); R_sq = resultsCell{end-1};
39+
eval(['current_results_out.' struct_string '_ROI_map = resultsCell{end};']); ROI = resultsCell{end};
40+
41+
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';
42+
43+
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);
44+
voxnum = numel(find(ROI));
45+
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)) ',' ...
46+
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)) ',' ...
47+
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)) ',' ...
48+
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)))];
49+
if ~exist(outcsv,'file')
50+
fid = fopen(outcsv,'w');
51+
fprintf(fid,'%s\n',outstringheader);
52+
else
53+
fid = fopen(outcsv,'a');
54+
end
55+
fprintf(fid,'%s\n',outstring);
56+
fclose(fid);
57+
eval(['current_results_out.' struct_string '_CSV = outstring ;']);

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

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
function [RSS,rms_val,chi,AIC,BIC,R_sq_mean] = qualityFit(Nb,Np,observed_arr,predicted_arr)
2+
3+
numVoxels = size(observed_arr,2);
4+
RSS = zeros(numVoxels,1);
5+
rms_val = zeros(numVoxels,1);
6+
chi = zeros(numVoxels,1);
7+
AIC = zeros(numVoxels,1);
8+
BIC = zeros(numVoxels,1);
9+
10+
for i = 1:numVoxels
11+
observed_signal = observed_arr(:,i);
12+
predicted_signal = predicted_arr(:,i);
13+
14+
signal_diff = observed_signal - predicted_signal;
15+
16+
obs_mean = mean(observed_signal);
17+
SST(i) = sum((observed_signal - obs_mean).^2);
18+
19+
RSS(i) = sum(signal_diff.^2); %SSE
20+
21+
R_sq(i) = 1 - (RSS(i)/SST(i));
22+
23+
rms_val(i) = sqrt(RSS(i) / Nb);
24+
chi(i) = (sum(signal_diff.^2 ./ predicted_signal)) / Nb;
25+
AIC(i) = 2 * Np + Nb * log(RSS(i) / Nb) + (2 * Np * (Np + 1)) / (Nb - Np - 1); % AICc
26+
BIC(i) = Nb * log(RSS(i)/Nb) + Np * log(Nb);
27+
end
28+
29+
disp(['mean Rsq ' num2str(mean(R_sq(find(R_sq > 0))))]);
30+
31+
R_sq_mean = mean(R_sq(find(R_sq > 0)));

‎src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/qa_fitting/returnFeatures.log

Whitespace-only changes.

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

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
function arrFeatures = returnFeatures(param_arr)
2+
3+
measure_set = {'Mean','Median','StDev','Skewness','Avg_Signal'};
4+
5+
arrFeatures = zeros(1,5);
6+
7+
arrFeatures(1) = mean(param_arr);
8+
arrFeatures(2) = median(param_arr);
9+
arrFeatures(3) = std(param_arr);
10+
arrFeatures(4) = skewness(param_arr);

0 commit comments

Comments
 (0)