Skip to content

Commit 37a2fba

Browse files
oliverchampionoliverchampion
authored
Merge pull request #20 from OSIPI/testing/add_ETP_and_add_configuration
Testing/add etp and add configuration
2 parents 05bed07 + b2ed0be commit 37a2fba

25 files changed

+757
-200
lines changed

‎.gitignore

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,6 @@
33
__pycache__/
44
*.nii.gz
55
bvals.txt
6-
*json
76
# Unit test / coverage reports
87
.tox/
98
.coverage

‎conftest.py

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
import pytest
2+
import pathlib
3+
import json
4+
import csv
5+
# import datetime
6+
7+
8+
def pytest_addoption(parser):
9+
parser.addoption(
10+
"--SNR",
11+
default=[100],
12+
nargs="+",
13+
type=int,
14+
help="Evaluation test SNR",
15+
)
16+
parser.addoption(
17+
"--ricianNoise",
18+
default=False,
19+
type=bool,
20+
help="Use Rician noise, non-rician is gaussian",
21+
)
22+
parser.addoption(
23+
"--algorithmFile",
24+
default="tests/IVIMmodels/unit_tests/algorithms.json",
25+
type=str,
26+
help="Algorithm file name",
27+
)
28+
parser.addoption(
29+
"--dataFile",
30+
default="tests/IVIMmodels/unit_tests/generic.json",
31+
type=str,
32+
help="Default data file name",
33+
)
34+
parser.addoption(
35+
"--saveFileName",
36+
default="",
37+
type=str,
38+
help="Saved results file name",
39+
)
40+
parser.addoption(
41+
"--rtol",
42+
default=1,
43+
type=float,
44+
help="Relative tolerance",
45+
)
46+
parser.addoption(
47+
"--atol",
48+
default=1,
49+
type=float,
50+
help="Absolute tolerance",
51+
)
52+
parser.addoption(
53+
"--fitCount",
54+
default=10,
55+
type=int,
56+
help="Number of fits to perform on the same parameters",
57+
)
58+
parser.addoption(
59+
"--saveDurationFileName",
60+
default="",
61+
type=str,
62+
help="Saved duration results file name",
63+
)
64+
65+
66+
@pytest.fixture(scope="session")
67+
def save_file(request):
68+
filename = request.config.getoption("--saveFileName")
69+
if filename:
70+
current_folder = pathlib.Path.cwd()
71+
filename = current_folder / filename
72+
# print(filename)
73+
# filename.unlink(missing_ok=True)
74+
filename = filename.as_posix()
75+
with open(filename, "w") as csv_file:
76+
writer = csv.writer(csv_file, delimiter=',')
77+
writer.writerow(("Algorithm", "Region", "SNR", "index", "f", "Dp", "D", "f_fitted", "Dp_fitted", "D_fitted"))
78+
# writer.writerow(["", datetime.datetime.now()])
79+
return filename
80+
81+
@pytest.fixture(scope="session")
82+
def save_duration_file(request):
83+
filename = request.config.getoption("--saveDurationFileName")
84+
if filename:
85+
current_folder = pathlib.Path.cwd()
86+
filename = current_folder / filename
87+
# print(filename)
88+
# filename.unlink(missing_ok=True)
89+
filename = filename.as_posix()
90+
with open(filename, "w") as csv_file:
91+
writer = csv.writer(csv_file, delimiter=',')
92+
writer.writerow(("Algorithm", "Region", "SNR", "Duration [us]", "Count"))
93+
# writer.writerow(["", datetime.datetime.now()])
94+
return filename
95+
96+
@pytest.fixture(scope="session")
97+
def rtol(request):
98+
return request.config.getoption("--rtol")
99+
100+
101+
@pytest.fixture(scope="session")
102+
def atol(request):
103+
return request.config.getoption("--atol")
104+
105+
@pytest.fixture(scope="session")
106+
def fit_count(request):
107+
return request.config.getoption("--fitCount")
108+
109+
@pytest.fixture(scope="session")
110+
def rician_noise(request):
111+
return request.config.getoption("--ricianNoise")
112+
113+
114+
def pytest_generate_tests(metafunc):
115+
if "SNR" in metafunc.fixturenames:
116+
metafunc.parametrize("SNR", metafunc.config.getoption("SNR"))
117+
if "ivim_algorithm" in metafunc.fixturenames:
118+
algorithms = algorithm_list(metafunc.config.getoption("algorithmFile"))
119+
metafunc.parametrize("ivim_algorithm", algorithms)
120+
if "algorithm_file" in metafunc.fixturenames:
121+
metafunc.parametrize("algorithm_file", metafunc.config.getoption("algorithmFile"))
122+
if "ivim_data" in metafunc.fixturenames:
123+
data = data_list(metafunc.config.getoption("dataFile"))
124+
metafunc.parametrize("ivim_data", data)
125+
if "data_file" in metafunc.fixturenames:
126+
metafunc.parametrize("data_file", metafunc.config.getoption("dataFile"))
127+
128+
129+
def algorithm_list(filename):
130+
current_folder = pathlib.Path.cwd()
131+
algorithm_path = current_folder / filename
132+
with algorithm_path.open() as f:
133+
algorithm_information = json.load(f)
134+
return algorithm_information["algorithms"]
135+
136+
def data_list(filename):
137+
current_folder = pathlib.Path.cwd()
138+
data_path = current_folder / filename
139+
with data_path.open() as f:
140+
all_data = json.load(f)
141+
142+
bvals = all_data.pop('config')
143+
bvals = bvals['bvalues']
144+
for name, data in all_data.items():
145+
yield name, bvals, data

‎phantoms/brain/sim_brain_phantom.py

Lines changed: 56 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -5,57 +5,59 @@
55
import nibabel as nib
66
from scipy.ndimage import zoom
77

8-
DIFFUSIVE_REGIME = 'diffusive'
9-
BALLISTIC_REGIME = 'ballistic'
10-
11-
folder = os.path.dirname(__file__)
12-
13-
###########################
14-
# Simulation parameters #
15-
regime = DIFFUSIVE_REGIME
16-
snr = 200
17-
resolution = [3,3,3]
18-
# #
19-
###########################
20-
21-
22-
# Ground truth
23-
nii = nib.load(os.path.join(folder,'ground_truth','hrgt_icbm_2009a_nls_3t.nii.gz'))
24-
segmentation = np.squeeze(nii.get_fdata()[...,-1])
25-
26-
with open(os.path.join(folder,'ground_truth',regime+'_groundtruth.json'), 'r') as f:
27-
ivim_pars = json.load(f)
28-
S0 = 1
29-
30-
# Sequence parameters
31-
bval_file = os.path.join(folder,'ground_truth',regime+'.bval')
32-
b = np.loadtxt(bval_file)
33-
if regime == BALLISTIC_REGIME:
34-
cval_file = bval_file.replace('bval','cval')
35-
c = np.loadtxt(cval_file)
36-
37-
# Calculate signal
38-
S = np.zeros(list(np.shape(segmentation))+[b.size])
39-
40-
if regime == BALLISTIC_REGIME:
41-
Db = ivim_pars["Db"]
42-
for i,(D,f,vd) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["vd"])):
43-
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Db-c**2*vd**2))
44-
else:
45-
for i,(D,f,Dstar) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"])):
46-
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Dstar))
47-
48-
# Resample to suitable resolution
49-
im = zoom(S,np.append(np.diag(nii.affine)[:3]/np.array(resolution),1),order=1)
50-
sz = im.shape
51-
52-
# Add noise
53-
im_noise = np.abs(im + S0/snr*(np.random.randn(sz[0],sz[1],sz[2],sz[3])+1j*np.random.randn(sz[0],sz[1],sz[2],sz[3])))
54-
55-
# Save as image and sequence parameters
56-
nii_out = nib.Nifti1Image(im_noise,np.eye(4))
57-
base_name = os.path.join(folder,'data','{}_snr{}'.format(regime,snr))
58-
nib.save(nii_out,base_name+'.nii.gz')
59-
shutil.copyfile(bval_file,base_name+'.bval')
60-
if regime == BALLISTIC_REGIME:
61-
shutil.copyfile(cval_file,base_name+'.cval')
8+
if __name__ == "__main__":
9+
10+
DIFFUSIVE_REGIME = 'diffusive'
11+
BALLISTIC_REGIME = 'ballistic'
12+
13+
folder = os.path.dirname(__file__)
14+
15+
###########################
16+
# Simulation parameters #
17+
regime = DIFFUSIVE_REGIME
18+
snr = 200
19+
resolution = [3,3,3]
20+
# #
21+
###########################
22+
23+
24+
# Ground truth
25+
nii = nib.load(os.path.join(folder,'ground_truth','hrgt_icbm_2009a_nls_3t.nii.gz'))
26+
segmentation = np.squeeze(nii.get_fdata()[...,-1])
27+
28+
with open(os.path.join(folder,'ground_truth',regime+'_groundtruth.json'), 'r') as f:
29+
ivim_pars = json.load(f)
30+
S0 = 1
31+
32+
# Sequence parameters
33+
bval_file = os.path.join(folder,'ground_truth',regime+'.bval')
34+
b = np.loadtxt(bval_file)
35+
if regime == BALLISTIC_REGIME:
36+
cval_file = bval_file.replace('bval','cval')
37+
c = np.loadtxt(cval_file)
38+
39+
# Calculate signal
40+
S = np.zeros(list(np.shape(segmentation))+[b.size])
41+
42+
if regime == BALLISTIC_REGIME:
43+
Db = ivim_pars["Db"]
44+
for i,(D,f,vd) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["vd"])):
45+
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Db-c**2*vd**2))
46+
else:
47+
for i,(D,f,Dstar) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"])):
48+
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Dstar))
49+
50+
# Resample to suitable resolution
51+
im = zoom(S,np.append(np.diag(nii.affine)[:3]/np.array(resolution),1),order=1)
52+
sz = im.shape
53+
54+
# Add noise
55+
im_noise = np.abs(im + S0/snr*(np.random.randn(sz[0],sz[1],sz[2],sz[3])+1j*np.random.randn(sz[0],sz[1],sz[2],sz[3])))
56+
57+
# Save as image and sequence parameters
58+
nii_out = nib.Nifti1Image(im_noise,np.eye(4))
59+
base_name = os.path.join(folder,'data','{}_snr{}'.format(regime,snr))
60+
nib.save(nii_out,base_name+'.nii.gz')
61+
shutil.copyfile(bval_file,base_name+'.bval')
62+
if regime == BALLISTIC_REGIME:
63+
shutil.copyfile(cval_file,base_name+'.cval')

‎pytest.ini

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
[pytest]
2+
markers =
3+
slow: marks tests as slow (deselect with '-m "not slow"')
4+
addopts =
5+
-m 'not slow'

‎src/original/ETP_SRI/LinearFitting.py

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,11 +66,20 @@ def ivim_fit(self, bvalues, signal):
6666
gt_cutoff = bvalues >= self.linear_cutoff
6767
linear_signal = np.log(signal)
6868
D = self.linear_fit(bvalues[gt_cutoff], linear_signal[gt_cutoff])
69+
# print(D)
70+
D[1] = max(D[1], 0) # constrain to positive values
6971

7072
if lt_cutoff.sum() > 0:
7173
signal_Dp = linear_signal[lt_cutoff] - gd.linear_signal(D[1], bvalues[lt_cutoff], np.log(D[0]))
72-
73-
Dp_prime = self.linear_fit(bvalues[lt_cutoff], np.log(signal_Dp))
74+
# print(signal_Dp)
75+
signal_valid = signal_Dp > 0
76+
lt_cutoff_dual = np.logical_and(lt_cutoff[:len(signal_Dp)], signal_valid)
77+
# print(lt_cutoff_dual)
78+
Dp_prime = [-1, -1]
79+
if lt_cutoff_dual.sum() > 0:
80+
# print(np.log(signal_Dp[lt_cutoff_dual]))
81+
Dp_prime = self.linear_fit(bvalues[:len(signal_Dp)][lt_cutoff_dual], np.log(signal_Dp[lt_cutoff_dual]))
82+
# print(Dp_prime)
7483

7584
if np.any(np.asarray(Dp_prime) < 0) or not np.all(np.isfinite(Dp_prime)):
7685
print('Perfusion fit failed')

‎src/original/PV_MUMC/two_step_IVIM_fit.py

Lines changed: 39 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -75,45 +75,48 @@ def fit_least_squares(bvalues, dw_data, IR=False, S0_output=False, fitS0=True,
7575
:return Dmv: scalar with Dmv of the specific voxel
7676
"""
7777

78-
try:
79-
def monofit(bvalues, Dpar):
80-
return np.exp(-bvalues * Dpar)
81-
82-
high_b = bvalues[bvalues >= cutoff]
83-
high_dw_data = dw_data[bvalues >= cutoff]
84-
boundspar = ([bounds[0][1]], [bounds[1][1]])
85-
params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2], bounds=boundspar)
86-
Dpar1 = params[0]
78+
#try:
79+
def monofit(bvalues, Dpar):
80+
return np.exp(-bvalues * Dpar)
81+
82+
high_b = bvalues[bvalues >= cutoff]
83+
high_dw_data = dw_data[bvalues >= cutoff]
84+
boundspar = ([bounds[0][1]], [bounds[1][1]])
85+
params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2], bounds=boundspar)
86+
Dpar = params[0]
8787

88-
if not fitS0:
89-
boundsupdated=([Dpar1 , bounds[0][2] , bounds[0][3] ],
90-
[Dpar1 , bounds[1][2] , bounds[1][3] ])
91-
params, _ = curve_fit(two_exp_noS0, bvalues, dw_data, p0=[Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
92-
Dpar, Fmv, Dmv = params[0], params[1], params[2]
93-
#when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
94-
if Fmv < 1e-4:
95-
Dmv = float("NaN")
88+
if not fitS0:
89+
boundsupdated=([Dpar1 , bounds[0][2] , bounds[0][3] ],
90+
[Dpar1 , bounds[1][2] , bounds[1][3] ])
91+
params, _ = curve_fit(two_exp_noS0, bvalues, dw_data, p0=[Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
92+
Dpar1, Fmv, Dmv = params[0], params[1], params[2]
93+
#when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
94+
if Fmv < 1e-4:
95+
Dmv = float("NaN")
96+
97+
else:
98+
#boundsupdated = ([bounds[0][0] , Dpar1 , bounds[0][2] , bounds[0][3] ],
99+
# [bounds[1][0] , Dpar1, bounds[1][2] , bounds[1][3] ])
100+
#params, _ = curve_fit(two_exp, bvalues, dw_data, p0=[1, Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
101+
boundsupdated = ([bounds[0][0] , bounds[0][2] , bounds[0][3] ],
102+
[bounds[1][0] , bounds[1][2] , bounds[1][3] ])
103+
params, _ = curve_fit(lambda bvalues, S0, Fmv, Dmv: two_exp(bvalues, S0, Dpar, Fmv, Dmv), bvalues, dw_data, p0=[1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
104+
S0 = params[0]
105+
Fmv, Dmv = params[1] , params[2]
106+
#when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
107+
if Fmv < 1e-4:
108+
Dmv = float("NaN")
96109

97-
else:
98-
boundsupdated = ([bounds[0][0] , Dpar1 , bounds[0][2] , bounds[0][3] ],
99-
[bounds[1][0] , Dpar1, bounds[1][2] , bounds[1][3] ])
100-
params, _ = curve_fit(two_exp, bvalues, dw_data, p0=[1, Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
101-
S0 = params[0]
102-
Dpar, Fmv, Dmv = params[1] , params[2] , params[3]
103-
#when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
104-
if Fmv < 1e-4:
105-
Dmv = float("NaN")
106-
107-
if S0_output:
108-
return Dpar, Fmv, Dmv, S0
109-
else:
110-
return Dpar, Fmv, Dmv
111-
except:
110+
if S0_output:
111+
return Dpar, Fmv, Dmv, S0
112+
else:
113+
return Dpar, Fmv, Dmv
114+
#except:
112115

113-
if S0_output:
114-
return 0, 0, 0, 0, 0, 0
115-
else:
116-
return 0, 0, 0, 0, 0
116+
if S0_output:
117+
return 0, 0, 0, 0, 0, 0
118+
else:
119+
return 0, 0, 0, 0, 0
117120

118121

119122

‎src/standardized wip/ETP_SRI_LinearFitting.py renamed to ‎src/standardized/ETP_SRI_LinearFitting.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
4646
# Check the inputs
4747

4848

49-
def ivim_fit(self, signals, bvalues=None, linear_fit_option=False):
49+
def ivim_fit(self, signals, bvalues=None, linear_fit_option=False, **kwargs):
5050
"""Perform the IVIM fit
5151
5252
Args:

0 commit comments

Comments
 (0)