Update to Photon calculation with TOPAS (new engine concept)

doseCalc/+DoseEngines/matRad_MonteCarloEngineAbstract.m

@@ -61,7 +61,7 @@ function setDefaults(this)
 
             % copy bixel weight vector into stf struct
             if nargin == 5
-                if sum([stf.totalNumOfBixels]) ~= numel(w)
+                if sum([stf.totalNumOfBixels]) ~= numel(w) && ~isfield([stf.ray],'shapes')
                     matRad_cfg.dispError('weighting does not match steering information')
                 end
                 counter = 0;
@@ -97,7 +97,7 @@ function setDefaults(this)
             % calculate cubes; use uniform weights here, weighting with actual fluence
             % already performed in dij construction
             resultGUI    = matRad_calcCubes(sum(w),dij);
-
+            
             % remember original fluence weights
             resultGUI.w  = w;
         end

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_NominalScenario(ct);
+% dose calculation settings
+pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
+
+%% Generate Beam Geometry STF
+stf = matRad_generateStf(ct,cst,pln);
+
+%% Dose Calculation
+
+dij = matRad_calcDoseInfluence(ct,cst,stf,pln);
+
+%% 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
+%resultGUI_MC = matRad_calcDoseInfluence(ct,cst,stf,pln);
+pln.propDoseCalc.engine = 'TOPAS';
+pln.propDoseCalc.beamProfile = 'phasespace';
+pln.propDoseCalc.externalCalculation =true;
+resultGUI_MC = matRad_calcDoseDirect(ct,stf,pln,cst,resultGUI.w);
+
+%% readout
+%foldername = 'FolderName';
+%pln = matRad_cfg.getDefaultClass(pln,'propMC');
+%resultGUI_MC = pln.propMC.readExternal(foldername);

matRad_aperture2collimation.m

@@ -122,6 +122,7 @@
 
     ray.shape = shapeTotalF;
     ray.weight = ones(1,nShapes);
+    ray.collimation = pln.propStf.collimation;
     stfTmp.ray = ray;
 
     stf(iBeam) = stfTmp;

matRad_calcCubes.m

@@ -165,7 +165,7 @@
 resultGUI = orderfields(resultGUI);
 
 % interpolation if dose grid does not match ct grid
-if any(dij.ctGrid.dimensions~=dij.doseGrid.dimensions)
+if isfield(dij,'ctGrid') && any(dij.ctGrid.dimensions~=dij.doseGrid.dimensions)
     myFields = fieldnames(resultGUI);
     for i = 1:numel(myFields)
         if numel(resultGUI.(myFields{i})) == dij.doseGrid.numOfVoxels

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"

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"