Description of matRad_calcParticleDoseMC

0001 function dij = matRad_calcParticleDoseMC(ct,stf,pln,cst,nCasePerBixel,calcDoseDirect)
0002 % matRad MCsqaure monte carlo photon dose calculation wrapper
0003 %
0004 % call
0005 %   dij = matRad_calcParticleDoseMc(ct,stf,pln,cst,calcDoseDirect)
0006 %
0007 % input
0008 %   ct:              matRad ct struct
0009 %   stf:             atRad steering information struct
0010 %   pln:            matRad plan meta information struct
0011 %   cst:            matRad cst struct
0012 %   nCasePerBixel:  number of histories per beamlet
0013 %   calcDoseDirect: binary switch to enable forward dose
0014 %                   calcualtion
0015 % output
0016 %   dij:            matRad dij struct
0017 %
0018 % References
0019 %
0020 %   https://aapm.onlinelibrary.wiley.com/doi/abs/10.1118/1.4943377
0021 %   http://www.openmcsquare.org/
0022 %
0023 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 %
0025 % Copyright 2019 the matRad development team.
0026 %
0027 % This file is part of the matRad project. It is subject to the license
0028 % terms in the LICENSE file found in the top-level directory of this
0029 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0030 % of the matRad project, including this file, may be copied, modified,
0031 % propagated, or distributed except according to the terms contained in the
0032 % LICENSE file.
0033 %
0034 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0035 
0036 
0037 matRad_cfg = MatRad_Config.instance();
0038 
0039 % check if valid machine
0040 if ~strcmp(pln.radiationMode,'protons') || ~strcmp(pln.machine,'generic_MCsquare')
0041     matRad_cfg.dispError('Wrong radiation modality and/or machine. For now MCsquare requires machine generic_MCsquare!');    
0042 end
0043 
0044 
0045 if nargin < 5
0046     % set number of particles simulated per pencil beam
0047     nCasePerBixel = matRad_cfg.propMC.MCsquare_defaultHistories;
0048     matRad_cfg.dispInfo('Using default number of Histories per Bixel: %d\n',nCasePerBixel);
0049 end
0050 
0051 if nargin < 6
0052     calcDoseDirect = false;
0053 end
0054 
0055 if isfield(pln,'propMC') && isfield(pln.propMC,'outputVariance')
0056     matRad_cfg.dispWarning('Variance scoring for MCsquare not yet supported.');
0057 end
0058 
0059 if ~strcmp(pln.radiationMode,'protons')
0060     errordlg('MCsquare is only supported for protons');
0061 end
0062 
0063 env = matRad_getEnvironment();
0064 
0065 %% check if binaries are available
0066 %Executables for simulation
0067 if ispc
0068     if exist('MCSquare_windows.exe','file') ~= 2
0069         matRad_cfg.dispError('Could not find MCsquare binary.\n');
0070     else
0071         mcSquareBinary = 'MCSquare_windows.exe';
0072     end
0073 elseif ismac
0074     if exist('MCsquare_mac','file') ~= 2
0075         matRad_cfg.dispError('Could not find MCsquare binary.\n');
0076     else
0077         mcSquareBinary = './MCsquare_mac';
0078     end
0079     %error('MCsquare binaries not available for mac OS.\n');
0080 elseif isunix
0081     if exist('MCsquare_linux','file') ~= 2
0082         matRad_cfg.dispError('Could not find MCsquare binary.\n');
0083     else
0084         mcSquareBinary = './MCsquare_linux';
0085     end
0086 end
0087 
0088 %Mex interface for import of sparse matrix
0089 if ~calcDoseDirect && ~matRad_checkMexFileExists('matRad_sparseBeamletsReaderMCsquare')
0090     matRad_cfg.dispWarning('Compiled sparse reader interface not found. Trying to compile it on the fly!');      
0091     try
0092         matRad_compileMCsquareSparseReader();        
0093     catch MException
0094         matRad_cfg.dispError('Could not find/generate mex interface for reading the sparse matrix. \nCause of error:\n%s\n Please compile it yourself.',MException.message);
0095     end
0096 end
0097 
0098 
0099 % set and change to MCsquare binary folder
0100 currFolder = pwd;
0101 fullfilename = mfilename('fullpath');
0102 MCsquareFolder = [fullfilename(1:find(fullfilename==filesep,1,'last')) 'MCsquare' filesep 'bin'];
0103 
0104 % cd to MCsquare folder (necessary for binary)
0105 cd(MCsquareFolder);
0106 
0107 %Check Materials
0108 if ~exist([MCsquareFolder filesep 'Materials'],'dir') || ~exist(fullfile(MCsquareFolder,'Materials','list.dat'),'file')
0109     matRad_cfg.dispInfo('First call of MCsquare: unzipping Materials...');    
0110     unzip('Materials.zip');
0111     matRad_cfg.dispInfo('Done');
0112 end
0113 
0114 % Since MCsquare 1.1 only allows similar resolution in x&y, we do some
0115 % extra checks on that before calling calcDoseInit. First, we make sure a
0116 % dose grid resolution is set in the pln struct
0117 if ~isfield(pln,'propDoseCalc') ...
0118         || ~isfield(pln.propDoseCalc,'doseGrid') ...
0119         || ~isfield(pln.propDoseCalc.doseGrid,'resolution') ...
0120         || ~all(isfield(pln.propDoseCalc.doseGrid.resolution,{'x','y','z'}))
0121     
0122     %Take default values
0123     pln.propDoseCalc.doseGrid.resolution = matRad_cfg.propDoseCalc.defaultResolution;
0124 end
0125 
0126 % Now we check for different x/y
0127 if pln.propDoseCalc.doseGrid.resolution.x ~= pln.propDoseCalc.doseGrid.resolution.y
0128     pln.propDoseCalc.doseGrid.resolution.x = mean([pln.propDoseCalc.doseGrid.resolution.x pln.propDoseCalc.doseGrid.resolution.y]);
0129     pln.propDoseCalc.doseGrid.resolution.y = pln.propDoseCalc.doseGrid.resolution.x;
0130     matRad_cfg.dispWarning('Anisotropic resolution in axial plane for dose calculation with MCsquare not possible\nUsing average x = y = %g mm\n',pln.propDoseCalc.doseGrid.resolution.x);
0131 end
0132 
0133 %Now we can run calcDoseInit as usual
0134 matRad_calcDoseInit;
0135 
0136 %Issue a warning when we have more than 1 scenario
0137 if dij.numOfScenarios ~= 1
0138     matRad_cfg.dispWarning('MCsquare is only implemented for single scenario use at the moment. Will only use the first Scenario for Monte Carlo calculation!');
0139 end
0140 
0141 % prefill ordering of MCsquare bixels
0142 dij.MCsquareCalcOrder = NaN*ones(dij.totalNumOfBixels,1);
0143 
0144 % We need to adjust the offset used in matRad_calcDoseInit
0145 mcSquareAddIsoCenterOffset = [dij.doseGrid.resolution.x/2 dij.doseGrid.resolution.y/2 dij.doseGrid.resolution.z/2] ...
0146                 - [dij.ctGrid.resolution.x   dij.ctGrid.resolution.y   dij.ctGrid.resolution.z];
0147 mcSquareAddIsoCenterOffset = mcSquareAddIsoCenterOffset - offset;
0148 
0149 % for MCsquare we explicitly downsample the ct to the dose grid (might not
0150 % be necessary in future MCsquare versions with separated grids)
0151 for s = 1:dij.numOfScenarios
0152     HUcube{s} =  matRad_interp3(dij.ctGrid.x,  dij.ctGrid.y',  dij.ctGrid.z,ct.cubeHU{s}, ...
0153                                 dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear');
0154 end
0155 
0156 % Explicitly setting the number of threads for MCsquare, 0 is all available
0157 nbThreads = 0;
0158 
0159 % set relative dose cutoff for storage in dose influence matrix, we use the
0160 % default value for the lateral cutoff here
0161 relDoseCutoff = 1 - matRad_cfg.propDoseCalc.defaultLateralCutOff;
0162 % set absolute calibration factor
0163 % convert from eV/g/primary to Gy 1e6 primaries
0164 absCalibrationFactorMC2 = 1.602176e-19 * 1.0e+9;
0165 
0166 if isequal(pln.propOpt.bioOptimization,'const_RBExD')
0167             dij.RBE = 1.1;
0168             matRad_cfg.dispInfo('matRad: Using a constant RBE of %g\n',dij.RBE);
0169 end
0170 
0171 % MCsquare settings
0172 MCsquareConfigFile = 'MCsquareConfig.txt';
0173 
0174 MCsquareConfig = MatRad_MCsquareConfig;
0175 
0176 MCsquareConfig.BDL_Plan_File = 'currBixels.txt';
0177 MCsquareConfig.CT_File       = 'MC2patientCT.mhd';
0178 MCsquareConfig.Num_Threads   = nbThreads;
0179 MCsquareConfig.RNG_Seed      = 1234;
0180 MCsquareConfig.Num_Primaries = nCasePerBixel;
0181 
0182 % turn simulation of individual beamlets
0183 MCsquareConfig.Beamlet_Mode = ~calcDoseDirect;
0184 % turn of writing of full dose cube
0185 MCsquareConfig.Dose_MHD_Output = calcDoseDirect;
0186 % turn on sparse output
0187 MCsquareConfig.Dose_Sparse_Output = ~calcDoseDirect;
0188 % set threshold of sparse matrix generation
0189 MCsquareConfig.Dose_Sparse_Threshold = relDoseCutoff;
0190 
0191 % write patient data
0192 MCsquareBinCubeResolution = [dij.doseGrid.resolution.x ...
0193                              dij.doseGrid.resolution.y ...
0194                              dij.doseGrid.resolution.z];   
0195 matRad_writeMhd(HUcube{1},MCsquareBinCubeResolution,MCsquareConfig.CT_File);
0196 
0197 
0198 
0199 counter = 0;             
0200 for i = 1:length(stf)
0201     stfMCsquare(i).gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry
0202     stfMCsquare(i).couchAngle  = stf(i).couchAngle;
0203     stfMCsquare(i).isoCenter   = stf(i).isoCenter + mcSquareAddIsoCenterOffset;
0204     stfMCsquare(i).energies    = unique([stf(i).ray.energy]);
0205     
0206     % allocate empty target point container
0207     for j = 1:numel(stfMCsquare(i).energies)
0208         stfMCsquare(i).energyLayer(j).targetPoints   = [];
0209         stfMCsquare(i).energyLayer(j).numOfPrimaries = [];
0210         stfMCsquare(i).energyLayer(j).rayNum         = [];
0211         stfMCsquare(i).energyLayer(j).bixelNum       = [];
0212     end
0213     
0214     for j = 1:stf(i).numOfRays
0215         for k = 1:stf(i).numOfBixelsPerRay(j)
0216             counter = counter + 1;
0217             dij.beamNum(counter)  = i;
0218             dij.rayNum(counter)   = j;
0219             dij.bixelNum(counter) = k;
0220         end
0221         
0222         for k = 1:numel(stfMCsquare(i).energies)
0223             if any(stf(i).ray(j).energy == stfMCsquare(i).energies(k))
0224                 stfMCsquare(i).energyLayer(k).rayNum   = [stfMCsquare(i).energyLayer(k).rayNum j];
0225                 stfMCsquare(i).energyLayer(k).bixelNum = [stfMCsquare(i).energyLayer(k).bixelNum ...
0226                     find(stf(i).ray(j).energy == stfMCsquare(i).energies(k))];
0227                 stfMCsquare(i).energyLayer(k).targetPoints = [stfMCsquare(i).energyLayer(k).targetPoints; ...
0228                                         -stf(i).ray(j).rayPos_bev(1) stf(i).ray(j).rayPos_bev(3)];
0229                 if calcDoseDirect
0230                     stfMCsquare(i).energyLayer(k).numOfPrimaries = [stfMCsquare(i).energyLayer(k).numOfPrimaries ...
0231                                          round(stf(i).ray(j).weight(stf(i).ray(j).energy == stfMCsquare(i).energies(k))*MCsquareConfig.Num_Primaries)];
0232                 else
0233                     stfMCsquare(i).energyLayer(k).numOfPrimaries = [stfMCsquare(i).energyLayer(k).numOfPrimaries ...
0234                         MCsquareConfig.Num_Primaries];
0235                 end
0236             end
0237         end
0238                
0239     end
0240     
0241 end
0242 
0243 % remember order
0244 counterMCsquare = 0;
0245 MCsquareOrder = NaN * ones(dij.totalNumOfBixels,1);
0246 for i = 1:length(stf)
0247     for j = 1:numel(stfMCsquare(i).energies)
0248         for k = 1:numel(stfMCsquare(i).energyLayer(j).numOfPrimaries)
0249             counterMCsquare = counterMCsquare + 1;
0250             ix = find(i                                         == dij.beamNum & ...
0251                       stfMCsquare(i).energyLayer(j).rayNum(k)   == dij.rayNum & ...
0252                       stfMCsquare(i).energyLayer(j).bixelNum(k) == dij.bixelNum);
0253         
0254             MCsquareOrder(ix) = counterMCsquare;
0255         end
0256     end
0257 end
0258 
0259 if any(isnan(MCsquareOrder))
0260     matRad_cfg.dispError('Invalid ordering of Beamlets for MCsquare computation!');
0261 end
0262 
0263 %% MC computation and dij filling
0264 matRad_writeMCsquareinputAllFiles(MCsquareConfigFile,MCsquareConfig,stfMCsquare);
0265 
0266 % run MCsquare
0267 [status,cmdout] = system([mcSquareBinary ' ' MCsquareConfigFile],'-echo');
0268     
0269 mask = false(dij.doseGrid.numOfVoxels,1);
0270 mask(VdoseGrid) = true;
0271 
0272 % read sparse matrix
0273 if ~calcDoseDirect
0274     dij.physicalDose{1} = absCalibrationFactorMC2 * matRad_sparseBeamletsReaderMCsquare ( ...
0275                     [MCsquareConfig.Output_Directory filesep 'Sparse_Dose.bin'], ...
0276                     dij.doseGrid.dimensions, ...
0277                     dij.totalNumOfBixels, ...
0278                     mask);
0279 else
0280     cube = matRad_readMhd(MCsquareConfig.Output_Directory,'Dose.mhd');
0281     dij.physicalDose{1} = sparse(VdoseGrid,ones(numel(VdoseGrid),1), ...
0282                                  absCalibrationFactorMC2 * cube(VdoseGrid), ...
0283                                  dij.doseGrid.numOfVoxels,1);
0284 end
0285 
0286 % reorder influence matrix to comply with matRad default ordering
0287 if MCsquareConfig.Beamlet_Mode
0288     dij.physicalDose{1} = dij.physicalDose{1}(:,MCsquareOrder);            
0289 end        
0290 
0291 matRad_cfg.dispInfo('matRad: done!\n');
0292 
0293 try
0294     % wait 0.1s for closing all waitbars
0295     allWaitBarFigures = findall(0,'type','figure','tag','TMWWaitbar');
0296     delete(allWaitBarFigures);
0297     pause(0.1);
0298 catch
0299 end
0300 
0301 %% clear all data
0302 delete([MCsquareConfig.CT_File(1:end-4) '.*']);
0303 delete('currBixels.txt');
0304 delete('MCsquareConfig.txt');
0305 
0306 %For Octave temporarily disable confirmation for recursive rmdir
0307 if strcmp(env,'OCTAVE')    
0308     rmdirConfirmState = confirm_recursive_rmdir(0);
0309 end
0310 rmdir(MCsquareConfig.Output_Directory,'s');
0311 
0312 %Reset to old confirmatoin state
0313 if strcmp(env,'OCTAVE')
0314     confirm_recursive_rmdir(rmdirConfirmState);
0315 end
0316 
0317 % cd back
0318 cd(currFolder);
0319 
0320 end
matRad_calcParticleDoseMC

Purpose

Synopsis

Description

Cross-reference information

Source code
