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)
Before running
Example of Parallel Computing using MDCS for COMKAT
After finishing the setting of MDCS, users could increase the number of workers to reduce the computational time of a pixel-wise parameter estimation. The maximum worker is the product of the number of computers and the number of CPUs per computer. To compare the computational time under different numbers of workers, we perform a pixel-wise estimation of FDG rate constants (k1~k4) for 128 x 128 pixels x 29 frames. This means that we estimate rate constants for 128 x 128 time-activity curves (TACs) and each TAC has 29 time frames. All data are generated from the FDG model (following the step 1 and step 2). The below example includes six different tests: 1, 2, 4, 8, 16 and 32 workers. In each test, we run five trials and calculate mean computational time from these five trials.
Step 1. This part is to create a 18F-FDG model by using COMKAT commands. 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 128 x 128 noisy TACs by adding noise to the noise-free TAC.
noise_level = 0.1; for i=1:128*128 noisy_data(:,i) = [addNoiseDefault(data,noise_level,scant)]; end
Step 3. This part is to perform parameter estimation for 128 x 128 TACs using MDCS. Here, there are six different tests: 1, 2, 4, 8, 16 and 32 workers. For each test, we run five trials (test_idx) and then calculate mean computational time from these five trials. Before fitting the noisy TACs, users must use the Matlab function 'matlabpool' to start parallel language worker pool. The value in function 'matlabpool' indicates the number of workers 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 Matlab command 'matlabpool close' to close matlabpool after each test.
%%%% Define settings specific for your environment %%%% nworkers = 16; % run the computation using 16 workers basepath = 'd:\COMKAT_R3.1'; % set this to the main comkat folder for your computer % IF YOU ARE USING FUNCTIONS IN ANY OTHER LOCATIONS THAN XXX, ENSURE TO ADD TO PATH BELOW %%%% End of environment-specific settings %%%% % run 5 trials for measuring mean and standard deviation of compute time for test_idx = 1:5 % specify the number of workers to use in this test matlabpool(nworkers); % add all required COMKAT m-files into the path % Distributed computing uses "headless" MATLAB sessions so path must be set from the command-line addpath(basepath); addpath([basepath '\utilities']); addpath([basepath '\validation']); t0 = clock; % use 'parfor' to perform parallel computing for 128x128 noisy data parfor i=1:128*128 % set the data, which will be fitted and use COMKAT's function 'fit' to fit the data cm2 = set(cm, 'ExperimentalData', noisy_data(:,i)); pfit(:,i) = fit(cm2, pinit, plb, pub); end time_consumed_parfor(test_idx) = etime(clock,t0); matlabpool close; end end