Feat: Upgrade Vanderbilt code to Python 3 and add unit tests (Fixes #102)

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')],
-)
-
+)

src/original/DS_BW_VanderBiltUMC_USA/dcemri.py

@@ -1,12 +1,12 @@
-
 # Pyhton module for DCE-MRI postprocessing 
 #
 # Copyright (C) 2014   David S. Smith
 # 
 
 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('converting DCE signal to effective R1')
     assert(flip > 0.0)
     assert(TR > 0.0 and TR < 1.0)
     S = S.T
@@ -78,7 +79,7 @@ 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('converting effective R1 to tracer tissue concentration')
     assert(relaxivity > 0.0)
     return (R1eff - R1map) / relaxivity
 
@@ -92,11 +93,13 @@ 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],
+            # FIXED: cumtrapz -> cumulative_trapezoid
+            tmp = cumulative_trapezoid(exp(-Kt*(t[k] - t[:k+1])/ve)*Cp[:k+1],
                            t[:k+1], initial=0.0) + vp * Cp[:k+1]
             Ct[k] = tmp[-1]
         else:
-            Ct[k] = simps(exp(-Kt*(t[k] - t[:k+1])/ve)*Cp[:k+1],
+            # FIXED: simps -> simpson
+            Ct[k] = simpson(exp(-Kt*(t[k] - t[:k+1])/ve)*Cp[:k+1],
                           t[:k+1]) + vp * Cp[:k+1]
     return Ct*Kt
 
@@ -108,13 +111,13 @@ 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 -> cumulative_trapezoid
+            tmp = cumulative_trapezoid(exp(-(Kt/ve)*(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(-(Kt/ve)*(t[k] - t[:k+1]))*Cp[:k+1], 
                           x=t[:k+1])
     return Kt*Ct
 
@@ -131,7 +134,7 @@ 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)
@@ -143,22 +146,22 @@ def fit_tofts_model(Ct, Cp, t, idxs=None, extended=False,
 
     # choose model and initialize fit parameters with reasonable values
     if extended:  # add vp if using Extended Tofts
-        print 'using Extended Tofts-Kety'
+        print('using Extended Tofts-Kety')
         fit_func = lambda t, Kt, ve, vp: \
                     ext_tofts_integral(t, Cp, Kt=Kt, ve=ve, vp=vp)
         coef0 = [0.01, 0.01, 0.01]
         popt_default = [-1,-1,-1]
         pcov_default = ones((3,3))
     else:
-        print 'using Standard Tofts-Kety'
+        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)
         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:
@@ -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

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, Kt=Kt, ve=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, Kt=Kt, ve=ve, vp=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