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 ROI (region-of-interest)-based method. This 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.

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