-
Notifications
You must be signed in to change notification settings - Fork 5
Expand file tree
/
Copy pathexample_script_ROI.m
More file actions
60 lines (50 loc) · 7.2 KB
/
Copy pathexample_script_ROI.m
File metadata and controls
60 lines (50 loc) · 7.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
%%Example script demonstrating how to fit a single timeseries of
%%DCE-MRI data using the Patlak model
%%Thanks to LK for generating and supplying test data
%%Parameters used to generate data were: PS = 0.05E-4 /s = 0.0003 /min, vP = 0.58e-2
clear; close all;
%% set up paths
procRoot='PathToMyCode'; %set this first
addpath([procRoot '/utilities']); %available at https://github.com/mjt320/utilities
addpath([procRoot '/DCE-functions']); %available at https://github.com/mjt320/DCE-functions
%% set acquisition parameters
TR_s=3e-3; %repetition time for spoiled gradient echo sequence
TE_s=0; %assume no transverse relaxation/dephasing
FA_deg=18; %flip angle
tRes_s=2.45; % time resolution
baselineScansIdx=1:14; % baseline scans to use for calculating enhancement
r1_permMperS=5.0; % relaxivity (/mMol/s) [gadovist at 3T]
r2s_permMperS=7.1; % relaxivity (/mMol/s) [gadovist at 3T]
T1_blood_s=1.68; % T1 of blood
T1_tissue_s=1; % T1 of tissue
Hct=0.4; % Hematocrit
%% Define column vector containing tissue and blood (VIF) signal intensity time series - this data
%% would normally be imported from a spreadsheet, read from a 4D image etc.
%Note that the DCE functions can operate on multiple time series if a
%multiple input column vectors are supplied as a 2D array
SI_tissue=[17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87254874 17.87313623 17.90735479 18.16327917 18.41037724 18.62693326 18.53816458 18.40988807 18.22729971 18.13460484 18.11104752 18.13418083 18.18191879 18.19378996 18.22781065 18.21045065 18.21368116 18.19097807 18.20445619 18.20006571 18.19076286 18.1914048 18.18212924 18.19524063 18.17929375 18.18731594 18.16815662 18.16377537 18.14605942 18.16009997 18.14144496 18.1516135 18.13505501 18.15478613 18.14569318 18.15467686 18.14911615 18.15564505 18.16484117 18.17029278 18.15896258 18.15086419 18.12529982 18.14009696 18.12978774 18.14426603 18.14114782 18.14763281 18.14935753 18.14507408 18.1360725 18.13299657 18.12992933 18.13395744 18.12509051 18.12674427 18.12602197 18.13819957 18.13632663 18.12863281 18.11969944 18.12250892 18.11722618 18.13065073 18.13111811 18.14080588 18.13533047 18.12531328 18.1269384 18.11576051 18.12904492 18.13421017 18.15590241 18.14934496 18.13816224 18.11856372 18.13939032 18.13530203 18.14389414 18.13704417 18.14904231 18.14485666 18.12994185 18.13749166 18.14261263 18.13963102 18.13771724 18.13817679 18.15388317 18.14385213 18.14549448 18.13534542 18.12111658 18.10993487 18.12697988 18.142087 18.13243097 18.11387361 18.10742936 18.12720187 18.1253138 18.12799876 18.10905239 18.11406403 18.11105294 18.12616817 18.10622044 18.11573599 18.11267627 18.12536401 18.12004225 18.11816224 18.11855833 18.12239207 18.12621666 18.1381937 18.13180288 18.1265215 18.11309795 18.1089967 18.13634598 18.13905339 18.15541581 18.13116958 18.12480547 18.11715598 18.10194511 18.1067912 18.09815531 18.10955838 18.12375356 18.13087385 18.12130958 18.1089864 18.10822049 18.12649817 18.13375045].';
SI_blood=[20,20,20,20,20,20,20,20,20,20,20,20,20,20,20.1156370399779,26.8351271371286,72.4441185804455,109.872368951734,138.392795717265,126.856781601193,109.223004849012,81.7976798687703,66.6046868109294,62.5353485112090,66.3609328557451,74.1575409695594,75.9769597194646,81.3096787165081,78.4339049637436,78.8398578248481,75.0738019174292,77.1447706207576,76.3294335758057,74.7177610526692,74.7177610526692,73.1011302646683,75.1336678596749,72.4300377025091,73.6426665907263,70.3955937509619,69.5747658178580,66.5283485045250,68.7803644531634,65.5670177290809,67.1837184176967,64.3102001228117,67.5391612320756,65.9246962602292,67.3401800176306,66.3166593715396,67.3204474142840,68.7607421887302,69.5708153367646,67.5928997286576,66.1460320259733,61.7287227100751,64.1625337381575,62.3245410863340,64.6990936411890,64.0851951911800,65.0944216088783,65.2951432227209,64.4874034261484,62.8747467145751,62.2684056859961,61.6626004920833,62.2664653877994,60.6710035846460,60.8732152675424,60.6703528374405,62.6695563094047,62.2684056859961,60.8707006849342,59.2563234529495,59.6621365575527,58.6735864773334,60.9064694196661,60.9064694196661,62.4803980364942,61.4642170724488,59.6659359328896,59.8669885905184,57.8590096285142,60.0775158759515,60.8859739334789,64.4939571253345,63.2970835079815,61.3091031118729,57.8551838019332,61.3648326216612,60.5831661068416,61.9735406178327,60.7210732191761,62.6884403193924,61.8932507581731,59.2537622208056,60.4732488629507,61.2727841324351,60.6800124514357,60.2718604920606,60.2718604920606,62.8760528898353,61.0817957750675,61.2830427965444,59.4602509515272,56.9205138751393,54.8962390104449,57.7990643095778,60.3303281524101,58.5861202295109,55.2865540827300,54.0855623885753,57.4694575659846,57.0692681746827,57.4631919982221,54.0792156989484,54.8906900290170,54.2947025244965,56.8678797567005,53.3035896029795,54.9084637727523,54.3030270548558,56.4548012203479,55.4540593018196,55.0545987988007,55.0545987988007,55.6545095375412,56.2530555306938,58.2605729910051,57.0768710447764,56.0847983068167,53.6614761156638,52.8715357651283,57.5929742779859,57.9861125960642,60.7305649056217,56.4611508835147,55.2805026216724,53.8676137466384,51.1116633350847,51.9093524891175,50.3115221203626,52.2765826692529,54.7125992259831,55.8884222492454,54.1418937445661,51.9013585036236,51.7012176176053,54.8612399475604,56.0586933127365].';
%% Set processing settings
NIgnore=0; % how many points to exclude from fitting (e.g. to reduce CBF effects, or avoid numerical instability in Patlak plots)
NFrames=size(SI_tissue,1); %number of time frames
t_s=(((1:NFrames)-1)*tRes_s).'; %time in seconds corresponding to each acquisition
%% DO THE PROCESSING...
%determine tissue Gd concentration
enh_tissue_pct=DCEFunc_Sig2Enh(SI_tissue,baselineScansIdx); %calculate tissue enhancement from signal intensity
conc_tissue_mM=DCEFunc_Enh2Conc_SPGR(enh_tissue_pct,T1_tissue_s,TR_s,TE_s,FA_deg,r1_permMperS,r2s_permMperS); %calculate tissue Gd concentration from signal enhancement
%determine blood plasma Gd concentration
enh_blood_pct=DCEFunc_Sig2Enh(SI_blood,baselineScansIdx); %calculate blood enhancement from signal intensity
conc_blood_mM=DCEFunc_Enh2Conc_SPGR(enh_blood_pct,T1_blood_s,TR_s,TE_s,FA_deg,r1_permMperS,r2s_permMperS); %calculate blood Gd concentration from signal enhancement
conc_plasma_mM=(1/(1-Hct))*conc_blood_mM; %use Hct to calculate plasma concentration
%fit to Patlak model
[Patlak_params, CtModelFitFast_mM]=... % Patlak fit using multiple linear regression (Patlak plot analysis also available)
DCEFunc_fitModel(tRes_s,conc_tissue_mM,conc_plasma_mM,'PatlakFast',struct('NIgnore',NIgnore));
%% Display results
figure(1)
subplot(3,2,1), plot(t_s,SI_tissue,'k.'); title('Tissue SI');
subplot(3,2,3), plot(t_s,enh_tissue_pct,'k.'); title('Tissue enh (%)');
subplot(3,2,2), plot(t_s,SI_blood,'k.'); title('Blood SI');
subplot(3,2,4), plot(t_s,enh_blood_pct,'k.'); title('Blood enhancement (%)');
subplot(3,2,6), plot(t_s,conc_plasma_mM,'k.'); title('Plasma conc (mM)'); xlabel('t (s)');
subplot(3,2,5), plot(t_s,conc_tissue_mM,'k.',t_s,CtModelFitFast_mM,'b-'); title('Tissue conc (points) + Patlak fit (line)'); xlabel('t (s)');
text(100,0.0075,{['vP = ' num2str(Patlak_params.vP)]; ['PS = ' num2str(Patlak_params.PS_perMin) ' per min']});