Difference between revisions of "Support:Documents:Examples: Reslicing using COMKAT image tool (advanced)"

From COMKAT wiki
Jump to navigation Jump to search
 
(18 intermediate revisions by 2 users not shown)
Line 1: Line 1:
==Estimation of Physiological Parameters Using a Physiologically Based Model==
+
==Reslicing 3D image volume using COMKAT image tool (advanced)==
  
 
===Overview===
 
===Overview===
 +
In the advanced 3D slicing using COMKAT image tool, we’re going to describe a slicing approach assisted by using the user-graphical-interface of COMKAT image tool.
  
To evaluate glucose transport and phosphorylation in skeletal muscle, [http://online.medphys.org/resource/1/mphya6/v38/i8/p4587_s1 a physiologically based model]  has been proposed by our group. Herein, we show an example of implementation of the proposed model using COMKAT software. In particular, the below example shows the implementation of the proposed model for a phosphorylatable glucose analog (e.g. 18F-labeled 2-fluoro-2-deoxy-D-glucose).  
+
A function for slicing is provided as an example to show the interaction method between the COMKAT GUI and slicing functions.
  
[[Image:Example-fig.jpg]]
 
  
The above figure shows the kinetics of [<sup>18</sup>F]2FDG in skeletal muscle. The model has three tissue compartments with five rate constants.  
+
====Background====
 +
If you are not familiar with the approach of using coordinateGen() and sliceVolume(), you may find more details in the page, [http://comkat.case.edu/index.php?title=3D_Reslicing_using_COMKAT_image_tool_(basic) 3D Slicing using COMKAT Image Tool (basic)]
  
The first compartment is the interstitial concentration of [<sup>18</sup>F]2FDG .
 
  
The second compartment is the intracellular concentration of [<sup>18</sup>F]2FDG .
 
  
The third compartment is the intracellular concentration of [<sup>18</sup>F]2FDG-6-P .  
+
----
 +
====Use COMKAT image tool to do the image processing====
 +
Load a file using COMKAT image tool and use the functions to translate and rotate the image volume as you desire.  
  
Rate constants '''k'''<sub>'''1'''</sub> and '''k'''<sub>'''2'''</sub> describe the passive exchange of the glucose analog between plasma and interstitial space.
 
  
Rate constants '''k'''<sub>'''3a'''</sub> and '''k'''<sub>'''4a'''</sub> describe the inward and backward transport of the glucose analog via glucose transporters.
+
====Send the handles to workspace====
 +
Send the GUI handles of COMKAT image tool to MATLAB workspace.
  
Rate constant '''k'''<sub>'''5a'''</sub> describe the phosphorylation of the intracellular glucose analog catalyzed by hexokinase.
+
[[Image: Example_Reslice_1.jpg‎]]
  
===COMKAT Implementation of the Physiologically Based Model===
 
  
 +
====Create the function for reslicing====
 +
Create following MATLAB function:
  
==== Create a new model ====
 
 
<pre>
 
<pre>
cm=compartmentModel;
+
function outputVolume = fcnReSlice(GUI_handles, a_row, a_column, a_plane, flagUseCurrIvd)
 +
%%*****************************************************************************%%
 +
% Example:  Using ImageVolumeData::sliceVolume() to reslice 3D images
 +
% e.g.
 +
%  outputVolume = fcnReSlice(GUI_handles [, a_row, a_column, a_plane, flagUseCurrIvd]);
 +
%
 +
%  Parameters -
 +
%      GUI_handles : GUI handles obtained from comkatimagetool
 +
%                a_row : Desired aspect ratio in row direction    (default: 1.0)
 +
%            a_column : Desired aspect ratio in colume direction  (default: 1.0)
 +
%              a_plane : Desired aspect ratio in plane direction  (default: 1.0)
 +
%  flagUseCurrIvd : A flag to control the data used to reslice
 +
%                                ( 1: Current COMKAT ImageVolumeData (default), 0: Original)
 +
%    outputVolume : Resliced 3D volume images
 +
%
 +
%
 +
%  Dylan Su    2012-07-30
 +
%  kuan-hao.su@case.edu
 +
%%*****************************************************************************%%
  
% Define parameters
+
%% Obtain the imageVolumeData object from GUI_handles
cm=addParameter(cm,'k1',0.1);                                            
+
ivd = GUI_handles.imageVolumeData{1}; % assume PET is image 1
cm=addParameter(cm,'k2','k1/Fis');                                       
 
cm=addParameter(cm,'k3','VG/(Fis*KGa+Fis*ISg*KGa/KGg)');
 
cm=addParameter(cm,'k4','VG/(Fic*KGa+Fic*ICg*KGa/KGg)');
 
cm=addParameter(cm,'k5','VH/(Fic*KHa+Fic*ICg*KHa/KHg)'); 
 
  
cm=addParameter(cm,'Fb',0.02);         % fraction of total space occupied by blood
+
if nargin == 1,
cm=addParameter(cm,'Fis',0.3);         % fraction of total space occupied by interstitial space
+
    a_row = 1;
cm=addParameter(cm,'Fic','1-Fis-Fb'); % fraction of total space occupied by intracellular space
+
    a_column = 1;
 +
    a_plane = 1;
 +
end
  
cm=addParameter(cm,'Pg',6);     % Plasma glucose concentration (mM)
+
if nargin < 5,
cm=addParameter(cm,'ISg',5.4); % Interstitial glucose concentration (mM)
+
     flagUseCurrIvd = 1; % 1: use current COMKAT image volume
cm=addParameter(cm,'ICg',0.2); % Intracellular glucose concentration (mM)
+
end
  
cm=addParameter(cm,'KGg',3.5); % Michaelis constant (mM) of glucose transporter (GLUT) for glucose
+
if (~flagUseCurrIvd),
cm=addParameter(cm,'KHg',0.13); % Michaelis constant (mM) of hexokinase (HK) for glucose
+
    %% Get the information from origianl data
 +
    positionInput = get(ivd,'ImagePositionPatient');
 +
    orientationInput = get(ivd,'ImageOrientationPatient');
 +
    pixelSpacing = get(ivd, 'PixelSpacing');
 +
    [nr, nc, np, ~] = get(ivd, 'VolumeDimension');
 +
    nrows = nr;
 +
    ncols = nc;
 +
    nplanes = np;
 +
   
 +
else
 +
    %% Get the information from current COMKAT imageVolumeData
 +
    idxSA = 1;  % get infor from short axis (SA) view
 +
   
 +
    positionInput = GUI_handles.view{idxSA}.position{1};
 +
    orientationInput = GUI_handles.view{idxSA}.orientation{1};
 +
    pixelSpacing = repmat( GUI_handles.view{idxSA}.pixelSpacing(1),[1,3]);
 +
    nrows = GUI_handles.rows;
 +
    ncols = GUI_handles.columns;
 +
   
 +
    [~, ~, np, ~] = get(ivd, 'VolumeDimension');
 +
   
 +
    pixelSpacingOrg = get(ivd, 'PixelSpacing');
 +
   
 +
    nplanes = ceil(np * pixelSpacingOrg(3) / pixelSpacing(3));
 +
end
  
cm=addParameter(cm,'KGa',14); % Michaelis constant (mM) of GLUT for glucose analog
+
%% Calculate the desired dimension
cm=addParameter(cm,'KHa',0.17); % Michaelis constant (mM) of HK for glucose analog
+
dnrows = ceil(nrows * a_row);
 +
dncols = ceil(ncols * a_column);
 +
dnplanes = ceil(nplanes * a_plane);
  
cm=addParameter(cm,'VG','(k1*Pg-k1*ISg)/(ISg/(KGg+ISg)-ICg/(KGg+ICg))'); % Maximal velocity of glucose transport for glucose and its analogs
+
% calcualte new pixelSpacing by preserving the original FOV size
cm=addParameter(cm,'VH','(k1*Pg-k1*ISg)/(ICg/(KHg+ICg))'); % Maximal velocity of glucose phosphorylation for glucose and its analogs
+
pixelSpacing = pixelSpacing ./ [dnrows/nrows, dncols/ncols, dnplanes/nplanes];
  
% Define compartment
+
%% Set the scale and offset the same as the original volume
cm=addCompartment(cm,'IS'); % interstitial
+
% use same scale and offset for subvolume as original volume
cm=addCompartment(cm,'IC'); % intracellular
+
% 0 = scaledPixel = rawPixel * s + o --> rawPixel = -o/s;
cm=addCompartment(cm,'IP'); % intracellular phosphorylated
+
s = get(ivd, 'VolumeFrameBufferScaleFactor');
cm=addCompartment(cm,'Junk');
+
o = get(ivd, 'VolumeFrameBufferRescaleIntercept');
  
% Define link
+
rawBackgroundPixelValue = -o/s;
cm=addLink(cm,'L','Cp','IS','k1');
 
cm=addLink(cm,'K','IS','Junk','k2');
 
cm=addLink(cm,'K','IS','IC','k3');
 
cm=addLink(cm,'K','IC','IS','k4');
 
cm=addLink(cm,'K','IC','IP','k5');
 
</pre>
 
 
 
==== Specify the input, output and scan time====
 
<pre>
 
% Specific activity (sa): if the unit of image data is the same with that of input function, the specific activity is 1.
 
cm=addParameter(cm,'sa',1);
 
% Usually, the decay time correction of image data is performed. So, dk is zero.
 
cm=addParameter(cm,'dk',0);
 
  
% Define scan time (minutes)
+
%% Determine location of first pixel in output volume
delay = 0.0;
+
posSA = positionInput;
scanduration = 120;  
 
  
t=[ones(5,1)/30;ones(10,1)/12;ones(12,1)*0.5;ones(8,1);ones(21,1)*5];
+
% xyz of the center of first pixel
scant = [[0;cumsum(t(1:(length(t)-1)))] cumsum(t)];
+
if (~flagUseCurrIvd),
scanTime  = [scant(:,1),scant(:,2)];
+
    % for subject's data, 'posSA' is the center of the first pixel
 +
    position = posSA;
 +
else
 +
    % for imageVolumeData, 'posSA' is the center of the folume
 +
    position = posSA - orientationInput(:,3) * (dnplanes * pixelSpacing(3)) / 2;
 +
end
  
cm = set(cm, 'ScanTime', scanTime);
+
planePosStep = orientationInput(:,3) * pixelSpacing(3); % calculate the step of new plane position
  
lambda = [-12.02 -2.57 -0.02];
+
%% Build output volume plane-by-plane
 +
outputVolume = zeros(dnrows, dncols, dnplanes); % initialize the output images
  
a = [1771.5 94.55 14.27];
+
for idxP = 1 : dnplanes,
 +
    planePos = position + (idxP - 1) * planePosStep; % position of plane to be interpolated
 +
   
 +
    % determine indicies in original volume corresponding to xyz physical (mm) location of plane
 +
    [u, v, w] = coordinateGen(ivd, ...
 +
        dncols, dnrows, pixelSpacing, planePos, orientationInput);
 +
   
 +
    % obtain a slice by interpolating from the ivd volumeFrameBuffer
 +
    fprintf('plane ==> %i/%i\n', idxP, dnplanes);
 +
    outputVolume(:, :, idxP) = sliceVolume(ivd, v , u , w, rawBackgroundPixelValue,'linear');
 +
   
 +
end
  
cm = addParameter(cm, 'pfeng', [delay a lambda]');
 
 
% Define input function
 
cm = addInput(cm, 'Cp', 'sa', 'dk', 'fengInputByPar','pfeng'); % Cp is the plasma input function
 
cm = addInput(cm, 'Ca',1,0, 'fengInputByPar', 'pfeng'); % Ca is the decay-corrected whole-blood input function. Herein, we assume that Cp=Ca.
 
 
% Define output obtained from each normalized compartment
 
wlistTotal={'IS',1;'IC',1;'IP',1};
 
xlistTotal={'Ca','Fb'};
 
cm=addOutput(cm,'TissueTotal',wlistTotal,xlistTotal);
 
 
</pre>
 
</pre>
  
Usually, the model outputs consist of compartments and weights with which they contribut to outputs (e.g. Fis*IS, Fic*IC, Fic*IP).
 
 
However, as described in detail in [http://online.medphys.org/resource/1/mphya6/v38/i8/p4587_s1 our published paper], tracer concentrations of all compartments (i.e. IS, IC and IP) in final model equations are normalized to total tissue space.
 
 
Accordingly, the weights of these compartment are 1.
 
 
==== Solve the model output and Generate synthetic data====
 
<pre>
 
[PET,PETIndex,Output,OutputIndex]=solve(cm);
 
 
noise_level=0.05;
 
  
sd=noise_level*sqrt(PET(:,3)./(PET(:,2)-PET(:,1)));
+
====Use the function to reslice the image volume====
 +
Then, you can reslice the volume data with this function.
 +
        e.g.   outputVolume = fcnReSlice( GUI_handles );
 +
                 
 +
                  The ’outputVolume‘ is the resliced 3D volume.
  
data=sd.*randn(size(PET(:,3)))+PET(:,3);
 
  
cm=set(cm,'ExperimentalData',data);
+
Or you may input the desired aspect ratios in each direction to vary the sampled resolutions of the resliced images.
 +
        e.g.  outputVolume = fcnReSlice( GUI_handles , 1.0, 1.0, 2.0);
 +
                 
 +
                  In this case, the number of sampling in the plane-direction would be doubled.
  
cm = set(cm,'ExperimentalDataSD',sd);
 
 
figure;
 
t = 0.5*(PET(:,1)+PET(:,2));
 
plot(t,PET(:,3),'g-',t,data,'ro','LineWidth',2);
 
xlabel('Time (minutes)');
 
ylabel('Concentrations (uCi/mL)');
 
legend('Noisy data','Noise-free data');
 
</pre>
 
 
[[Image:Example11.jpg]]
 
 
==== Fit Data ====
 
<pre>
 
% Define model parameters to be estimated
 
cm=addSensitivity(cm,'k1','ISg','ICg','Fis','Fb');
 
 
optsIRLS = setopt('SDModelFunction', @IRLSnoiseModel);
 
cm = set(cm, 'IRLSOptions', optsIRLS);
 
 
oo = optimset('TolFun', 1e-8, 'TolX', 1e-4,'Algorithm','interior-point');
 
cm = set(cm, 'OptimizerOptions', oo);
 
 
% Set initial conditions and bounds
 
%        k1,    ISg,  ICg,    Fs,    Fb
 
pinit = [0.01  ; 4.4 ; 0.1    ;0.15 ; 0.01];
 
plb  = [0.001; 1    ; 0.001;0.10 ; 0    ];
 
pub  = [0.5  ; 6    ; 1      ;0.60 ; 0.04];
 
 
[pfit, qfit, modelfit, exitflag, output, lambda, grad, hessian, objfunval] = fitGen(cm, pinit, plb, pub, 'IRLS');
 
</pre>
 
 
====Results====
 
<pre>
 
figure;
 
t = 0.5*(PET(:,1)+PET(:,2));
 
plot(t,data,'ro',t,modelfit,'b-','LineWidth',2);
 
xlabel('Time (minutes)');
 
ylabel('Concentrations (uCi/mL)');
 
legend('Noisy data', 'Model fit');
 
 
Pg=6;KGa=14;KHa=0.17;KGg=3.5;KHg=0.13; % Units are mM
 
 
k1=pfit(1);k2=pfit(1)/pfit(4);ISg=pfit(2);ICg=pfit(3);Fis=pfit(4);Fb=pfit(5);
 
 
% Physiological parameters
 
 
VG=(k1*Pg-k1*ISg)/(ISg/(KGg+ISg)-ICg/(KGg+ICg))     
 
 
VH=(k1*Pg-k1*ISg)/(ICg/(KHg+ICg))           
 
         
 
CI=VG*ISg/(KGg+ISg)         
 
                         
 
CE=VG*ICg/(KGg+ICg)   
 
                             
 
PR=VH*ICg/(KHg+ICg)                 
 
</pre>
 
  
[[Image:Example12.jpg]]
+
If you’d like to reslice the image volume from original image data without pre-sampling, translation and rotation, you may set the fifth parameter to zero.
 +
        e.g.  outputVolume = fcnReSlice( GUI_handles , 1.0, 1.0, 2.0, 0);
 +
                 
 +
                  So now the ‘outputVolume’ is the resliced volume of the original data.

Latest revision as of 20:45, 8 August 2012

Reslicing 3D image volume using COMKAT image tool (advanced)

Overview

In the advanced 3D slicing using COMKAT image tool, we’re going to describe a slicing approach assisted by using the user-graphical-interface of COMKAT image tool.

A function for slicing is provided as an example to show the interaction method between the COMKAT GUI and slicing functions.


Background

If you are not familiar with the approach of using coordinateGen() and sliceVolume(), you may find more details in the page, 3D Slicing using COMKAT Image Tool (basic)



Use COMKAT image tool to do the image processing

Load a file using COMKAT image tool and use the functions to translate and rotate the image volume as you desire.


Send the handles to workspace

Send the GUI handles of COMKAT image tool to MATLAB workspace.

Example Reslice 1.jpg


Create the function for reslicing

Create following MATLAB function:

function outputVolume = fcnReSlice(GUI_handles, a_row, a_column, a_plane, flagUseCurrIvd)
%%*****************************************************************************%%
% Example:  Using ImageVolumeData::sliceVolume() to reslice 3D images
% e.g. 
%  outputVolume = fcnReSlice(GUI_handles [, a_row, a_column, a_plane, flagUseCurrIvd]);
% 
%  Parameters -
%       GUI_handles : GUI handles obtained from comkatimagetool
%                 a_row : Desired aspect ratio in row direction     (default: 1.0)
%            a_column : Desired aspect ratio in colume direction  (default: 1.0)
%              a_plane : Desired aspect ratio in plane direction   (default: 1.0)
%   flagUseCurrIvd : A flag to control the data used to reslice 
%                                 ( 1: Current COMKAT ImageVolumeData (default), 0: Original)
%     outputVolume : Resliced 3D volume images
%
%
%   Dylan Su     2012-07-30
%   kuan-hao.su@case.edu
%%*****************************************************************************%%

%% Obtain the imageVolumeData object from GUI_handles
ivd = GUI_handles.imageVolumeData{1}; % assume PET is image 1

if nargin == 1,
    a_row = 1;
    a_column = 1;
    a_plane = 1;
end

if nargin < 5,
    flagUseCurrIvd = 1;  % 1: use current COMKAT image volume
end

if (~flagUseCurrIvd),
    %% Get the information from origianl data
    positionInput = get(ivd,'ImagePositionPatient');
    orientationInput = get(ivd,'ImageOrientationPatient');
    pixelSpacing = get(ivd, 'PixelSpacing');
    [nr, nc, np, ~] = get(ivd, 'VolumeDimension');
    nrows = nr;
    ncols = nc;
    nplanes = np;
    
else
    %% Get the information from current COMKAT imageVolumeData
    idxSA = 1;  % get infor from short axis (SA) view
    
    positionInput = GUI_handles.view{idxSA}.position{1};
    orientationInput = GUI_handles.view{idxSA}.orientation{1};
    pixelSpacing = repmat( GUI_handles.view{idxSA}.pixelSpacing(1),[1,3]);
    nrows = GUI_handles.rows;
    ncols = GUI_handles.columns;
    
    [~, ~, np, ~] = get(ivd, 'VolumeDimension');
    
    pixelSpacingOrg = get(ivd, 'PixelSpacing');
    
    nplanes = ceil(np * pixelSpacingOrg(3) / pixelSpacing(3)); 
end

%% Calculate the desired dimension
dnrows = ceil(nrows * a_row);
dncols = ceil(ncols * a_column);
dnplanes = ceil(nplanes * a_plane);

% calcualte new pixelSpacing by preserving the original FOV size
pixelSpacing = pixelSpacing ./ [dnrows/nrows, dncols/ncols, dnplanes/nplanes];

%% Set the scale and offset the same as the original volume
% use same scale and offset for subvolume as original volume
% 0 = scaledPixel = rawPixel * s + o --> rawPixel = -o/s;
s = get(ivd, 'VolumeFrameBufferScaleFactor');
o = get(ivd, 'VolumeFrameBufferRescaleIntercept');

rawBackgroundPixelValue = -o/s;

%% Determine location of first pixel in output volume
posSA = positionInput;

% xyz of the center of first pixel
if (~flagUseCurrIvd),
    % for subject's data, 'posSA' is the center of the first pixel
    position = posSA;
else
    % for imageVolumeData, 'posSA' is the center of the folume
    position = posSA - orientationInput(:,3) * (dnplanes * pixelSpacing(3)) / 2;
end

planePosStep = orientationInput(:,3) * pixelSpacing(3); % calculate the step of new plane position

%% Build output volume plane-by-plane
outputVolume = zeros(dnrows, dncols, dnplanes);  % initialize the output images

for idxP = 1 : dnplanes,
    planePos = position + (idxP - 1) * planePosStep; % position of plane to be interpolated
    
    % determine indicies in original volume corresponding to xyz physical (mm) location of plane
    [u, v, w] = coordinateGen(ivd, ...
        dncols, dnrows, pixelSpacing, planePos, orientationInput);
    
    % obtain a slice by interpolating from the ivd volumeFrameBuffer
    fprintf('plane ==> %i/%i\n', idxP, dnplanes);
    outputVolume(:, :, idxP) = sliceVolume(ivd, v , u , w, rawBackgroundPixelValue,'linear');
    
end


Use the function to reslice the image volume

Then, you can reslice the volume data with this function.

        e.g.   outputVolume = fcnReSlice( GUI_handles );
                 
                 The ’outputVolume‘ is the resliced 3D volume. 


Or you may input the desired aspect ratios in each direction to vary the sampled resolutions of the resliced images.

        e.g.   outputVolume = fcnReSlice( GUI_handles , 1.0, 1.0, 2.0);
                 
                 In this case, the number of sampling in the plane-direction would be doubled.


If you’d like to reslice the image volume from original image data without pre-sampling, translation and rotation, you may set the fifth parameter to zero.

        e.g.   outputVolume = fcnReSlice( GUI_handles , 1.0, 1.0, 2.0, 0);
                 
                 So now the ‘outputVolume’ is the resliced volume of the original data.