diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m new file mode 100644 index 000000000..7b09a466a --- /dev/null +++ b/examples/matRad_example15_photonMC_MLC.m @@ -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); diff --git a/matRad_calcDoseInit.m b/matRad_calcDoseInit.m index e0d41e560..91bb30277 100644 --- a/matRad_calcDoseInit.m +++ b/matRad_calcDoseInit.m @@ -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; diff --git a/matRad_calcPhotonDoseMC.m b/matRad_calcPhotonDoseMC.m index 2c6342060..9f98802fc 100644 --- a/matRad_calcPhotonDoseMC.m +++ b/matRad_calcPhotonDoseMC.m @@ -37,7 +37,7 @@ matRad_cfg = MatRad_Config.instance(); -if nargin < 6 +if nargin < 5 calcDoseDirect = false; end diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m index ed478a595..d3996ce0f 100644 --- a/matRad_calcPhotonDoseMCtopas.m +++ b/matRad_calcPhotonDoseMCtopas.m @@ -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 diff --git a/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in b/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in index 16ee09185..9077118d8 100644 --- a/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in +++ b/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in @@ -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" diff --git a/topas/beamSetup/TOPAS_beamSetup_phasespace.txt.in b/topas/beamSetup/TOPAS_beamSetup_phasespace.txt.in index fd0da7aa0..62143c902 100644 --- a/topas/beamSetup/TOPAS_beamSetup_phasespace.txt.in +++ b/topas/beamSetup/TOPAS_beamSetup_phasespace.txt.in @@ -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" diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 250c21590..23726e3c7 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -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 @@ -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 @@ -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 @@ -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 @@ -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