Skip to content

Use per-fraction reference doses in quality indicators #904

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
acsevillam wants to merge 1 commit into
base: dev_varRBErobOpt
Choose a base branch
from
+486 −145
Open

Use per-fraction reference doses in quality indicators #904

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
150 changes: 108 additions & 42 deletions matRad/matRad_planAnalysis.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,29 +10,22 @@
% stf: matRad stf struct with beam information
% pln: matRad pln struct with plan information
% name / value pairs: Optional parameters for analysis customization
% refGy: (optional) Dose values for V_XGy calculation (default: [40 50 60])
% quantity: (optional) resultGUI dose quantity to analyse
% evaluationMode:(optional) 'perFraction' or 'total' for figures/tables
% doseWindow: (optional) dose axis window for DVH display
% refGy: (optional) Per-fraction dose values for V_XGy calculation
% refVol:(optional) Volume percentages for D_X calculation (default: [2 5 95 98])
%
% output
% resultGUI: Updated resultGUI with analysis data

% Initialize input parser for function arguments
p = inputParser();

% Define required inputs
p.addRequired('ct',@isstruct);
p.addRequired('cst',@iscell);
p.addRequired('stf',@isstruct);
p.addRequired('pln',@isstruct);

% Copyright 2024 the matRad development team.
%
% Copyright 2024 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/LICENSE.md. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% 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/LICENSE.md. 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expand All @@ -48,10 +41,13 @@
p.addRequired('pln',@isstruct);

% Define optional parameters with default values
p.addParameter('refGy',[40 50 60],@isnumeric); % Reference dose values for V_XGy calculation
p.addParameter('refGy',[],@isnumeric); % Per-fraction reference dose values for V_XGy calculation
p.addParameter('refVol',[2 5 95 98],@isnumeric); % Reference volume percentages for D_X calculation
p.addParameter('showDVH',true,@islogical); % Flag to show or hide the DVH plot
p.addParameter('showQI',true,@islogical); % Flag to show or hide the Quality Indicators plot
p.addParameter('quantity','',@(x) ischar(x) || isstring(x)); % Dose quantity to analyze
p.addParameter('evaluationMode','perFraction',@(x) ischar(x) || isstring(x)); % Figure/table evaluation mode
p.addParameter('doseWindow',[],@(x) isempty(x) || (isnumeric(x) && numel(x) == 2)); % DVH dose axis window

% Parse input arguments to extract values
p.parse(ct,cst,stf,pln,varargin{:});
Expand All @@ -65,9 +61,26 @@
refVol = p.Results.refVol;
showDVH = p.Results.showDVH;
showQI = p.Results.showQI;
quantity = p.Results.quantity;
doseWindow = p.Results.doseWindow;
evaluationMode = p.Results.evaluationMode;

if isstring(quantity)
quantity = char(quantity);
end

[~,evaluationMode,evaluationScale] = matRad_convertToEvaluationMode([],pln,evaluationMode);

if ~isempty(doseWindow)
doseWindow = doseWindow(:)';
end

% Determine which dose cube to use based on resultGUI structure
if ~isfield(pln,'bioParam')
if ~isempty(quantity)
visQ = quantity;
elseif isfield(pln,'bioParam') && isobject(pln.bioParam) && isprop(pln.bioParam,'quantityVis')
visQ = pln.bioParam.quantityVis;
elseif isfield(pln,'bioParam') && isstruct(pln.bioParam) && isfield(pln.bioParam,'quantityVis')
visQ = pln.bioParam.quantityVis;
elseif isfield(resultGUI,'RBExD')
visQ = 'RBExD';
Expand All @@ -85,52 +98,105 @@
% Calculate DVH and quality indicators
resultGUI.dvh = matRad_calcDVH(cst,doseCube,'cum'); % Calculate cumulative DVH
resultGUI.qi = matRad_calcQualityIndicators(cst,pln,doseCube,refGy,refVol); % Calculate quality indicators
resultGUI.analysisQuantity = visQ;
resultGUI.evaluationModeBase = 'perFraction';
resultGUI.evaluationMode = evaluationMode;
resultGUI.evaluationScale = evaluationScale;
resultGUI.displayDvh = convertDvhForEvaluation(resultGUI.dvh,pln,evaluationMode);
resultGUI.displayQi = convertQiForEvaluation(resultGUI.qi,pln,evaluationMode);

dvhScen = {};
if isfield(pln,'multScen') && pln.multScen.totNumScen > 1
for i = 1:pln.multScen.totNumScen
scenFieldName = sprintf('%s_scen%d',visQ,i);
if isfield(resultGUI,scenFieldName)
dvhScen{i} = matRad_calcDVH(cst,resultGUI.(scenFieldName),'cum'); % Calculate cumulative scenario DVH
dvhScen{i} = convertDvhForEvaluation( ...
matRad_calcDVH(cst,resultGUI.(scenFieldName),'cum'),pln,evaluationMode); % Calculate cumulative scenario DVH
end
end
end


% Configuration for GUI appearance
matRad_cfg = MatRad_Config.instance();

% Create figure for plots with background color from configuration
hF = figure('Color',matRad_cfg.gui.backgroundColor);

colorSpec = {'Color',matRad_cfg.gui.elementColor,...
'XColor',matRad_cfg.gui.textColor,...
'YColor',matRad_cfg.gui.textColor,...
'GridColor',matRad_cfg.gui.textColor,...
'MinorGridColor',matRad_cfg.gui.backgroundColor};
if showDVH || showQI
% Configuration for GUI appearance
matRad_cfg = MatRad_Config.instance();

% Determine subplot layout based on flags
if showDVH && showQI
hDVHax = subplot(2,1,1,colorSpec{:}); % DVH plot area
hQIax = subplot(2,1,2,colorSpec{:}); % Quality Indicators plot area
elseif showDVH
hDVHax = subplot(1,1,1,colorSpec{:}); % Only DVH plot
elseif showQI
hQIax = subplot(1,1,1,colorSpec{:}); % Only Quality Indicators plot
% Create figure for plots with background color from configuration
hF = figure('Color',matRad_cfg.gui.backgroundColor);

colorSpec = {'Color',matRad_cfg.gui.elementColor,...
'XColor',matRad_cfg.gui.textColor,...
'YColor',matRad_cfg.gui.textColor,...
'GridColor',matRad_cfg.gui.textColor,...
'MinorGridColor',matRad_cfg.gui.backgroundColor};

% Determine subplot layout based on flags
if showDVH && showQI
hDVHax = subplot(2,1,1,colorSpec{:}); % DVH plot area
hQIax = subplot(2,1,2,colorSpec{:}); % Quality Indicators plot area
elseif showDVH
hDVHax = subplot(1,1,1,colorSpec{:}); % Only DVH plot
elseif showQI
hQIax = subplot(1,1,1,colorSpec{:}); % Only Quality Indicators plot
end
end

% Display DVH if enabled
if showDVH
matRad_showDVH(resultGUI.dvh,cst,pln,'axesHandle',hDVHax,'LineWidth',3); % Show DVH plot
matRad_showDVH(resultGUI.displayDvh,cst,pln,'axesHandle',hDVHax,'LineWidth',3); % Show DVH plot

for i = 1:numel(dvhScen)
matRad_showDVH(dvhScen{i},cst,pln,'axesHandle',hDVHax,'LineWidth',0.5,'plotLegend',false,'LineStyle','--'); % Show DVH plot
end

if ~isempty(doseWindow)
xlim(hDVHax,doseWindow);
end
end

% Display Quality Indicators if enabled
if showQI
matRad_showQualityIndicators(hQIax,resultGUI.qi); % Show Quality Indicators plot
matRad_showQualityIndicators(hQIax,resultGUI.displayQi); % Show Quality Indicators plot
end

end
end

function dvh = convertDvhForEvaluation(dvh,pln,evaluationMode)
if isempty(dvh)
return;
end

for i = 1:numel(dvh)
if isfield(dvh(i),'doseGrid') && ~isempty(dvh(i).doseGrid)
dvh(i).doseGrid = matRad_convertToEvaluationMode( ...
dvh(i).doseGrid,pln,evaluationMode);
end
end
end

function qi = convertQiForEvaluation(qi,pln,evaluationMode)
if isempty(qi)
return;
end

doseFields = {'mean','std','max','min','referenceDose'};
for i = 1:numel(qi)
for j = 1:numel(doseFields)
fieldName = doseFields{j};
if isfield(qi(i),fieldName) && isnumeric(qi(i).(fieldName))
qi(i).(fieldName) = matRad_convertToEvaluationMode( ...
qi(i).(fieldName),pln,evaluationMode);
end
end

fields = fieldnames(qi(i));
for j = 1:numel(fields)
fieldName = fields{j};
if strncmp(fieldName,'D_',2) && isnumeric(qi(i).(fieldName))
qi(i).(fieldName) = matRad_convertToEvaluationMode( ...
qi(i).(fieldName),pln,evaluationMode);
end
end

end
end
65 changes: 35 additions & 30 deletions matRad/planAnalysis/matRad_EQD2accumulation.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
function result = matRad_EQD2accumulation(pln1,ct1,cst1,dose1,prescribedDose1, ...
pln2,ct2,cst2,dose2,prescribedDose2)

% matRad function to accumulate and compare dose and EQD2 for two treatment
% plans
% matRad function to accumulate and compare dose and EQD2 for two treatment
% plans
%
% call
% result = matRad_EQD2accumulation(pln1,ct1,cst1,dose1,prescribedDose1, ...
Expand All @@ -17,22 +16,22 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2019 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/LICENSE.md. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% Copyright 2019 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/LICENSE.md. 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% parameters that can be changed to fit to the RT plans, also alpha/beta
% for the EQD_2 calculation can be changed
checkImageRegis = true; % true or false, shows pictures of the particle and photon CT to check if the image registration went wrong (true shows pictures)

checkImageRegis = true; % show CT registration check figures

alphaBetaRatio = 2; % alpha/beta value

Expand Down Expand Up @@ -76,30 +75,30 @@

% deleting empty structs in reference cst
cst2(cellfun(@isempty,cst2(:,4)),:) = [];

% adding target structs with the help of the geometrical information
% from the image registration
for i = 1 : size(cst1)

if strcmp(cst1{i,3} ,'TARGET')

cube1 = zeros(ct1.cubeDim);

structIndices = cst1{i,4}{1};
cube1(structIndices) = 1;
newPhotonCube = imwarp(cube1,Rmoving,geomtform,'OutputView',Rfixed);

% ist das hier zul�ssig? vergr�ssern wir hier durch interpolation
% bei imwarp nicht die volumina?
newStructIndices = find(newPhotonCube > 0);

cst2{end+1,1} = size(cst2) + 1;
cst2{end ,2} = [cst1{i,2} '_Photon'];
cst2{end ,3} = 'TARGET';
cst2{end ,4}{1} = newStructIndices;
cst2{end ,5} = cst1{i,5};
cst2{end ,6} = cst1{i,6};

end
end

Expand All @@ -116,10 +115,13 @@


%% DVH calculation

aimedDose = numOfFractions1 * prescribedDose1 + numOfFractions2 * prescribedDose2;
aimedEQD2 = numOfFractions1 * prescribedDose1 .* ((prescribedDose1 + alphaBetaRatio) / (2 + alphaBetaRatio)) + numOfFractions2 * prescribedDose2 .* ((prescribedDose2 + alphaBetaRatio) / (2 + alphaBetaRatio));

aimedEQD2 = numOfFractions1 * prescribedDose1 .* ...
((prescribedDose1 + alphaBetaRatio) / (2 + alphaBetaRatio)) + ...
numOfFractions2 * prescribedDose2 .* ...
((prescribedDose2 + alphaBetaRatio) / (2 + alphaBetaRatio));

dvh_dose = matRad_calcDVH(cst2,result.totalDose ./ aimedDose);
dvh_EQD2 = matRad_calcDVH(cst2,result.totalEQD2 ./ aimedEQD2);
dvh_EQD2ratio = matRad_calcDVH(cst2,result.EQD2ratio);
Expand All @@ -128,7 +130,9 @@
matRad_showDVH(dvh_dose,cst2,pln2,'plotLegend',false,'axesHandle',gca);
hold on
matRad_showDVH(dvh_EQD2,cst2,pln2,'plotLegend',false,'axesHandle',gca,'lineStyle','--');
title(['Added dose cubes divided by aimed dose, straight line RBExD added [aimed: ' num2str(aimedDose) 'Gy], dashed line EQD_2 added [aimed : ' num2str(aimedEQD2) 'Gy]']);
title(['Added dose cubes divided by aimed dose, straight line RBExD added ', ...
'[aimed: ' num2str(aimedDose) 'Gy], dashed line EQD_2 added ', ...
'[aimed : ' num2str(aimedEQD2) 'Gy]']);
xlabel('relative Dose [%]');
legend(cst2{:,2});
hold off
Expand All @@ -140,17 +144,18 @@
xlabel('relative Dose [%]');

%% Qualitiy indicators
fieldNamesResult = fieldnames(result);
doseResultFields = {'totalDose','EQD2_1','EQD2_2','totalEQD2'};

for resultGUInumber = 1 : numel(fieldNamesResult)
qualityIndicators.(fieldNamesResult{resultGUInumber}) = matRad_calcQualityIndicators(cst2,pln2,result.(fieldNamesResult{resultGUInumber}));
for resultGUInumber = 1 : numel(doseResultFields)
qualityIndicators.(doseResultFields{resultGUInumber}) = matRad_calcQualityIndicators( ...
cst2,pln2,result.(doseResultFields{resultGUInumber}) ./ pln2.numOfFractions,[],[]);
end


%% image registration check (shows ct pictures)
if checkImageRegis == true
movedCT = imwarp(moving,Rmoving,geomtform,'OutputView',Rfixed);

figure,imshowpair(squeeze( fixed(:,50,:)),squeeze( moving(:,50,:)),'Scaling','joint');
title('fixed vs moving');
figure,imshowpair(squeeze( fixed(:,50,:)),squeeze( movedCT(:,50,:)),'scaling','joint');
Expand All @@ -160,4 +165,4 @@
title('fixed vs moving');
figure,imshowpair( fixed(:,:,50),movedCT(:,:,50),'scaling','joint');
title('fixed vs moved');
end
end
Loading
Loading