Skip to content
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

Update to the Photon MLC for TOPAS #699

Closed
wants to merge 14 commits into from
Closed

Update to the Photon MLC for TOPAS #699

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
f9b6294
Create matRad_example15_photonMC_MLC.m
JenHardt Sep 1, 2023
dd82e79
small bug fix
JenHardt Sep 1, 2023
7c21e98
Update matRad_calcPhotonDoseMCtopas.m
JenHardt Sep 1, 2023
cdbcca5
MLC updates
JenHardt Sep 1, 2023
18f33fd
Update matRad_example15_photonMC_MLC.m
JenHardt Sep 1, 2023
2ece06e
Update matRad_TopasConfig.m
JenHardt Sep 1, 2023
2184a10
small bugfix
JenHardt Sep 4, 2023
9865942
Update matRad_TopasConfig.m
JenHardt Sep 4, 2023
146ddd3
Update matRad_TopasConfig.m
JenHardt Sep 6, 2023
da650b1
update
JenHardt Oct 27, 2023
a6e679a
Update matRad_TopasConfig.m
JenHardt Oct 30, 2023
cdf7cd3
use other phase space file and adjustments to that
JenHardt Jan 26, 2024
cc42adf
small changes
JenHardt Feb 1, 2024
c2c293d
Update matRad_TopasConfig.m
JenHardt Feb 28, 2024
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
97 changes: 97 additions & 0 deletions examples/matRad_example15_photonMC_MLC.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
%% Example: Photon Treatment Plan using VMC++ dose calculation
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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/LICENSES.txt. 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% In this example we will show
% (i) how to load patient data into matRad
% (ii) how to setup a photon dose calculation based on the VMC++ Monte Carlo algorithm
% (iii) how to inversely optimize the beamlet intensities directly from command window in MATLAB.
% (iv) how to visualize the result

%% set matRad runtime configuration
matRad_rc %If this throws an error, run it from the parent directory first to set the paths

%% Patient Data Import
% Let's begin with a clear Matlab environment and import the boxphantom
% into your workspace.
load('BOXPHANTOM.mat');

%% Treatment Plan
% The next step is to define your treatment plan labeled as 'pln'. This
% structure requires input from the treatment planner and defines the most
% important cornerstones of your treatment plan.

pln.radiationMode = 'photons';
pln.machine = 'Generic';
pln.numOfFractions = 30;
pln.propStf.gantryAngles = [0:72:359];
pln.propStf.couchAngles = [0 0 0 0 0];
%pln.propStf.gantryAngles = [0];
%pln.propStf.couchAngles = [0];
pln.propStf.bixelWidth = 10;
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
% Enable sequencing and direct aperture optimization (DAO).
pln.propOpt.runSequencing = 1;
pln.propOpt.runDAO = 1;

quantityOpt = 'physicalDose';
modelName = 'none';

% retrieve bio model parameters
pln.bioParam = matRad_bioModel(pln.radiationMode,quantityOpt, modelName);

% retrieve scenarios for dose calculation and optimziation
pln.multScen = matRad_multScen(ct,'nomScen');
% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]

%%
pln.propMC.engine = 'TOPAS';
pln.propMC.beamProfile = 'phasespace';
pln.propMC.externalCalculation =true;
%% Generate Beam Geometry STF
stf = matRad_generateStf(ct,cst,pln);

%% Dose Calculation
dij = matRad_calcPhotonDose(ct,stf,pln,cst);

%% Inverse Optimization for IMRT
resultGUI = matRad_fluenceOptimization(dij,cst,pln);
%% Sequencing
% This is a multileaf collimator leaf sequencing algorithm that is used in
% order to modulate the intensity of the beams with multiple static
% segments, so that translates each intensity map into a set of deliverable
% aperture shapes.
resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,5,1);
[pln,stf] = matRad_aperture2collimation(pln,stf,resultGUI.sequencing,resultGUI.apertureInfo);
%% Aperture visualization
% Use a matrad function to visualize the resulting aperture shapes
matRad_visApertureInfo(resultGUI.apertureInfo)
%% Plot the Resulting Dose Slice
% Just let's plot the transversal iso-center dose slice
slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z);
figure,
imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet)

%% Dose Calculation
pln.propMC.numHistories = 1e8;
resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1);

%% readout
%foldername = 'FolderName';
%pln = matRad_cfg.getDefaultClass(pln,'propMC');
%resultGUI_MC = pln.propMC.readExternal(foldername);
2 changes: 1 addition & 1 deletion matRad_calcDoseInit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
~isfield(pln.propDoseCalc,'doseGrid') || ...
~isfield(pln.propDoseCalc.doseGrid,'resolution')
% default values
dij.doseGrid.resolution = matRad_cfg.propDoseCalc.defaultResolution;
dij.doseGrid.resolution = matRad_cfg.propDoseCalc.doseGrid.defaultResolution;
else
% take values from pln strcut
dij.doseGrid.resolution.x = pln.propDoseCalc.doseGrid.resolution.x;
Expand Down
2 changes: 1 addition & 1 deletion matRad_calcPhotonDoseMC.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@

matRad_cfg = MatRad_Config.instance();

if nargin < 6
if nargin < 5
calcDoseDirect = false;
end

Expand Down
21 changes: 17 additions & 4 deletions matRad_calcPhotonDoseMCtopas.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,21 +79,34 @@

load([pln.radiationMode,'_',pln.machine]);

for i = 1:size(stf,2)
stf(i).SCD = machine.meta.SCD;
end


%Collect weights
if calcDoseDirect
if calcDoseDirect
w = zeros(sum([stf(:).totalNumOfBixels]),1);
ct = 1;
counter = 1;
for i = 1:length(stf)
for j = 1:stf(i).numOfRays
rayBix = stf(i).numOfBixelsPerRay(j);
w(ct:ct+rayBix-1) = stf(i).ray(j).weight;
ct = ct + rayBix;
if isfield(stf(1).ray, 'shapes')
w(counter:counter+rayBix-1) = [stf(i).ray.shapes.weight];
else
w(counter:counter+rayBix-1) = stf(i).ray(j).weight;
end
counter = counter + rayBix;
end
end
else
w = ones(sum([stf(:).totalNumOfBixels]),1);
end

for i = 1:size(stf,2)
stf(i).ray.energy = stf(i).ray.energy.*ones(size(w));
end

currDir = cd;

for shiftScen = 1:pln.multScen.totNumShiftScen
Expand Down
5 changes: 3 additions & 2 deletions topas/beamSetup/TOPAS_beamSetup_mlc.txt.in
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
s:Ge/MultiLeafCollimatorA/Type = "TsMultiLeafCollimator"
s:Ge/MultiLeafCollimatorA/Parent = "Nozzle"
s:Ge/MultiLeafCollimatorA/Material = "Aluminum"
s:Ge/MultiLeafCollimatorA/Material = "G4_W"
d:Ge/MultiLeafCollimatorA/TransX = 0.0 cm
d:Ge/MultiLeafCollimatorA/TransY = 0.0 cm
d:Ge/MultiLeafCollimatorA/TransZ = 0.5 * Ge/MultiLeafCollimatorA/Thickness cm
d:Ge/MultiLeafCollimatorA/TransZ = Sim/Ge/MultiLeafCollimatorA/TransZ cm
d:Ge/MultiLeafCollimatorA/RotX = 0.0 deg
d:Ge/MultiLeafCollimatorA/RotY = 0.0 deg
d:Ge/MultiLeafCollimatorA/RotZ = 90.0 deg
s:Ge/MultiLeafCollimatorA/DrawingStyle = "Solid"
2 changes: 1 addition & 1 deletion topas/beamSetup/TOPAS_beamSetup_phasespace.txt.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,6 @@ s:So/Phasespace/Component = "Nozzle"
#b:So/Example/LimitedAssumeEveryParticleIsNewHistory = "true"
#b:So/Example/LimitedAssumePhotonIsNewHistory = "true"
i:So/Phasespace/PhaseSpaceMultipleUse = 0
i:So/Phasespace/NumberOfHistoriesInRun = 5
i:So/Phasespace/NumberOfHistoriesInRun = Tf/Phasespace/NumberOfHistoriesInRun/Value
i:So/Phasespace/PreCheckShowParticleCountAtInterval = 100000
b:So/Phasespace/PhaseSpacePreCheck = "False"
89 changes: 63 additions & 26 deletions topas/matRad_TopasConfig.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,10 @@
'Scorer_RBE_libamtrack','TOPAS_scorer_doseRBE_libamtrack.txt.in',...
'Scorer_RBE_LEM1','TOPAS_scorer_doseRBE_LEM1.txt.in',...
'Scorer_RBE_WED','TOPAS_scorer_doseRBE_Wedenberg.txt.in',...
'Scorer_RBE_MCN','TOPAS_scorer_doseRBE_McNamara.txt.in');
'Scorer_RBE_MCN','TOPAS_scorer_doseRBE_McNamara.txt.in', ...
... %PhaseSpace Source
'phaseSpaceSourcePhotons' ,'VarianClinaciX_6MV_20x20_aboveMLC_w2' );


end

Expand Down Expand Up @@ -689,6 +692,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w)
% Allocate possible scored quantities
processedQuantities = {'','_std','_batchStd'};
topasCubesTallies = unique(erase(topasCubesTallies,processedQuantities(2:end)));


% Loop through 4D scenarios
for ctScen = 1:dij.numOfScenarios
Expand Down Expand Up @@ -1380,8 +1384,8 @@ function writeStfFields(obj,ct,stf,pln,w,baseData)
end

% NozzleAxialDistance
if isPhoton
fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', 1000 + ct.cubeDim(3)*ct.resolution.z);%Not sure if this is correct,100 cm is SSD and probably distance from surface to isocenter needs to be added
if isPhoton
fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', stf(beamIx).SCD+40); %Phasespace hardcorded infront of MLC at SSD 46 cm
else
fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', nozzleToAxisDistance);
end
Expand Down Expand Up @@ -1569,20 +1573,20 @@ function writeStfFields(obj,ct,stf,pln,w,baseData)
end

case 'phasespace'

fprintf(fileID,'d:Sim/GantryAngle = %f deg\n',stf(beamIx).gantryAngle); %just one beam angle for now
fprintf(fileID,'d:Sim/CouchAngle = %f deg\n',stf(beamIx).couchAngle);
% Here the phasespace file is loaded and referenced in the beamSetup file
phaseSpaceFileName = 'SIEMENS_PRIMUS_6.0_0.10_15.0x15.0';
% Here the phasespace file is loaded and referenced in the beamSetup file
if obj.externalCalculation
matRad_cfg.dispWarning(['External calculation and phaseSpace selected, manually place ' phaseSpaceFileName '.header and ' phaseSpaceFileName '.phsp into your simulation directory.']);
matRad_cfg.dispWarning(['External calculation and phaseSpace selected, manually place ' obj.infilenames.phaseSpaceSourcePhotons '.header and ' obj.infilenames.phaseSpaceSourcePhotons '.phsp into your simulation directory.']);
else
if length(dir([obj.thisFolder filesep 'beamSetup' filesep 'phasespace' filesep phaseSpaceFileName '*'])) < 2
if length(dir([obj.thisFolder filesep 'beamSetup' filesep 'phasespace' filesep obj.infilenames.phaseSpaceSourcePhotons '*'])) < 2
matRad_cfg.dispError([phaseSpaceFileName ' header or phsp file could not be found in beamSetup/phasespace folder.']);
end
end
phasespaceStr = ['..' filesep 'beamSetup' filesep 'phasespace' filesep phaseSpaceFileName];
phasespaceStr = replace(phasespaceStr, '\', '/');
fprintf(fileID,'s:So/Phasespace/PhaseSpaceFileName = "%s"\n', phasespaceStr);
%phasespaceStr = ['..' filesep 'beamSetup' filesep 'phasespace' filesep phaseSpaceFileName];
%&phasespaceStr = replace(phasespaceStr, '\', '/');
fprintf(fileID,'s:So/Phasespace/PhaseSpaceFileName = "%s"\n', obj.infilenames.phaseSpaceSourcePhotons );

end

Expand Down Expand Up @@ -1644,49 +1648,82 @@ function writeStfFields(obj,ct,stf,pln,w,baseData)

% Write MLC if available
if isfield(stf(beamIx).ray, 'shapes')
SCD = stf(beamIx).SCD;
fname = fullfile(obj.thisFolder,obj.infilenames.beam_mlc);
TOPAS_mlcSetup = fileread(fname);

fprintf(fileID,'%s\n',TOPAS_mlcSetup);
% For now only for one ray
[numOfLeaves,leafTimes]=size([stf(beamIx).ray.shapes(:).leftLeafPos]); %there are #numOfLeaves leaves and #leafTimes times/shapes
leftLeafPos = [stf(beamIx).ray.shapes(:).leftLeafPos];
rightLeafPos = [stf(beamIx).ray.shapes(:).rightLeafPos];
leftLeafPos = [stf(beamIx).ray.shapes(:).leftLeafPos]*SCD./SAD;
rightLeafPos = [stf(beamIx).ray.shapes(:).rightLeafPos]*SCD./SAD;
% Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator
fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n', 4);
fprintf(fileID,'d:Ge/MultiLeafCollimatorA/MaximumLeafOpen = %f cm\n',15);
fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Thickness = %f cm\n',15);
fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Length = %f cm\n',6);
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves);
fprintf(fileID,num2str(pln.propStf.collimation.leafWidth*ones(1,numOfLeaves),' % 2d'));
fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Thickness = %f cm\n',8);
fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Length = %f cm\n',15);
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves+2);
fprintf(fileID, '%f ', [200, pln.propStf.collimation.leafWidth*ones(1,numOfLeaves)*SCD./SAD , 200]);
fprintf(fileID,' mm \n');
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',numOfLeaves);
for i = 1:numOfLeaves
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',numOfLeaves+2);
for i = 0:numOfLeaves+1
fprintf( fileID,'Tf/LeafXPlus%i/Value ',i);
end
fprintf(fileID,'mm \n');
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',numOfLeaves);
for i = 1:numOfLeaves
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',numOfLeaves+2);
for i = 0:numOfLeaves+1
fprintf( fileID,'Tf/LeafXMinus%i/Value ',i);
end
fprintf(fileID,'mm \n');

%initilization of time features
fprintf(fileID,'d:Tf/TimelineStart = 0 ms\n');
fprintf(fileID,'d:Tf/TimelineEnd = %f ms\n',leafTimes*10);
fprintf(fileID,'i:Tf/NumberOfSequentialTimes = %i \n',leafTimes);

for i = 1:numOfLeaves
fprintf(fileID,'s:Tf/LeafXMinus%i/Function = "Step"\n',i);
fprintf(fileID,'dv:Tf/LeafXMinus%i/Times = %i ', i,leafTimes);
fprintf(fileID,num2str([1:leafTimes]*10,' % 2d'));
fprintf(fileID,'%i ', [1:leafTimes]*10);
fprintf(fileID,' ms\n');
fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,leafTimes);
fprintf(fileID,num2str(leftLeafPos(i,:),' % 2d'));
fprintf(fileID,'%f ', leftLeafPos(i,:));
fprintf(fileID,' mm\n\n');

fprintf(fileID,'s:Tf/LeafXPlus%i/Function = "Step"\n',i);
fprintf(fileID,'dv:Tf/LeafXPlus%i/Times = %i ',i,leafTimes);
fprintf(fileID,num2str([1:leafTimes]*10,' % 2d'));
fprintf(fileID,'%i ',[1:leafTimes]*10);
fprintf(fileID,' ms\n');
fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes);
fprintf(fileID,num2str(rightLeafPos(i,:),' % 2d'));
fprintf(fileID,'%f ', rightLeafPos(i,:));
fprintf(fileID,' mm\n\n');
end
%Add aditional Leaf at the top and bottom to catch
%scattering
for i = [0,numOfLeaves+1]
fprintf(fileID,'s:Tf/LeafXMinus%i/Function = "Step"\n',i);
fprintf(fileID,'dv:Tf/LeafXMinus%i/Times = %i ', i,leafTimes);
fprintf(fileID,'%i ',[1:leafTimes]*10);
fprintf(fileID,' ms\n');
fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,leafTimes);
fprintf(fileID,'%f ', zeros(size([1:leafTimes])));
fprintf(fileID,' mm\n\n');

fprintf(fileID,'s:Tf/LeafXPlus%i/Function = "Step"\n',i);
fprintf(fileID,'dv:Tf/LeafXPlus%i/Times = %i ',i,leafTimes);
fprintf(fileID,'%i ',[1:leafTimes]*10);
fprintf(fileID,' ms\n');
fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes);
fprintf(fileID,'%f ', zeros(size([1:leafTimes])));
fprintf(fileID,' mm\n\n');
end

fprintf(fileID, 's:Tf/Phasespace/NumberOfHistoriesInRun/Function = "Step" \n');
fprintf(fileID, 'dv:Tf/Phasespace/NumberOfHistoriesInRun/Times = %i ', leafTimes);
fprintf(fileID,'%i ',[1:leafTimes]*10);
fprintf(fileID,' ms\n');
fprintf(fileID, 'iv:Tf/Phasespace/NumberOfHistoriesInRun/Values = %i ', leafTimes);
fprintf(fileID,'%i ',[dataTOPAS(:).current]);
fprintf(fileID,' \n');

end

% Translate patient according to beam isocenter
Expand Down
Loading