Support:Documents:Examples:Estimate Parametric Image with Matlab Distributed Computing Server
Estimation of Parametric Image using Matlab Distributed Computing Server (MDCS)
Overview
Generally, the estimation of kinetic parameters is performed by a ROI (region-of-interest)-based method. This means that time-activity curves (TACs) are generated by calculating mean activity from a user-defined ROI at each time point in a dynamic image data set. This method is simple and robust. However, it cannot represent physiological properties of the tissue, which is heterogeneous. For example, drawing a ROI in a tumor may include different tissue types, which have different biological properties. Therefore, a ROI-based method may fail to represent some significant characteristics of the tumor. One solution is to estimate kinetic parameters pixel-by-pixel method. This approach may generate several parametric images, and pixels value in each parametric image represents the value of a kinetic parameter. The generation of parametric images is time-consuming and the accuracy of parameter estimation is easily affected by image noise. To reduce the computational time, one alternate approach is to use parallel computing, which speeds up process of parameter estimation by using several computers or multiple CPUs. This example demonstrates how the MATLAB Distributed Computing Server (MDCS) can be used to reduce the time required for parameter estimation.
Setting Matlab Distributed Computing Server (MDCS)
To start the parallel computing, user should install MDCS as described in this document. Note: In order to use MDCS for COMKAT, users must have COMKAT folders with the same pathway in both the client and the cluster. Users could use funciton 'addpath' to add pathway for COMKAT folders.
Example of Parallel Computing using MDCS for COMKAT
Commonly, the image size of a PET image is 128x128. Instead of performing parameter estimation for a PET image, we use the 18F-FDG model to generate 128x128 noise-free TACs and then estimate kinetic parameters using MDCS. To compare the reduction of computational time at different numbers of workers, the number of workers is changed from 1, 2, 4, 8, 16 and 32. Also, there are 5 trials for each test.
Step 1. Create a 18F-FDG model. The basic commands can be found in the user manual and the overview of the 18F-FDG model can be found in the example.
cm = compartmentModel; % start with a new, empty model
% k1 k2 k3 k4
ktrueA=[0.1 ; 0.13 ; 0.06 ; 0.0068];
% define the parameters
cm = addParameter(cm, 'sa', 1); % specific activity of injection, kBq/pmol
cm = addParameter(cm, 'dk', log(2)/109.8); % radioactive decay
cm = addParameter(cm, 'PV', 1); % (none)
cm = addParameter(cm, 'k1', 0.1); % 1/min
cm = addParameter(cm, 'k2', 0.13); % 1/min
cm = addParameter(cm, 'k3', 0.06); % ml/(pmol*min)
cm = addParameter(cm, 'k4', 0.0068); % 1/min
% define input function parameter vector
cm = addParameter(cm, 'pin', [28; 0.75; 0.70; 4.134; 0.1191; 0.01043]);
% define compartments
cm = addCompartment(cm, 'Junk');
cm = addCompartment(cm, 'Ce' );
cm = addCompartment(cm, 'Cm' );
% define plasma input function
% specifying function as refCp with parameters pin
cm = addInput(cm, 'Cp', 'sa', 'dk', 'refCp', 'pin'); % plamsa pmol/ml
% connect inputs and compartments
cm = addLink(cm, 'L', 'Cp', 'Ce', 'k1');
cm = addLink(cm, 'K', 'Ce', 'Junk','k2');
cm = addLink(cm, 'K', 'Ce', 'Cm', 'k3');
cm = addLink(cm, 'K', 'Cm', 'Ce', 'k4');
% specify scan begin and end times
ttt=[ ones(6,1)*5/60; ... % 6 frames x 5 sec
ones(2,1)*15/60; ... % 2 frames x 15 sec
ones(6,1)*0.5;... % 6 frames x 0.5 min
ones(3,1)*2;... % 3 frames x 2 min
ones(2,1)*5;... % 2 frames x 5 min
ones(10,1)*10]; % 10 frames x 10 min
scant = [[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
cm = set(cm, 'ScanTime', scant);
% define an outputs, one for each region
cm = addOutput(cm, 'RegA', {'Ce', 'PV'; 'Cm', 'PV'}, {});
% solve model and generate example output
[PET, PETindex]=solve(cm);
data = PET(:,3); % data will have 3 columns, one for each region
% specify parameters to be adjusted in fitting
cm = addSensitivity(cm, 'pin', 'k1', 'k2', 'k3', 'k4');
% set parameter values initial guess, lower and upper bounds. values are in same order as sensitivities
% _____________pin_________________ ______Reg______
pinit = [ 10; 0.4; 0.4; 3; 0.05; 0.01; 0.1; 0.1; 0.05; 0.001; ];
plb = [ 10; 0.1; 0.1; 1; 0.05; 0.001; 1e-3; 1e-3; 1e-3 ; 1e-5];
pub = [100; 2. ; 2. ; 10; 1. ; 0.05; 1.; 1.; 1.; 1.;];
Step 2. Generate 128x128 noisy TACs by adding noise to noise-free TACs.
noise_level = 0.1;
for i=1:128*128
noisy_data(:,i) = [addNoiseDefault(data,noise_level,scant)];
end
Step 3. Perform parameter estimation for 128x128 TACs with different workers (matl_pool_n). For each test, the number of trials is 5 (test_idx). Before fitting the noisy TACs, users must use the Matlab function 'matlabpool' to start parallel language worker pool. This indicates that the number of workers is used in the test. To use parallel computing for fitting all TACs , the 'for' loop is replaced by the 'parfor' loop. Note: users should use 'command matlabpool close' to close matlabpool after each test.
for pool_idx =1:6
mat_pool_n = [32 16 8 4 2 1];
for test_idx = 1:5 % the number of trials
matlabpool(mat_pool_n(pool_idx));
addpath('D:\COMKAT_R3.1');
addpath('D:\COMKAT_R3.1\utilities');
addpath('D:\COMKAT_R3.1\validation');
t0 = clock;
for_times = 128*128;
parfor i=1:for_times
cm2 = set(cm, 'ExperimentalData', noisy_data(:,i));
pfit(:,i) = fit(cm2, pinit, plb, pub);
end
time_consumed_parfor(pool_idx,test_idx) = etime(clock,t0);
matlabpool close;
end
end
Results: The below figure is the computational time (minute) versus the number of workers. As we expected, the computational time reduces with the increase of the number of workers.
