diff --git a/matRad/matRad_planAnalysis.m b/matRad/matRad_planAnalysis.m index e9f9231ff..015340fd7 100644 --- a/matRad/matRad_planAnalysis.m +++ b/matRad/matRad_planAnalysis.m @@ -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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -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{:}); @@ -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'; @@ -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 \ No newline at end of file +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 diff --git a/matRad/planAnalysis/matRad_EQD2accumulation.m b/matRad/planAnalysis/matRad_EQD2accumulation.m index ff4133745..b175f40d7 100644 --- a/matRad/planAnalysis/matRad_EQD2accumulation.m +++ b/matRad/planAnalysis/matRad_EQD2accumulation.m @@ -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, ... @@ -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 @@ -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 @@ -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); @@ -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 @@ -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'); @@ -160,4 +165,4 @@ title('fixed vs moving'); figure,imshowpair( fixed(:,:,50),movedCT(:,:,50),'scaling','joint'); title('fixed vs moved'); -end \ No newline at end of file +end diff --git a/matRad/planAnalysis/matRad_calcQualityIndicators.m b/matRad/planAnalysis/matRad_calcQualityIndicators.m index 8088b7ee3..7bbb7be9e 100644 --- a/matRad/planAnalysis/matRad_calcQualityIndicators.m +++ b/matRad/planAnalysis/matRad_calcQualityIndicators.m @@ -1,6 +1,6 @@ function qi = matRad_calcQualityIndicators(cst,pln,doseCube,refGy,refVol) % matRad QI calculation -% +% % call % qi = matRad_calcQualityIndicators(cst,pln,doseCube) % qi = matRad_calcQualityIndicators(cst,pln,doseCube,refGy,refVol) @@ -8,16 +8,19 @@ % input % cst: matRad cst struct % pln: matRad pln struct -% doseCube: arbitrary doseCube (e.g. physicalDose) +% doseCube: per-fraction doseCube (e.g. physicalDose) % refGy: (optional) array of dose values used for V_XGy calculation -% default is [40 50 60] +% defaults are derived from doseCube % refVol:(optional) array of volumes (0-100) used for D_X calculation -% default is [2 5 95 98] +% default is [2 5 50 95 98] % NOTE: Call either both or none! +% Dose values are evaluated per fraction. Use display +% scaling outside this function for total-dose plots. % % output -% qi various quality indicators like CI, HI (for -% targets) and DX, VX within a structure set +% qi: various quality indicators like CI, HI (for +% targets), coverage indicators, and DX, VX within a +% structure set % % References % van't Riet et. al., IJROBP, 1997 Feb 1;37(3):731-6. @@ -25,13 +28,13 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Copyright 2016 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 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 % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -39,30 +42,44 @@ matRad_cfg = MatRad_Config.instance(); +if ~exist('refGy', 'var') + refGy = []; +end + +if ~exist('refVol', 'var') + refVol = []; +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; + 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; +targetDoseInfo = matRad_getTargetReferenceDoses(cst,pln); +targetDoseCstIndex = [targetDoseInfo.cstIndex]; + for runVoi = 1:size(cst,1) - + indices = cst{runVoi,4}{1}; - numOfVoxels = numel(indices); + numOfVoxels = numel(indices); voiPrint = sprintf('%3d %20s',cst{runVoi,1},cst{runVoi,2}); %String that will print quality indicators - + qi(runVoi).name = cst{runVoi,2}; + % get Dose, dose is sorted to simplify calculations 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); @@ -89,54 +106,44 @@ voiPrint = sprintf('%s\n%27s',voiPrint,' '); % if current voi is a target -> calculate homogeneity and conformity - if strcmp(cst{runVoi,3},'TARGET') > 0 + if strcmp(cst{runVoi,3},'TARGET') > 0 - % 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 + targetDoseIx = find(targetDoseCstIndex == runVoi,1,'first'); + if isempty(targetDoseIx) + referenceDose = inf; + else + referenceDose = targetDoseInfo(targetDoseIx).refDose; end - if referenceDose == inf + 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); + qi(runVoi).(['CI_' StringReferenceDose 'Gy']) = VTarget95^2/(numOfVoxels * VTreated95); - % Homogeneity Index (one out of many), fieldname contains reference dose + % 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); + 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 @@ -154,4 +161,3 @@ end end - diff --git a/matRad/planAnalysis/matRad_convertFromEvaluationMode.m b/matRad/planAnalysis/matRad_convertFromEvaluationMode.m new file mode 100644 index 000000000..e732d2b04 --- /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..5c47f2860 --- /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 = 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 = 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..188d4a68a --- /dev/null +++ b/matRad/planAnalysis/matRad_getTargetReferenceDoses.m @@ -0,0 +1,91 @@ +function targetDoseInfo = matRad_getTargetReferenceDoses(cst,pln) +% matRad_getTargetReferenceDoses resolves target reference doses from the cst +% +% 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 = getReferenceDoseFromObjectives(cst{cstIdx,6}); + if isfinite(refDose) + refDose = refDose/pln.numOfFractions; + 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 = 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_indicatorWrapper.m b/matRad/planAnalysis/matRad_indicatorWrapper.m index 9e3bf18f2..938fc11b4 100644 --- a/matRad/planAnalysis/matRad_indicatorWrapper.m +++ b/matRad/planAnalysis/matRad_indicatorWrapper.m @@ -1,6 +1,6 @@ function [dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI,refGy,refVol) % matRad indictor wrapper -% +% % call % [dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI) % [dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI,refGy,refVol) @@ -9,8 +9,8 @@ % cst: matRad cst struct % pln: matRad pln struct % resultGUI: matRad resultGUI struct -% refGy: (optional) array of dose values used for V_XGy calculation -% default is [40 50 60] +% refGy: (optional) per-fraction dose values used for V_XGy +% calculation. Defaults are derived from doseCube. % refVol:(optional) array of volumes (0-100) used for D_X calculation % default is [2 5 95 98] % NOTE: Call either both or none! @@ -25,13 +25,13 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Copyright 2017 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 2017 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -39,18 +39,19 @@ % Initialize the matRad configuration instance matRad_cfg = MatRad_Config.instance(); % Display a deprecation warning for the current function -matRad_cfg.dispDeprecationWarning('The matRad_indicatorWrapper function will be deprecated soon!\nPlan analysis is now handled by matRad_planAnalysis!'); +matRad_cfg.dispDeprecationWarning(['The matRad_indicatorWrapper function will be deprecated soon!\n', ... + 'Plan analysis is now handled by matRad_planAnalysis!']); % Initialize an empty cell array for optional arguments to translate the into key-value pairs args = {}; % Check if 'refVol' variable exists and add it to the arguments if it does -if exist('refVol', 'var') - args{end+1,end+2} = {'refVol',refVol}; +if exist('refVol', 'var') && ~isempty(refVol) + args = [args {'refVol',refVol}]; end % Check if 'refGy' variable exists and add it to the arguments if it does -if exist('refGy', 'var') - args{end+1,end+2} = {'refGy',refGy}; +if exist('refGy', 'var') && ~isempty(refGy) + args = [args {'refGy',refGy}]; end % Initialize empty structures for ct and stf, required for matRad_planAnalysis @@ -66,5 +67,3 @@ end - - diff --git a/test/planAnalysis/test_calcQualityIndicators.m b/test/planAnalysis/test_calcQualityIndicators.m new file mode 100644 index 000000000..741e56a9f --- /dev/null +++ b/test/planAnalysis/test_calcQualityIndicators.m @@ -0,0 +1,64 @@ +function test_suite = test_calcQualityIndicators + +test_functions = localfunctions(); + +initTestSuite; + +function test_cov_uses_per_fraction_reference_dose + [cst,pln] = 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_cov_thresholds_are_fractional + [cst,pln] = 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_total_display_does_not_change_canonical_qi + [cst,pln] = qiFixture(78,39,10); + resultGUI.physicalDose = 2 * ones(10,1); + + resultGUI = matRad_planAnalysis(resultGUI,struct(),cst,struct(),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_default_reference_doses_ignore_nan_outside_voi + [cst,pln] = 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] = qiFixture(prescriptionDose,numOfFractions,numOfVoxels) + pln = struct(); + pln.numOfFractions = numOfFractions; + pln.bioParam = struct('model','none'); + + 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;