Skip to content

Commit bd82813

Browse files
locastrelocastre
committed
include IVIM code from Dr. Amita Shukla-Dave Lab at MSKCC
1 parent 1f4cbf4 commit bd82813

Some content is hidden

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

75 files changed

+16020
-0
lines changed

‎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.asv

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
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)
6+
Institution: Memorial Sloan Kettering Cancer Center
7+
Department: Medical Physics
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+
21+
22+
qamper_path = 'path:\to\MRI-QAMPER_IVIM';
23+
addpath(genpath(qamper_path));
24+
25+
img_nii = 'dwi.nii';
26+
bval_file = fullfile(qamper_path,'test_data','702-HN401D-D2019_10_08.bval');
27+
roi_nii = fullfile(qamper_path,'test_data','702-HN401D-D2019_10_08_BD_REDO_SV.nii.gz');
28+
29+
save_folder = fullfile(qamper_path,'test_data');
30+
31+
batchResultsFolder = run_QAMPER_IVIM(img_nii,bval_file,roi_nii,save_folder);
32+
33+
```

‎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/buildBatchFromCSV.m

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
function batchStruct = buildBatchFromCSV(img_nii,roi_nii,qamper_mode,optionS)
2+
3+
batchStruct = [];
4+
batchStruct.batch_id = num2str(now);
5+
6+
if strcmp(qamper_mode,'DWI')
7+
if isfield(optionS,'DWIstruct')
8+
DWIstruct = optionS.DWIstruct;
9+
end
10+
11+
DWIstruct.runDWI = 1;
12+
[droot,fbase,ext] = fileparts(img_nii);
13+
batchStruct.currentNiiPath = fullfile(droot,'dwi_maps');
14+
15+
16+
if strcmp(ext,'.gz')
17+
[~,fbase,~] = fileparts(fbase);
18+
end
19+
20+
bval_mat = fullfile(droot,[fbase '.mat']);
21+
22+
T = load(bval_mat);
23+
24+
DWIstruct.bval_arr = T.bval_arr;
25+
26+
img = nii_load(img_nii,1);
27+
vol_num = size(img.img,4);
28+
29+
DWIstruct.files{1} = img_nii; DWIstruct.vols = vol_num;
30+
DWIstruct.files{2} = roi_nii;
31+
32+
if ~isfield(DWIstruct,'ADC')
33+
DWIstruct.ADC = 0;
34+
end
35+
if ~isfield(DWIstruct,'IVIM')
36+
DWIstruct.IVIM = 0;
37+
end
38+
if ~isfield(DWIstruct,'NGIVIM')
39+
DWIstruct.NGIVIM = 0;
40+
end
41+
if ~isfield(DWIstruct,'Kurtosis')
42+
DWIstruct.Kurtosis = 0;
43+
end
44+
45+
if ~isfield(DWIstruct,'roiType')
46+
DWIstruct.roiType = 'node';
47+
end
48+
49+
if ~isfield(DWIstruct,'avgRepeatBvals')
50+
DWIstruct.avgRepeatBvals = 1;
51+
end
52+
53+
if ~isfield(DWIstruct,'channels')
54+
DWIstruct.channels = 16; %
55+
end
56+
57+
if ~isfield(DWIstruct, 'avgChannels')
58+
DWIstruct.avgChannels = 2; %
59+
end
60+
61+
if DWIstruct.IVIM
62+
if ~isfield(DWIstruct,'LB_IVIM')
63+
IVIMBoundsTable = [{0} {0} {0} {0}; {0.5} {0.005} {0.5} {1.0}; {0.05} {0.0005} {0.01} {0.5}];
64+
65+
DWIstruct.LB_IVIM = cell2mat(IVIMBoundsTable(1,:));
66+
DWIstruct.UB_IVIM = cell2mat(IVIMBoundsTable(2,:));
67+
DWIstruct.x0_IVIM = cell2mat(IVIMBoundsTable(3,:));
68+
end
69+
end
70+
if DWIstruct.NGIVIM
71+
if ~isfield(DWIstruct,'LB_NG')
72+
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}];
73+
74+
DWIstruct.LB_NG = cell2mat(NGIVIMDefaultBoundsTable(1,:));
75+
DWIstruct.UB_NG = cell2mat(NGIVIMDefaultBoundsTable(2,:));
76+
DWIstruct.x0_NG = cell2mat(NGIVIMDefaultBoundsTable(3,:));
77+
end
78+
end
79+
80+
if ~isfield(DWIstruct,'Parallel')
81+
DWIstruct.Parallel = 1;
82+
end
83+
84+
DWIstruct.bval_tf = ones(size(DWIstruct.bval_arr));
85+
%
86+
% DWIstruct.previewAxes = handles.fittingAxes;
87+
%
88+
DWIstruct.nanopt = 1;
89+
90+
if ~isfield(DWIstruct,'kernelsize')
91+
DWIstruct.kernelsize = 1;
92+
end
93+
if ~isfield(DWIstruct,'smoothopt')
94+
DWIstruct.smoothopt = 1;
95+
end
96+
97+
DWIstruct.previewMode = 0;
98+
DWIstruct.previewAxes = 0;
99+
100+
else
101+
DWIstruct.runDWI = 0;
102+
end
103+
104+
batchStruct.DWIstruct = DWIstruct;
105+
batchStruct.T1struct.runT1 = 0;
106+
batchStruct.DCEstruct.runDCE = 0;
107+
batchStruct.T2struct.runT2 = 0;

0 commit comments

Comments
 (0)