Difference between revisions of "Support:Documents:Examples:Estimate Parametric Image with Matlab Distributed Computing Server"

From COMKAT wiki
Jump to navigation Jump to search
 
(61 intermediate revisions by 4 users not shown)
Line 3: Line 3:
 
===Overview===
 
===Overview===
  
Generally, the estimation of kinetic parameters is performed by a ROI (region-of-interest)-based method. It means that time-activity curves are generated by calculating mean activity from a user-defined ROI at each dynamic image data. 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 by a pixel-by-pixel method. It generates several parametric images, and pixel value in each parametric image represents a kinetic parameter. However, the generation of parametric images is time-consuming and the accuracy of parameter estimation is easily affected by image noise. To reduce the compuational time, one alternate approach is to use parallel computing, which speeds up processof parameter estimation by using several computers or multiple CPUs. Here, we propose an example to reduce the computational load of parameter estimation of the pixel-by-pixel method by Matlab Distributed Computing Server (MDCS).
+
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)===
 
===Setting Matlab Distributed Computing Server (MDCS)===
  
 +
Before running the example, settings for MDCS must be finished. The introduction of setting MDCS can be found in the [[Support:Documents:Manual:Distributed Computing with COMKAT]].
  
 +
===Quick test===
 +
After finishing the setting of MDCS, users could run a quick test for performing parallel computing using MDCS with COMKAT functions.
 +
Note: Since MDCS uses headless MATLAB sessions, COMKAT must be copied to all workers (or nodes), and complete path name must be the same as 'local' client. For example, in the client, we put the COMKAT_R4.1b folder in 'd:\COMKAT_R4.1b'. In all wokers, the COMKAT_R4.1b folder must be putted in the same path (i.e., 'd:\COMKAT_R4.1b'). Also, setting path for 'all required COMKAT functions' must be performed in the command-line (i.e., addpath).
  
===Example===
+
<pre>
 +
nworkers = 16;                      % run the computation using 16 workers
 +
parpool(nworkers);            % specify the number of workers to use in this test
 +
basepath = 'd:\COMKAT_R4.1b'; % set this to the main comkat folder for your computer
 +
addpath(basepath);
 +
parfor i=1:100                        % use 'parfor' to perform parallel computing
 +
cm = compartmentModel;
 +
end
 +
delete(gcp);
 +
</pre>
 +
 
 +
==Example of Parallel Computing using MDCS for COMKAT==
 +
 
 +
After finishing the above test, users could change the number of workers for reducing the computational time of a pixel-wise parameter estimation. The below example is to 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 each TAC has 29 time frames. All data are generated from the FDG model (following the step 1 and step 2).
 +
 
 +
'''Step 1.''' Create a <sup>18</sup>F-FDG model. The introduction of basic commands can be found in the [http://comkat.case.edu/comkat/comkat_wiki/index.php?title=Support:Documents:User_manual user manual] and the overview of the <sup>18</sup>F-FDG model can be found in the [http://comkat.case.edu/comkat/comkat_wiki/index.php?title=Support:Documents:Examples:FDG_with_Time-varying_Rate_Constants example].
  
 
<pre>
 
<pre>
 
cm = compartmentModel;  % start with a new, empty model
 
cm = compartmentModel;  % start with a new, empty model
  
%       k1     k2      k3     k4
+
%         k1     k2      k3     k4
 
ktrueA=[0.1 ;  0.13 ; 0.06 ; 0.0068];
 
ktrueA=[0.1 ;  0.13 ; 0.06 ; 0.0068];
 
   
 
   
 
% define the parameters
 
% define the parameters
cm = addParameter(cm, 'sa',    1);           % specific activity of injection, kBq/pmol
+
cm = addParameter(cm, 'sa',    1);                 % specific activity of injection, kBq/pmol
 
cm = addParameter(cm, 'dk',    log(2)/109.8); % radioactive decay
 
cm = addParameter(cm, 'dk',    log(2)/109.8); % radioactive decay
cm = addParameter(cm, 'PV',    1);           % (none)
+
cm = addParameter(cm, 'PV',    1);                 % (none)
 
   
 
   
% region A
+
cm = addParameter(cm, 'k1',    0.1);             % 1/min
cm = addParameter(cm, 'k1',    0.1);     % 1/min
+
cm = addParameter(cm, 'k2',    0.13);           % 1/min
cm = addParameter(cm, 'k2',    0.13);     % 1/min
+
cm = addParameter(cm, 'k3',    0.06);           % ml/(pmol*min)
cm = addParameter(cm, 'k3',    0.06);     % ml/(pmol*min)
+
cm = addParameter(cm, 'k4',    0.0068);       % 1/min
cm = addParameter(cm, 'k4',    0.0068);   % 1/min
 
 
   
 
   
 
% define input function parameter vector
 
% define input function parameter vector
Line 41: Line 59:
  
 
% connect inputs and compartments
 
% connect inputs and compartments
% region A
 
 
cm = addLink(cm, 'L', 'Cp',  'Ce', 'k1');
 
cm = addLink(cm, 'L', 'Cp',  'Ce', 'k1');
 
cm = addLink(cm, 'K', 'Ce', 'Junk','k2');
 
cm = addLink(cm, 'K', 'Ce', 'Junk','k2');
Line 49: Line 66:
 
% specify scan begin and end times
 
% specify scan begin and end times
 
ttt=[ ones(6,1)*5/60; ...    %  6 frames x  5  sec
 
ttt=[ ones(6,1)*5/60; ...    %  6 frames x  5  sec
       ones(2,1)*15/60; ...   %  2 frames x 15  sec
+
       ones(2,1)*15/60; ...   %  2 frames x 15  sec
       ones(6,1)*0.5;...     %  6 frames x  0.5 min
+
       ones(6,1)*0.5;...         %  6 frames x  0.5 min
       ones(3,1)*2;...       %  3 frames x  2  min
+
       ones(3,1)*2;...           %  3 frames x  2  min
       ones(2,1)*5;...       %  2 frames x  5  min
+
       ones(2,1)*5;...           %  2 frames x  5  min
       ones(10,1)*10];       % 10 frames x 10  min
+
       ones(10,1)*10];         % 10 frames x 10  min
  
 
scant = [[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
 
scant = [[0;cumsum(ttt(1:(length(ttt)-1)))] cumsum(ttt)];
Line 69: Line 86:
 
cm = addSensitivity(cm, 'pin', 'k1', 'k2', 'k3', 'k4');
 
cm = addSensitivity(cm, 'pin', 'k1', 'k2', 'k3', 'k4');
 
   
 
   
% set parameter values initial guess, lower and upper bounds.  values are in same order as ensitivities
+
% set parameter values initial guess, lower and upper bounds.  values are in same order as sensitivities
 
%        _____________pin_________________  ______Reg______     
 
%        _____________pin_________________  ______Reg______     
 
pinit = [ 10; 0.4;  0.4;  3;  0.05; 0.01;  0.1;  0.1; 0.05; 0.001; ];
 
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];
+
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.;];
+
pub = [100; 2. ;  2. ; 10;  1.  ; 0.05;  1.;  1.;  1.;    1.;];
 +
</pre>
  
 +
'''Step 2.''' Generate 128 x 128 noisy TACs by adding noise to the noise-free TAC.
 +
 +
<pre>
 
noise_level = 0.1;
 
noise_level = 0.1;
 
for i=1:128*128
 
for i=1:128*128
 
     noisy_data(:,i) = [addNoiseDefault(data,noise_level,scant)];
 
     noisy_data(:,i) = [addNoiseDefault(data,noise_level,scant)];
 
end
 
end
 +
</pre>
 +
 +
'''-------------------------------------------The below step is to use MDCS for estimating FDG rate constants ----------------------------------------------------'''
 +
 +
'''Step 3.''' Estimate FDG rate constants (k1~k4) for 128 x 128 TACs (generated in the step 2) using MDCS.
 +
 +
<pre>
 +
 +
%%%% Define specific settings for your environment %%%%
 +
 +
nworkers = 16;  % run the computation using 16 workers
 +
 +
basepath = 'd:\COMKAT_R4.1b'; % 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
  
for pool_idx =1
+
%%%% End of environment-specific settings %%%%
   
+
 
    mat_pool_n = [31 16 8 4 2 1];
+
% run 5 trials for measuring mean and standard deviation of compute time
   
+
for test_idx = 1:5        
    for test_idx = 1:5
+
 
          
+
         % specify the number of workers to use in this test
         matlabpool(mat_pool_n(pool_idx));
+
         parpool(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;
 
         t0 = clock;
        for_times = 128*128;
 
  
         parfor i=1:for_times
+
      % use 'parfor' to perform parallel computing for 128x128 noisy data
 +
         parfor i=1:128*128
  
             cm2 = set(cm, 'ExperimentalData', noisy_data(:,i));
+
             % set the data, which will be fitted and use COMKAT's function 'fit' to fit the data
 
+
            cm = set(cm, 'ExperimentalData', noisy_data(:,i));
             pfit(:,i) = fit(cm2, pinit, plb, pub);
+
             pfit(:,i) = fit(cm, pinit, plb, pub);
  
 
         end
 
         end
  
         time_consumed_parfor(pool_idx,test_idx) = etime(clock,t0);
+
         time_consumed_parfor(test_idx) = etime(clock,t0);
 
 
        matlabpool close;
 
  
    end
+
        delete(gcp);
 
end
 
end
 
</pre>
 
</pre>

Latest revision as of 15:34, 2 March 2018

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 the example, settings for MDCS must be finished. The introduction of setting MDCS can be found in the Support:Documents:Manual:Distributed Computing with COMKAT.

Quick test

After finishing the setting of MDCS, users could run a quick test for performing parallel computing using MDCS with COMKAT functions. Note: Since MDCS uses headless MATLAB sessions, COMKAT must be copied to all workers (or nodes), and complete path name must be the same as 'local' client. For example, in the client, we put the COMKAT_R4.1b folder in 'd:\COMKAT_R4.1b'. In all wokers, the COMKAT_R4.1b folder must be putted in the same path (i.e., 'd:\COMKAT_R4.1b'). Also, setting path for 'all required COMKAT functions' must be performed in the command-line (i.e., addpath).

nworkers = 16;                       % run the computation using 16 workers
parpool(nworkers);             % specify the number of workers to use in this test
basepath = 'd:\COMKAT_R4.1b'; % set this to the main comkat folder for your computer
addpath(basepath);
parfor i=1:100                         % use 'parfor' to perform parallel computing 
cm = compartmentModel;
end
delete(gcp);

Example of Parallel Computing using MDCS for COMKAT

After finishing the above test, users could change the number of workers for reducing the computational time of a pixel-wise parameter estimation. The below example is to 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 each TAC has 29 time frames. All data are generated from the FDG model (following the step 1 and step 2).

Step 1. Create a 18F-FDG model. The introduction of 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

-------------------------------------------The below step is to use MDCS for estimating FDG rate constants ----------------------------------------------------

Step 3. Estimate FDG rate constants (k1~k4) for 128 x 128 TACs (generated in the step 2) using MDCS.


%%%% Define specific settings for your environment %%%%

nworkers = 16;   % run the computation using 16 workers

basepath = 'd:\COMKAT_R4.1b'; % 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
        parpool(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
            cm = set(cm, 'ExperimentalData', noisy_data(:,i));
            pfit(:,i) = fit(cm, pinit, plb, pub);

        end

        time_consumed_parfor(test_idx) = etime(clock,t0);

        delete(gcp);
end