From 6f7361e4ed1a5d16fc959dccbe48e16cf04b6183 Mon Sep 17 00:00:00 2001 From: Camilo Sevilla Date: Tue, 19 May 2026 11:08:37 -0500 Subject: [PATCH] Use per-fraction reference doses in quality indicators --- matRad/matRad_planAnalysis.m | 224 ++++++++++-------- matRad/planAnalysis/matRad_EQD2accumulation.m | 66 +++--- .../matRad_calcQualityIndicators.m | 204 ++++++++-------- .../matRad_convertFromEvaluationMode.m | 39 +++ .../matRad_convertToEvaluationMode.m | 71 ++++++ .../matRad_getTargetReferenceDoses.m | 91 +++++++ .../matRad_resolveDoseAnalysisQuantity.m | 79 ++++++ .../planAnalysis/test_calcQualityIndicators.m | 67 ++++++ 8 files changed, 611 insertions(+), 230 deletions(-) create mode 100644 matRad/planAnalysis/matRad_convertFromEvaluationMode.m create mode 100644 matRad/planAnalysis/matRad_convertToEvaluationMode.m create mode 100644 matRad/planAnalysis/matRad_getTargetReferenceDoses.m create mode 100644 matRad/planAnalysis/matRad_resolveDoseAnalysisQuantity.m create mode 100644 test/planAnalysis/test_calcQualityIndicators.m diff --git a/matRad/matRad_planAnalysis.m b/matRad/matRad_planAnalysis.m index 0b709383f..90f6a4a06 100644 --- a/matRad/matRad_planAnalysis.m +++ b/matRad/matRad_planAnalysis.m @@ -1,144 +1,176 @@ -function resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln,varargin) -% matRad plan analysis function -% This function performs analysis on radiation therapy plans, including DVH (Dose-Volume Histogram) and quality indicators. -% It optionally displays these analyses based on input parameters. +function resultGUI = matRad_planAnalysis(resultGUI, ct, cst, stf, pln, varargin) +% matRad_planAnalysis calculates DVH and quality indicators for a plan. % -% input: -% resultGUI: matRad resultGUI struct containing the analysis results +% call +% resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln) +% resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln,varargin) +% +% input +% resultGUI: matRad resultGUI struct containing dose cubes % ct: matRad ct struct with computed tomography data % cst: matRad cst cell array with structure definitions % 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]) -% refVol:(optional) Volume percentages for D_X calculation (default: [2 5 95 98]) +% 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 % -% output: +% 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. -% -% 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 2024-2026 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Initialize input parser for optional parameters p = inputParser(); +p.addRequired('ct', @isstruct); +p.addRequired('cst', @iscell); +p.addRequired('stf', @isstruct); +p.addRequired('pln', @isstruct); +p.addParameter('refGy', [], @isnumeric); +p.addParameter('refVol', [2 5 95 98], @isnumeric); +p.addParameter('showDVH', true, @islogical); +p.addParameter('showQI', true, @islogical); +p.addParameter('quantity', '', @(x) ischar(x) || isstring(x)); +p.addParameter('evaluationMode', 'perFraction', @(x) ischar(x) || isstring(x)); +p.addParameter('doseWindow', [], @(x) isempty(x) || (isnumeric(x) && numel(x) == 2)); +p.parse(ct, cst, stf, pln, varargin{:}); -% Define required inputs again for clarity -p.addRequired('ct',@isstruct); -p.addRequired('cst',@iscell); -p.addRequired('stf',@isstruct); -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('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 - -% Parse input arguments to extract values -p.parse(ct,cst,stf,pln,varargin{:}); - -% Assign parsed values to variables ct = p.Results.ct; cst = p.Results.cst; -stf = p.Results.stf; pln = p.Results.pln; refGy = p.Results.refGy; refVol = p.Results.refVol; showDVH = p.Results.showDVH; showQI = p.Results.showQI; +quantity = p.Results.quantity; +doseWindow = p.Results.doseWindow; +evaluationMode = p.Results.evaluationMode; -% Determine which dose cube to use based on resultGUI structure -if isfield(resultGUI,'RBExDose') - visQ = 'RBExDose'; -else - visQ = 'physicalDose'; -end +[~, evaluationMode, evaluationScale] = matRad_convertToEvaluationMode([], pln, evaluationMode); -if ~isfield(resultGUI,visQ) - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Unknown quantity ''%s'' to analyse!',visQ); +if ~isempty(doseWindow) + doseWindow = doseWindow(:)'; end -% Validate / Create Scenario model -if ~isfield(pln,'multScen') +visQ = matRad_resolveDoseAnalysisQuantity(resultGUI, pln, quantity); + +if ~isfield(pln, 'multScen') pln.multScen = 'nomScen'; end -if ~isa(pln.multScen,'matRad_ScenarioModel') - pln.multScen = matRad_ScenarioModel.create(pln.multScen,ct); -end +if ~isa(pln.multScen, 'matRad_ScenarioModel') + pln.multScen = matRad_ScenarioModel.create(pln.multScen, ct); +end doseCube = resultGUI.(visQ); - -% 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.dvh = matRad_calcDVH(cst, doseCube, 'cum'); +resultGUI.qi = matRad_calcQualityIndicators(cst, pln, doseCube, refGy, refVol); +resultGUI.analysisQuantity = visQ; +resultGUI.evaluationModeBase = 'perFraction'; +resultGUI.evaluationMode = evaluationMode; +resultGUI.evaluationScale = evaluationScale; +resultGUI.displayDvh = matRad_convertDvhForEvaluation(resultGUI.dvh, pln, evaluationMode); +resultGUI.displayQi = matRad_convertQiForEvaluation(resultGUI.qi, pln, evaluationMode); dvhScen = {}; - - -if 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 +numScenarios = pln.multScen.totNumScen; +if numScenarios > 1 + for i = 1:numScenarios + scenFieldName = sprintf('%s_scen%d', visQ, i); + if isfield(resultGUI, scenFieldName) + dvhScen{i} = matRad_convertDvhForEvaluation( ... + matRad_calcDVH(cst, resultGUI.(scenFieldName), 'cum'), pln, evaluationMode); %#ok 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}; - -% 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 +if showDVH || showQI + matRad_cfg = MatRad_Config.instance(); + 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 + hDVHax = subplot(2, 1, 1, colorSpec{:}); + hQIax = subplot(2, 1, 2, colorSpec{:}); + elseif showDVH + hDVHax = subplot(1, 1, 1, colorSpec{:}); + elseif showQI + hQIax = subplot(1, 1, 1, colorSpec{:}); + 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); for i = 1:numel(dvhScen) - matRad_showDVH(dvhScen{i},cst,pln,'axesHandle',hDVHax,'LineWidth',0.5,'plotLegend',false,'LineStyle','--'); % Show DVH plot + matRad_showDVH(dvhScen{i}, cst, pln, 'axesHandle', hDVHax, ... + 'LineWidth', 0.5, 'plotLegend', false, 'LineStyle', '--'); + 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); +end + +end + +function dvh = matRad_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 = matRad_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 -end \ No newline at end of file + 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 diff --git a/matRad/planAnalysis/matRad_EQD2accumulation.m b/matRad/planAnalysis/matRad_EQD2accumulation.m index 7603d8c3e..010c4760d 100644 --- a/matRad/planAnalysis/matRad_EQD2accumulation.m +++ b/matRad/planAnalysis/matRad_EQD2accumulation.m @@ -1,8 +1,8 @@ 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, ... @@ -18,21 +18,21 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Copyright 2019-2026 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - + + % 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 @@ -76,34 +76,36 @@ % 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 -%% EQD_2 calculations (calculates RBExDose added dose, photon EQD_2, particle EQD_2, added EQD_2 and RBExDose added divided by EQD_2 added [and invers ^-1]) +%% EQD_2 calculations +% Calculates RBExDose added dose, photon EQD_2, particle EQD_2, added EQD_2, +% and RBExDose added divided by EQD_2 added [and invers ^-1]. result.totalDose = numOfFractions1 * warpDose1 + numOfFractions2 * dose2; @@ -116,10 +118,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); @@ -128,7 +133,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 RBExDose added [aimed: ' num2str(aimedDose) 'Gy], dashed line EQD_2 added [aimed : ' num2str(aimedEQD2) 'Gy]']); +title(['Added dose cubes divided by aimed dose, straight line RBExDose added ', ... + '[aimed: ' num2str(aimedDose) 'Gy], dashed line EQD_2 added ', ... + '[aimed : ' num2str(aimedEQD2) 'Gy]']); xlabel('relative Dose [%]'); legend(cst2{:,2}); hold off @@ -140,17 +147,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'); diff --git a/matRad/planAnalysis/matRad_calcQualityIndicators.m b/matRad/planAnalysis/matRad_calcQualityIndicators.m index ebc724618..d90b92443 100644 --- a/matRad/planAnalysis/matRad_calcQualityIndicators.m +++ b/matRad/planAnalysis/matRad_calcQualityIndicators.m @@ -1,23 +1,19 @@ -function qi = matRad_calcQualityIndicators(cst,pln,doseCube,refGy,refVol) -% matRad QI calculation -% -% call: +function qi = matRad_calcQualityIndicators(cst, pln, doseCube, refGy, refVol) +% matRad_calcQualityIndicators calculates DVH quality indicators. +% +% call % qi = matRad_calcQualityIndicators(cst,pln,doseCube) % qi = matRad_calcQualityIndicators(cst,pln,doseCube,refGy,refVol) % -% input: -% cst: matRad cst struct +% input +% cst: matRad cst cell array % pln: matRad pln struct -% doseCube: arbitrary doseCube (e.g. physicalDose) -% refGy: (optional) array of dose values used for V_XGy calculation -% default is [40 50 60] -% refVol:(optional) array of volumes (0-100) used for D_X calculation -% default is [2 5 95 98] -% NOTE: Call either both or none! +% doseCube: per-fraction dose cube +% refGy: (optional) per-fraction dose values for V_XGy calculation +% refVol:(optional) volume percentages for D_X calculation % -% output: -% qi various quality indicators like CI, HI (for -% targets) and DX, VX within a structure set +% output +% qi: quality indicators in per-fraction dose units % % References % van't Riet et. al., IJROBP, 1997 Feb 1;37(3):731-6. @@ -26,132 +22,130 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Copyright 2016-2026 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - matRad_cfg = MatRad_Config.instance(); +if ~exist('refGy', 'var') + refGy = []; +end + if ~exist('refVol', 'var') || isempty(refVol) refVol = [2 5 50 95 98]; end -if ~exist('refGy', 'var') || isempty(refGy) - refGy = floor(linspace(0,max(doseCube(:)),6)*10)/10; +if isempty(refGy) + finiteDose = doseCube(isfinite(doseCube)); + if isempty(finiteDose) + refGy = 0; + else + refGy = floor(linspace(0, max(finiteDose(:)), 6) * 10) / 10; + end end - -% calculate QIs per VOI qi = struct; -for runVoi = 1:size(cst,1) - - indices = cst{runVoi,4}{1}; - numOfVoxels = numel(indices); - voiPrint = sprintf('%3d %20s',cst{runVoi,1},cst{runVoi,2}); %String that will print quality indicators - - % get Dose, dose is sorted to simplify calculations - doseInVoi = sort(doseCube(indices)); - +targetDoseInfo = matRad_getTargetReferenceDoses(cst, pln); +targetDoseCstIndex = [targetDoseInfo.cstIndex]; + +for runVoi = 1:size(cst, 1) + indices = cst{runVoi, 4}{1}; + numOfVoxels = numel(indices); + voiPrint = sprintf('%3d %20s', cst{runVoi, 1}, cst{runVoi, 2}); + qi(runVoi).name = cst{runVoi, 2}; + + doseInVoi = sort(doseCube(indices)); + if ~isempty(doseInVoi) - - qi(runVoi).name = cst{runVoi,2}; - - % easy stats qi(runVoi).mean = mean(doseInVoi); qi(runVoi).std = std(doseInVoi); qi(runVoi).max = doseInVoi(end); qi(runVoi).min = doseInVoi(1); - voiPrint = sprintf('%s - Mean dose = %5.2f Gy +/- %5.2f Gy (Max dose = %5.2f Gy, Min dose = %5.2f Gy)\n%27s', ... - voiPrint,qi(runVoi).mean,qi(runVoi).std,qi(runVoi).max,qi(runVoi).min,' '); + voiPrint = sprintf(['%s - Mean dose = %5.2f Gy +/- %5.2f Gy ', ... + '(Max dose = %5.2f Gy, Min dose = %5.2f Gy)\n%27s'], ... + voiPrint, qi(runVoi).mean, qi(runVoi).std, ... + qi(runVoi).max, qi(runVoi).min, ' '); - DX = @(x) matRad_interp1(linspace(0,1,numOfVoxels),doseInVoi,(100-x)*0.01); + DX = @(x) matRad_interp1(linspace(0, 1, numOfVoxels), doseInVoi, (100 - x) * 0.01); VX = @(x) numel(doseInVoi(doseInVoi >= x)) / numOfVoxels; - % create VX and DX struct fieldnames at runtime and fill for runDX = 1:numel(refVol) - qi(runVoi).(strcat('D_',num2str(refVol(runDX)))) = DX(refVol(runDX)); - voiPrint = sprintf('%sD%d%% = %5.2f Gy, ',voiPrint,refVol(runDX),DX(refVol(runDX))); + qi(runVoi).(strcat('D_', num2str(refVol(runDX)))) = DX(refVol(runDX)); + voiPrint = sprintf('%sD%d%% = %5.2f Gy, ', voiPrint, ... + refVol(runDX), DX(refVol(runDX))); end - voiPrint = sprintf('%s\n%27s',voiPrint,' '); + voiPrint = sprintf('%s\n%27s', voiPrint, ' '); + for runVX = 1:numel(refGy) - sRefGy = num2str(refGy(runVX),3); - qi(runVoi).(['V_' strrep(sRefGy,'.','_') 'Gy']) = VX(refGy(runVX)); - voiPrint = sprintf(['%sV' sRefGy 'Gy = %6.2f%%, '],voiPrint,VX(refGy(runVX))*100); + sRefGy = num2str(refGy(runVX), 3); + qi(runVoi).(['V_' strrep(sRefGy, '.', '_') 'Gy']) = VX(refGy(runVX)); + voiPrint = sprintf(['%sV' sRefGy 'Gy = %6.2f%%, '], ... + voiPrint, VX(refGy(runVX)) * 100); end - voiPrint = sprintf('%s\n%27s',voiPrint,' '); - - % if current voi is a target -> calculate homogeneity and conformity - if strcmp(cst{runVoi,3},'TARGET') > 0 + voiPrint = sprintf('%s\n%27s', voiPrint, ' '); - % loop over target objectives and get the lowest dose objective - referenceDose = inf; - - if isstruct(cst{runVoi,6}) - cst{runVoi,6} = num2cell(arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{runVoi,6})); - end - - for runObjective = 1:numel(cst{runVoi,6}) - % check if this is an objective that penalizes underdosing - obj = cst{runVoi,6}{runObjective}; - if ~isa(obj,'matRad_DoseOptimizationFunction') - try - obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); - catch ME - matRad_cfg.dispWarning('Objective/Constraint not valid!\n%s',ME.message) - continue; - end - end - - %if strcmp(cst{runVoi,6}(runObjective).type,'square deviation') > 0 || strcmp(cst{runVoi,6}(runObjective).type,'square underdosing') > 0 - if isa(obj,'DoseObjectives.matRad_SquaredDeviation') || isa(obj,'DoseObjectives.matRad_SquaredUnderdosing') - referenceDose = (min(obj.getDoseParameters(),referenceDose))/pln.numOfFractions; - end + if strcmp(cst{runVoi, 3}, 'TARGET') > 0 + targetDoseIx = find(targetDoseCstIndex == runVoi, 1, 'first'); + if isempty(targetDoseIx) + referenceDose = inf; + else + referenceDose = targetDoseInfo(targetDoseIx).refDose; end - if referenceDose == inf - voiPrint = sprintf('%s%s',voiPrint,'Warning: target has no objective that penalizes underdosage, '); + if ~isfinite(referenceDose) + voiPrint = sprintf('%s%s', voiPrint, ... + 'Warning: target has no objective that penalizes underdosage, '); else - - StringReferenceDose = regexprep(num2str(round(referenceDose*100)/100),'\D','_'); - % Conformity Index, fieldname contains reference dose - VTarget95 = sum(doseInVoi >= 0.95*referenceDose); % number of target voxels recieving dose >= 0.95 dPres - VTreated95 = sum(doseCube(:) >= 0.95*referenceDose); %number of all voxels recieving dose >= 0.95 dPres ("treated volume") - qi(runVoi).(['CI_' StringReferenceDose 'Gy']) = VTarget95^2/(numOfVoxels * VTreated95); - - % Homogeneity Index (one out of many), fieldname contains reference dose - qi(runVoi).(['HI_' StringReferenceDose 'Gy']) = (DX(5) - DX(95))/referenceDose * 100; - - voiPrint = sprintf('%sCI = %6.4f, HI = %5.2f for reference dose of %3.1f Gy\n',voiPrint,... - qi(runVoi).(['CI_' StringReferenceDose 'Gy']),qi(runVoi).(['HI_' StringReferenceDose 'Gy']),referenceDose); + StringReferenceDose = regexprep(num2str(round(referenceDose * 100) / 100), '\D', '_'); + VTarget95 = sum(doseInVoi >= 0.95 * referenceDose); + VTreated95 = sum(doseCube(:) >= 0.95 * referenceDose); + qi(runVoi).(['CI_' StringReferenceDose 'Gy']) = ... + VTarget95^2 / (numOfVoxels * VTreated95); + + qi(runVoi).(['HI_' StringReferenceDose 'Gy']) = ... + (DX(5) - DX(95)) / referenceDose * 100; + qi(runVoi).referenceDose = referenceDose; + qi(runVoi).COV_95 = VX(0.95 * referenceDose); + qi(runVoi).COV_98 = VX(0.98 * referenceDose); + qi(runVoi).COV_99 = VX(0.99 * referenceDose); + qi(runVoi).COV1 = VX(referenceDose); + + voiPrint = sprintf(['%sCI = %6.4f, HI = %5.2f for ', ... + 'reference dose of %3.1f Gy\n%27s'], ... + voiPrint, qi(runVoi).(['CI_' StringReferenceDose 'Gy']), ... + qi(runVoi).(['HI_' StringReferenceDose 'Gy']), ... + referenceDose, ' '); + voiPrint = sprintf(['%sCOV95 = %6.2f%%, COV98 = %6.2f%%, ', ... + 'COV99 = %6.2f%%, COV1 = %6.2f%%\n'], ... + voiPrint, 100 * qi(runVoi).COV_95, 100 * qi(runVoi).COV_98, ... + 100 * qi(runVoi).COV_99, 100 * qi(runVoi).COV1); end end - %We do it this way so the percentages in the string are not interpreted as format specifiers - matRad_cfg.dispInfo('%s\n',voiPrint); - else - matRad_cfg.dispInfo('%d %s - No dose information.',cst{runVoi,1},cst{runVoi,2}); + + matRad_cfg.dispInfo('%s\n', voiPrint); + else + matRad_cfg.dispInfo('%d %s - No dose information.', cst{runVoi, 1}, cst{runVoi, 2}); end end -% assign VOI names which could be corrupted due to empty structures listOfFields = fieldnames(qi); -for i = 1:size(cst,1) - indices = cst{i,4}{1}; - doseInVoi = sort(doseCube(indices)); - if isempty(doseInVoi) - for j = 1:numel(listOfFields) - qi(i).(listOfFields{j}) = NaN; - end - qi(i).name = cst{i,2}; - end +for i = 1:size(cst, 1) + indices = cst{i, 4}{1}; + doseInVoi = sort(doseCube(indices)); + if isempty(doseInVoi) + for j = 1:numel(listOfFields) + qi(i).(listOfFields{j}) = NaN; + end + qi(i).name = cst{i, 2}; + end end end - diff --git a/matRad/planAnalysis/matRad_convertFromEvaluationMode.m b/matRad/planAnalysis/matRad_convertFromEvaluationMode.m new file mode 100644 index 000000000..50f57c584 --- /dev/null +++ b/matRad/planAnalysis/matRad_convertFromEvaluationMode.m @@ -0,0 +1,39 @@ +function [baseValue, evaluationMode, evaluationScale] = matRad_convertFromEvaluationMode(value, pln, evaluationMode) +% matRad_convertFromEvaluationMode Convert evaluation values to per-fraction. +% +% call +% baseValue = matRad_convertFromEvaluationMode(value,pln,evaluationMode) +% [baseValue,evaluationMode,evaluationScale] = matRad_convertFromEvaluationMode(value,pln,evaluationMode) +% +% input +% value: numeric value, array, or empty value expressed in +% evaluationMode +% pln: matRad pln struct +% evaluationMode: 'perFraction' or 'total' +% +% output +% baseValue: value converted to the per-fraction evaluation base +% evaluationMode: normalized evaluation mode string +% evaluationScale: scalar conversion factor from base to evaluationMode +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2026 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +[~, evaluationMode, evaluationScale] = matRad_convertToEvaluationMode( ... + 1, pln, evaluationMode); +baseValue = value ./ evaluationScale; + +end diff --git a/matRad/planAnalysis/matRad_convertToEvaluationMode.m b/matRad/planAnalysis/matRad_convertToEvaluationMode.m new file mode 100644 index 000000000..decd78d44 --- /dev/null +++ b/matRad/planAnalysis/matRad_convertToEvaluationMode.m @@ -0,0 +1,71 @@ +function [convertedValue, evaluationMode, evaluationScale] = matRad_convertToEvaluationMode(rawValue, pln, evaluationMode) +% matRad_convertToEvaluationMode Convert per-fraction values for evaluation. +% +% call +% convertedValue = matRad_convertToEvaluationMode(rawValue,pln,evaluationMode) +% [convertedValue,evaluationMode,evaluationScale] = matRad_convertToEvaluationMode(rawValue,pln,evaluationMode) +% +% input +% rawValue: numeric value, array, or empty per-fraction value +% pln: matRad pln struct +% evaluationMode: 'perFraction' or 'total' +% +% output +% convertedValue: rawValue converted for the requested evaluation mode +% evaluationMode: normalized evaluation mode string +% evaluationScale: scalar conversion factor applied to rawValue +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2026 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin < 3 || isempty(evaluationMode) + evaluationMode = 'perFraction'; +end + +if isstring(evaluationMode) + evaluationMode = char(evaluationMode); +end + +switch lower(evaluationMode) + case 'perfraction' + evaluationMode = 'perFraction'; + evaluationScale = 1; + case 'total' + evaluationMode = 'total'; + evaluationScale = matRad_getNumOfFractions(pln); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Unknown evaluationMode "%s". Use "perFraction" or "total".', evaluationMode); +end + +convertedValue = rawValue .* evaluationScale; + +end + +function numOfFractions = matRad_getNumOfFractions(pln) + +numOfFractions = 1; +if isfield(pln, 'numOfFractions') && ~isempty(pln.numOfFractions) + numOfFractions = pln.numOfFractions; +end + +if ~(isnumeric(numOfFractions) && isscalar(numOfFractions) && ... + isfinite(numOfFractions) && numOfFractions > 0) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('pln.numOfFractions must be a positive finite scalar.'); +end + +end diff --git a/matRad/planAnalysis/matRad_getTargetReferenceDoses.m b/matRad/planAnalysis/matRad_getTargetReferenceDoses.m new file mode 100644 index 000000000..695f5b465 --- /dev/null +++ b/matRad/planAnalysis/matRad_getTargetReferenceDoses.m @@ -0,0 +1,91 @@ +function targetDoseInfo = matRad_getTargetReferenceDoses(cst, pln) +% matRad_getTargetReferenceDoses resolves per-fraction target doses. +% +% call +% targetDoseInfo = matRad_getTargetReferenceDoses(cst,pln) +% +% input +% cst: matRad cst cell array +% pln: matRad pln struct +% +% output +% targetDoseInfo: struct array with one entry per TARGET VOI. Each +% entry contains cstIndex, name, and refDose in dose +% per fraction. +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2026 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRad_cfg = MatRad_Config.instance(); + +isTarget = cellfun(@(voiType) isequal(voiType, 'TARGET'), cst(:, 3)); +targetRows = find(isTarget); +targetDoseInfo = repmat(struct('cstIndex', [], 'name', [], 'refDose', []), ... + numel(targetRows), 1); + +for refTarget = 1:numel(targetRows) + cstIdx = targetRows(refTarget); + + refDose = matRad_getReferenceDoseFromObjectives(cst{cstIdx, 6}); + if isfinite(refDose) + refDose = matRad_convertFromEvaluationMode(refDose, pln, 'total'); + else + matRad_cfg.dispWarning('Target %s has no objective that penalizes underdosage. Reference dose unavailable.\n', cst{cstIdx, 2}); + end + + targetDoseInfo(refTarget).cstIndex = cstIdx; + targetDoseInfo(refTarget).name = cst{cstIdx, 2}; + targetDoseInfo(refTarget).refDose = refDose; +end + +end + +function refDose = matRad_getReferenceDoseFromObjectives(objectives) +matRad_cfg = MatRad_Config.instance(); +refDose = inf; + +if isempty(objectives) + return +end + +if isstruct(objectives) + objectives = num2cell(arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct, objectives)); +end + +for runObjective = 1:numel(objectives) + obj = objectives{runObjective}; + if ~isa(obj, 'matRad_DoseOptimizationFunction') + try + obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); + catch ME + matRad_cfg.dispWarning('Objective/Constraint not valid!\n%s', ME.message); + continue + end + end + + isUnderdoseObjective = isa(obj, 'DoseObjectives.matRad_SquaredDeviation') || ... + isa(obj, 'DoseObjectives.matRad_SquaredUnderdosing') || ... + isa(obj, 'DoseObjectives.matRad_MinDVH'); + + if isUnderdoseObjective + doseParameters = obj.getDoseParameters(); + doseParameters = doseParameters(isfinite(doseParameters)); + if ~isempty(doseParameters) + refDose = min([refDose; doseParameters(:)]); + end + end +end +end diff --git a/matRad/planAnalysis/matRad_resolveDoseAnalysisQuantity.m b/matRad/planAnalysis/matRad_resolveDoseAnalysisQuantity.m new file mode 100644 index 000000000..92443cb1e --- /dev/null +++ b/matRad/planAnalysis/matRad_resolveDoseAnalysisQuantity.m @@ -0,0 +1,79 @@ +function quantity = matRad_resolveDoseAnalysisQuantity(resultGUI, pln, requestedQuantity) +% matRad_resolveDoseAnalysisQuantity resolves the resultGUI dose field. +% +% call +% quantity = matRad_resolveDoseAnalysisQuantity(resultGUI,pln) +% quantity = matRad_resolveDoseAnalysisQuantity(resultGUI,pln,requestedQuantity) +% +% input +% resultGUI: matRad resultGUI struct +% pln: matRad pln struct +% requestedQuantity: optional explicit resultGUI dose field +% +% output +% quantity: selected resultGUI dose field +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2026 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin < 3 + requestedQuantity = ''; +end + +quantity = matRad_charIfString(requestedQuantity); +if ~isempty(quantity) + matRad_validateQuantity(resultGUI, quantity); + return +end + +quantity = matRad_getPlanQuantity(pln, 'quantityOpt'); +if ~isempty(quantity) && isfield(resultGUI, quantity) + return +end + +if isfield(resultGUI, 'analysisQuantity') && ... + isfield(resultGUI, matRad_charIfString(resultGUI.analysisQuantity)) + quantity = matRad_charIfString(resultGUI.analysisQuantity); +elseif isfield(resultGUI, 'RBExDose') + quantity = 'RBExDose'; +else + quantity = 'physicalDose'; +end + +matRad_validateQuantity(resultGUI, quantity); + +end + +function quantity = matRad_getPlanQuantity(pln, fieldName) +quantity = ''; +if isfield(pln, 'propOpt') && isfield(pln.propOpt, fieldName) && ... + ~isempty(pln.propOpt.(fieldName)) + quantity = matRad_charIfString(pln.propOpt.(fieldName)); +end +end + +function value = matRad_charIfString(value) +if isstring(value) + value = char(value); +end +end + +function matRad_validateQuantity(resultGUI, quantity) +if ~isfield(resultGUI, quantity) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Unknown quantity ''%s'' to analyse!', quantity); +end +end diff --git a/test/planAnalysis/test_calcQualityIndicators.m b/test/planAnalysis/test_calcQualityIndicators.m new file mode 100644 index 000000000..249667a8d --- /dev/null +++ b/test/planAnalysis/test_calcQualityIndicators.m @@ -0,0 +1,67 @@ +function test_suite = test_calcQualityIndicators + +test_functions = localfunctions(); + +initTestSuite; + +function test_covUsesPerFractionReferenceDose +[cst, pln] = helper_qiFixture(78, 39, 10); +doseCube = 2 * ones(10, 1); + +qi = matRad_calcQualityIndicators(cst, pln, doseCube, [], []); + +assertElementsAlmostEqual(qi.COV1, 1); +assertElementsAlmostEqual(qi.referenceDose, 2); + +function test_covThresholdsAreFractional +[cst, pln] = helper_qiFixture(78, 39, 10); +doseCube = 0.99 * 2 * ones(10, 1); + +qi = matRad_calcQualityIndicators(cst, pln, doseCube, [], []); + +assertElementsAlmostEqual(qi.COV_99, 1); +assertElementsAlmostEqual(qi.COV1, 0); + +function test_totalDisplayDoesNotChangeCanonicalQi +[cst, pln] = helper_qiFixture(78, 39, 10); +ct.numOfCtScen = 1; +stf = struct(); +resultGUI.physicalDose = 2 * ones(10, 1); + +resultGUI = matRad_planAnalysis(resultGUI, ct, cst, stf, pln, ... + 'showDVH', false, 'showQI', false, 'evaluationMode', 'total'); + +assertElementsAlmostEqual(resultGUI.qi.COV1, 1); +assertElementsAlmostEqual(resultGUI.qi.referenceDose, 2); +assertEqual(resultGUI.evaluationModeBase, 'perFraction'); +assertEqual(resultGUI.evaluationMode, 'total'); +assertElementsAlmostEqual(resultGUI.displayQi.referenceDose, 78); +assertTrue(isfield(resultGUI.qi, 'V_2Gy')); +assertFalse(isfield(resultGUI.qi, 'V_40Gy')); +assertElementsAlmostEqual(resultGUI.displayDvh(1).doseGrid, ... + resultGUI.dvh(1).doseGrid * pln.numOfFractions); + +function test_defaultReferenceDosesIgnoreNanOutsideVoi +[cst, pln] = helper_qiFixture(78, 39, 2); +doseCube = [1; 2; NaN]; + +qi = matRad_calcQualityIndicators(cst, pln, doseCube, [], []); + +assertTrue(isfield(qi, 'V_2Gy')); +assertFalse(isfield(qi, 'V_NaNGy')); + +function [cst, pln] = helper_qiFixture(prescriptionDose, numOfFractions, numOfVoxels) +pln = struct(); +pln.numOfFractions = numOfFractions; +pln.bioModel = 'none'; +pln.propOpt.quantityOpt = 'physicalDose'; + +objective = DoseObjectives.matRad_SquaredDeviation(1, prescriptionDose); + +cst = cell(1, 6); +cst{1, 1} = 1; +cst{1, 2} = 'TARGET'; +cst{1, 3} = 'TARGET'; +cst{1, 4}{1} = 1:numOfVoxels; +cst{1, 5} = struct('Visible', true, 'Priority', 1); +cst{1, 6}{1} = objective;