diff --git a/setup.py b/setup.py index 513933f0..8d7c0cee 100644 --- a/setup.py +++ b/setup.py @@ -3,12 +3,24 @@ setup( name='osipi_code_collection', version='0.0.2', - install_requires=['dicom==0.9.9.post1', 'imageio==2.19.3', 'joblib==1.1.0', 'lmfit==1.0.3', - 'mat73==0.59', 'matplotlib==3.5.2', 'nibabel==3.2.2', 'numpy==1.23.0rc2', - 'opencv-python==4.5.5.64', 'openpyxl==3.2.0b1', 'pandas==1.4.2', - 'progressbar==2.5', 'pytest==7.1.2', 'scipy==1.8.1'], + install_requires=[ + 'dicom==0.9.9.post1', + 'imageio>=2.19.3', + 'joblib>=1.1.0', + 'lmfit>=1.0.3', + 'mat73>=0.59', + 'matplotlib>=3.5.2', + 'nibabel>=3.2.2', + 'numpy>=1.23.0', + 'opencv-python>=4.5.5.64', + 'openpyxl>=3.2.0b1', + 'pandas>=1.4.2', + 'progressbar>=2.5', + 'pytest>=7.1.2', + 'scipy>=1.8.1', + 'scikit-image>=0.19.0' # Added this dependency + ], include_package_data=True, package_dir={'osipi_code_collection': 'src'}, packages=[f'osipi_code_collection.{module}' for module in find_packages('src')], -) - +) \ No newline at end of file diff --git a/src/original/DS_BW_VanderBiltUMC_USA/dcemri.py b/src/original/DS_BW_VanderBiltUMC_USA/dcemri.py index ddb7febc..9f24d157 100644 --- a/src/original/DS_BW_VanderBiltUMC_USA/dcemri.py +++ b/src/original/DS_BW_VanderBiltUMC_USA/dcemri.py @@ -1,4 +1,3 @@ - # Pyhton module for DCE-MRI postprocessing # # Copyright (C) 2014 David S. Smith @@ -6,7 +5,8 @@ import time from pylab import * -from scipy.integrate import cumtrapz, simps +# FIXED: Updated deprecated imports for SciPy 1.14+ +from scipy.integrate import cumulative_trapezoid, simpson from scipy.optimize import curve_fit @@ -17,17 +17,18 @@ def status_check(k, N, tstart, nupdates=10): telapsed = time.time() - tstart ttotal = telapsed * 100.0 / pct_complete trem = ttotal - telapsed - print '%.0f%% complete, %d of %d s remain' % \ - (pct_complete, trem, ttotal) + print('%.0f%% complete, %d of %d s remain' % \ + (pct_complete, trem, ttotal)) if k == N - 1: - print '%d s elapsed' % (time.time() - tstart) + print('%d s elapsed' % (time.time() - tstart)) def signal_to_noise_ratio(im1, im2, mask=None, thresh=None): ''' Compute SNR of two images (see Dietrich et al. 2007, JMRI, 26, 375) ''' - print 'computing signal-to-noise ratio' - from skimage.filter import threshold_otsu + print('computing signal-to-noise ratio') + # FIXED: updated from skimage.filter to skimage.filters for Python 3 + from skimage.filters import threshold_otsu if mask is None: if thresh is None: thresh = threshold_otsu(im1) @@ -38,7 +39,7 @@ def signal_to_noise_ratio(im1, im2, mask=None, thresh=None): def signal_enhancement_ratio(data, thresh=0.01): ''' Compute max signal enhancement ratio for dynamic data ''' - print 'computing signal enhancement ratios' + print('computing signal enhancement ratios') assert(thresh > 0.0) ndyn = data.shape[-1] image_shape = data.shape[:-1] @@ -54,7 +55,7 @@ def signal_enhancement_ratio(data, thresh=0.01): def dce_to_r1eff(S, S0, R1, TR, flip): - print 'converting DCE signal to effective R1' + print('conv_erting DCE signal to effectiv_e R1') assert(flip > 0.0) assert(TR > 0.0 and TR < 1.0) S = S.T @@ -69,7 +70,7 @@ def dce_to_r1eff(S, S0, R1, TR, flip): def dce_to_r1eff_old(S, S0map, idxs, TR, flip): - ''' Convert DCE signal to effective R1, based on the FLASH signal equation ''' + ''' Conv_ert DCE signal to effectiv_e R1, based on the FLASH signal equation ''' T = zeros_like(S) T[idxs,:] = (S[idxs,:].T / S0map.flat[idxs] / sin(flip)).T # normalize by pre-contrast signal R1 = zeros_like(T) @@ -78,12 +79,12 @@ def dce_to_r1eff_old(S, S0map, idxs, TR, flip): def r1eff_to_conc(R1eff, R1map, relaxivity): - print 'converting effective R1 to tracer tissue concentration' + print('conv_erting effectiv_e R1 to tracer tissue concentration') assert(relaxivity > 0.0) return (R1eff - R1map) / relaxivity -def ext_tofts_integral(t, Cp, Kt=0.1, ve=0.2, vp=0.1, +def ext_tofts_integral(t, Cp, Ktrans=0.1, v_e=0.2, v_p=0.1, uniform_sampling=True): """ Extended Tofts Model, with time t in min. Works when t_dce = t_aif only and t is uniformly spaced. @@ -92,15 +93,17 @@ def ext_tofts_integral(t, Cp, Kt=0.1, ve=0.2, vp=0.1, Ct = zeros(nt) for k in range(nt): if uniform_sampling: - tmp = cumtrapz(exp(-Kt*(t[k] - t[:k+1])/ve)*Cp[:k+1], - t[:k+1], initial=0.0) + vp * Cp[:k+1] + # FIXED: cumtrapz -> cumulativ_e_trapezoid + tmp = cumulative_trapezoid(exp(-Ktrans*(t[k] - t[:k+1])/v_e)*Cp[:k+1], + t[:k+1], initial=0.0) + v_p * Cp[:k+1] Ct[k] = tmp[-1] else: - Ct[k] = simps(exp(-Kt*(t[k] - t[:k+1])/ve)*Cp[:k+1], - t[:k+1]) + vp * Cp[:k+1] - return Ct*Kt + # FIXED: simps -> simpson + Ct[k] = simpson(exp(-Ktrans*(t[k] - t[:k+1])/v_e)*Cp[:k+1], + t[:k+1]) + v_p * Cp[:k+1] + return Ct*Ktrans -def tofts_integral(t, Cp, Kt=0.1, ve=0.2, uniform_sampling=True): +def tofts_integral(t, Cp, Ktrans=0.1, v_e=0.2, uniform_sampling=True): ''' Standard Tofts Model, with time t in min. Current works only when AIF and DCE data are sampled on same grid. ''' @@ -108,20 +111,20 @@ def tofts_integral(t, Cp, Kt=0.1, ve=0.2, uniform_sampling=True): Ct = zeros(nt) for k in range(nt): if uniform_sampling: - tmp = cumtrapz(exp(-(Kt/ve)*(t[k] - t[:k+1]))*Cp[:k+1], + # FIXED: cumtrapz -> cumulativ_e_trapezoid + tmp = cumulative_trapezoid(exp(-(Ktrans/v_e)*(t[k] - t[:k+1]))*Cp[:k+1], t[:k+1], initial=0.0) Ct[k] = tmp[-1] - #Ct[k] = simps(exp(-(Kt/ve)*(t[k] - t[:k+1]))*Cp[:k+1], - # dx=t[1]-t[0]) else: - Ct[k] = simps(exp(-(Kt/ve)*(t[k] - t[:k+1]))*Cp[:k+1], + # FIXED: simps -> simpson + Ct[k] = simpson(exp(-(Ktrans/v_e)*(t[k] - t[:k+1]))*Cp[:k+1], x=t[:k+1]) - return Kt*Ct + return Ktrans*Ct def fit_tofts_model(Ct, Cp, t, idxs=None, extended=False, plot_each_fit=False): - ''' Solve tissue model for each voxel and return parameter maps. + ''' Solv_e tissue model for each voxel and return parameter maps. Ct: tissue concentration of CA, expected to be N x Ndyn @@ -131,51 +134,51 @@ def fit_tofts_model(Ct, Cp, t, idxs=None, extended=False, idxs: indices of ROI to fit ''' - print 'fitting perfusion parameters' + print('fitting perfusion parameters') N, ndyn = Ct.shape - Kt = zeros(N) - ve = zeros(N) - Kt_cov = zeros(N) - ve_cov = zeros(N) + Ktrans = zeros(N) + v_e = zeros(N) + Ktrans_cov = zeros(N) + v_e_cov = zeros(N) if idxs is None: idxs = range(N) # choose model and initialize fit parameters with reasonable values - if extended: # add vp if using Extended Tofts - print 'using Extended Tofts-Kety' - fit_func = lambda t, Kt, ve, vp: \ - ext_tofts_integral(t, Cp, Kt=Kt, ve=ve, vp=vp) + if extended: # add v_p if using Extended Tofts + print('using Extended Tofts-Kety') + fit_func = lambda t, Ktrans, v_e, v_p: \ + ext_tofts_integral(t, Cp, Ktrans=Ktrans, v_e=v_e, v_p=v_p) coef0 = [0.01, 0.01, 0.01] popt_default = [-1,-1,-1] pcov_default = ones((3,3)) else: - print 'using Standard Tofts-Kety' - vp = zeros(N) - vp_cov= zeros(N) - fit_func = lambda t, Kt, ve: tofts_integral(t, Cp, Kt=Kt, ve=ve) + print('using Standard Tofts-Kety') + v_p = zeros(N) + v_p_cov= zeros(N) + fit_func = lambda t, Ktrans, v_e: tofts_integral(t, Cp, Ktrans=Ktrans, v_e=v_e) coef0 = [0.01, 0.01] popt_default = [-1,-1] pcov_default = ones((2,2)) - print 'fitting %d voxels' % len(idxs) + print('fitting %d voxels' % len(idxs)) tstart = time.time() for k, idx in enumerate(idxs): try: - popt, pcov = curve_fit(fit_func, t, Ct[idx,:], p0=coef0) + popt, pcov = curv_e_fit(fit_func, t, Ct[idx,:], p0=coef0) except RuntimeError: popt = popt_default pcov = pcov_default - Kt[idx] = popt[0] - ve[idx] = popt[1] + Ktrans[idx] = popt[0] + v_e[idx] = popt[1] try: - Kt_cov[idx] = pcov[0,0] - ve_cov[idx] = pcov[1,1] + Ktrans_cov[idx] = pcov[0,0] + v_e_cov[idx] = pcov[1,1] except TypeError: None #print idx, popt, pcov if extended: - vp[idx] = popt[2] - vp_cov[idx] = pcov[2,2] + v_p[idx] = popt[2] + v_p_cov[idx] = pcov[2,2] if plot_each_fit: figure(1) clf() @@ -186,11 +189,11 @@ def fit_tofts_model(Ct, Cp, t, idxs=None, extended=False, status_check(k, len(idxs), tstart=tstart) # bundle parameters for return - params = [Kt, ve] - stds = [sqrt(Kt_cov), sqrt(ve_cov)] + params = [Ktrans, v_e] + stds = [sqrt(Ktrans_cov), sqrt(v_e_cov)] if extended: - params.append(vp) - stds.append(sqrt(vp_cov)) + params.append(v_p) + stds.append(sqrt(v_p_cov)) return (params, stds) @@ -215,7 +218,7 @@ def t1_signal_eqn(x, M0, R1): for j in range(n): if images[j,:].mean() > 0.1: try: - popt, pcov = curve_fit(t1_signal_eqn, flip_angles, + popt, pcov = curv_e_fit(t1_signal_eqn, flip_angles, images[j,:].copy()) except RuntimeError: popt = [0, 0] @@ -236,5 +239,4 @@ def process(dcefile, t1file, t1_flip, R, TE, TR, dce_flip, extended=False, plotting=False): ''' Compute perfusion parameters for a DCE-MRI data set. ''' - return None - + return None \ No newline at end of file diff --git a/test/test_DS_BW_VanderBilt.py b/test/test_DS_BW_VanderBilt.py new file mode 100644 index 00000000..49138a96 --- /dev/null +++ b/test/test_DS_BW_VanderBilt.py @@ -0,0 +1,43 @@ +import pytest +import numpy as np +from src.original.DS_BW_VanderBiltUMC_USA.dcemri import tofts_integral, ext_tofts_integral, signal_to_noise_ratio + +def test_tofts_integral_basics(): + # Define time points (0 to 10 minutes) + t = np.linspace(0, 10, 20) + # Define an idealized AIF (exponential decay) + Cp = np.exp(-t) + + # Run the model + Kt = 0.2 + ve = 0.1 + Ct = tofts_integral(t, Cp, Ktrans=Kt, v_e=ve) + + # Checks + assert Ct.shape == t.shape + assert np.all(np.isfinite(Ct)) # No NaNs or Infs + assert np.max(Ct) > 0 # Should have some signal + +def test_extended_tofts_integral_basics(): + t = np.linspace(0, 10, 20) + Cp = np.exp(-t) + + # Run extended model (adds vp) + Kt = 0.2 + ve = 0.1 + vp = 0.05 + Ct = ext_tofts_integral(t, Cp, Ktrans=Kt, v_e=ve, v_p=vp) + + assert Ct.shape == t.shape + assert np.all(np.isfinite(Ct)) + +def test_signal_to_noise_ratio_basics(): + # Create two identical images with slight noise difference + img1 = np.ones((10, 10)) + img2 = np.ones((10, 10)) + 0.1 * np.random.rand(10, 10) + + # This function prints to stdout, so we just check it runs without crashing + snr, mask = signal_to_noise_ratio(img1, img2) + + assert isinstance(snr, float) + assert mask.shape == img1.shape \ No newline at end of file