Skip to content

Commit

Permalink
PhotonTOPAS Update
Browse files Browse the repository at this point in the history 
- export leafWidth from aperture2collimation
- small bug with additional "end" fix in ParticleDoseMCtopas
- added switch for uniform weights in case of ~calcDoseDirect
- changed materialConversion to 'RSP'
- renamed numOfLeaves and leafTimes parameters for easier understanding
- Translated MLC so that the Nozzle lies at beginning of MLC, not centered
  • Loading branch information
@HomolkaN
HomolkaN committed Sep 16, 2022
1 parent d43b80c commit 5ee5467
Show file tree
Hide file tree
Showing 6 changed files with 72 additions and 70 deletions.
2 changes: 1 addition & 1 deletion matRad_aperture2collimation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@
pln.propStf.bixelWidth = 'field';
pln.propStf.collimation.convResolution = 0.5; %[mm]
pln.propStf.collimation.fieldWidth = fieldWidth;

pln.propStf.collimation.leafWidth = leafWidth;

%
%[bixelFieldX,bixelFieldY] = ndgrid(-fieldWidth/2:bixelWidth:fieldWidth/2,-fieldWidth/2:leafWidth:fieldWidth/2);
Expand Down
83 changes: 41 additions & 42 deletions matRad_calcParticleDoseMCtopas.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,13 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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/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
% 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/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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expand All @@ -51,48 +51,47 @@
matRad_cfg.dispWarning('Variance scoring for TOPAS not yet supported.');
end

if isfield(pln,'propMC') && isfield(pln.propMC,'config')
if isfield(pln,'propMC') && isfield(pln.propMC,'config')
if isa(pln.propMC.config,'MatRad_TopasConfig')
matRad_cfg.dispInfo('Using given Topas Configuration in pln.propMC.config!\n');
topasConfig = pln.propMC.config;
else
else
%Create a default instance of the configuration
topasConfig = MatRad_TopasConfig();
end
%Overwrite parameters
%mc = metaclass(topasConfig); %get metaclass information to check if we can overwrite properties

if isstruct(pln.propMC.config)
props = fieldnames(pln.propMC.config);
for fIx = 1:numel(props)
fName = props{fIx};
if isprop(topasConfig,fName)
%We use a try catch block to catch errors when trying
%to overwrite protected/private properties instead of a
%metaclass approach
try
topasConfig.(fName) = pln.propMC.config.(fName);
catch
matRad_cfg.dispWarning('Property ''%s'' for MatRad_TopasConfig will be omitted due to protected/private access or invalid value.',fName);
end
else
matRad_cfg.dispWarning('Unkown property ''%s'' for MatRad_TopasConfig will be omitted.',fName);
topasConfig = MatRad_TopasConfig();
end
%Overwrite parameters
%mc = metaclass(topasConfig); %get metaclass information to check if we can overwrite properties

if isstruct(pln.propMC.config)
props = fieldnames(pln.propMC.config);
for fIx = 1:numel(props)
fName = props{fIx};
if isprop(topasConfig,fName)
%We use a try catch block to catch errors when trying
%to overwrite protected/private properties instead of a
%metaclass approach
try
topasConfig.(fName) = pln.propMC.config.(fName);
catch
matRad_cfg.dispWarning('Property ''%s'' for MatRad_TopasConfig will be omitted due to protected/private access or invalid value.',fName);
end
else
matRad_cfg.dispWarning('Unkown property ''%s'' for MatRad_TopasConfig will be omitted.',fName);
end
else
matRad_cfg.dispError('Invalid Configuration in pln.propMC.config');
end
else
matRad_cfg.dispError('Invalid Configuration in pln.propMC.config');
end
else
topasConfig = MatRad_TopasConfig();
end


if ~calcDoseDirect
matRad_cfg.dispError('matRad so far only supports direct dose calculation for TOPAS!\n');
end

if ~isfield(pln.propStf,'useRangeShifter')
if ~isfield(pln.propStf,'useRangeShifter')
pln.propStf.useRangeShifter = false;
end

Expand All @@ -107,9 +106,9 @@

for s = 1:ct.numOfCtScen
cubeHUresampled{s} = matRad_interp3(dij.ctGrid.x, dij.ctGrid.y', dij.ctGrid.z,ct.cubeHU{s}, ...
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear');
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear');
cubeResampled{s} = matRad_interp3(dij.ctGrid.x, dij.ctGrid.y', dij.ctGrid.z,ct.cube{s}, ...
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear');
dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear');
end

%Allocate temporary resampled CT
Expand Down Expand Up @@ -150,7 +149,7 @@
% manipulate isocenter
for k = 1:length(stf)
stf(k).isoCenter = stf(k).isoCenter + pln.multScen.isoShift(shiftScen,:);
end
end

for ctScen = 1:pln.multScen.numOfCtScen
for rangeShiftScen = 1:pln.multScen.totNumRangeScen
Expand All @@ -160,14 +159,14 @@
%TOPAS here)
ctR.cubeHU = cubeHUresampled(ctScen);
ctR.cube = cubeResampled(ctScen);
topasConfig.writeAllFiles(ctR,pln,stf,topasBaseData,w);

topasConfig.writeAllFiles(ctR,pln,stf,topasBaseData,w);

% Run simulation for current scenario
cd(topasConfig.workingDir);

for beamIx = 1:numel(stf)
for runIx = 1:topasConfig.numOfRuns
for runIx = 1:topasConfig.numOfRuns
fname = sprintf('%s_field%d_run%d',topasConfig.label,beamIx,runIx);
topasCall = sprintf('%s %s.txt > %s.out 2> %s.err',topasConfig.topasExecCommand,fname,fname,fname);
if topasConfig.parallelRuns
Expand All @@ -183,7 +182,7 @@
matRad_cfg.dispError('TOPAS simulation exited with error code %d\n',status);
end
end

if topasConfig.parallelRuns
runsFinished = false;
pause('on');
Expand All @@ -193,15 +192,15 @@
runsFinished = all(fin);
end
end

end

cd(currDir);

%% Simulation finished - read out volume scorers from topas simulation

topasCubes = matRad_readTopasData(topasConfig.workingDir);

fnames = fieldnames(topasCubes);
dij.MC_tallies = fnames;
for f = 1:numel(fnames)
Expand Down
18 changes: 11 additions & 7 deletions matRad_calcPhotonDoseMCtopas.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,14 +126,18 @@
topasConfig.numOfRuns = 1; %matRad_cfg.propMC.topas_defaultNumBatches;
topasConfig.beamProfile = "uniform"; %"virtualGaussian";
%Collect weights
w = zeros(sum([stf(:).totalNumOfBixels]),1);
ct = 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 calcDoseDirect
w = zeros(sum([stf(:).totalNumOfBixels]),1);
ct = 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;
end
end
else
w = ones(sum([stf(:).totalNumOfBixels]),1);
end

currDir = cd;
Expand Down
5 changes: 2 additions & 3 deletions topas/MatRad_TopasConfig.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,8 +110,7 @@
%Let's set some default commands taken from topas installation
%instructions for mac & debain/ubuntu
if ispc %We assume topas is installed in wsl (since no windows version)
%obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; ~/topas/bin/topas';
obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; export LD_LIBRARY_PATH=~/topas/libexternal/:$LD_LIBRARY_PATH; ~/topas/topas'; %<-- for Pia's laptop uncomment
obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; ~/topas/bin/topas';
elseif ismac
obj.topasExecCommand = 'export TOPAS_G4_DATA_DIR=/Applications/G4Data; export QT_QPA_PLATFORM_PLUGIN_PATH=/Applications/topas/Frameworks; /Applications/topas/bin/topas';
elseif isunix
Expand Down Expand Up @@ -673,7 +672,7 @@ function writePatient(obj,ct,pln)
end

%cubeExport = obj.materialConversion; %'RSP'; %'HUSchneiderToWater'; 'CustomWaterRSP'
cubeExport = 'CustomWaterRSP';
cubeExport = 'RSP';
rspCubeMethod = 2;
checkMaterial = false;

Expand Down
2 changes: 1 addition & 1 deletion topas/TOPAS_beamSetup_mlc.txt.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ s:Ge/MultiLeafCollimatorA/Parent = "Nozzle"
s:Ge/MultiLeafCollimatorA/Material = "Aluminum"
d:Ge/MultiLeafCollimatorA/TransX = 0.0 cm
d:Ge/MultiLeafCollimatorA/TransY = 0.0 cm
d:Ge/MultiLeafCollimatorA/TransZ = 0.0 cm
d:Ge/MultiLeafCollimatorA/TransZ = 7.5 cm
d:Ge/MultiLeafCollimatorA/RotX = 0.0 deg
d:Ge/MultiLeafCollimatorA/RotY = 0.0 deg
d:Ge/MultiLeafCollimatorA/RotZ = 90.0 deg
32 changes: 16 additions & 16 deletions topas/matRad_TopasConfigPhotons.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,42 +430,42 @@ function writeStfFields(obj,ct,stf,pln,baseData,w)

if exist('TOPAS_mlcSetup')
fprintf(fileID,'%s\n',TOPAS_mlcSetup);
%Erstmal nur für einen ray!
[l,t]=size([stf(1).ray.shapes(:).leftLeafPos]); %there are l leaves and t times/shapes
%Erstmal nur f??r einen ray!
[numOfLeaves,leafTimes]=size([stf(1).ray.shapes(:).leftLeafPos]); %there are #numOfLeaves leaves and #leafTimes times/shapes
leftLeafPos = [stf(1).ray.shapes(:).leftLeafPos];
rightLeafPos = [stf(1).ray.shapes(:).rightLeafPos];
%Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator
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 ', l)
fprintf(fileID,num2str(ones(1,l),' % 2d'));
fprintf(fileID,' cm \n');
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',l);
for i=1:l
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves);
fprintf(fileID,num2str(pln.propStf.collimation.leafWidth*ones(1,numOfLeaves),' % 2d'));
fprintf(fileID,' mm \n');
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',numOfLeaves);
for i=1:numOfLeaves
fprintf( fileID,'Tf/LeafXPlus%i/Value ',i);
end
fprintf(fileID,'mm \n');
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',l);
for i=1:l
fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',numOfLeaves);
for i=1:numOfLeaves
fprintf( fileID,'Tf/LeafXMinus%i/Value ',i);
end
fprintf(fileID,'mm \n');

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

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

Expand Down

0 comments on commit 5ee5467

Please sign in to comment.