Comparison of two dose cubes in terms of gamma index, absolute and visual difference
function [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)
0001 function [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069 matRad_cfg = MatRad_Config.instance();
0070
0071
0072 if ~isequal(size(cube1),size(cube2))
0073 matRad_cfg.dispError('dose cubes must be the same size\n');
0074 end
0075
0076 if ~exist('localglobal','var')
0077 localglobal = 'global';
0078 end
0079 if ~exist('n','var')
0080 n = 0;
0081 end
0082 if ~exist('criteria','var')
0083 criteria = [3 3];
0084 end
0085 if ~exist('cst','var') || isempty(cst)
0086 cst = [];
0087 skip = 1;
0088 matRad_cfg.dispWarning('cst not specified, skipping DVH, isocenter in the center of the doseCube');
0089 else
0090 skip = 0;
0091 end
0092 if ~exist('pln','var') || isempty(pln)
0093 pln = [];
0094 end
0095 if ~exist('contours','var') || isempty(contours)
0096 contours = true;
0097 end
0098 if exist('contours','var') || ~isempty(contours)
0099 if strcmp(contours,'on')
0100 contours = true;
0101 else
0102 contours = false;
0103 end
0104 end
0105
0106 if (~exist('enable','var') || isempty(enable)) || strcmp(enable,'all')
0107 enable = [1 1 1];
0108 end
0109 if enable(1)==0
0110 gammaCube = [];
0111 gammaPassRate = [];
0112 end
0113
0114 diffCMap = matRad_getColormap('diffMap');
0115
0116
0117 if isempty(cst)
0118 [~,s(1)] = max(sum(sum(cube1,1),3));
0119 [~,s(2)] = max(sum(sum(cube1,2),3));
0120 [~,s(3)] = max(sum(sum(cube1,1),2));
0121 isoCenter = [ct.resolution.y*s(1) ct.resolution.x*s(2) ct.resolution.z*s(3)];
0122 else
0123 isoCenter = matRad_getIsoCenter(cst,ct,0);
0124 end
0125
0126 resolution = [ct.resolution.x ct.resolution.y ct.resolution.z];
0127
0128 slicename = {round(isoCenter(2)./resolution(2)),round(isoCenter(1)./resolution(1)),round(isoCenter(3)./resolution(3))};
0129 doseWindow = [0 max([cube1(:); cube2(:)])];
0130 planename = {'coronal','sagittal','axial'};
0131
0132
0133 intEnergy1 = matRad_calcIntEnergy(cube1,ct,pln);
0134 intEnergy2 = matRad_calcIntEnergy(cube2,ct,pln);
0135
0136 matRad_cfg.dispInfo('Integral energy comparison: Cube 1 = %1.4g MeV, Cube 2 = %1.4g MeV, difference = %1.4g Mev\n',intEnergy1,intEnergy2,intEnergy1-intEnergy2);
0137
0138
0139 if enable(1) == 1
0140
0141
0142 matRad_cfg.dispInfo('Calculating gamma index cube...\n');
0143 if exist('criteria','var')
0144 relDoseThreshold = criteria(1);
0145 dist2AgreeMm = criteria(2);
0146 else
0147 dist2AgreeMm = 3;
0148 relDoseThreshold = 3;
0149 end
0150
0151 [gammaCube,gammaPassRate] = matRad_gammaIndex(cube1,cube2,resolution,criteria,[],n,localglobal,cst);
0152
0153
0154
0155 differenceCube = cube1-cube2;
0156 doseDiffWindow = [-max(differenceCube(:)) max(differenceCube(:))];
0157
0158 doseGammaWindow = [0 2];
0159
0160
0161
0162
0163 if contours == false
0164 cstHandle = [];
0165 else
0166 cstHandle = cst;
0167 end
0168
0169 for plane = 1:3
0170 disp(['Plotting ',planename{plane},' plane...']);
0171
0172
0173 hfig.(planename{plane}).('fig') = figure('Renderer', 'painters', 'Position', [10 50 800 800]);
0174 set(gcf,'Color',[1 1 1]);
0175
0176
0177 hfig.(planename{plane}).('cube1').Axes = subplot(2,2,1);
0178 [hfig.(planename{plane}).('cube1').CMap,...
0179 hfig.(planename{plane}).('cube1').Dose,...
0180 hfig.(planename{plane}).('cube1').Ct,...
0181 hfig.(planename{plane}).('cube1').Contour,...
0182 hfig.(planename{plane}).('cube1').IsoDose] = ...
0183 matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube1,plane,slicename{plane},[],[],colorcube,jet,doseWindow,[],100);
0184
0185
0186 hfig.(planename{plane}).('cube2').Axes = subplot(2,2,2);
0187 [hfig.(planename{plane}).('cube2').CMap,...
0188 hfig.(planename{plane}).('cube2').Dose,...
0189 hfig.(planename{plane}).('cube2').Ct,...
0190 hfig.(planename{plane}).('cube2').Contour,...
0191 hfig.(planename{plane}).('cube2').IsoDose] = ...
0192 matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube2,plane,slicename{plane},[],[],colorcube,jet,doseWindow,[],100);
0193
0194
0195 hfig.(planename{plane}).('diff').Axes = subplot(2,2,3);
0196 [hfig.(planename{plane}).('diff').CMap,...
0197 hfig.(planename{plane}).('diff').Dose,...
0198 hfig.(planename{plane}).('diff').Ct,...
0199 hfig.(planename{plane}).('diff').Contour,...
0200 hfig.(planename{plane}).('diff').IsoDose] = ...
0201 matRad_plotSliceWrapper(gca,ct,cstHandle,1,differenceCube,plane,slicename{plane},[],[],colorcube,diffCMap,doseDiffWindow,[],100);
0202
0203
0204 hfig.(planename{plane}).('gamma').Axes = subplot(2,2,4);
0205 gammaCMap = matRad_getColormap('gammaIndex');
0206 [hfig.(planename{plane}).('gamma').CMap,...
0207 hfig.(planename{plane}).('gamma').Dose,...
0208 hfig.(planename{plane}).('gamma').Ct,...
0209 hfig.(planename{plane}).('gamma').Contour,...
0210 hfig.(planename{plane}).('gamma').IsoDose]=...
0211 matRad_plotSliceWrapper(gca,ct,cstHandle,1,gammaCube,plane,slicename{plane},[],[],colorcube,gammaCMap,doseGammaWindow,[],100);
0212
0213
0214 matRad_plotAxisLabels(hfig.(planename{plane}).('cube1').Axes,ct,plane,slicename{plane},[],100);
0215 set(get(hfig.(planename{plane}).('cube1').Axes, 'title'), 'string', 'Dose 1');
0216 matRad_plotAxisLabels(hfig.(planename{plane}).('cube2').Axes,ct,plane,slicename{plane},[],100);
0217 set(get(hfig.(planename{plane}).('cube2').Axes, 'title'), 'string', 'Dose 2');
0218 matRad_plotAxisLabels(hfig.(planename{plane}).('diff').Axes,ct,plane,slicename{plane},[],100);
0219 set(get(hfig.(planename{plane}).('diff').Axes, 'title'), 'string', 'Absolute difference');
0220 matRad_plotAxisLabels(hfig.(planename{plane}).('gamma').Axes,ct,plane,slicename{plane},[],100);
0221 set(get(hfig.(planename{plane}).('gamma').Axes, 'title'), 'string', {[num2str(gammaPassRate{1,2},5) '% of points > ' num2str(relDoseThreshold) '% pass gamma criterion (' num2str(relDoseThreshold) '% / ' num2str(dist2AgreeMm) 'mm)']; ['with ' num2str(2^n-1) ' interpolation points']});
0222
0223 end
0224 end
0225
0226
0227 centerAtIsocenter = false;
0228 if enable(2) == 1
0229 matRad_cfg.dispInfo('Plotting profiles...\n');
0230 fontsize = 12;
0231 profilex{1} = squeeze(cube1(slicename{1},:,slicename{3}));
0232 profiley{1} = squeeze(cube1(:,slicename{2},slicename{3}));
0233 profilez{1} = squeeze(cube1(slicename{1},slicename{2},:));
0234
0235 profilex{2} = squeeze(cube2(slicename{1},:,slicename{3}));
0236 profiley{2} = squeeze(cube2(:,slicename{2},slicename{3}));
0237 profilez{2} = squeeze(cube2(slicename{1},slicename{2},:));
0238
0239 posX = resolution(1)*(1:length(profilex{1}));
0240 posY = resolution(2)*(1:length(profiley{1}));
0241 posZ = resolution(3)*(1:length(profilez{1}));
0242 if centerAtIsocenter
0243 posX = posX - isoCenter(1);
0244 posY = posY - isoCenter(2);
0245 posZ = posZ - isoCenter(3);
0246 end
0247
0248 if exist('pln','var') && ~isempty(pln)
0249 if strcmp(pln.propOpt.bioOptimization,'none')
0250 yLabelString = 'Dose [Gy]';
0251 else
0252 yLabelString = 'RBE x Dose [Gy(RBE)]';
0253 end
0254 else
0255 yLabelString = 'Dose [Gy]';
0256 end
0257
0258 hfig.profiles.fig = figure('Renderer', 'painters', 'Position', [10 50 800 800]);
0259 set(gcf,'Color',[1 1 1]);
0260
0261 hfig.profiles.x = subplot(2,2,1);
0262 plot(posX,profilex{1},'r')
0263 hold on
0264 plot(posX,profilex{2},'r--')
0265 xlabel('X [mm]','FontSize',fontsize)
0266 ylabel(yLabelString,'FontSize',fontsize);
0267 title('x-Profiles');
0268 legend({'Dose 1','Dose 2'},'Location','southeast')
0269 legend boxoff
0270
0271 hfig.profiles.y = subplot(2,2,2);
0272 plot(posY,profiley{1},'r')
0273 hold on
0274 plot(posY,profiley{2},'r--')
0275 xlabel('Y [mm]','FontSize',fontsize)
0276 ylabel(yLabelString,'FontSize',fontsize);
0277 title('y-Profiles');
0278 legend({'Dose 1','Dose 2'},'Location','southeast')
0279 legend boxoff
0280
0281 hfig.profiles.z = subplot(2,2,3);
0282 plot(posZ,profilez{1},'r')
0283 hold on
0284 plot(posZ,profilez{2},'r--')
0285 xlabel('Z [mm]','FontSize',fontsize)
0286 ylabel(yLabelString,'FontSize',fontsize);
0287 title('z-Profiles');
0288 legend({'Dose 1','Dose 2'},'Location','southeast')
0289 legend boxoff
0290
0291 set(hfig.profiles.fig,'name',['Profiles:, x=',num2str(slicename{1}),'mm, y=',num2str(slicename{2}),'mm, z=',num2str(slicename{3}),'mm']);
0292
0293 end
0294
0295
0296 if enable(3) == 1 && ~isempty(cst)
0297 matRad_cfg.dispInfo('Calculating DVH...\n');
0298 dvh1 = matRad_calcDVH(cst,cube1);
0299 dvh2 = matRad_calcDVH(cst,cube2);
0300 dvhWindow = max([dvh1(1).doseGrid dvh2(1).doseGrid]);
0301
0302 disp('Plotting DVH...');
0303
0304 hfig.dvh.fig = figure('Renderer', 'painters', 'Position', [10 100 1000 700]);
0305 set(gcf,'Color',[1 1 1]);
0306 matRad_showDVH(dvh1,cst,pln);
0307 hold on
0308 matRad_showDVH(dvh2,cst,pln,2);
0309 xlim([0 dvhWindow*1.2])
0310 title('Dose Volume Histrogram, Dose 1: solid, Dose 2: dashed')
0311 end
0312
0313 matRad_cfg.dispInfo('Done!');
0314
0315 end