matRad_compareDose

Purpose ^

Comparison of two dose cubes in terms of gamma index, absolute and visual difference

Synopsis ^

function [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)

Description ^

 Comparison of two dose cubes in terms of gamma index, absolute and visual difference

 call
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, enable)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, contours)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, pln)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, criteria)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, n)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, localglobal)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst, contours)
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst, pln)
               ...
   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)

 input
   cube1:         dose cube 1 as an M x N x O array
   cube2:         dose cube 2 as an M x N x O array
   ct:            ct struct with ct cube
   cst:           list of interesting volumes inside the patient as matRad
                  struct (optional, does not calculate gamma and DVH)
   enable         (optional) specify if all sections are evaluated
                  boolean 3x1 array
                  [1 0 0]: evaluate only basic plots
                  [0 1 0]: evaluate only line profiles
                  [0 0 1]: evaluate only DVH
   contours       (optional) specify if contours are plotted,
                  'on' or 'off'
   pln            (optional) specify BioModel for DVH plot
   criteria:      (optional)[1x2] vector specifying the distance to agreement
                  criterion; first element is percentage difference,
                  second element is distance [mm], default [3 3]
   n:             (optional) number of interpolations. there will be 2^n-1
                  interpolation points. The maximum suggested value is 3.
                  default n=0
   localglobal:   (optional) parameter to choose between 'global' and 'local'
                  normalization 


 output

   gammaCube:      result of gamma index calculation
   gammaPassRate:  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.
   hfig:           Figure handle struct for all 3 figures and subplots,
                   indexed by plane names

 References gamma analysis:
   [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.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cross-reference information ^

This function calls: This function is called by:

Source code ^

0001 function [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)
0002 % Comparison of two dose cubes in terms of gamma index, absolute and visual difference
0003 %
0004 % call
0005 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct)
0006 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst)
0007 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, enable)
0008 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, contours)
0009 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, pln)
0010 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, criteria)
0011 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, n)
0012 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, localglobal)
0013 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable)
0014 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst, contours)
0015 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst, pln)
0016 %               ...
0017 %   [gammaCube,gammaPassRate,hfig] = matRad_compareDose(cube1, cube2, ct, cst,enable , contours, pln, criteria, n,localglobal)
0018 %
0019 % input
0020 %   cube1:         dose cube 1 as an M x N x O array
0021 %   cube2:         dose cube 2 as an M x N x O array
0022 %   ct:            ct struct with ct cube
0023 %   cst:           list of interesting volumes inside the patient as matRad
0024 %                  struct (optional, does not calculate gamma and DVH)
0025 %   enable         (optional) specify if all sections are evaluated
0026 %                  boolean 3x1 array
0027 %                  [1 0 0]: evaluate only basic plots
0028 %                  [0 1 0]: evaluate only line profiles
0029 %                  [0 0 1]: evaluate only DVH
0030 %   contours       (optional) specify if contours are plotted,
0031 %                  'on' or 'off'
0032 %   pln            (optional) specify BioModel for DVH plot
0033 %   criteria:      (optional)[1x2] vector specifying the distance to agreement
0034 %                  criterion; first element is percentage difference,
0035 %                  second element is distance [mm], default [3 3]
0036 %   n:             (optional) number of interpolations. there will be 2^n-1
0037 %                  interpolation points. The maximum suggested value is 3.
0038 %                  default n=0
0039 %   localglobal:   (optional) parameter to choose between 'global' and 'local'
0040 %                  normalization
0041 %
0042 %
0043 % output
0044 %
0045 %   gammaCube:      result of gamma index calculation
0046 %   gammaPassRate:  rate of voxels passing the specified gamma criterion
0047 %                   evaluated for every structure listed in 'cst'.
0048 %                   note that only voxels exceeding the dose threshold are
0049 %                   considered.
0050 %   hfig:           Figure handle struct for all 3 figures and subplots,
0051 %                   indexed by plane names
0052 %
0053 % References gamma analysis:
0054 %   [1]  http://www.ncbi.nlm.nih.gov/pubmed/9608475
0055 %
0056 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 %
0058 % Copyright 2015 the matRad development team.
0059 %
0060 % This file is part of the matRad project. It is subject to the license
0061 % terms in the LICENSE file found in the top-level directory of this
0062 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0063 % of the matRad project, including this file, may be copied, modified,
0064 % propagated, or distributed except according to the terms contained in the
0065 % LICENSE file.
0066 %
0067 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0068 
0069 matRad_cfg = MatRad_Config.instance();
0070 
0071 %% check if cubes consistent
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 % Load colormap for difference
0114 diffCMap = matRad_getColormap('diffMap');
0115 
0116 %% Calculate iso-center slices and resolution
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 %% Integral Energy Output
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 %% Colorwash images
0139 if enable(1) == 1
0140     
0141     % Get the gamma cube
0142     matRad_cfg.dispInfo('Calculating gamma index cube...\n');
0143     if exist('criteria','var')
0144         relDoseThreshold = criteria(1); % in [%]
0145         dist2AgreeMm     = criteria(2); % in [mm]
0146     else
0147         dist2AgreeMm     = 3; % in [mm]
0148         relDoseThreshold = 3; % in [%]
0149     end
0150     
0151     [gammaCube,gammaPassRate] = matRad_gammaIndex(cube1,cube2,resolution,criteria,[],n,localglobal,cst);
0152     
0153     
0154     % Calculate absolute difference cube and dose windows for plots
0155     differenceCube  = cube1-cube2;
0156     doseDiffWindow  = [-max(differenceCube(:)) max(differenceCube(:))];
0157     %doseGammaWindow = [0 max(gammaCube(:))];
0158     doseGammaWindow = [0 2]; %We choose 2 as maximum value since the gamma colormap has a sharp cut in the middle
0159     
0160     
0161     % Plot everything
0162     % Plot dose slices
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         % Initialize Figure
0173         hfig.(planename{plane}).('fig') = figure('Renderer', 'painters', 'Position', [10 50 800 800]);
0174         set(gcf,'Color',[1 1 1]);
0175         
0176         % Plot Dose 1
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         % Plot Dose 2
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         % Plot absolute difference
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         % Plot gamma analysis
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         % Adjusting axes
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 %% Plot profiles through isoCenter
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 %% Calculate and plot DVH
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     % Plot DVH
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

| Generated by m2html © 2005