From 141ffb3665a59bf1b990da8b9ee4d2a99e5b5135 Mon Sep 17 00:00:00 2001 From: Camilo Sevilla Date: Tue, 19 May 2026 12:49:15 -0500 Subject: [PATCH] Fix biological warm start for cell ax bx --- matRad/matRad_fluenceOptimization.m | 55 +---- .../matRad_prepareBiologicalOptimizationDij.m | 194 ++++++++++++++++++ .../test_biologicalOptimizationDij.m | 116 +++++++++++ 3 files changed, 312 insertions(+), 53 deletions(-) create mode 100644 matRad/optimization/matRad_prepareBiologicalOptimizationDij.m create mode 100644 test/optimization/test_biologicalOptimizationDij.m diff --git a/matRad/matRad_fluenceOptimization.m b/matRad/matRad_fluenceOptimization.m index 081f1ac20..2cc2c9988 100644 --- a/matRad/matRad_fluenceOptimization.m +++ b/matRad/matRad_fluenceOptimization.m @@ -182,32 +182,8 @@ backProjection = matRad_DoseProjection; end -% Check minimum biological quantities available -if isa(backProjection,'matRad_EffectProjection') && ~all(isfield(dij,{'ax','bx'})) - matRad_cfg.dispWarning('Biological optimization requested, but no ax & bx provided in dij. Getting from cst...'); - - %First get the voxels where we need it - validScen = ~cellfun(@isempty,dij.physicalDose); - d = cellfun(@(D) D*ones(dij.totalNumOfBixels,1),dij.physicalDose(validScen),'UniformOutput',false); - d = sum(cell2mat(d'),2); - ixZeroDose = d == 0; - - numOfCtScenarios = numel(cst{1,4}); - for i = 1:numOfCtScenarios - dij.ax{i} = zeros(dij.doseGrid.numOfVoxels,1); - dij.bx{i} = zeros(dij.doseGrid.numOfVoxels,1); - - for v = 1:size(cst,1) - %We already did the overlap stuff so we do not need to care for - %overlaps here - dij.ax{i}(cst{v,4}{i}) = cst{v,5}.alphaX; - dij.bx{i}(cst{v,4}{i}) = cst{v,5}.betaX; - end - - dij.ax{i}(ixZeroDose) = 0; - dij.bx{i}(ixZeroDose) = 0; - end -end +% Check and prepare biological quantities for all initialization paths. +dij = matRad_prepareBiologicalOptimizationDij(dij,cst,backProjection); % calculate initial beam intensities wInit @@ -217,15 +193,6 @@ %do nothing as wInit was passed to the function matRad_cfg.dispInfo('chosen provided wInit!\n'); - % Write ixDose which is needed for the optimizer - if isa(backProjection, 'matRad_EffectProjection') - dij.ixDose = dij.bx~=0; - - %pre-calculations - dij.gamma = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); - dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); - end - elseif isa(backProjection, 'matRad_ConstantRBEProjection') && strcmp(pln.radiationMode,'protons') % check if a constant RBE is defined - if not use 1.1 if ~isfield(dij,'RBE') @@ -238,13 +205,6 @@ matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); elseif isa(backProjection, 'matRad_EffectProjection') - % retrieve photon LQM parameter - [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels); - checkAxBx = cellfun(@(ax1,bx1,ax2,bx2) isequal(ax1(ax1~=0),ax2(ax1~=0)) && isequal(bx1(bx1~=0),bx2(bx1~=0)),dij.ax,dij.bx,ax,bx); - if ~all(checkAxBx) - matRad_cfg.dispError('Inconsistent biological parameters in dij.ax and/or dij.bx - please recalculate dose influence matrix before optimization!\n'); - end - for i = 1:size(cst,1) for j = 1:size(cst{i,6},2) @@ -256,10 +216,6 @@ end end - for s = 1:numel(dij.bx) - dij.ixDose{s} = dij.bx{s}~=0; - end - doseTmp = dij.physicalDose{1}*wOnes; if all(isfield(dij,{'mAlphaDose','mSqrtBetaDose'})) aTmp = dij.mAlphaDose{1}*wOnes; @@ -279,13 +235,6 @@ elseif isequal(pln.propOpt.quantityOpt,'RBExDose') - %pre-calculations - for s = 1:numel(dij.ixDose) - dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); - dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); - end - - % calculate current effect in target CurrEffectTarget = aTmp(V) + bTmp(V).^2; % ensure a underestimated biological effective dose diff --git a/matRad/optimization/matRad_prepareBiologicalOptimizationDij.m b/matRad/optimization/matRad_prepareBiologicalOptimizationDij.m new file mode 100644 index 000000000..056f8e6b8 --- /dev/null +++ b/matRad/optimization/matRad_prepareBiologicalOptimizationDij.m @@ -0,0 +1,194 @@ +function dij = matRad_prepareBiologicalOptimizationDij(dij, cst, backProjection) +% matRad_prepareBiologicalOptimizationDij prepares biological dij metadata +% +% call +% dij = matRad_prepareBiologicalOptimizationDij(dij,cst,backProjection) +% +% input +% dij: matRad dij struct +% cst: matRad cst struct on the dose grid +% backProjection: matRad backprojection object +% +% output +% dij: dij with biological ax/bx, ixDose and gamma metadata +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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 ~isa(backProjection, 'matRad_EffectProjection') + return +end + +matRad_cfg = MatRad_Config.instance(); +numOfVoxels = matRad_getDoseGridNumOfVoxels(dij); +numOfCtScen = matRad_getNumOfCtScenarios(dij, cst); + +axBxFromCst = ~all(isfield(dij, {'ax', 'bx'})); +if axBxFromCst + matRad_cfg.dispWarning(['Biological optimization requested, but no ' ... + 'ax & bx provided in dij. Getting from cst...']); + [dij.ax, dij.bx] = matRad_getPhotonLQMParameters(cst, numOfVoxels); + ixZeroDose = matRad_getZeroDoseVoxelMask(dij, numOfVoxels); + for scenIx = 1:numel(dij.ax) + dij.ax{scenIx}(ixZeroDose) = 0; + dij.bx{scenIx}(ixZeroDose) = 0; + end +else + dij.ax = matRad_asCtScenarioCell(dij.ax, numOfCtScen, ... + numOfVoxels, 'dij.ax'); + dij.bx = matRad_asCtScenarioCell(dij.bx, numOfCtScen, ... + numOfVoxels, 'dij.bx'); +end + +[dij.ax, dij.bx] = matRad_validateAxBxCells(dij.ax, dij.bx, ... + numOfCtScen, numOfVoxels); +if ~axBxFromCst + matRad_validateAxBxAgainstCst(dij.ax, dij.bx, cst, numOfVoxels); +end +dij.ixDose = cell(numOfCtScen, 1); +for scenIx = 1:numOfCtScen + dij.ixDose{scenIx} = dij.bx{scenIx} ~= 0; +end + +if isa(backProjection, 'matRad_VariableRBEProjection') + dij.gamma = cell(numOfCtScen, 1); + for scenIx = 1:numOfCtScen + dij.gamma{scenIx} = zeros(numOfVoxels, 1); + dij.gamma{scenIx}(dij.ixDose{scenIx}) = ... + dij.ax{scenIx}(dij.ixDose{scenIx}) ./ ... + (2 * dij.bx{scenIx}(dij.ixDose{scenIx})); + end +end + +end + +function numOfVoxels = matRad_getDoseGridNumOfVoxels(dij) +matRad_cfg = MatRad_Config.instance(); +if ~isstruct(dij) || ~isfield(dij, 'doseGrid') || ... + ~isfield(dij.doseGrid, 'numOfVoxels') || ... + ~isnumeric(dij.doseGrid.numOfVoxels) || ... + ~isscalar(dij.doseGrid.numOfVoxels) || ... + dij.doseGrid.numOfVoxels < 1 + matRad_cfg.dispError(['dij.doseGrid.numOfVoxels must be a positive ' ... + 'numeric scalar for biological optimization!']); +end +numOfVoxels = dij.doseGrid.numOfVoxels; +end + +function numOfCtScen = matRad_getNumOfCtScenarios(dij, cst) +matRad_cfg = MatRad_Config.instance(); +if isfield(dij, 'physicalDose') && iscell(dij.physicalDose) && ... + ~isempty(dij.physicalDose) + numOfCtScen = size(dij.physicalDose, 1); +elseif ~isempty(cst) && size(cst, 2) >= 4 && iscell(cst{1, 4}) + numOfCtScen = numel(cst{1, 4}); +else + matRad_cfg.dispError(['Unable to determine number of CT scenarios ' ... + 'for biological optimization!']); +end + +if numOfCtScen < 1 + matRad_cfg.dispError(['Biological optimization requires at least one ' ... + 'CT scenario!']); +end +end + +function values = matRad_asCtScenarioCell(values, numOfCtScen, numOfVoxels, fieldName) +matRad_cfg = MatRad_Config.instance(); +if iscell(values) + values = values(:); +elseif isnumeric(values) + if numOfCtScen == 1 && isvector(values) + values = {values(:)}; + elseif size(values, 1) == numOfVoxels && size(values, 2) == numOfCtScen + values = mat2cell(values, numOfVoxels, ones(1, numOfCtScen)); + values = values(:); + else + matRad_cfg.dispError(['%s must be a cell array by CT scenario or ' ... + 'a numeric array with one column per CT scenario!'], fieldName); + end +else + matRad_cfg.dispError(['%s must be numeric or a cell array of numeric ' ... + 'vectors!'], fieldName); +end +end + +function [ax, bx] = matRad_validateAxBxCells(ax, bx, numOfCtScen, numOfVoxels) +matRad_cfg = MatRad_Config.instance(); +if ~iscell(ax) || ~iscell(bx) || numel(ax) ~= numOfCtScen || ... + numel(bx) ~= numOfCtScen + matRad_cfg.dispError(['dij.ax and dij.bx must have one cell entry per ' ... + 'CT scenario!']); +end + +for scenIx = 1:numOfCtScen + matRad_validateBiologicalParameterVector(ax{scenIx}, numOfVoxels, ... + sprintf('dij.ax{%d}', scenIx)); + matRad_validateBiologicalParameterVector(bx{scenIx}, numOfVoxels, ... + sprintf('dij.bx{%d}', scenIx)); + ax{scenIx} = ax{scenIx}(:); + bx{scenIx} = bx{scenIx}(:); +end +end + +function matRad_validateBiologicalParameterVector(value, numOfVoxels, fieldName) +matRad_cfg = MatRad_Config.instance(); +if ~isnumeric(value) || ~isvector(value) || numel(value) ~= numOfVoxels + matRad_cfg.dispError(['%s must be a numeric vector with one entry per ' ... + 'dose-grid voxel!'], fieldName); +end +end + +function ixZeroDose = matRad_getZeroDoseVoxelMask(dij, numOfVoxels) +matRad_cfg = MatRad_Config.instance(); +if ~isfield(dij, 'physicalDose') || ~iscell(dij.physicalDose) + matRad_cfg.dispError(['dij.physicalDose is required to prepare ' ... + 'biological optimization quantities!']); +end + +validScen = find(~cellfun(@isempty, dij.physicalDose)); +if isempty(validScen) + matRad_cfg.dispError(['At least one physical dose influence matrix is ' ... + 'required to prepare biological optimization quantities!']); +end + +doseSum = zeros(numOfVoxels, 1); +for scenIx = validScen(:)' + doseMatrix = dij.physicalDose{scenIx}; + if size(doseMatrix, 1) ~= numOfVoxels + matRad_cfg.dispError(['dij.physicalDose{%d} has an inconsistent ' ... + 'number of dose-grid voxels!'], scenIx); + end + doseSum = doseSum + doseMatrix * ones(size(doseMatrix, 2), 1); +end +ixZeroDose = doseSum == 0; +end + +function matRad_validateAxBxAgainstCst(ax, bx, cst, numOfVoxels) +matRad_cfg = MatRad_Config.instance(); +[refAx, refBx] = matRad_getPhotonLQMParameters(cst, numOfVoxels); +if numel(ax) ~= numel(refAx) || numel(bx) ~= numel(refBx) + matRad_cfg.dispError(['Inconsistent number of biological CT scenarios ' ... + 'in dij.ax and/or dij.bx!']); +end + +isConsistent = cellfun(@(ax1, bx1, ax2, bx2) ... + isequal(ax1(ax1 ~= 0), ax2(ax1 ~= 0)) && ... + isequal(bx1(bx1 ~= 0), bx2(bx1 ~= 0)), ... + ax, bx, refAx, refBx); +if ~all(isConsistent) + matRad_cfg.dispError(['Inconsistent biological parameters in dij.ax ' ... + 'and/or dij.bx - please recalculate dose influence matrix before ' ... + 'optimization!']); +end +end diff --git a/test/optimization/test_biologicalOptimizationDij.m b/test/optimization/test_biologicalOptimizationDij.m new file mode 100644 index 000000000..e015bf84f --- /dev/null +++ b/test/optimization/test_biologicalOptimizationDij.m @@ -0,0 +1,116 @@ +function test_suite = test_biologicalOptimizationDij + +test_functions = localfunctions(); + +initTestSuite; + +function test_variableRBECellAxBxComputesIxDoseAndGamma +[dij, cst] = helper_biologicalDijAndCst(2); + +dij = matRad_prepareBiologicalOptimizationDij( ... + dij, cst, matRad_VariableRBEProjection); + +assertTrue(iscell(dij.ax)); +assertTrue(iscell(dij.bx)); +assertTrue(iscell(dij.ixDose)); +assertTrue(iscell(dij.gamma)); +assertEqual(find(dij.ixDose{1}), [1; 2; 3]); +assertEqual(find(dij.ixDose{2}), [1; 2; 3]); +assertElementsAlmostEqual(dij.gamma{1}(1:3), ... + [2; 5 / 3; 2], 'absolute', 1e-12); +assertElementsAlmostEqual(dij.gamma{2}(1:3), ... + [5 / 3; 2; 2], 'absolute', 1e-12); + +function test_effectCellAxBxComputesIxDoseWithoutGamma +[dij, cst] = helper_biologicalDijAndCst(1); + +dij = matRad_prepareBiologicalOptimizationDij( ... + dij, cst, matRad_EffectProjection); + +assertTrue(iscell(dij.ixDose)); +assertEqual(find(dij.ixDose{1}), [1; 2; 3]); +assertFalse(isfield(dij, 'gamma')); + +function test_missingAxBxIsRebuiltFromCst +[dij, cst] = helper_biologicalDijAndCst(2); +dij = rmfield(dij, {'ax', 'bx'}); + +dij = matRad_prepareBiologicalOptimizationDij( ... + dij, cst, matRad_BEDProjection); + +assertTrue(iscell(dij.ax)); +assertTrue(iscell(dij.bx)); +assertElementsAlmostEqual(dij.ax{1}, [0.2; 0.1; 0.2; 0], ... + 'absolute', 1e-12); +assertElementsAlmostEqual(dij.bx{2}, [0.03; 0.05; 0.05; 0], ... + 'absolute', 1e-12); +assertEqual(find(dij.ixDose{1}), [1; 2; 3]); + +function test_numericSingleCtAxBxIsNormalizedToCell +[dij, cst] = helper_biologicalDijAndCst(1); +dij.ax = dij.ax{1}; +dij.bx = dij.bx{1}'; + +dij = matRad_prepareBiologicalOptimizationDij( ... + dij, cst, matRad_EffectProjection); + +assertTrue(iscell(dij.ax)); +assertTrue(iscell(dij.bx)); +assertEqual(size(dij.ax{1}), [4 1]); +assertEqual(size(dij.bx{1}), [4 1]); + +function test_invalidAxBxSizeFailsClearly +[dij, cst] = helper_biologicalDijAndCst(1); +dij.ax{1} = dij.ax{1}(1:3); + +assertExceptionThrown(@() matRad_prepareBiologicalOptimizationDij( ... + dij, cst, matRad_EffectProjection), 'matRad:Error'); + +function [dij, cst] = helper_biologicalDijAndCst(numOfCtScen) +numOfVoxels = 4; +dij = struct(); +dij.doseGrid.numOfVoxels = numOfVoxels; +dij.totalNumOfBixels = 2; +dij.numOfScenarios = numOfCtScen; +dij.physicalDose = cell(numOfCtScen, 1, 1); +dij.ax = cell(numOfCtScen, 1); +dij.bx = cell(numOfCtScen, 1); + +cst = helper_biologicalCst(numOfCtScen); +for ctScen = 1:numOfCtScen + doseMatrix = sparse(numOfVoxels, 2); + doseMatrix(cst{1, 4}{ctScen}, 1) = 1; + doseMatrix(cst{2, 4}{ctScen}, 2) = 1; + dij.physicalDose{ctScen, 1, 1} = doseMatrix; + dij.ax{ctScen} = zeros(numOfVoxels, 1); + dij.bx{ctScen} = zeros(numOfVoxels, 1); + dij.ax{ctScen}(cst{1, 4}{ctScen}) = cst{1, 5}.alphaX; + dij.bx{ctScen}(cst{1, 4}{ctScen}) = cst{1, 5}.betaX; + dij.ax{ctScen}(cst{2, 4}{ctScen}) = cst{2, 5}.alphaX; + dij.bx{ctScen}(cst{2, 4}{ctScen}) = cst{2, 5}.betaX; +end + +function cst = helper_biologicalCst(numOfCtScen) +cst = cell(2, 6); +cst{1, 1} = 1; +cst{1, 2} = 'TARGET'; +cst{1, 3} = 'TARGET'; +cst{1, 4} = cell(1, numOfCtScen); +cst{1, 5} = struct('alphaX', 0.2, 'betaX', 0.05); +cst{1, 6} = {}; +cst{2, 1} = 2; +cst{2, 2} = 'OAR'; +cst{2, 3} = 'OAR'; +cst{2, 4} = cell(1, numOfCtScen); +cst{2, 5} = struct('alphaX', 0.1, 'betaX', 0.03); +cst{2, 6} = {}; + +for ctScen = 1:numOfCtScen + if ctScen == 1 + cst{1, 4}{ctScen} = [1; 3]; + cst{2, 4}{ctScen} = 2; + else + cst{1, 4}{ctScen} = [2; 3]; + cst{2, 4}{ctScen} = 1; + end +end