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. 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 many different tissue types, which have different biological properties. Therefore, a ROI-based method may fail to represent some significant characteristics of the tumor. One proposed method 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).

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)
 
% region A
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
% region A
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 ensitivities
%        _____________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.;];

noise_level = 0.1;
for i=1:128*128
    noisy_data(:,i) = [addNoiseDefault(data,noise_level,scant)];
end

for pool_idx =1
    
    mat_pool_n = [31 16 8 4 2 1];
    
    for test_idx = 1:5
        
        matlabpool(mat_pool_n(pool_idx));
        
        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