matRad_calcPhotonDoseMC

 matRad ompMC monte carlo photon dose calculation wrapper

 call
   dij = matRad_calcPhotonDoseMc(ct,stf,pln,cst,visBool)

 input
   ct:                         matRad ct struct
   stf:                        matRad steering information struct
   pln:                        matRad plan meta information struct
   cst:                        matRad cst struct
   visBool:                    binary switch to enable visualization
 output
   dij:                        matRad dij struct

 References
   -

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Copyright 2018 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.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,nCasePerBixel,visBool)
0002 % matRad ompMC monte carlo photon dose calculation wrapper
0003 %
0004 % call
0005 %   dij = matRad_calcPhotonDoseMc(ct,stf,pln,cst,visBool)
0006 %
0007 % input
0008 %   ct:                         matRad ct struct
0009 %   stf:                        matRad steering information struct
0010 %   pln:                        matRad plan meta information struct
0011 %   cst:                        matRad cst struct
0012 %   visBool:                    binary switch to enable visualization
0013 % output
0014 %   dij:                        matRad dij struct
0015 %
0016 % References
0017 %   -
0018 %
0019 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 %
0021 % Copyright 2018 the matRad development team.
0022 %
0023 % This file is part of the matRad project. It is subject to the license
0024 % terms in the LICENSE file found in the top-level directory of this
0025 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0026 % of the matRad project, including this file, may be copied, modified,
0027 % propagated, or distributed except according to the terms contained in the
0028 % LICENSE file.
0029 %
0030 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0031 
0032 
0033 matRad_cfg =  MatRad_Config.instance();
0034 
0035 tic
0036 
0037 % disable visualiazation by default
0038 if nargin < 6
0039     visBool = false;
0040 end
0041 
0042 if nargin < 5
0043     nCasePerBixel = matRad_cfg.propMC.ompMC_defaultHistories;
0044     matRad_cfg.dispInfo('Using default number of Histories per Bixel: %d\n',nCasePerBixel);
0045 end
0046 
0047 fileFolder = fileparts(mfilename('fullpath'));
0048 
0049 if ~matRad_checkMexFileExists('omc_matrad') %exist('matRad_ompInterface','file') ~= 3
0050     matRad_cfg.dispWarning('Compiled mex interface not found. Trying to compile the ompMC interface on the fly!');    
0051     try        
0052         matRad_compileOmpMCInterface();
0053     catch MException       
0054         matRad_cfg.dispError('Could not find/generate mex interface for MC dose calculation.\nCause of error:\n%s\n Please compile it yourself (preferably with OpenMP support).',MException.message);
0055     end
0056 end
0057 
0058 matRad_calcDoseInit;
0059 
0060 % set up arrays for book keeping
0061 dij.bixelNum = NaN*ones(dij.totalNumOfBixels,1);
0062 dij.rayNum   = NaN*ones(dij.totalNumOfBixels,1);
0063 dij.beamNum  = NaN*ones(dij.totalNumOfBixels,1);
0064 
0065 dij.numHistoriesPerBeamlet = nCasePerBixel;
0066 
0067 omcFolder = [matRad_cfg.matRadRoot filesep 'ompMC'];
0068 %omcFolder = [matRad_cfg.matRadRoot filesep 'submodules' filesep 'ompMC'];
0069 
0070 %% Setup OmpMC options / parameters
0071 
0072 %display options
0073 ompMCoptions.verbose = matRad_cfg.logLevel - 1;
0074 
0075 % start MC control
0076 ompMCoptions.nHistories = nCasePerBixel;
0077 ompMCoptions.nSplit = 20;
0078 ompMCoptions.nBatches = 10;
0079 ompMCoptions.randomSeeds = [97 33];
0080 
0081 %start source definition
0082 ompMCoptions.spectrumFile = [omcFolder filesep 'spectra' filesep 'mohan6.spectrum'];
0083 ompMCoptions.monoEnergy   = 0.1; 
0084 ompMCoptions.charge       = 0;
0085                                                                     
0086 % start MC transport
0087 ompMCoptions.dataFolder   = [omcFolder filesep 'data' filesep];
0088 ompMCoptions.pegsFile     = [omcFolder filesep 'pegs4' filesep '700icru.pegs4dat'];
0089 ompMCoptions.pgs4formFile = [omcFolder filesep 'pegs4' filesep 'pgs4form.dat'];
0090 
0091 ompMCoptions.global_ecut = 0.7;
0092 ompMCoptions.global_pcut = 0.010; 
0093 
0094 % Relative Threshold for dose
0095 ompMCoptions.relDoseThreshold = 1 - matRad_cfg.propDoseCalc.defaultLateralCutOff;
0096 
0097 % Output folders
0098 ompMCoptions.outputFolder = [omcFolder filesep 'output' filesep];
0099 
0100 % Create Material Density Cube
0101 material = cell(4,5);
0102 material{1,1} = 'AIR700ICRU';
0103 material{1,2} = -1024; 
0104 material{1,3} = -974;
0105 material{1,4} = 0.001;
0106 material{1,5} = 0.044;
0107 material{2,1} = 'LUNG700ICRU';
0108 material{2,2} = -974; 
0109 material{2,3} = -724;
0110 material{2,4} = 0.044; 
0111 material{2,5} = 0.302;
0112 material{3,1} = 'ICRUTISSUE700ICRU';
0113 material{3,2} = -724; 
0114 material{3,3} = 101;
0115 material{3,4} = 0.302; 
0116 material{3,5} = 1.101;
0117 material{4,1} = 'ICRPBONE700ICRU';
0118 material{4,2} = 101; 
0119 material{4,3} = 1976;
0120 material{4,4} = 1.101; 
0121 material{4,5} = 2.088;
0122 
0123 % conversion from HU to densities & materials
0124 for s = 1:dij.numOfScenarios
0125 
0126     HUcube{s} =  matRad_interp3(dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z,ct.cubeHU{s}, ...
0127                             dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest');
0128     
0129     % projecting out of bounds HU values where necessary
0130     if max(HUcube{s}(:)) > material{end,3}
0131         matRad_cfg.dispWarning('Projecting out of range HU values');
0132         HUcube{s}(HUcube{s}(:) > material{end,3}) = material{end,3};
0133     end
0134     if min(HUcube{s}(:)) < material{1,2}
0135         matRad_cfg.dispWarning('Projecting out of range HU values');
0136         HUcube{s}(HUcube{s}(:) < material{1,2}) = material{1,2};
0137     end
0138 
0139     % find material index
0140     cubeMatIx{s} = NaN*ones(dij.doseGrid.dimensions,'int32');
0141     for i = size(material,1):-1:1
0142         cubeMatIx{s}(HUcube{s} <= material{i,3}) = i;
0143     end
0144     
0145     % create an artificial HU lookup table
0146     hlut = [];
0147     for i = 1:size(material,1)       
0148         hlut = [hlut;material{i,2} material{i,4};material{i,3}-1e-10 material{i,5}]; % add eps for interpolation
0149     end
0150     
0151     cubeRho{s} = interp1(hlut(:,1),hlut(:,2),HUcube{s});
0152 
0153 end
0154 
0155 ompMCgeo.material = material;
0156 
0157 scale = 10; % to convert to cm
0158 
0159 ompMCgeo.xBounds = (dij.doseGrid.resolution.y * (0.5 + [0:dij.doseGrid.dimensions(1)])) ./ scale;
0160 ompMCgeo.yBounds = (dij.doseGrid.resolution.x * (0.5 + [0:dij.doseGrid.dimensions(2)])) ./ scale;
0161 ompMCgeo.zBounds = (dij.doseGrid.resolution.z * (0.5 + [0:dij.doseGrid.dimensions(3)])) ./ scale;
0162 
0163 %% debug visualization
0164 if visBool
0165     
0166     figure
0167     hold on
0168 
0169     axis equal
0170     
0171     % ct box
0172     ctCorner1 = [ompMCgeo.xBounds(1) ompMCgeo.yBounds(1) ompMCgeo.zBounds(1)];
0173     ctCorner2 = [ompMCgeo.xBounds(end) ompMCgeo.yBounds(end) ompMCgeo.zBounds(end)];
0174     plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner1(3)],'k' )
0175     plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
0176     plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
0177     plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
0178     plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner2(3) ctCorner2(3)],'k' )
0179     plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
0180     plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
0181     plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
0182     plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' )
0183     plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' )
0184     plot3([ctCorner1(1) ctCorner1(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' )
0185     plot3([ctCorner2(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' )
0186         
0187     xlabel('x [cm]')
0188     ylabel('y [cm]')
0189     zlabel('z [cm]')
0190 
0191     rotate3d on
0192     
0193 end
0194 
0195 %% Create beamlet source
0196 useCornersSCD = false; %false -> use ISO corners
0197 
0198 numOfBixels = [stf(:).numOfRays];
0199 beamSource = zeros(dij.numOfBeams, 3);
0200 
0201 bixelCorner = zeros(dij.totalNumOfBixels,3);
0202 bixelSide1 = zeros(dij.totalNumOfBixels,3);
0203 bixelSide2 = zeros(dij.totalNumOfBixels,3);
0204 
0205 counter = 0;
0206 
0207 for i = 1:dij.numOfBeams % loop over all beams
0208    
0209     % define beam source in physical coordinate system in cm
0210     beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10;
0211 
0212     for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray!
0213         
0214         counter = counter + 1;
0215         
0216         dij.beamNum(counter)  = i;
0217         dij.rayNum(counter)   = j;
0218         dij.bixelNum(counter) = j;
0219         
0220         if useCornersSCD
0221             beamletCorners = stf(i).ray(j).rayCorners_SCD;
0222         else    
0223             beamletCorners = stf(i).ray(j).beamletCornersAtIso;
0224         end
0225         
0226         % get bixel corner and delimiting vectors.
0227         % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm
0228         currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ scale;
0229         bixelCorner(counter,:) = currCorner;
0230         bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ scale - currCorner;
0231         bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ scale - currCorner;
0232         
0233         if visBool
0234             for k = 1:4
0235                 currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10;
0236                 % rays connecting source and ray corner
0237                 plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b')
0238                 % connection between corners
0239                 lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10;
0240                 plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r')
0241             end
0242         end
0243         
0244     end
0245         
0246 end
0247 
0248 ompMCsource.nBeams = dij.numOfBeams;
0249 ompMCsource.iBeam = dij.beamNum(:);
0250 
0251 % Switch x and y directions to match ompMC cs.
0252 ompMCsource.xSource = beamSource(:,2);
0253 ompMCsource.ySource = beamSource(:,1);
0254 ompMCsource.zSource = beamSource(:,3);
0255 
0256 ompMCsource.nBixels = sum(numOfBixels(:));
0257 ompMCsource.xCorner = bixelCorner(:,2);
0258 ompMCsource.yCorner = bixelCorner(:,1);
0259 ompMCsource.zCorner = bixelCorner(:,3);
0260 
0261 ompMCsource.xSide1 = bixelSide1(:,2);
0262 ompMCsource.ySide1 = bixelSide1(:,1);
0263 ompMCsource.zSide1 = bixelSide1(:,3);
0264 
0265 ompMCsource.xSide2 = bixelSide2(:,2);
0266 ompMCsource.ySide2 = bixelSide2(:,1);
0267 ompMCsource.zSide2 = bixelSide2(:,3);
0268 
0269 if visBool
0270     plot3(ompMCsource.ySource,ompMCsource.xSource,ompMCsource.zSource,'rx')
0271 end
0272 
0273 %% Call the OmpMC interface
0274 
0275 %ompMC for matRad returns dose/history * nHistories.
0276 % This factor calibrates to 1 Gy in a %(5x5)cm^2 open field (1 bixel) at
0277 % 5cm depth for SSD = 900 which corresponds to the calibration for the
0278 % analytical base data.
0279 absCalibrationFactor = 3.49056 * 1e12; %Approximate!
0280 
0281 %Now we have to calibrate to the the beamlet width.
0282 absCalibrationFactor = absCalibrationFactor * (pln.propStf.bixelWidth/50)^2;
0283 
0284 matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation... \n');
0285 
0286 outputVariance = matRad_cfg.propMC.ompMC_defaultOutputVariance;
0287 
0288 if isfield(pln,'propMC') && isfield(pln.propMC,'outputVariance')
0289     outputVariance = pln.propMC.outputVariance;
0290 end
0291 
0292 
0293 %run over all scenarios
0294 for s = 1:dij.numOfScenarios
0295     ompMCgeo.isoCenter = [stf(:).isoCenter];
0296     
0297     %Run the Monte Carlo simulation and catch  possible mex-interface
0298     %issues
0299     try
0300         %If we ask for variance, a field in the dij will be filled
0301         if outputVariance
0302             [dij.physicalDose{s},dij.physicalDose_MCvar{s}] = omc_matrad(cubeRho{s},cubeMatIx{s},ompMCgeo,ompMCsource,ompMCoptions);
0303         else
0304             [dij.physicalDose{s}] = omc_matrad(cubeRho{s},cubeMatIx{s},ompMCgeo,ompMCsource,ompMCoptions);
0305         end
0306     catch ME
0307         errorString = [ME.message '\nThis error was thrown by the MEX-interface of ompMC.\nMex interfaces can raise compatability issues which may be resolved by compiling them by hand directly on your particular system.'];
0308         matRad_cfg.dispError(ME.identifier,errorString);
0309     end
0310     
0311     %Calibrate the dose with above factor
0312     dij.physicalDose{s} = dij.physicalDose{s} * absCalibrationFactor;
0313     if isfield(dij,'physicalDose_MCvar')
0314         dij.physicalDose_MCvar{s} = dij.physicalDose_MCvar{s} * absCalibrationFactor^2;
0315     end
0316 end
0317 
0318 matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n');
0319 matRad_cfg.dispInfo(evalc('toc'));
0320 
0321 try
0322     % wait 0.1s for closing all waitbars
0323     allWaitBarFigures = findall(0,'type','figure','tag','TMWWaitbar');
0324     delete(allWaitBarFigures);
0325     pause(0.1);
0326 catch
0327 end
0328 
0329 end