Skip to content

Implementation of multi alpha/beta for proton plans with TOPAS #884

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
SimonaFa wants to merge 4 commits into
base: dev
Choose a base branch
from
+143 −29
Open

Implementation of multi alpha/beta for proton plans with TOPAS #884

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f0fc4bb
Changes to topas interface to enable multiple radiosensitivities
SimonaFa Nov 3, 2025
42d0e7c
Writing and reading modification to TOPAS engine for multi alpha beta…
SimonaFa Dec 1, 2025
dba9f96
Multi Alpha/Beta implementation concept for proton plans with TOPAS
SimonaFa Feb 2, 2026
6d649e8
Merge branch 'dev' into dev_TOPASinterface_multiAlphaBeta
SimonaFa Feb 2, 2026
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
169 changes: 140 additions & 29 deletions matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,14 +246,27 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
end
end

% Get alpha beta parameters from bioParam struct
if isfield(obj.bioParameters, 'tissuseAlphaX')
obj.bioParameters.AlphaX = obj.bioModel.tissueAlphaX(1);
obj.bioParameters.BetaX = obj.bioModel.tissueBetaX(1);
tmpAlphaX = [];
tmpBetaX = [];
for idx = 1:size(cst, 1)
if ~isempty(cst{idx,5}) && isfield(cst{idx,5}, 'alphaX')
tmpAlphaX = [tmpAlphaX cst{idx,5}.alphaX];
tmpBetaX = [tmpBetaX cst{idx,5}.betaX];
end
end
if numel(obj.bioParameters.AlphaX)>1
matRad_cfg.dispWarning('!!! Only a unique alpha/beta ratio supported at the moment. Found multiple, only the first one will be used !!!!');
abX = [tmpAlphaX(:) tmpBetaX(:)];
unique_abX = unique(abX, 'rows');
obj.bioParameters.AlphaX = unique_abX(:, 1);
obj.bioParameters.BetaX = unique_abX(:, 2);

% Get alpha beta parameters from bioParam struct
if isfield(obj.bioParameters, 'tissueAlphaX')
Comment on lines +262 to +263
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The condition checks isfield(obj.bioParameters,'tissueAlphaX'), but the values are then read from obj.bioModel.tissueAlphaX/BetaX (and bioModel is a class, not a struct). As written, this branch will never run, so tabulated-model tissue arrays won’t be used. Consider checking isprop(obj.bioModel,'tissueAlphaX') (or isa(obj.bioModel,'matRad_LQRBETabulatedModel')) and then copying obj.bioModel.tissueAlphaX/tissueBetaX into obj.bioParameters.

Suggested change
% Get alpha beta parameters from bioParam struct
if isfield(obj.bioParameters, 'tissueAlphaX')
% Get alpha beta parameters from bioModel if it provides tissue-wise values
if isprop(obj.bioModel, 'tissueAlphaX') && isprop(obj.bioModel, 'tissueBetaX')

Copilot uses AI. Check for mistakes.
obj.bioParameters.AlphaX = obj.bioModel.tissueAlphaX;
obj.bioParameters.BetaX = obj.bioModel.tissueBetaX;
end
%if numel(obj.bioParameters.AlphaX)>1
% matRad_cfg.dispWarning('!!! Only a unique alpha/beta ratio supported at the moment. Found multiple, only the first one will be used !!!!');
%end

end
if obj.scorer.LET
Expand Down Expand Up @@ -508,10 +521,15 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
end

% Get photon parameters for RBExDose calculation
if this.calcBioDose
if this.calcBioDose || this.scorer.RBE
this.scorer.RBE = true;
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);
dij.abx(dij.bx>0) = dij.ax(dij.bx>0)./dij.bx(dij.bx>0);
this.calcBioDose = true;
[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,this.VdoseGrid);
dij.abx = zeros(size(dij.ax{1}));
ax = dij.ax{1};
bx = dij.bx{1};
dij.abx(bx>0) = ax(bx>0)./bx(bx>0);
dij.abx = {dij.abx};
Comment on lines +528 to +532
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

dij.abx is computed only for scenario 1 (dij.ax{1}/dij.bx{1}) and then wrapped into a single-element cell. If multiple CT scenarios exist, dij.abx won’t match the size/ordering of dij.ax/dij.bx. Consider computing abx for each scenario (similar to how obj.MCparam.abx is built earlier) and storing it as a cell array with numOfCtScen entries.

Suggested change
dij.abx = zeros(size(dij.ax{1}));
ax = dij.ax{1};
bx = dij.bx{1};
dij.abx(bx>0) = ax(bx>0)./bx(bx>0);
dij.abx = {dij.abx};
% Compute alpha/beta (abx) for each CT scenario, matching dij.ax/dij.bx
numCtScen = numel(dij.ax);
dij.abx = cell(numCtScen,1);
for ctScenIdx = 1:numCtScen
ax = dij.ax{ctScenIdx};
bx = dij.bx{ctScenIdx};
abx = zeros(size(ax));
mask = bx > 0;
abx(mask) = ax(mask) ./ bx(mask);
dij.abx{ctScenIdx} = abx;
end

Copilot uses AI. Check for mistakes.
end

% save current directory to revert back to later
Expand Down Expand Up @@ -1171,16 +1189,39 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
end
% Handle RBE-related quantities (not multiplied by sum(w)!)
elseif ~isempty(strfind(lower(topasCubesTallies{j}),'alpha'))
modelName = strsplit(topasCubesTallies{j},'_');
modelName = modelName{end};
talliesFlags = strsplit(topasCubesTallies{j},'_');
modelName = talliesFlags{end};
if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)]))
dij.(['mAlphaDose_' modelName]){ctScen}(:,d) = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1) .* dij.physicalDose{ctScen}(:,d);

if contains(topasCubesTallies{j}, 'CellType')
ab_idx = str2num(talliesFlags{2});
organAlpha = obj.bioParameters.AlphaX(ab_idx);
organBeta = obj.bioParameters.BetaX(ab_idx);
mask = find( (dij.ax{1} == organAlpha) & (dij.bx{1} == organBeta));
topasCube_values = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1);
topasCube_values = topasCube_values(mask,d);
dij.(['mAlphaDose_' modelName]){ctScen}(mask,d) = topasCube_values .* dij.physicalDose{ctScen}(mask,d);
Comment on lines +1200 to +1203
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the CellType branch, topasCube_values is reshaped to a vector, but then indexed as topasCube_values(mask,d). This will error because the reshaped array is 1-D. Use 1-D indexing (topasCube_values(mask)) and keep the assignment to dij...(...)(mask,d) for the beam column.

Copilot uses AI. Check for mistakes.
else

dij.(['mAlphaDose_' modelName]){ctScen}(:,d) = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1) .* dij.physicalDose{ctScen}(:,d);
end
end
elseif ~isempty(strfind(lower(topasCubesTallies{j}),'beta'))
modelName = strsplit(topasCubesTallies{j},'_');
modelName = modelName{end};
if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)]))
if contains(topasCubesTallies{j}, 'CellType')
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the beta-tally CellType branch, talliesFlags{2} is used to derive the CellType index, but talliesFlags is not defined in this branch (it’s only created in the alpha branch). This can fail depending on tally iteration order. Define talliesFlags = strsplit(topasCubesTallies{j},'_') locally before accessing talliesFlags{2}.

Suggested change
if contains(topasCubesTallies{j}, 'CellType')
if contains(topasCubesTallies{j}, 'CellType')
talliesFlags = strsplit(topasCubesTallies{j},'_');

Copilot uses AI. Check for mistakes.
ab_idx = str2num(talliesFlags{2});
organAlpha = obj.bioParameters.AlphaX(ab_idx);
Comment on lines +1192 to +1215
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new use of str2num(talliesFlags{2}) on topasCubesTallies introduces a code‑injection risk, because talliesFlags is derived from TOPAS tally names that ultimately come from file names in an external simulation folder and may be attacker‑controlled (e.g. when using readExternal on untrusted output). Since str2num evaluates its input via eval, a crafted tally/file name like alpha_CellType=1;system('malicious_command');1 can cause arbitrary MATLAB commands to execute when this code parses the cell‑type index. To fix this, avoid str2num entirely here (e.g. use str2double or explicit numeric parsing on a constrained pattern) and ensure that only strictly numeric indices are accepted from tally names.

Copilot uses AI. Check for mistakes.
organBeta = obj.bioParameters.BetaX(ab_idx);
mask = find( (dij.ax{1} == organAlpha) & (dij.bx{1} == organBeta));
topasCube_values = reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1);
topasCube_values = topasCube_values(mask,d);
dij.(['mBetaDose_' modelName]){ctScen}(mask,d) = topasCube_values .* dij.physicalDose{ctScen}(mask,d);
else

dij.(['mSqrtBetaDose_' modelName]){ctScen}(:,d) = sqrt(reshape(topasCubes.([topasCubesTallies{j} '_beam',num2str(d)]){ctScen},[],1)) .* dij.physicalDose{ctScen}(:,d);
Comment on lines +1219 to 1223
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this CellType beta branch, topasCube_values is reshaped to a vector but then indexed as topasCube_values(mask,d), which will error for 1-D data. Additionally, the assignment targets mBetaDose_..., but prepareDij() only allocates mSqrtBetaDose_... and other paths store sqrt(beta) in that field. Use 1-D indexing and align with the existing mSqrtBetaDose_... + sqrt(...) convention.

Copilot uses AI. Check for mistakes.
end
end
elseif ~isempty(strfind(topasCubesTallies{j},'LET'))
if isfield(topasCubes,[topasCubesTallies{j} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} '_beam' num2str(d)]))
Expand Down Expand Up @@ -1355,7 +1396,31 @@ function writeScorers(obj,fID,beamIx)
% Read appropriate scorer from file and write to config file
matRad_cfg.dispDebug('Reading RBE Scorer from %s\n',fname);
scorerName = fileread(fname);
fprintf(fID,'\n%s\n\n',scorerName);
generalScorer = scorerName;

if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'\n%s\n\n',scorerName);
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
scorerName = generalScorer;
insertText = ['_CellType_' num2str(idxCell)];
scorerName = strrep(scorerName, ...
'Alpha/', ['Alpha' insertText '/']);
scorerName = strrep(scorerName, ...
'Beta/', ['Beta' insertText '/']);
scorerName = strrep(scorerName, ...
'Sc/PrescribedDose Gy', ['Sc/PrescribedDose' insertText ' Gy']);
scorerName = strrep(scorerName, ...
'Sc/CellLines', ['Sc/CellLines' insertText]);
scorerName = strrep(scorerName, ...
'Sc/SimultaneousExposure', ['Sc/SimultaneousExposure' insertText]);
scorerName = strrep(scorerName, ...
'Sim/ScoreLabel + "', ['Sim/ScoreLabel + "' insertText]);

fprintf(fID,'\n%s\n\n',scorerName);
end

end

if obj.calc4DInterplay
for PhaseNum = obj.MCparam.Phases{beamIx}'
Expand Down Expand Up @@ -1402,28 +1467,60 @@ function writeScorers(obj,fID,beamIx)
% Begin writing biological scorer components: cell lines
switch obj.radiationMode
case 'protons'
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID,'sv:Sc/CellLines = 1 "CellLineGeneric"\n');
fprintf(fID,'d:Sc/CellLineGeneric/Alphax = Sc/AlphaX /Gy\n');
fprintf(fID,'d:Sc/CellLineGeneric/Betax = Sc/BetaX /Gy2\n');
fprintf(fID,'d:Sc/CellLineGeneric/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n\n');
if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID,'sv:Sc/CellLines = 1 "CellLineGeneric"\n');
fprintf(fID,'d:Sc/CellLineGeneric/Alphax = Sc/AlphaX /Gy\n');
fprintf(fID,'d:Sc/CellLineGeneric/Betax = Sc/BetaX /Gy2\n');
fprintf(fID,'d:Sc/CellLineGeneric/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n\n');
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
insertText = ['_CellType_' num2str(idxCell)];
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID, ['sv:Sc/CellLines' insertText ' = 1 "CellLineGeneric' insertText '"\n']);
fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/Alphax = Sc/AlphaX' insertText ' /Gy\n']);
fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/Betax = Sc/BetaX' insertText ' /Gy2\n']);
fprintf(fID, ['d:Sc/CellLineGeneric' insertText '/AlphaBetaRatiox = Sc/AlphaBetaX' insertText ' Gy\n\n']);
end
end
case {'carbon','helium'}
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID,'sv:Sc/CellLines = 1 "CellGeneric_abR2"\n');
fprintf(fID,'d:Sc/CellGeneric_abR2/Alphax = Sc/AlphaX /Gy\n');
fprintf(fID,'d:Sc/CellGeneric_abR2/Betax = Sc/BetaX /Gy2\n\n');
% fprintf(fID,'d:Sc/CellGeneric_abR2/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n');
if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID,'sv:Sc/CellLines = 1 "CellGeneric_abR2"\n');
fprintf(fID,'d:Sc/CellGeneric_abR2/Alphax = Sc/AlphaX /Gy\n');
fprintf(fID,'d:Sc/CellGeneric_abR2/Betax = Sc/BetaX /Gy2\n\n');
% fprintf(fID,'d:Sc/CellGeneric_abR2/AlphaBetaRatiox = Sc/AlphaBetaX Gy\n');
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
insertText = ['_CellType_' num2str(idxCell)];
fprintf(fID,'\n### Biological Parameters ###\n');
fprintf(fID, ['sv:Sc/CellLines' insertText ' = 1 "CellLineGeneric_abR2' insertText '"\n']);
fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Alphax = Sc/AlphaX' inserText ' /Gy\n']);
Copy link

Copilot AI Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo: inserText is undefined (should be insertText). This will throw at runtime when writing the carbon/helium multi-CellType scorer section.

Suggested change
fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Alphax = Sc/AlphaX' inserText ' /Gy\n']);
fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Alphax = Sc/AlphaX' insertText ' /Gy\n']);

Copilot uses AI. Check for mistakes.
fprintf(fID, ['d:Sc/CellLineGeneric_abR2' insertText '/Betax = Sc/BetaX' insertText ' /Gy2\n\n']);
end
end
otherwise
matRad_cfg.dispError([obj.radiationMode ' not implemented']);
end

% write biological scorer components: dose parameters
matRad_cfg.dispDebug('Writing Biologial Scorer components.\n');
fprintf(fID,'d:Sc/PrescribedDose = %.4f Gy\n',obj.bioParameters.PrescribedDose);
fprintf(fID,'b:Sc/SimultaneousExposure = %s\n',obj.bioParameters.SimultaneousExposure);
fprintf(fID,'d:Sc/AlphaX = %.4f /Gy\n',obj.bioParameters.AlphaX);
fprintf(fID,'d:Sc/BetaX = %.4f /Gy2\n',obj.bioParameters.BetaX);
fprintf(fID,'d:Sc/AlphaBetaX = %.4f Gy\n',obj.bioParameters.AlphaX/obj.bioParameters.BetaX);
if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'d:Sc/PrescribedDose = %.4f Gy\n',obj.bioParameters.PrescribedDose);
fprintf(fID,'b:Sc/SimultaneousExposure = %s\n',obj.bioParameters.SimultaneousExposure);
fprintf(fID,'d:Sc/AlphaX = %.4f /Gy\n',obj.bioParameters.AlphaX);
fprintf(fID,'d:Sc/BetaX = %.4f /Gy2\n',obj.bioParameters.BetaX);
fprintf(fID,'d:Sc/AlphaBetaX = %.4f Gy\n',obj.bioParameters.AlphaX/obj.bioParameters.BetaX);
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
insertText = ['_CellType_' num2str(idxCell)];
fprintf(fID, ['d:Sc/PrescribedDose' insertText ' = %.4f Gy\n'],obj.bioParameters.PrescribedDose);
fprintf(fID, ['b:Sc/SimultaneousExposure' insertText ' = %s\n'],obj.bioParameters.SimultaneousExposure);
fprintf(fID, ['d:Sc/AlphaX' insertText ' = %.4f /Gy\n'],obj.bioParameters.AlphaX(idxCell));
fprintf(fID, ['d:Sc/BetaX' insertText ' = %.4f /Gy2\n'],obj.bioParameters.BetaX(idxCell));
fprintf(fID, ['d:Sc/AlphaBetaX' insertText ' = %.4f Gy\n\n'],obj.bioParameters.AlphaX(idxCell)/obj.bioParameters.BetaX(idxCell));
end
end

% Update MCparam.tallies with processed scorer
for i = 1:length(obj.scorer.RBE_model)
Expand Down Expand Up @@ -1451,9 +1548,23 @@ function writeScorers(obj,fID,beamIx)
% Write subscorer to config files
for s = 1:length(scorerNames)
if strcmp(obj.radiationMode,'protons')
fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_LET = "ProtonLET"\n',scorerPrefix,scorerNames{s});
if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_LET = "ProtonLET"\n',scorerPrefix,scorerNames{s});
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
insertText = ['_CellType_' num2str(idxCell)];
fprintf(fID,['s:Sc/%s%s' insertText '/ReferencedSubScorer_LET = "ProtonLET"\n'],scorerPrefix,scorerNames{s});
end
end
end
if length(obj.bioParameters.AlphaX) ==1
fprintf(fID,'s:Sc/%s%s/ReferencedSubScorer_Dose = "Tally_DoseToWater"\n',scorerPrefix,scorerNames{s});
else
for idxCell = 1:length(obj.bioParameters.AlphaX)
insertText = ['_CellType_' num2str(idxCell)];
fprintf(fID,['s:Sc/%s%s' insertText '/ReferencedSubScorer_Dose = "Tally_DoseToWater"\n'],scorerPrefix,scorerNames{s});
end
end
if obj.calc4DInterplay
for PhaseNum = obj.MCparam.Phases{beamIx}'
if strcmp(obj.radiationMode,'protons')
Expand Down
3 changes: 3 additions & 0 deletions matRad/planAnalysis/matRad_compareDose.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,9 @@
% Calculate absolute difference cube and dose windows for plots
differenceCube = cube1-cube2;
doseDiffWindow = [-max(abs(differenceCube(:))) max(abs(differenceCube(:)))];
if doseDiffWindow(1)==0 && doseDiffWindow(2)==0
doseDiffWindow(2) = 1;
end
%doseGammaWindow = [0 max(gammaCube(:))];
doseGammaWindow = [0 2]; %We choose 2 as maximum value since the gamma colormap has a sharp cut in the middle

Expand Down
Loading