create ETM non_linear fitting

Added fitting for 1 voxel and for the whole image

src/osipi/DRO/Model.py

@@ -1,5 +1,7 @@
 import numpy as np
+from scipy import integrate
 from scipy.integrate import cumtrapz, trapz
+from scipy.optimize import curve_fit
 
 
 def modifiedToftsMurase(Cp, Ctiss, dt, datashape):
@@ -95,55 +97,66 @@ def tofts(cp, c_tiss, dt, datashape):
     return ktrans, kep
 
 
-def exntended_tofts(cp, c_tiss, dt, datashape):
+# Curve fitting function for a single voxel's time series data
+
+
+def Extended_Tofts_Integral(t, Cp, Kt=0.1, ve=0.1, vp=0.2, uniform_sampling=True):
+    nt = len(t)
+    Ct = np.zeros(nt)
+    for k in range(nt):
+        tmp = vp * Cp[: k + 1] + integrate.cumtrapz(
+            np.exp(-Kt * (t[k] - t[: k + 1]) / ve) * Cp[: k + 1], t[: k + 1], initial=0.0
+        )
+        Ct[k] = tmp[-1]
+    return Ct
+
+
+def FIT_single_voxel(ct, ca, time):
+    def fit_func(t, kt, ve, vp):
+        return Extended_Tofts_Integral(t, ca, Kt=kt, ve=ve, vp=vp)
+
+    ini = [0.1, 0.1, 0.2]  # Initial guess for [Kt, ve, vp]
+    popt, pcov = curve_fit(fit_func, time, ct, p0=ini)
+    return popt, pcov
+
+
+def extended_tofts_model(ca, c_tiss, t, datashape):
     """
     Extended Tofts model for DCE-MRI DRO
-
+    ca -- arterial input function
+    c_tiss -- 4D array of tissue concentration data (x, y, z, time)
+    dt -- time interval between samples
+    datashape -- shape of the data
     """
     ktrans = np.zeros(c_tiss.shape[:-1])
-    kep = np.zeros(c_tiss.shape[:-1])
+    ve = np.zeros(c_tiss.shape[:-1])
     vp = np.zeros(c_tiss.shape[:-1])
 
     for k in range(0, datashape[2]):
+        print(f"Processing slice {k+1}/{datashape[2]}")
         for j in range(0, datashape[0]):
+            print(f"Processing row {j+1}/{datashape[0]}")
             for i in range(0, datashape[1]):
-                c = c_tiss[j, i, k, :]
-                cp_integrated = cumtrapz(cp, dx=dt, initial=0)
-                c_tiss_integrated = -cumtrapz(c_tiss[i, j, k, :], dx=dt, initial=0)
-                A = np.column_stack((cp_integrated, c_tiss_integrated, cp))
-                ktrans_voxel, kep_voxel, vp_voxel = np.linalg.lstsq(A, c, rcond=None)[0]
-                ktrans[j, i, k] = ktrans_voxel
-                kep[j, i, k] = kep_voxel
-                vp[j, i, k] = vp_voxel
+                ct = c_tiss[j, i, k, :]
+                popt, _ = FIT_single_voxel(ct, ca, t)
+                ktrans[j, i, k], ve[j, i, k], vp[j, i, k] = popt
 
-    return ktrans, kep, vp
+    return ktrans, ve, vp
 
 
-def forward_tofts(ktrans, kep, cp, vp, dt):
+def extended_tofts_model_1vox(ca, c_tiss, t):
+    """
+    Extended Tofts model for DCE-MRI DRO
+    ca -- arterial input function
+    c_tiss -- 1D array of tissue concentration data (time)
+    dt -- time interval between samples
     """
-    Forward Tofts model for DCE-MRI DRO
 
-    Parameters:
-    ktrans (numpy.ndarray): The transfer constant Ktrans.
-    kep (numpy.ndarray): The rate constant kep.
-    cp (numpy.ndarray): The plasma concentration C_p(t).
-    vp (numpy.ndarray): The plasma volume fraction v_p.
-    dt (float): The time step between measurements.
+    ct = c_tiss[:]
+    popt, _ = FIT_single_voxel(ct, ca, t)
+    ktrans, ve, vp = popt
 
-    Returns:
-    numpy.ndarray: The tissue concentration C_tiss(t).
-    """
-    time_points = cp.shape[-1]
-    c_tiss = np.zeros(ktrans.shape)
-    for t in range(time_points):
-        if t == 0:
-            c_tiss[..., t] = vp * cp[..., t]
-        else:
-            exp = np.exp(-kep * np.arange(t + 1)[::-1] * dt)
-            integral = trapz(cp[..., : t + 1] * exp, dx=dt, axis=-1)
-            c_tiss[..., t] = vp * cp[..., t] + ktrans * integral
-
-    return c_tiss
+    return ktrans, ve, vp
 
 
 def ForwardsModTofts(K1, k2, Vp, Cp, dt):

src/osipi/DRO/main.py

@@ -1,16 +1,27 @@
+import time
+
 import numpy as np
 from DICOM_processing import (
     SignalEnhancementExtract,
-    read_dicom_slices_as_signal,
+    read_dicom_slices_as_4d_signal,
 )
 from Display import animate_mri
+from matplotlib import pyplot as plt
 from roi_selection import ICfromROI
 
-from osipi.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref
+import osipi
+from osipi.DRO.Conc2DROSignal import Conc2Sig_Linear, STDmap, addnoise, calcR1_R2, createR10_withref
 from osipi.DRO.filters_and_noise import median_filter
-from osipi.DRO.Model import ForwardsModTofts, modifiedToftsMurase, modifiedToftsMurase1Vox
+from osipi.DRO.Model import (
+    ForwardsModTofts,
+    extended_tofts_model_1vox,
+    modifiedToftsMurase,
+    modifiedToftsMurase1Vox,
+)
 
-signal, slices, dicom_ref = read_dicom_slices_as_signal("data/subject2/12.000000-perfusion-17557")
+signal, slices, dicom_ref = read_dicom_slices_as_4d_signal(
+    "data/subject2/12.000000-perfusion-17557"
+)
 # anim = animate_mri(slices, mode="time", slice_index=7, time_index=5)
 data_shape = signal.shape
 
@@ -34,16 +45,49 @@
 Ev = ICfromROI(E, sagmask, roivoxv, sagmask.ndim - 1)
 S0ref = ICfromROI(S0[:, :, z], sagmask[:, :, z, 0], roivoxv, sagmask.ndim - 2)
 
+max_index = np.unravel_index(np.argmax(E * aifmask, axis=None), E.shape)
+print(max_index)
+
 Hct = 0.45
 
-K1, k2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape)
+E_vox = E[max_index[0], max_index[1], max_index[2], :]
+
+start_time = time.time()
+
+k1, ve, vp = extended_tofts_model_1vox(Ea, E_vox, t)
+
+end_time = time.time()
+
+execution_time = end_time - start_time
+
+print(f"The execution time of the line is: {execution_time} seconds")
+
+K1, K2, Vp = modifiedToftsMurase(Ea / (1 - Hct), E, dt, data_shape)
+
+k1_vox = K1[max_index[0], max_index[1], max_index[2]]
+k2_vox = K2[max_index[0], max_index[1], max_index[2]]
+Vp_vox = Vp[max_index[0], max_index[1], max_index[2]]
+
+ct = osipi.extended_tofts(t, Ea, k1_vox, k1_vox / k2_vox, Vp_vox)
+ct_real = E[max_index[0], max_index[1], max_index[2], :]
+ct_extended = osipi.extended_tofts(t, Ea, k1, ve, vp)
+
+plt.figure(figsize=(10, 6))
+plt.plot(t, ct, label="ct_Mudified_Tofts")
+plt.plot(t, ct_real, label="ct_Raw")
+plt.plot(t, ct_extended, label="ct_Extended_Tofts")
+plt.xlabel("Time")
+plt.ylabel("Concentration")
+plt.title("ct_raw vs model")
+plt.legend()
+plt.show()
 
 pvc_K1, pvc_k2, pvc_Vp = modifiedToftsMurase1Vox(Ea / (1 - Hct), Ev, dt, data_shape)
 pvc = abs((1 - Hct) / pvc_Vp)  # sagittal sinus Vp should be 0.55, so calc correction factor if not
 
 # Apply correction factor to fitted parameters
 cor_K1 = K1 * pvc
-cor_k2 = k2 * pvc
+cor_k2 = K2 * pvc
 cor_Vp = Vp * pvc
 # Apply Median Filter to parameters all with footprint (3,3)