From 5c684abf81818205ab9440f68738744545c70d21 Mon Sep 17 00:00:00 2001 From: mjt320 Date: Sat, 19 Mar 2022 01:02:53 +0000 Subject: [PATCH 1/6] add EnhToConcSPGR (analytical) --- .idea/SEPAL.iml | 2 +- .idea/misc.xml | 2 +- demo/demo_fit_dce.ipynb | 4 ++-- src/dce_fit.py | 53 +++++++++++++++++++++++++++++++++++++++++ src/signal_models.py | 4 ++-- 5 files changed, 59 insertions(+), 6 deletions(-) diff --git a/.idea/SEPAL.iml b/.idea/SEPAL.iml index 06b0085..d49a298 100644 --- a/.idea/SEPAL.iml +++ b/.idea/SEPAL.iml @@ -2,7 +2,7 @@ - + diff --git a/.idea/misc.xml b/.idea/misc.xml index 8161a60..c3334de 100644 --- a/.idea/misc.xml +++ b/.idea/misc.xml @@ -1,4 +1,4 @@ - + \ No newline at end of file diff --git a/demo/demo_fit_dce.ipynb b/demo/demo_fit_dce.ipynb index 44ce0c3..3f59d90 100644 --- a/demo/demo_fit_dce.ipynb +++ b/demo/demo_fit_dce.ipynb @@ -584,7 +584,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -598,7 +598,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.8" } }, "nbformat": 4, diff --git a/src/dce_fit.py b/src/dce_fit.py index e93a905..56ae32e 100644 --- a/src/dce_fit.py +++ b/src/dce_fit.py @@ -130,6 +130,59 @@ def proc(self, enh, t10, k_fa=1): return C_func(enh) +class EnhToConcSPGR(Fitter): + """Convert enhancement to concentration. + + Subclass of Fitter. Calculates points on the enh vs. conc curve, + interpolates and uses this to "look up" concentration values given the + enhancement values. It assumes the fast water exchange limit. + """ + + def __init__(self, tr, fa, r1): + """ + + Args: + c_to_r_model (CRModel): concentration to relaxation + relationship + signal_model (SignalModel): relaxation to signal relationship + C_min (float, optional): minimum value of concentration to look for + C_max (float, optional): maximum value of concentration to look for + n_samples (int, optional): number of points to sample the enh-conc + function, prior to interpolation + """ + self.tr = tr + self.fa = fa + self.r1 = r1 + + def output_info(self): + """Get output info. Overrides superclass method. + """ + return ('C_t', True), + + def proc(self, enh, t10, k_fa=1): + """Calculate concentration time series. Overrides superclass method. + + Args: + enh (ndarray): 1D array of enhancements (%) + t10 (float): tissue T10 (s) + k_fa (float, optional): B1 correction factor (actual/nominal flip + angle). Defaults to 1. + + Returns: + ndarray: 1D array of tissue concentrations (mM) + """ + if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa): + raise ValueError( + f'Unable to calculate concentration: nan arguments received.') + cos_fa_true = np.cos(k_fa * self.fa * np.pi/180) + exp_r10_tr = np.exp(self.tr/t10) + C_t = -np.log((exp_r10_tr * (enh-100*cos_fa_true-enh*exp_r10_tr+100)) / + (100 * exp_r10_tr + enh * cos_fa_true - 100 * exp_r10_tr * + cos_fa_true - enh * exp_r10_tr * cos_fa_true) + ) / (self.tr * self.r1) + return C_t + + class ConcToPKP(Fitter): """Fit tissue concentrations using pharmacokinetic model. diff --git a/src/signal_models.py b/src/signal_models.py index 1b16755..a900b83 100644 --- a/src/signal_models.py +++ b/src/signal_models.py @@ -64,12 +64,12 @@ def __init__(self, tr, fa, te): te (float): echo time (s) """ self.tr = tr - self.fa = fa + self.fa = fa * np.pi / 180 self.te = te def R_to_s(self, s0, R1, R2=None, R2s=0, k_fa=1): """Get signal for this model. Overrides superclass method.""" - fa = k_fa * self.fa * np.pi / 180 + fa = k_fa * self.fa s = s0 * (((1.0-np.exp(-self.tr*R1))*np.sin(fa)) / (1.0-np.exp(-self.tr*R1)*np.cos(fa)) ) * np.exp(-self.te*R2s) From 1a8f006fb7f906352025161a96f66fe5e6851b52 Mon Sep 17 00:00:00 2001 From: mjt320 Date: Sat, 19 Mar 2022 01:03:10 +0000 Subject: [PATCH 2/6] add EnhToConcSPGR (analytical) --- demo/demo_conc_to_enh.ipynb | 100 ++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 demo/demo_conc_to_enh.ipynb diff --git a/demo/demo_conc_to_enh.ipynb b/demo/demo_conc_to_enh.ipynb new file mode 100644 index 0000000..0d467fc --- /dev/null +++ b/demo/demo_conc_to_enh.ipynb @@ -0,0 +1,100 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 34, + "id": "3b2d4aa7-5b11-44fc-ad48-069caaeab47a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "import sys\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "sys.path.append('../src')\n", + "import dce_fit, relaxivity, signal_models, water_ex_models, aifs, pk_models\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "9439e538-cdd8-4849-bedf-1ead532f0811", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgsUlEQVR4nO3deXxU9b3/8dcnJIBsBkiCQAgEZJWAhIGwJqlIFdS63HtbKHqpK7YutbZWudKLV0ur7c+ttdbiWn8q6NUrUqt1QUNYIiRAhEBYhUAgECCXJYQAYb6/PzL6o8gSMpOcmeT9fDx4ZHIyk/P+csibk+85c4455xARkcgT5XUAERGpHRW4iEiEUoGLiEQoFbiISIRSgYuIRKjo+lxZXFyc69atW32uUkQk4i1btmyPcy7+5OX1WuDdunUjLy+vPlcpIhLxzKzoVMs1hSIiEqFU4CIiEUoFLiISoVTgIiIR6qwFbmYvmVmpmRWcsOz3ZrbWzFaa2btmFlunKUVE5Ftqsgf+CnD5Scs+Afo75wYA64GpIc4lIiJncdYCd85lA2UnLfvYOVcV+PQLILEOsomIyBmEYg78JuDDEHwfEZEG51BpKfMHDaLos89C/r2DKnAzexCoAl4/w3NuM7M8M8vbvXt3MKsTEYk4y++/n4z8fPavXRvy7201uaGDmXUD3nfO9T9h2WTgdmCMc66iJivz+XxO78QUkcZkTcuWxBw/zoUVFVhU7faZzWyZc8538vJafTczuxy4H/heTctbRKSxWfPqq/SrqGDHVVfVurzPpCanEc4CcoDeZlZsZjcDzwCtgU/MLN/Mngt5MhGRCLfnN7/hEHDx44/Xyfc/68WsnHMTT7H4xTrIIiLSYOzfupXB69axrHdv0pOS6mQdeiemiEgdyP/5z2kJxE+bVmfrUIGLiISY8/vp9Le/saZFC/pef32drUcFLiISYquee46eR46w9/vfr9P1qMBFRELs4O9+x34g9bHH6nQ9KnARkRAqyctjaFER+YMG0TIhoU7XpQIXEQmhdffcQxMg+Ykn6nxdKnARkRA5XFZGyuLF5F5wAUmZmXW+PhW4iEiI5N17L+2do9n999fL+lTgIiIh4Px+Lpg9m3XNmzPw7rvrZZ0qcBGREMh/8kl6HjlC6YQJdXLdk1NRgYuIhMDR3/+ePWYMefLJelunClxEJEhF8+YxZNcuCkaOpHlsbL2tVwUuIhKkop/+lGNAn6efrtf1qsBFRIKwe/VqhqxezdJevbggNbVe160CFxEJwuopU2gGdH7qqXpftwpcRKSWDu7YwcWLF7O0Y0e6jxtX7+tXgYuI1NKyH/+YWOdo9cgjnqxfBS4iUgtHy8vp/f775J9/Pv1vvtmTDCpwEZFaWHrPPXT0+zl2772eZTDnXL2tzOfzuby8vHpbn4hIXfBXVbG5ZUuqoqLodehQnb/z0syWOed8Jy/XHriIyDlaMnUqPY4eZe8tt9Tb2+ZPRQUuInIO/FVVtHvmGb5q2pS0xx/3NIsKXETkHCx98EF6V1ZSctNNNGna1NMsKnARkRpyfj+xf/wjm2NiSKvHi1adzlkL3MxeMrNSMys4YVk7M/vEzDYEPrat25giIt5bOm0afQ4fZvuPfkR08+Zex6nRHvgrwOUnLXsAmOec6wnMC3wuItJgOb+fNk89xZaYGIb94Q9exwFqUODOuWyg7KTFVwN/DTz+K3BNaGOJiISX3OnT6Xv4MMX//u9hsfcNtZ8D7+CcKwEIfEwIXSQRkfDi/H5aP/44RdHRDHvmGa/jfKPOD2Ka2W1mlmdmebt3767r1YmIhNwX991H38OH2RYmc99fq22B7zKzjgCBj6Wne6JzbqZzzuec88XHx9dydSIi3qiqrKTDM8+wsVkzhv/pT17H+Se1LfC5wOTA48nAe6GJIyISXnJuv53uR4+y92c/8/y875Od9VooZjYLyATigF3AdGAO8BaQBGwF/s05d/KBzm/RtVBEJJJU7ttHWVwcZc2bc9GBA569bf5010KJPtsLnXMTT/OlMUGnEhEJY19Mnkzm8ePsevhhT695cjrhl0hEJAwcKC4m5W9/I699ewZ5eMnYM1GBi4icwvJJk2jvHC09uNdlTanARUROUpKXx5DsbBZ36ULf66/3Os5pqcBFRE6yacIEooEur73mdZQzUoGLiJyg8LXXGLVpEzlDh9IlPd3rOGekAhcRCXB+P0fuvJPdZgz67//2Os5ZqcBFRAKWTJ3Kxfv3UzhhAucnJXkd56xU4CIiwNHycjo+8QQbmzVjxEsveR2nRlTgIiLA4kmT6FpVxb5p08LqglVnogIXkUZvV34+qXPnkhsfj2/aNK/j1JgKXEQavQ3XXUczIP71172Ock5U4CLSqH35xz8yavNmckaNotvYsV7HOScqcBFptI5VVNDivvsobtKEoe++63Wcc6YCF5FGa/GkSfQ8coTt991Hi7g4r+OcMxW4iDRKu/LzGTRnDrnx8QydMcPrOLWiAheRRmnjNdfQDEiYNSssr/VdE5GZWkQkCLkPP8zIoiJyMjLoOiZy702jAheRRqV85046Pvwwm5o2ZfjcuV7HCYoKXEQalWXjx9Pp+HEq/vAHmrVp43WcoKjARaTRWP3yy4xesYIF/fuTMmWK13GCpgIXkUbhWEUFMT/5Cbuiohj04YdexwkJFbiINAqLrrmGXpWVbH3gAdokJnodJyRU4CLS4G187z2Gf/IJOYmJpEXoOd+nogIXkQbtWEUFxyZO5IAZPT/+2Os4IaUCF5EGbdH48fQ9fJiv7r+fuL59vY4TUkEVuJn9zMxWm1mBmc0ys8i4CrqINAprXn2VUfPns7B7d9J++1uv44RcrQvczDoDdwM+51x/oAkwIVTBRESCUblvH81uvZXdUVH0//RTr+PUiWCnUKKB88wsGmgB7Ag+kohI8L4YO5YeR4+y/eGHiU1O9jpOnah1gTvntgP/B9gKlAD7nXPfOkJgZreZWZ6Z5e3evbv2SUVEaij/qadIz8sju29ffA8+6HWcOhPMFEpb4GogGegEtDSz609+nnNupnPO55zzxcfH1z6piEgNlG3YQIef/5yimBhSP/vM6zh1KpgplEuBzc653c65Y8D/ACNCE0tE5Nw5v5/1GRm09/upfPllWl1wgdeR6lQwBb4VGGZmLczMgDFAYWhiiYicuwXXX8+wkhIWf+979J00yes4dS6YOfAlwNvAcmBV4HvNDFEuEZFzsvG99xgyaxZ57duT/s47XsepF9HBvNg5Nx2YHqIsIiK1crisDPeDH1BuRtesLKKig6q2iKF3YopIxMsdPZqeR45Q9MgjxPfv73WceqMCF5GItvD220lfs4astLQGfcrgqajARSRirX/nHVL/8hdWxMYyKivL6zj1TgUuIhHpQHExMRMnciAqis7z5xPdvPFdikkFLiIRx/n9rB42jC7HjrHr6adJGDDA60ieUIGLSMSZf911DN++nYVXXMHAO+/0Oo5nVOAiElFWPPEEI997jy86diRj7lyv43hKBS4iEWNbdjZJv/gFW5s2pe8XX2BRjbvCGvfoRSRiHNyxg8rvfpcoIPqDDzg/KcnrSJ5TgYtI2PNXVbHa56P7kSN89Zvf0HXMGK8jhQUVuIiEvexLLmFYSQkLr72WwQ884HWcsKECF5GwlnPvvWQuWMCCnj1Jf/ttr+OEFRW4iIStVTNncvGTT7KydWuG5uU1+oOWJ9PfhoiEpaJ58+h4++2URkfTackSmrVp43WksKMCF5Gws3fdOvzjxmGA++AD4vr29TpSWFKBi0hYqdy3j+1DhtDx2DF2/PnPdBs71utIYUsFLiJhw19VxYqUFAYcPMiKn/2MlClTvI4U1lTgIhIWnN9Pdloaw4uLybriCoY/8YTXkcKeClxEwsL8yy4jc/ly5g8c2OivcVJTKnAR8Vz2xIlkfvopC5OTGa3TBWtMf0si4qnFd93FqNmzWdKhA2kFBY3mhsShoAIXEc/kPvwwQ555hpXnn8+ANWuIadHC60gRRQUuIp5Y+eyzXDR9OhvPO4/klSs5r107ryNFHBW4iNS7VTNnknzHHeyKiSEuN1eXhq2loArczGLN7G0zW2tmhWY2PFTBRKRhKnjxRbpOmcLemBiaL15M/EUXeR0pYgV7tOBp4B/OuX81s6aAJrBE5LTWvPoqXW65hf+NjqbZokV09Pm8jhTRal3gZtYGSAd+BOCcOwocDU0sEWloCl9/nU4/+hEHoqOJzs6m45AhXkeKeMFMoXQHdgMvm9kKM3vBzFqe/CQzu83M8swsb/fu3UGsTkQi1bo336TDDTdQHhWFff45nYdrtjUUginwaCAV+LNzbhBwCPjWrTKcczOdcz7nnC8+Pj6I1YlIJCp48UU6TJzI4ago/J9+SuKoUV5HajCCKfBioNg5tyTw+dtUF7qICABf/vGPdL3lFvY3aYLLyiIpM9PrSA1KrQvcObcT2GZmvQOLxgBrQpJKRCJe3owZ9Lz7bkqbNqXZkiXa864DwZ6FchfweuAMlK+AG4OPJCKRbsnUqVz86KNsPu884pYt0w0Z6khQBe6cywd0HpCIfGPRHXeQ9uyzFLZqRVJBAed37ep1pAZLV40RkZCZf+21ZMyZw4rYWC5cvZrWnTp5HalBU4GLSND8VVVkDx9OZl4eOZ07M6iggOaxsV7HavBU4CISlKPl5eSmpJC5ZQvzU1IYlZdHk6ZNvY7VKOhiViJSawd37GBV166M3LKFrLFjSc/PV3nXIxW4iNTKrvx8ii+8kIFlZSy48UYyP/5Yd9KpZ/rbFpFztu7NNznu85F0+DArpk9n9EsveR2pUVKBi8g5WfLgg3SeMAFzjm1vvMGQhx7yOlKjpYOYIlIjzu9n/tVXk/7++6xt0YJ2CxbQJ1VXz/CSClxEzupoeTlf+HxkrltHTufODMzPp0VcnNexGj1NoYjIGe0pLGR1ly6kr1tH1siRpG3ZovIOEypwETmtNa++ypGUFPrs28fCKVPIXLiQqGj94h4uVOAickoLbrqJ7pMn44Atr73GqOee8zqSnET/lYrIPzlaXk7O8OFkFBSwvG1buubkkNi799lfKPVOe+Ai8o1d+fms7dSJjIICsnw+BuzYQXuVd9hSgYsIAMsefZSo1FR6HDzI4rvvJjM3l+jmzb2OJWegKRSRRq6qspKFl15K+qJFfNWsGfvffJMRV1/tdSypARW4SCNWkpfHrksvJXP/fhb07MngxYt1imAE0RSKSCOVN2MGMUOHcuH+/SycMoXR69ervCOM9sBFGpkjBw6Qc+mlZObmsr55c6LfeYdR48d7HUtqQQUu0ohsfO89jk+YQGZlJdl9+zJk4ULOa9fO61hSS5pCEWkE/FVVzL/uOhKvuYZ2R46wdNo00tesUXlHOO2BizRwJXl57LjsMjLKyliakEDyvHkM7d/f61gSAtoDF2nAcu69l2ZDh9KnrIzsSZMYUlJCvMq7wdAeuEgDtLuggI2XX87w7dspaNmSlu+8Q/pll3kdS0Is6D1wM2tiZivM7P1QBBKR2nN+P4t+8hOaDBjAoO3bybr8cvrs2UOyyrtBCsUe+E+BQqBNCL6XiNTSzuXLKRo/npG7drGqVStazJpF5pVXeh1L6lBQe+BmlghcAbwQmjgicq6c38/CW26huc9Hyq5dZF19Nf327qWHyrvBC3YP/Cngl0Dr0z3BzG4DbgNISkoKcnUicqKvPvyQ/T/8IaP27ePLNm04/+23yRw71utYUk9qvQduZlcCpc65ZWd6nnNupnPO55zzxcfH13Z1InKCyn37yMrIoPP48XTft4/siRNJ2buXbirvRiWYKZSRwPfMbAswG7jEzF4LSSoROa1ljz7KzoQEMrOzyevWjaOrVpH+xhu61VkjVOsCd85Ndc4lOue6AROAz5xz14csmYj8k9KVK1nUrRuDp07FAcsfe4yRmzfrvO5GTP9li4S5o+XlLJ40iUFz5+IDsjIyGDZnDsmxsV5HE4+F5J2Yzrks55wOeYuEWO7DD1Pcvj2Zc+eyPi6OHf/4B5lZWTRXeQvaAxcJS1s++YTS669naGkpm2NiyJ0+nSEPPeR1LAkzKnCRMHKguJjl113HiNxc2gFZV1zBiNmzSW7VyutoEoZ0MSuRMHCsooLsiRM5kpREem4uS3r2pPLLL8l8/32aqrzlNFTgIh5yfj9f/PKXFMfGkj57NtvbtGHtX//K6PXrSRgwwOt4EuY0hSLikZXPPovdfz/DysvZ2KwZS3/1K4Y89BAWpf0qqRkVuEg92/T+++y59VbSdu6kJCqKBZMnM/y557iweXOvo0mEUYGL1JOtWVkU3XorIzZuJB7IGjuWoW+8wWjdCV5qSQUuUse25+Sw6aabGL52LfHAgtRU+r/2Gpl9+3odTSKcClykjpTk5bH+ppsYtmoVccDilBR6v/IKmampXkeTBkJHS0RCbHdBAfNTU2k7ZAgjVq1iSZ8+7Fm8mIyVK7lA5S0hpAIXCZEdS5Ywf+BAWqWkMHLFCnJ79mTn/PmkFxbSefhwr+NJA6QpFJEgbf7oI4rvuou0DRuIB77o2ZMuf/oTo3VtbqljKnCRWlr31lvs/cUvSNu2jQuAnAEDuPC55xitvW2pJypwkXO06i9/4fB//idDS0s5ACwYNoyLnn+eDF2XW+qZClykBqoqK8n9j/+g1fPPk1Jezl4zssaM4eLnnyczOdnreNJI6SCmyBkc3LGD+ddeS0nr1gx/8knaVFYy/9/+jeY7d5L56afEqrzFQypwkVPYsWQJWUOH4u/cmYw5cyhr0YIlDzxA4qFDZLz1Fi0TEryOKKIpFJGvOb+fghde4MBvfsPQoiISgKVdunD+f/0XA2+80et4It+iApdGr2LPHpbddx9xb75JyuHD7AcWDR7MhU8/zYiRI72OJ3JamkKRRqvos8/I8vk4kpDA6FdeIdrvJ3viRJqUlJCZl0eiylvCnPbApVHxV1Wx7Ne/hmefZfDu3XQC8hITafHLXzLgjjvoqWtxSwRRgUujUJKby7qpU+kxfz5DqqrYFRVFdkYGfZ54guG6PolEKBW4NFjHKipY/utfYy++yODSUjoCy9u2ZdsNN+CbMYNM3WtSIpwKXBqconnz2DxtGn2XLiXN72dnVBQLRowg+ZFHSL3kEq/jiYRMrQvczLoArwIXAH5gpnPu6VAFEzkXh0pLyX/oIVrMmsWgfftIBPI6dGDLzTcz+Fe/IlO3K5MGKJg98Crg58655WbWGlhmZp8459aEKJvIGfmrqvjy6ac59Oc/M3DTJkYCW6Ojybr0Unr/9rek+XxeRxSpU7UucOdcCVASeHzQzAqBzoAKXOrUpvffZ9uMGfTMzWXQ8ePsB/J79aLNnXeS8uMfkxStmUFpHELyL93MugGDgCWn+NptwG0ASUlJoVidNEJ7CgtZPW0a8f/4B/0qKugKrIiPp2jiRAZNn87odu28jihS78w5F9w3MGsFzAdmOOf+50zP9fl8Li8vL6j1SeOxf+tWVj3yCM3mzOHiPXuIAQrPO4/Syy6j3yOPEK/Lt0ojYWbLnHPfmhMMag/czGKAd4DXz1beIjVRvnMnK2fMoMnbb3Pxzp2MArZFR7MoLY3O999P32uvRfdyF6kWzFkoBrwIFDrnnghdJGlsDpeVkf/b38Ls2QwsLmYEUBIVRU5qKnE/+QkX3XgjXfQOSZFvCWYPfCRwA7DKzPIDy/7DOfdB0KmkwavYs4eVv/89x996iwFbtjAc2G1Gbv/+xE6ZQsrtt9NRByNFziiYs1AWAhbCLNLA7du8mYLHHiN67lwGlJQwDCgzY0Xv3rS6+WYG3HUXGTpfW6TGtIsjdWpXfj7rfvc7Wn70EQPKyhhF9fRIbkoKrW+4gZQ77iC9RQuvY4pEJBW4hFzRvHlsfvJJ2mVn0//gQToAm2NiWJSWRtytt9Jv8mRNj4iEgH6KJGjHKipYPXMm+157jaRVq+h+9ChdqT7lL/uSS0i8+256XHUVyToQKRJSKnCplb3r1lH45JNEffABF23bxsXAEWBV+/Zsy8yk+9130zc9Xaf8idQhFbjUiPP72fDuu+yYOZN2OTn0P3iQUUBpVBQre/Ui5ppruOinP8XXqZPXUUUaDRW4nNb/btpE4TPPcPyDD+ixaRO9jh+nF7CmRQuyMzJIuOkm+vzwhyRoPlvEE/rJk29UVVay5uWXKZs1i7jly+l36BAjgP1AYadObBwzhl733EO/1FT6eR1WRFTgjd227Gw2P/ccTbOy6FtSwgDgOLCmVSuyMzJoP3EifSdPZpjOzxYJOyrwRmZPYSEbnn+eYx99RJcNG0g+dowuQHGTJnzZpw8x48fT9847SUlO9jqqiJyFCryB219UxNq//IXKDz+k49q19KqsJA44AKxNSGDrqFF0ueUWki+7jESd5icSUVTgDUz5zp2sfeEFyufOJWH1anpXVJAGVACF7dqRlZ5O3Pe/T59JkxiqaRGRiKYCj3D7t25l/SuvcOijj2i3ahV9Dx7EBxwF1px/PgsyMmj7L/9Cn8mTGdymjddxRSSEVOARZld+Pl+9+ipH582jw4YN9Dx8mCFU36B0bcuWLBo+nNbf+x59b7mFi+PivI4rInVIBR7GnN9P0bx5bHvjDVi4kKQtW+haVUUH4BCwrm1bFgwdSpsrrqD35Mn0T0jwOrKI1CMVeBg5cuAAG956i71//ztNc3PpUVJCN7+fbsBeMzZ26MBmn4+4a66h98SJpOoqfiKNmgrcI87vZ8eSJRS99RbH5s+n3YYN9Cov5+u7PBY3acKGpCTWjhhBpx/8gO7jx5OmdzyKyAnUCPXkcFkZG2bPpuzvf6d5fj7ddu6ks99PZ+AwsL5NG3J8Ppp/5zt0/f73SfT5SPQ6tIiENRV4HTh+9ChbPv6YkrlzOf7FF8Rv2kTPigoGBL5eFB3NpqQk1vl8xF91FT3/9V8ZqOkQETlHKvAgOb+frZ9/zvY5cziak0Psxo302L+fHkAPoBzYEBvLooEDaTFmDMkTJtD1oovo6nFuEYl8KvBz4Px+ti9ezLZ33+XIokW0Wb+e7vv20dU5ugKVwIZWrViekkKTtDQuuPJKkseNY1DTpl5HF5EGSAV+Gv6qKrZlZ7Pjww85kpNDy3XrSN67l0TnSKT6jTIbW7RgZZ8+2JAhJIwfT/erriJFUyEiUk9U4FSfvvfV3/7Gnk8/xb98ObFbtpB84ABdga5Uv0lmU/PmFF54IW7wYOIuv5we115LP72zUUQ81OgKfN/mzWx+9132Z2cTvWoV8Tt20L2y8ptbfx0ENrduzYqUFKIGDaL9mDEkX3klvdu1o7eXwUVETtJgC/xYRQVFn35K6eefczQvj/M2bqTz7t0kHj/OoMBzSqKiKG7fnkWpqTRLS6PjuHEkfec7DND51iISAYJqKjO7HHgaaAK84Jx7NCSpzoHz+9mek8OOjz6iYulSYtauJW7XLrpVVnIhcCHVUyBFTZtSlJjIxn79aD16NElXXUXH/v3pWN+BRURCpNYFbmZNgD8BY4FiINfM5jrn1oQq3Mn2FBay7cMPObB4MbZ6NW2Li+laXk4ifPOml+ImTShp146cgQOJGTSIuO98h67f/S49YmPpUVfBREQ8EMwe+FBgo3PuKwAzmw1cDYS8wLMuuYR+8+eT4Pfz9fX1yszY2qYNK1JSsJQU2qankzRuHIlJSXoHo4g0CsEUeGdg2wmfFwNpJz/JzG4DbgNISkqq1YqaJCWxrnt31vTrR+thw0gcN46EAQNopzvIiEgjFkyB2ymWuW8tcG4mMBPA5/N96+s1MfqVV2rzMhGRBi2YXdhioMsJnycCO4KLIyIiNRVMgecCPc0s2cyaAhOAuaGJJSIiZ1PrKRTnXJWZ3Ql8RPVphC8551aHLJmIiJxRUOeBO+c+AD4IURYRETkHOo1DRCRCqcBFRCKUClxEJEKpwEVEIpQ5V6v31tRuZWa7gaJavjwO2BPCOF7SWMJPQxkHaCzhKpixdHXOxZ+8sF4LPBhmluec83mdIxQ0lvDTUMYBGku4qouxaApFRCRCqcBFRCJUJBX4TK8DhJDGEn4ayjhAYwlXIR9LxMyBi4jIP4ukPXARETmBClxEJEJFRIGb2eVmts7MNprZA17nORsz22Jmq8ws38zyAsvamdknZrYh8LHtCc+fGhjbOjO7zLvkYGYvmVmpmRWcsOycs5vZ4MDfwUYz+4OZneoGIF6M5SEz2x7YNvlmNj7cx2JmXczsczMrNLPVZvbTwPKI2y5nGEskbpfmZrbUzL4MjOW/Asvrb7s458L6D9WXqt0EdAeaAl8C/bzOdZbMW4C4k5b9Dngg8PgB4LHA436BMTUDkgNjbeJh9nQgFSgIJjuwFBhO9Z2bPgTGhclYHgJ+cYrnhu1YgI5AauBxa2B9IG/EbZczjCUSt4sBrQKPY4AlwLD63C6RsAf+zc2TnXNHga9vnhxprgb+Gnj8V+CaE5bPds4dcc5tBjZSPWZPOOeygbKTFp9TdjPrCLRxzuW46n+dr57wmnpzmrGcTtiOxTlX4pxbHnh8ECik+p60EbddzjCW0wnnsTjnXHng05jAH0c9bpdIKPBT3Tz5TBs8HDjgYzNbZtU3dQbo4Jwrgep/xEBCYHkkjO9cs3cOPD55ebi408xWBqZYvv71NiLGYmbdgEFU7+1F9HY5aSwQgdvFzJqYWT5QCnzinKvX7RIJBV6jmyeHmZHOuVRgHHCHmaWf4bmROL6vnS57OI/pz0AP4GKgBHg8sDzsx2JmrYB3gHuccwfO9NRTLAv3sUTkdnHOHXfOXUz1PYGHmln/Mzw95GOJhAKPuJsnO+d2BD6WAu9SPSWyK/CrEoGPpYGnR8L4zjV7ceDxycs955zbFfih8wPP8/+nq8J6LGYWQ3Xhve6c+5/A4ojcLqcaS6Rul6855/YBWcDl1ON2iYQCj6ibJ5tZSzNr/fVj4LtAAdWZJweeNhl4L/B4LjDBzJqZWTLQk+oDGuHknLIHfm08aGbDAkfT//2E13jq6x+sgGup3jYQxmMJrPdFoNA598QJX4q47XK6sUTodok3s9jA4/OAS4G11Od2qc+jtkEc7R1P9dHqTcCDXuc5S9buVB9p/hJY/XVeoD0wD9gQ+NjuhNc8GBjbOjw4W+Ok/LOo/hX2GNV7BjfXJjvgo/qHcBPwDIF3/YbBWP4vsApYGfiB6hjuYwFGUf0r9UogP/BnfCRulzOMJRK3ywBgRSBzAfCfgeX1tl30VnoRkQgVCVMoIiJyCipwEZEIpQIXEYlQKnARkQilAhcRiVAqcBGRCKUCFxGJUP8P3RUpr/2wS80AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "r1, r2 = 5.0, 0\n", + "tr, fa, te = 3.4e-3, 15, 1.7e-3\n", + "t10 = 2\n", + "k_fa = 1.4\n", + "\n", + "enh = np.arange(3000)\n", + "\n", + "c_to_r_model = relaxivity.CRLinear(r1, r2)\n", + "signal_model = signal_models.SPGR(tr, fa, te)\n", + "e_to_c_num = dce_fit.EnhToConc(c_to_r_model, signal_model)\n", + "C_t_num = e_to_c_num.proc(enh, t10, k_fa)\n", + "\n", + "e_to_c_ana = dce_fit.EnhToConcSPGR(tr, fa, r1)\n", + "C_t_ana = e_to_c_ana.proc(enh, t10, k_fa)\n", + "\n", + "plt.plot(enh, C_t_num, 'k-',\n", + " enh, C_t_ana, 'r-')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 277c0a6f528bf36ffae79278cdaad9b1a3eb2ae6 Mon Sep 17 00:00:00 2001 From: mjt320 Date: Mon, 21 Mar 2022 10:33:31 +0000 Subject: [PATCH 3/6] Linear Patlak fitting implemented --- README.md | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index bcbbd5a..4047315 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ **Please note: This library is also hosted in the [OSIPI DCE-DSC-MRI_CodeCollection repository](https://github.com/OSIPI/DCE-DSC-MRI_CodeCollection), where unit tests and perfusion code by other authors can also be found.** -Python library for simulating and fitting DCE-MRI data. It permits arbitrary combinations of pulse sequence, pharmacokinetic model, water exchange model, etc. The code is a work-in-progress, has not been extensively tested and is not recommended or approved for use. +Python library for simulating and fitting DCE-MRI data. It permits arbitrary combinations of pulse sequence, pharmacokinetic model, water exchange model, etc. The code is a work-in-progress, has not been extensively tested and is not recommended or approved for clinical use. Created 28 September 2020 @authors: Michael Thrippleton @@ -14,7 +14,7 @@ Created 28 September 2020 - Fit tissue concentration using pharmacokinetic model - Fit signal enhancement using pharmacokinetic model - Pharmacokinetic models: steady-state, Patlak, extended Tofts, Tofts, 2CXM, 2CUM -- AIFs: patient-specific (measured), Parker, bi-exponential Parker +- AIFs: including patient-specific (measured), Parker, bi-exponential Parker - Fitting free AIF time delay parameter - Relaxivity models: linear - Signal models: spoiled gradient echo @@ -22,7 +22,6 @@ Created 28 September 2020 - T1 fitting using variable flip angle method, IR-SPGR and DESPOT1-HIFI ### Not yet implemented/limitations: -- Generally untested. Not optimised for speed or robustness. - Additional pharmacokinetic models (add by inheriting from PkModel class) - Additional relaxivity models (add by inheriting from CRModel class) - Additional water exchange models, e.g. 3S2X, 2S1X (add by inheriting from WaterExModel class) @@ -30,8 +29,6 @@ Created 28 September 2020 - R2/R2* effects not included in fitting of enhancement curves (but is included for enhancement-to-concentration conversion) - Compartment-specific relaxivity parameters/models - Fitting free water exchange parameters -- Special model implementations, e.g. linear and graphical versions of Patlak model ### TODO: - fast C calculation for SPGR with r2=0 -- inversion recovery T1 measurement From db737263fc5950f92945b892b125eac08d748e09 Mon Sep 17 00:00:00 2001 From: mjt320 Date: Mon, 21 Mar 2022 10:52:28 +0000 Subject: [PATCH 4/6] EnhToConcSPGR working and documented --- src/dce_fit.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/dce_fit.py b/src/dce_fit.py index 56ae32e..0f01790 100644 --- a/src/dce_fit.py +++ b/src/dce_fit.py @@ -9,6 +9,7 @@ SigToEnh EnhToConc ConcToPKP + EnhToConcSPGR EnhToPKP PatlakLinear @@ -133,25 +134,21 @@ def proc(self, enh, t10, k_fa=1): class EnhToConcSPGR(Fitter): """Convert enhancement to concentration. - Subclass of Fitter. Calculates points on the enh vs. conc curve, - interpolates and uses this to "look up" concentration values given the - enhancement values. It assumes the fast water exchange limit. + Subclass of Fitter. Uses analytical formula for SPGR signal, + excluding T2* effects and assuming the fast water exchange limit. This + approach is faster than EnhToConc. """ def __init__(self, tr, fa, r1): """ Args: - c_to_r_model (CRModel): concentration to relaxation - relationship - signal_model (SignalModel): relaxation to signal relationship - C_min (float, optional): minimum value of concentration to look for - C_max (float, optional): maximum value of concentration to look for - n_samples (int, optional): number of points to sample the enh-conc - function, prior to interpolation + tr (float): repetition time (s) + fa (float): flip angle (deg) + r1 (float): R1 relaxivity (s^-1 mM^-1) """ self.tr = tr - self.fa = fa + self.fa = fa * np.pi/180 self.r1 = r1 def output_info(self): @@ -174,7 +171,7 @@ def proc(self, enh, t10, k_fa=1): if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa): raise ValueError( f'Unable to calculate concentration: nan arguments received.') - cos_fa_true = np.cos(k_fa * self.fa * np.pi/180) + cos_fa_true = np.cos(k_fa * self.fa) exp_r10_tr = np.exp(self.tr/t10) C_t = -np.log((exp_r10_tr * (enh-100*cos_fa_true-enh*exp_r10_tr+100)) / (100 * exp_r10_tr + enh * cos_fa_true - 100 * exp_r10_tr * From 9e76149b31ba2eee7852b6d5047a9e72a8e0e1a1 Mon Sep 17 00:00:00 2001 From: mjt320 Date: Mon, 21 Mar 2022 11:02:46 +0000 Subject: [PATCH 5/6] update --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 4047315..737efcf 100644 --- a/README.md +++ b/README.md @@ -31,4 +31,4 @@ Created 28 September 2020 - Fitting free water exchange parameters ### TODO: -- fast C calculation for SPGR with r2=0 +- inversion recovery T1 measurment From e73df8f9277a59429e29b5807df0987aba12b1fd Mon Sep 17 00:00:00 2001 From: mjt320 Date: Thu, 14 Apr 2022 14:39:25 +0100 Subject: [PATCH 6/6] document patlak linear update readme --- README.md | 3 ++- src/dce_fit.py | 31 +++++++++++++++++-------------- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index 737efcf..97e88c8 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Python library for simulating and fitting DCE-MRI data. It permits arbitrary combinations of pulse sequence, pharmacokinetic model, water exchange model, etc. The code is a work-in-progress, has not been extensively tested and is not recommended or approved for clinical use. -Created 28 September 2020 +Created 28 September 2020 @authors: Michael Thrippleton @email: m.j.thrippleton@ed.ac.uk @institution: University of Edinburgh, UK @@ -14,6 +14,7 @@ Created 28 September 2020 - Fit tissue concentration using pharmacokinetic model - Fit signal enhancement using pharmacokinetic model - Pharmacokinetic models: steady-state, Patlak, extended Tofts, Tofts, 2CXM, 2CUM +- Patlak fitting with multiple linear regression - AIFs: including patient-specific (measured), Parker, bi-exponential Parker - Fitting free AIF time delay parameter - Relaxivity models: linear diff --git a/src/dce_fit.py b/src/dce_fit.py index 0f01790..4db5286 100644 --- a/src/dce_fit.py +++ b/src/dce_fit.py @@ -374,16 +374,18 @@ class PatlakLinear(Fitter): def __init__(self, t, aif, upsample_factor=1, include=None): """ Args: - pk_model (PkModel): Pharmacokinetic model used to predict tracer - distribution. - pk_pars_0 (list, optional): list of dicts containing starting values - of pharmacokinetic parameters. If there are >1 dicts then the - optimisation will be run multiple times and the global minimum - used. - Example: [{'vp': 0.1, 'ps': 1e-3, 've': 0.5}] - Defaults to values in PkModel.typical_vals. - include (ndarray, optional): 1D float array of true/false or 1/0 - indicating which points to include in the linear regression + t (ndarray): 1D float array of times (s) at which concentration + should be calculated. Normally these are the times at which + data points were measured. The sequence of times does not + have to start at zero. + aif (aifs.AIF): AIF object to use. + upsample_factor (int, optional): The IRF and AIF are upsampled by + this factor when calculating concentration. For non-uniform + temporal resolution, the smallest time difference between time + points is divided by this number. The default is 1. + include (ndarray, optional): 1D float array of True/False or 1/0 + indicating which points to include in the linear regression. + Defaults to None, in which case all points are included. """ self.t = t self.aif = aif @@ -415,10 +417,11 @@ def proc(self, C_t): volume. Returns: - tuple: (pk_par_1, pk_par_2, ..., Ct_fit) - pk_par_i (float): fitted parameters (in the order given in - self.PkModel.parameter_names) - Ct_fit (ndarray): best-fit tissue concentration (mM). + tuple: vp, ps, Ct_fit + vp (float): blood plasma volume fraction (fraction) + ps (float): permeability-surface area product (min^-1) + Ct_fit (ndarray): 1D array of floats containing fitted tissue + concentrations (mM) """ if any(np.isnan(C_t[self.include])): raise ValueError(f'Unable to fit model: nan arguments received.')