gamma index calculation
according to http://www.ncbi.nlm.nih.gov/pubmed/9608475
call
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,cst)
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,cst)
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,n,cst)
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,localglobal,cst)
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,cst)
...
[gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,localglobal,cst)
input
cube1: dose cube as an M x N x O array
cube2: dose cube as an M x N x O array
resolution: resolution of the cubes [mm/voxel]
criteria: [1x2] vector specifying the distance to agreement
criterion; first element is percentage difference,
second element is distance [mm]
slice: (optional) slice in cube1/2 that will be visualized
n: (optional) number of interpolations. there will be 2^n-1
interpolation points. The maximum suggested value is 3.
localglobal: (optional) parameter to choose between 'global' and 'local'
normalization
cst: list of interessing volumes inside the patient
output
gammaCube: result of gamma index calculation
gammaPassRateCell: rate of voxels passing the specified gamma criterion
evaluated for every structure listed in 'cst'.
note that only voxels exceeding the dose threshold are
considered.
References
[1] http://www.ncbi.nlm.nih.gov/pubmed/9608475
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Copyright 2015 the matRad development team.
This file is part of the matRad project. It is subject to the license
terms in the LICENSE file found in the top-level directory of this
distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
of the matRad project, including this file, may be copied, modified,
propagated, or distributed except according to the terms contained in the
LICENSE file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0001 function [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,localglobal,cst) 0002 % gamma index calculation 0003 % according to http://www.ncbi.nlm.nih.gov/pubmed/9608475 0004 % 0005 % call 0006 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,cst) 0007 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,cst) 0008 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,n,cst) 0009 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,localglobal,cst) 0010 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,cst) 0011 % ... 0012 % [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,localglobal,cst) 0013 % 0014 % input 0015 % cube1: dose cube as an M x N x O array 0016 % cube2: dose cube as an M x N x O array 0017 % resolution: resolution of the cubes [mm/voxel] 0018 % criteria: [1x2] vector specifying the distance to agreement 0019 % criterion; first element is percentage difference, 0020 % second element is distance [mm] 0021 % slice: (optional) slice in cube1/2 that will be visualized 0022 % n: (optional) number of interpolations. there will be 2^n-1 0023 % interpolation points. The maximum suggested value is 3. 0024 % localglobal: (optional) parameter to choose between 'global' and 'local' 0025 % normalization 0026 % cst: list of interessing volumes inside the patient 0027 % 0028 % output 0029 % 0030 % gammaCube: result of gamma index calculation 0031 % gammaPassRateCell: rate of voxels passing the specified gamma criterion 0032 % evaluated for every structure listed in 'cst'. 0033 % note that only voxels exceeding the dose threshold are 0034 % considered. 0035 % 0036 % References 0037 % [1] http://www.ncbi.nlm.nih.gov/pubmed/9608475 0038 % 0039 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0040 % 0041 % Copyright 2015 the matRad development team. 0042 % 0043 % This file is part of the matRad project. It is subject to the license 0044 % terms in the LICENSE file found in the top-level directory of this 0045 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 0046 % of the matRad project, including this file, may be copied, modified, 0047 % propagated, or distributed except according to the terms contained in the 0048 % LICENSE file. 0049 % 0050 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0051 0052 [env, ~] = matRad_getEnvironment(); 0053 0054 % set parameters for gamma index calculation 0055 if exist('criteria','var') 0056 relDoseThreshold = criteria(1); % in [%] 0057 dist2AgreeMm = criteria(2); % in [mm] 0058 else 0059 dist2AgreeMm = 3; % in [mm] 0060 relDoseThreshold = 3; % in [%] 0061 end 0062 0063 % set parameters for gamma index calculation 0064 if ~exist('n','var') 0065 n = 0; 0066 end 0067 0068 % set global or local gamma evaluation 0069 if ~exist('localglobal','var') 0070 localglobal = 'global'; 0071 end 0072 0073 % check if cubes consistent 0074 if ~isequal(size(cube1),size(cube2)) 0075 error('dose cubes must be the same size\n'); 0076 end 0077 0078 % define search nighborhood and resolution 0079 if round(n) ~= n 0080 error('n must be an integer value'); 0081 else 0082 resolution = resolution./(2.^n); 0083 end 0084 0085 neighborhood = 1.5* dist2AgreeMm; % [mm] 0086 searchX = ceil(neighborhood./resolution(1)); 0087 searchY = ceil(neighborhood./resolution(2)); 0088 searchZ = ceil(neighborhood./resolution(3)); 0089 0090 % init cube 0091 % cut all the zeros surrounding the interesting part 0092 cut1 = any(any(cube1,2),3); 0093 cut2 = any(any(cube1,1),3); 0094 cut3 = any(any(cube1,1),2); 0095 0096 % avoids that "zero-slides" are deleted between two interest regions 0097 k = find(cut1); cut1(k(1):k(end)) = 1; 0098 k = find(cut2); cut2(k(1):k(end)) = 1; 0099 k = find(cut2); cut2(k(1):k(end)) = 1; 0100 0101 % copy cubes 0102 cubex1c = cube1; 0103 cubex2c = cube2; 0104 0105 % cut cubes 0106 cubex1c( ~cut1, :, :) = []; % rows 0107 cubex1c( :, ~cut2, :) = []; % columns 0108 cubex1c( :, :, ~cut3) = []; % slices 0109 cubex2c( ~cut1, :, :) = []; % rows 0110 cubex2c( :, ~cut2, :) = []; % columns 0111 cubex2c( :, :, ~cut3) = []; % slices 0112 0113 % pad cubes 0114 cubex2 = zeros([size(cubex2c,1)+2*searchX size(cubex2c,2)+2*searchY size(cubex2c,3)+2*searchZ]); 0115 cubex2((1+searchX):(end-searchX), ... 0116 (1+searchY):(end-searchY), ... 0117 (1+searchZ):(end-searchZ)) = cubex2c; 0118 0119 % interpolate if necessary 0120 if n > 0 0121 switch env 0122 case 'MATLAB' 0123 cubex2 = interp3(cubex2,n,'cubic'); 0124 case 'OCTAVE' 0125 cubex2 = interp3(cubex2,n,'linear'); 0126 end 0127 end 0128 0129 % set up temporary cubes required for calculation 0130 gammaCubeSq = inf*ones(size(cubex1c)); 0131 0132 % adjust dose threshold 0133 if strcmp(localglobal,'local') 0134 doseThreshold = cubex1c .* relDoseThreshold/100; 0135 elseif strcmp(localglobal,'global') 0136 doseThreshold = relDoseThreshold/100 * max(cube1(:)); 0137 end 0138 0139 % search for min 0140 for i = -searchX:searchX 0141 for j = -searchY:searchY 0142 for k = -searchZ:searchZ 0143 0144 delta_sq = ((i*resolution(1))^2 + ... 0145 (j*resolution(2))^2 + ... 0146 (k*resolution(3))^2) / dist2AgreeMm^2; 0147 0148 tmpCube = cubex1c - cubex2((1+((2^n)*searchX)+i) : 2^n : (end-((2^n)*searchX)+i), ... 0149 (1+((2^n)*searchY)+j) : 2^n : (end-((2^n)*searchY)+j), ... 0150 (1+((2^n)*searchZ)+k) : 2^n : (end-((2^n)*searchZ)+k)); 0151 0152 tmpCube = tmpCube.^2 ./ doseThreshold.^2 + delta_sq; 0153 0154 gammaCubeSq = min(gammaCubeSq,tmpCube); 0155 0156 0157 end 0158 end 0159 0160 % display '.'; 0161 0162 end 0163 0164 % evaluate gamma cube and set to zero all the voxel that contain an 0165 % infinite value 0166 gammaCubeSqx = zeros(size(cube1)); 0167 gammaCubeSqx(cut1,cut2,cut3) = gammaCubeSq; 0168 gammaCube = sqrt(gammaCubeSqx); 0169 0170 % set values where we did not compute gamma keeping inf in the cube to zero 0171 gammaCube(cube1<=0 & cube2<=0) = 0; 0172 0173 % compute gamma pass rate 0174 doseIx = cube1 > relDoseThreshold/100*max(cube1(:)) | cube2 > relDoseThreshold/100*max(cube2(:)); 0175 numOfPassGamma = sum(gammaCube(doseIx) < 1); 0176 gammaPassRate = 100 * numOfPassGamma / sum(doseIx(:)); 0177 0178 % compute stats for all segmented sturtures 0179 if exist('cst','var') 0180 gammaPassRateCell = cell(1,2); 0181 gammaPassRateCell{1,1} = 'Whole CT'; 0182 gammaPassRateCell{1,2} = gammaPassRate; 0183 0184 for i = 1:size(cst, 1) 0185 volume = cst{i,4}{1,1}; % indices of voxels of the interesting volume 0186 doseIxVol = false(size(doseIx)); 0187 doseIxVol(volume) = doseIx(volume); 0188 numOfPassGammaVol = sum(gammaCube(doseIxVol) < 1); 0189 gammaPassRateVol = 100 * numOfPassGammaVol / sum(doseIxVol(:)); 0190 if isnan(gammaPassRateVol) % remove organs not receiving any dose (resulting in NaN) 0191 continue; 0192 end 0193 gammaPassRateCell{end+1, 1} = cst{i,2}; 0194 gammaPassRateCell{end, 2} = gammaPassRateVol; 0195 0196 end 0197 end 0198 0199 % visualize if applicable 0200 if exist('slice','var') && ~isempty(slice) 0201 figure 0202 set(gcf,'Color',[1 1 1]); 0203 imagesc(gammaCube(:,:,slice),[0 2]) 0204 myColormap = matRad_getColormap('gammaIndex'); 0205 0206 colormap(gca,myColormap); 0207 colorbar 0208 0209 title({[num2str(gammaPassRate,5) '% of points > ' num2str(relDoseThreshold) ... 0210 '% pass gamma criterion (' num2str(relDoseThreshold) '% / ' ... 0211 num2str(dist2AgreeMm) 'mm)']; ['with ' num2str(2^n-1) ' interpolation points']}); 0212 end