From f9b6294f6f154bc365962ed3b26d16379e7fcfda Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:35:52 +0200 Subject: [PATCH 01/14] Create matRad_example15_photonMC_MLC.m --- examples/matRad_example15_photonMC_MLC.m | 97 ++++++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 examples/matRad_example15_photonMC_MLC.m diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m new file mode 100644 index 000000000..5127a17b5 --- /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]; +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; +pln.propMC.infilenames.phaseSpaceSourcePhotons = 'SIEMENS_PRIMUS_6mv_15x15'; +%% 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 +stf.ray.energy = [6,6,6]; +stf.ray.weight = [stf.ray.shapes.weight]; +resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1); + +%% readout +foldername = 'E:\Code\matRad\topas\MCrun\photons_Generic_01-09-23'; +pln = matRad_cfg.getDefaultClass(pln,'propMC'); +resultGUI_MC = pln.propMC.readExternal(foldername); From dd82e795a19dc737828dbe77cf3eb5f63daf9835 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:36:09 +0200 Subject: [PATCH 02/14] small bug fix --- matRad_calcPhotonDoseMC.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 7c21e988e9ec031048cfd971522424dadec9fce2 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:36:41 +0200 Subject: [PATCH 03/14] Update matRad_calcPhotonDoseMCtopas.m small bug fix, and read in of MLC weights and energy for TOPAS simulation --- matRad_calcPhotonDoseMCtopas.m | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m index ed478a595..04bb53abb 100644 --- a/matRad_calcPhotonDoseMCtopas.m +++ b/matRad_calcPhotonDoseMCtopas.m @@ -79,21 +79,28 @@ load([pln.radiationMode,'_',pln.machine]); +stf.SCD = machine.meta.SCD; + %Collect weights 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; + w(counter:counter+rayBix-1) = stf(i).ray(j).weight; + counter = counter + rayBix; end end else w = ones(sum([stf(:).totalNumOfBixels]),1); end +if isfield(stf(1).ray, 'shapes') %weight is in stf.ray.shapes + stf.ray.weight = [stf.ray.shapes.weight]; %how for multiple fields ? +end +stf.ray.energy = stf.ray.energy*ones(size(w)); %how for multiple fields ? + currDir = cd; for shiftScen = 1:pln.multScen.totNumShiftScen From cdbcca58b9ec73129707e293108624a42581e6fe Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:41:47 +0200 Subject: [PATCH 04/14] MLC updates Changed Positioning of MLC so that is SCD away, Made MLC bigger to prevent photons flying past it Initilise Time Features. and us to weight different MLC shapes Cchanged Material of MLC to thungsten To Do: Check out multiple Photon Fields, Improve histories per run for the different shapes - do it like it is done for the particle bixels --- topas/beamSetup/TOPAS_beamSetup_mlc.txt.in | 5 +- .../TOPAS_beamSetup_phasespace.txt.in | 2 +- topas/matRad_TopasConfig.m | 68 +++++++++++++++---- 3 files changed, 57 insertions(+), 18 deletions(-) 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..cb6d395ea 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' ,'SIEMENS_PRIMUS_6.0_0.10_15.0x15.0' ); + end @@ -1569,20 +1572,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 @@ -1653,23 +1656,30 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) leftLeafPos = [stf(beamIx).ray.shapes(:).leftLeafPos]; rightLeafPos = [stf(beamIx).ray.shapes(:).rightLeafPos]; % Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator + MLCshift = 0.5*15 + stf.SCD; %MLC thickness + SCD in cm + fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n',MLCshift); 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/Length = %f cm\n',15); + fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves+2); + fprintf(fileID,num2str([200, pln.propStf.collimation.leafWidth*ones(1,numOfLeaves) , 200],' % 2d')); 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); @@ -1687,6 +1697,34 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) fprintf(fileID,num2str(rightLeafPos(i,:),' % 2d')); fprintf(fileID,' mm\n\n'); end + %Add aditional Leaf at the top and bottom to catch + %scattering + for i = [0,10] + 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,' ms\n'); + fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,leafTimes); + fprintf(fileID,num2str([0,0,0],' % 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,leafTimes); + fprintf(fileID,num2str([1:leafTimes]*10,' % 2d')); + fprintf(fileID,' ms\n'); + fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes); + fprintf(fileID,num2str([0,0,0],' % 2d')); + 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,num2str([1:leafTimes]*10,' % 2d')); + fprintf(fileID,' ms\n'); + fprintf(fileID, 'iv:Tf/Phasespace/NumberOfHistoriesInRun/Values = %i ', leafTimes); + fprintf(fileID,num2str(round(historyCount(beamIx) * w'./sum(w)),' % 2d')); + fprintf(fileID,' \n'); + end % Translate patient according to beam isocenter From 18f33fd0be800fe30b803334a17ea4c2e88e0e6b Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:52:48 +0200 Subject: [PATCH 05/14] Update matRad_example15_photonMC_MLC.m --- examples/matRad_example15_photonMC_MLC.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m index 5127a17b5..1ae649b60 100644 --- a/examples/matRad_example15_photonMC_MLC.m +++ b/examples/matRad_example15_photonMC_MLC.m @@ -87,8 +87,9 @@ imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet) %% Dose Calculation -stf.ray.energy = [6,6,6]; -stf.ray.weight = [stf.ray.shapes.weight]; +%stf.ray.energy = [6,6,6]; +%stf.ray.weight = [stf.ray.shapes.weight]; +pln.propMC.numHistories = 1e8; resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1); %% readout From 2ece06e7e417c5e1c1fae6f467cebbb6cb89f104 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 1 Sep 2023 16:58:10 +0200 Subject: [PATCH 06/14] Update matRad_TopasConfig.m --- topas/matRad_TopasConfig.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index cb6d395ea..2cd7446b6 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -1656,7 +1656,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) leftLeafPos = [stf(beamIx).ray.shapes(:).leftLeafPos]; rightLeafPos = [stf(beamIx).ray.shapes(:).rightLeafPos]; % Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator - MLCshift = 0.5*15 + stf.SCD; %MLC thickness + SCD in cm + MLCshift = 0.5*15 + stf.SCD/10; %MLC thickness + SCD in cm fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n',MLCshift); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/MaximumLeafOpen = %f cm\n',15); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Thickness = %f cm\n',15); From 2184a10e7adbca02c3e1e5a9fbdfa00bca5d5485 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Mon, 4 Sep 2023 09:14:44 +0200 Subject: [PATCH 07/14] small bugfix --- examples/matRad_example15_photonMC_MLC.m | 8 ++++---- matRad_calcPhotonDoseMCtopas.m | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m index 1ae649b60..6147b0a30 100644 --- a/examples/matRad_example15_photonMC_MLC.m +++ b/examples/matRad_example15_photonMC_MLC.m @@ -89,10 +89,10 @@ %% Dose Calculation %stf.ray.energy = [6,6,6]; %stf.ray.weight = [stf.ray.shapes.weight]; -pln.propMC.numHistories = 1e8; +pln.propMC.numHistories = 1e7; resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1); %% readout -foldername = 'E:\Code\matRad\topas\MCrun\photons_Generic_01-09-23'; -pln = matRad_cfg.getDefaultClass(pln,'propMC'); -resultGUI_MC = pln.propMC.readExternal(foldername); +%foldername = 'E:\Code\matRad\topas\MCrun\photons_Generic_01-09-23'; +%pln = matRad_cfg.getDefaultClass(pln,'propMC'); +%resultGUI_MC = pln.propMC.readExternal(foldername); diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m index 04bb53abb..9cc0e986a 100644 --- a/matRad_calcPhotonDoseMCtopas.m +++ b/matRad_calcPhotonDoseMCtopas.m @@ -97,9 +97,9 @@ end if isfield(stf(1).ray, 'shapes') %weight is in stf.ray.shapes - stf.ray.weight = [stf.ray.shapes.weight]; %how for multiple fields ? + w = [stf.ray.shapes.weight]; %how for multiple fields ? end -stf.ray.energy = stf.ray.energy*ones(size(w)); %how for multiple fields ? +stf.ray.energy = stf.ray.energy.*ones(size(w)); %how for multiple fields ? currDir = cd; From 98659424efa7b08c2ccf61cbd63753b463a67074 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Mon, 4 Sep 2023 09:24:26 +0200 Subject: [PATCH 08/14] Update matRad_TopasConfig.m --- topas/matRad_TopasConfig.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 2cd7446b6..899247f6c 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -1722,7 +1722,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) fprintf(fileID,num2str([1:leafTimes]*10,' % 2d')); fprintf(fileID,' ms\n'); fprintf(fileID, 'iv:Tf/Phasespace/NumberOfHistoriesInRun/Values = %i ', leafTimes); - fprintf(fileID,num2str(round(historyCount(beamIx) * w'./sum(w)),' % 2d')); + fprintf(fileID,num2str(round(historyCount(beamIx) * w./sum(w)),' % 2d')); fprintf(fileID,' \n'); end From 146ddd39191bce81c4b08ea662f8ea2c73a312ac Mon Sep 17 00:00:00 2001 From: JenHardt Date: Wed, 6 Sep 2023 15:57:05 +0200 Subject: [PATCH 09/14] Update matRad_TopasConfig.m --- topas/matRad_TopasConfig.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 899247f6c..d5ad00486 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -1656,7 +1656,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) leftLeafPos = [stf(beamIx).ray.shapes(:).leftLeafPos]; rightLeafPos = [stf(beamIx).ray.shapes(:).rightLeafPos]; % Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator - MLCshift = 0.5*15 + stf.SCD/10; %MLC thickness + SCD in cm + MLCshift = stf.SCD/10 - 0.5*15; %MLC thickness + SCD in cm SCD defined at the edge facing pacient? fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n',MLCshift); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/MaximumLeafOpen = %f cm\n',15); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Thickness = %f cm\n',15); From da650b1b1e66c81e7c3c300f55b94b29e97d2f62 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 27 Oct 2023 16:48:22 +0200 Subject: [PATCH 10/14] update extension to multiple field, corrected leaf position to be at collimator not as isocenter --- examples/matRad_example15_photonMC_MLC.m | 12 +++---- matRad_calcPhotonDoseMCtopas.m | 18 +++++++---- topas/matRad_TopasConfig.m | 40 +++++++++++------------- 3 files changed, 37 insertions(+), 33 deletions(-) diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m index 6147b0a30..589991f12 100644 --- a/examples/matRad_example15_photonMC_MLC.m +++ b/examples/matRad_example15_photonMC_MLC.m @@ -35,8 +35,10 @@ pln.radiationMode = 'photons'; pln.machine = 'Generic'; pln.numOfFractions = 30; -pln.propStf.gantryAngles = [0]; -pln.propStf.couchAngles = [0]; +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); @@ -87,12 +89,10 @@ imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet) %% Dose Calculation -%stf.ray.energy = [6,6,6]; -%stf.ray.weight = [stf.ray.shapes.weight]; -pln.propMC.numHistories = 1e7; +pln.propMC.numHistories = 1e9; resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1); %% readout -%foldername = 'E:\Code\matRad\topas\MCrun\photons_Generic_01-09-23'; +%foldername = 'FolderName'; %pln = matRad_cfg.getDefaultClass(pln,'propMC'); %resultGUI_MC = pln.propMC.readExternal(foldername); diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m index 9cc0e986a..d3996ce0f 100644 --- a/matRad_calcPhotonDoseMCtopas.m +++ b/matRad_calcPhotonDoseMCtopas.m @@ -79,16 +79,23 @@ load([pln.radiationMode,'_',pln.machine]); -stf.SCD = machine.meta.SCD; +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); counter = 1; for i = 1:length(stf) for j = 1:stf(i).numOfRays rayBix = stf(i).numOfBixelsPerRay(j); - w(counter:counter+rayBix-1) = stf(i).ray(j).weight; + 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 @@ -96,10 +103,9 @@ w = ones(sum([stf(:).totalNumOfBixels]),1); end -if isfield(stf(1).ray, 'shapes') %weight is in stf.ray.shapes - w = [stf.ray.shapes.weight]; %how for multiple fields ? +for i = 1:size(stf,2) + stf(i).ray.energy = stf(i).ray.energy.*ones(size(w)); end -stf.ray.energy = stf.ray.energy.*ones(size(w)); %how for multiple fields ? currDir = cd; diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index d5ad00486..9b414aec4 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -1383,8 +1383,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', (SAD-stf(beamIx).SCD) +50); %Phasespace behind MLC else fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', nozzleToAxisDistance); end @@ -1647,22 +1647,20 @@ 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 - MLCshift = stf.SCD/10 - 0.5*15; %MLC thickness + SCD in cm SCD defined at the edge facing pacient? - fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n',MLCshift); + fprintf(fileID,'d:Sim/Ge/MultiLeafCollimatorA/TransZ = %f cm\n', 0.5*5); 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/Thickness = %f cm\n',5); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Length = %f cm\n',15); fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves+2); - fprintf(fileID,num2str([200, pln.propStf.collimation.leafWidth*ones(1,numOfLeaves) , 200],' % 2d')); + 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+2); for i = 0:numOfLeaves+1 @@ -1683,46 +1681,46 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) 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,10] + 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,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([0,0,0],' % 2d')); + fprintf(fileID,'%f ', [0,0,0]); 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([0,0,0],' % 2d')); + fprintf(fileID,'%f ',[0,0,0]); 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,num2str([1:leafTimes]*10,' % 2d')); + fprintf(fileID,'%i ',[1:leafTimes]*10); fprintf(fileID,' ms\n'); fprintf(fileID, 'iv:Tf/Phasespace/NumberOfHistoriesInRun/Values = %i ', leafTimes); - fprintf(fileID,num2str(round(historyCount(beamIx) * w./sum(w)),' % 2d')); + fprintf(fileID,'%i ',[dataTOPAS(:).current]); fprintf(fileID,' \n'); end From a6e679a36f71328a195cb48a07f5ec2e8b326fb8 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Mon, 30 Oct 2023 12:56:22 +0100 Subject: [PATCH 11/14] Update matRad_TopasConfig.m --- topas/matRad_TopasConfig.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 9b414aec4..6c5ebd4d9 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -1703,7 +1703,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) fprintf(fileID,'%i ',[1:leafTimes]*10); fprintf(fileID,' ms\n'); fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,leafTimes); - fprintf(fileID,'%f ', [0,0,0]); + fprintf(fileID,'%f ', zeros(size([1:leafTimes]))); fprintf(fileID,' mm\n\n'); fprintf(fileID,'s:Tf/LeafXPlus%i/Function = "Step"\n',i); @@ -1711,7 +1711,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) fprintf(fileID,'%i ',[1:leafTimes]*10); fprintf(fileID,' ms\n'); fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes); - fprintf(fileID,'%f ',[0,0,0]); + fprintf(fileID,'%f ', zeros(size([1:leafTimes]))); fprintf(fileID,' mm\n\n'); end From cdf7cd3eeff7374b3a505e97f0d11b61a1645883 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 26 Jan 2024 17:09:28 +0100 Subject: [PATCH 12/14] use other phase space file and adjustments to that --- examples/matRad_example15_photonMC_MLC.m | 3 +-- topas/matRad_TopasConfig.m | 19 +++++++++++-------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m index 589991f12..787794101 100644 --- a/examples/matRad_example15_photonMC_MLC.m +++ b/examples/matRad_example15_photonMC_MLC.m @@ -62,8 +62,7 @@ %% pln.propMC.engine = 'TOPAS'; pln.propMC.beamProfile = 'phasespace'; -pln.propMC.externalCalculation = true; -pln.propMC.infilenames.phaseSpaceSourcePhotons = 'SIEMENS_PRIMUS_6mv_15x15'; +pln.propMC.externalCalculation =true; %% Generate Beam Geometry STF stf = matRad_generateStf(ct,cst,pln); diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 6c5ebd4d9..4f62849ff 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -136,7 +136,7 @@ 'Scorer_RBE_WED','TOPAS_scorer_doseRBE_Wedenberg.txt.in',... 'Scorer_RBE_MCN','TOPAS_scorer_doseRBE_McNamara.txt.in', ... ... %PhaseSpace Source - 'phaseSpaceSourcePhotons' ,'SIEMENS_PRIMUS_6.0_0.10_15.0x15.0' ); + 'phaseSpaceSourcePhotons' ,'VarianClinaciX_6MV_20x20_aboveMLC_w2' ); end @@ -693,6 +693,9 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) processedQuantities = {'','_std','_batchStd'}; topasCubesTallies = unique(erase(topasCubesTallies,processedQuantities(2:end))); + normFac = sum(w); + + % Loop through 4D scenarios for ctScen = 1:dij.numOfScenarios @@ -708,7 +711,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]) ... && iscell(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); end end end @@ -721,7 +724,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) for p = 1:length(processedQuantities) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); end end end @@ -742,7 +745,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]) ... && iscell(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); end end % Handle RBE-related quantities (not multiplied by sum(w)!) @@ -779,7 +782,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) for p = 1:length(processedQuantities) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); end end % Handle RBE-related quantities (not multiplied by sum(w)!) @@ -1384,7 +1387,7 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) % NozzleAxialDistance if isPhoton - fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', (SAD-stf(beamIx).SCD) +50); %Phasespace behind MLC + 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 @@ -1655,9 +1658,9 @@ function writeStfFields(obj,ct,stf,pln,w,baseData) 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', 0.5*5); + 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',5); + 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]); From cc42adfeb59941430a314d63c367e83731890aab Mon Sep 17 00:00:00 2001 From: JenHardt Date: Thu, 1 Feb 2024 17:37:32 +0100 Subject: [PATCH 13/14] small changes --- examples/matRad_example15_photonMC_MLC.m | 2 +- matRad_calcDoseInit.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/matRad_example15_photonMC_MLC.m b/examples/matRad_example15_photonMC_MLC.m index 787794101..7b09a466a 100644 --- a/examples/matRad_example15_photonMC_MLC.m +++ b/examples/matRad_example15_photonMC_MLC.m @@ -88,7 +88,7 @@ imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet) %% Dose Calculation -pln.propMC.numHistories = 1e9; +pln.propMC.numHistories = 1e8; resultGUI_MC = matRad_calcPhotonDoseMC(ct,stf,pln,cst,1); %% readout 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; From c2c293d2ea98f3333a95b88461f9a641adbd0d5a Mon Sep 17 00:00:00 2001 From: JenHardt Date: Wed, 28 Feb 2024 09:29:30 +0100 Subject: [PATCH 14/14] Update matRad_TopasConfig.m --- topas/matRad_TopasConfig.m | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/topas/matRad_TopasConfig.m b/topas/matRad_TopasConfig.m index 4f62849ff..23726e3c7 100644 --- a/topas/matRad_TopasConfig.m +++ b/topas/matRad_TopasConfig.m @@ -692,8 +692,6 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Allocate possible scored quantities processedQuantities = {'','_std','_batchStd'}; topasCubesTallies = unique(erase(topasCubesTallies,processedQuantities(2:end))); - - normFac = sum(w); % Loop through 4D scenarios @@ -711,7 +709,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]) ... && iscell(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); end end end @@ -724,7 +722,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) for p = 1:length(processedQuantities) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); end end end @@ -745,7 +743,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]) ... && iscell(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen,1}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} '_ray' num2str(dij.rayNum(d)) '_bixel' num2str(dij.bixelNum(d)) processedQuantities{p} '_beam' num2str(dij.beamNum(d))]){ctScen},[],1); end end % Handle RBE-related quantities (not multiplied by sum(w)!) @@ -782,7 +780,7 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) for p = 1:length(processedQuantities) % Check if current quantity is available and write to dij if isfield(topasCubes,[topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)]) && iscell(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam' num2str(d)])) - dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = normFac*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); + dij.([topasCubesTallies{j} processedQuantities{p}]){ctScen}(:,d) = sum(w)*reshape(topasCubes.([topasCubesTallies{j} processedQuantities{p} '_beam',num2str(d)]){ctScen},[],1); end end % Handle RBE-related quantities (not multiplied by sum(w)!)