matRad_calcParticleDose

 matRad particle dose calculation wrapper
 
 call
   dij = matRad_calcParticleDose(ct,stf,pln,cst,calcDoseDirect)

 input
   ct:             ct cube
   stf:            matRad steering information struct
   pln:            matRad plan meta information struct
   cst:            matRad cst struct
   calcDoseDirect: boolian switch to bypass dose influence matrix
                   computation and directly calculate dose; only makes
                   sense in combination with matRad_calcDoseDirect.m

 output
   dij:            matRad dij struct

 References
   [1] http://iopscience.iop.org/0031-9155/41/8/005

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

 Copyright 2015 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_calcParticleDose(ct,stf,pln,cst,calcDoseDirect)
0002 % matRad particle dose calculation wrapper
0003 %
0004 % call
0005 %   dij = matRad_calcParticleDose(ct,stf,pln,cst,calcDoseDirect)
0006 %
0007 % input
0008 %   ct:             ct cube
0009 %   stf:            matRad steering information struct
0010 %   pln:            matRad plan meta information struct
0011 %   cst:            matRad cst struct
0012 %   calcDoseDirect: boolian switch to bypass dose influence matrix
0013 %                   computation and directly calculate dose; only makes
0014 %                   sense in combination with matRad_calcDoseDirect.m
0015 %
0016 % output
0017 %   dij:            matRad dij struct
0018 %
0019 % References
0020 %   [1] http://iopscience.iop.org/0031-9155/41/8/005
0021 %
0022 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 %
0024 % Copyright 2015 the matRad development team.
0025 %
0026 % This file is part of the matRad project. It is subject to the license
0027 % terms in the LICENSE file found in the top-level directory of this
0028 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0029 % of the matRad project, including this file, may be copied, modified,
0030 % propagated, or distributed except according to the terms contained in the
0031 % LICENSE file.
0032 %
0033 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034 
0035 
0036 matRad_cfg =  MatRad_Config.instance();
0037 
0038 % init dose calc
0039 matRad_calcDoseInit;
0040 
0041 % initialize waitbar
0042 figureWait = waitbar(0,'calculate dose influence matrix for particles...');
0043 % prevent closure of waitbar and show busy state
0044 set(figureWait,'pointer','watch');
0045 
0046 % helper function for energy selection
0047 round2 = @(a,b)round(a*10^b)/10^b;
0048 
0049 if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
0050         && strcmp(pln.radiationMode,'carbon')
0051    
0052         alphaDoseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios);
0053         betaDoseTmpContainer  = cell(numOfBixelsContainer,dij.numOfScenarios);
0054         for i = 1:dij.numOfScenarios
0055             dij.mAlphaDose{i}    = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1);
0056             dij.mSqrtBetaDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1);
0057         end
0058         
0059 elseif isequal(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons')
0060             dij.RBE = 1.1;
0061             matRad_cfg.dispInfo('matRad: Using a constant RBE of %g\n',dij.RBE);   
0062 end
0063 
0064 if isfield(pln,'propDoseCalc') && ...
0065    isfield(pln.propDoseCalc,'calcLET') && ...
0066    pln.propDoseCalc.calcLET
0067   if isfield(machine.data,'LET')
0068     letDoseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios);
0069     % Allocate space for dij.dosexLET sparse matrix
0070     for i = 1:dij.numOfScenarios
0071         dij.mLETDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1);
0072     end
0073   else
0074     matRad_cfg.dispWarning('LET not available in the machine data. LET will not be calculated.');
0075   end
0076 end
0077 
0078 % generates tissue class matrix for biological optimization
0079 if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
0080         && strcmp(pln.radiationMode,'carbon')
0081     
0082     if   isfield(machine.data,'alphaX') && isfield(machine.data,'betaX')
0083             
0084         matRad_cfg.dispInfo('matRad: loading biological base data... ');
0085         vTissueIndex = zeros(size(VdoseGrid,1),1);
0086         dij.ax       = zeros(dij.doseGrid.numOfVoxels,1);
0087         dij.bx       = zeros(dij.doseGrid.numOfVoxels,1);
0088 
0089         cst = matRad_setOverlapPriorities(cst);
0090     
0091         % resizing cst to dose cube resolution
0092         cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
0093                                          dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
0094         % retrieve photon LQM parameter for the current dose grid voxels
0095         [dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid);
0096 
0097         for i = 1:size(cst,1)
0098 
0099             % check if cst is compatiable
0100             if ~isempty(cst{i,5}) && isfield(cst{i,5},'alphaX') && isfield(cst{i,5},'betaX') 
0101 
0102                 % check if base data contains alphaX and betaX
0103                 IdxTissue = find(ismember(machine.data(1).alphaX,cst{i,5}.alphaX) & ...
0104                                  ismember(machine.data(1).betaX,cst{i,5}.betaX));
0105 
0106                 % check consitency of biological baseData and cst settings
0107                 if ~isempty(IdxTissue)
0108                     isInVdoseGrid = ismember(VdoseGrid,cst{i,4}{1});
0109                     vTissueIndex(isInVdoseGrid) = IdxTissue;
0110                 else
0111                     matRad_cfg.dispError('biological base data and cst inconsistent\n');
0112                 end
0113                     
0114             else
0115                     vTissueIndex(row) = 1;
0116                     matRad_cfg.dispInfo(['matRad: tissue type of ' cst{i,2} ' was set to 1\n']);          
0117             end
0118         end
0119         matRad_cfg.dispInfo('done.\n');
0120 
0121     else
0122         
0123         matRad_cfg.dispError('base data is incomplement - alphaX and/or betaX is missing');
0124         
0125     end
0126     
0127 % issue warning if biological optimization not possible
0128 elseif sum(strcmp(pln.propOpt.bioOptimization,{'LEMIV_effect','LEMIV_RBExD'}))>0 && ~strcmp(pln.radiationMode,'carbon') ||...
0129        ~strcmp(pln.radiationMode,'protons') && strcmp(pln.propOpt.bioOptimization,'const_RBExD')
0130     warndlg([pln.propOpt.bioOptimization ' optimization not possible with ' pln.radiationMode '- physical optimization is carried out instead.']);
0131     pln.propOpt.bioOptimization = 'none';      
0132 end
0133 
0134 % lateral cutoff for raytracing and geo calculations
0135 effectiveLateralCutoff = matRad_cfg.propDoseCalc.defaultGeometricCutOff;
0136 
0137 matRad_cfg.dispInfo('matRad: Particle dose calculation...\n');
0138 counter = 0;
0139 
0140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0141 for i = 1:length(stf) % loop over all beams
0142   
0143     % init beam
0144     matRad_calcDoseInitBeam;
0145   
0146     % Determine lateral cutoff
0147     matRad_cfg.dispInfo('matRad: calculate lateral cutoff...');
0148     cutOffLevel = matRad_cfg.propDoseCalc.defaultLateralCutOff;
0149     visBoolLateralCutOff = 0;
0150     machine = matRad_calcLateralParticleCutOff(machine,cutOffLevel,stf(i),visBoolLateralCutOff);
0151     matRad_cfg.dispInfo('done.\n');    
0152 
0153     for j = 1:stf(i).numOfRays % loop over all rays
0154 
0155         if ~isempty(stf(i).ray(j).energy)
0156 
0157             % find index of maximum used energy (round to keV for numerical
0158             % reasons
0159             energyIx = max(round2(stf(i).ray(j).energy,4)) == round2([machine.data.energy],4);
0160 
0161             maxLateralCutoffDoseCalc = max(machine.data(energyIx).LatCutOff.CutOff);
0162 
0163             % Ray tracing for beam i and ray j
0164             [ix,radialDist_sq] = matRad_calcGeoDists(rot_coordsVdoseGrid, ...
0165                                                      stf(i).sourcePoint_bev, ...
0166                                                      stf(i).ray(j).targetPoint_bev, ...
0167                                                      machine.meta.SAD, ...
0168                                                      find(~isnan(radDepthVdoseGrid{1})), ...
0169                                                      maxLateralCutoffDoseCalc);
0170              
0171             radDepths = radDepthVdoseGrid{1}(ix);   
0172                        
0173             % just use tissue classes of voxels found by ray tracer
0174             if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
0175                 && strcmp(pln.radiationMode,'carbon')
0176                     vTissueIndex_j = vTissueIndex(ix,:);
0177             end
0178 
0179             for k = 1:stf(i).numOfBixelsPerRay(j) % loop over all bixels per ray
0180 
0181                 counter = counter + 1;
0182                 bixelsPerBeam = bixelsPerBeam + 1;
0183                 
0184                 % Display progress and update text only 200 times
0185                 if mod(bixelsPerBeam,max(1,round(stf(i).totalNumOfBixels/200))) == 0
0186                         matRad_progress(bixelsPerBeam/max(1,round(stf(i).totalNumOfBixels/200)),...
0187                                         floor(stf(i).totalNumOfBixels/max(1,round(stf(i).totalNumOfBixels/200))));
0188                 end
0189                 
0190                 % update waitbar only 100 times if it is not closed
0191                 if mod(counter,round(dij.totalNumOfBixels/100)) == 0 && ishandle(figureWait)
0192                     waitbar(counter/dij.totalNumOfBixels,figureWait);
0193                 end
0194 
0195                 % remember beam and bixel number
0196                 if ~calcDoseDirect
0197                    dij.beamNum(counter)  = i;
0198                    dij.rayNum(counter)   = j;
0199                    dij.bixelNum(counter) = k;
0200                 end
0201                 
0202                 % find energy index in base data
0203                 energyIx = find(round2(stf(i).ray(j).energy(k),4) == round2([machine.data.energy],4));
0204                 
0205                 % create offset vector to account for additional offsets modelled in the base data and a potential
0206                 % range shifter. In the following, we only perform dose calculation for voxels having a radiological depth
0207                 % that is within the limits of the base data set (-> machine.data(i).dephts). By this means, we only allow
0208                 % interpolations in matRad_calcParticleDoseBixel() and avoid extrapolations.
0209                 offsetRadDepth = machine.data(energyIx).offset - stf(i).ray(j).rangeShifter(k).eqThickness;
0210                 
0211                 % find depth depended lateral cut off
0212                 if cutOffLevel >= 1
0213                     currIx = radDepths <= machine.data(energyIx).depths(end) + offsetRadDepth;
0214                 elseif cutOffLevel < 1 && cutOffLevel > 0
0215                     % perform rough 2D clipping
0216                     currIx = radDepths <= machine.data(energyIx).depths(end) + offsetRadDepth & ...
0217                          radialDist_sq <= max(machine.data(energyIx).LatCutOff.CutOff.^2);
0218 
0219                     % peform fine 2D clipping
0220                     if length(machine.data(energyIx).LatCutOff.CutOff) > 1
0221                         currIx(currIx) = matRad_interp1((machine.data(energyIx).LatCutOff.depths + offsetRadDepth)',...
0222                             (machine.data(energyIx).LatCutOff.CutOff.^2)', radDepths(currIx)) >= radialDist_sq(currIx);
0223                     end
0224                 else
0225                     matRad_cfg.dispError('cutoff must be a value between 0 and 1')
0226                 end
0227                 
0228                 % empty bixels may happen during recalculation of error
0229                 % scenarios -> skip to next bixel
0230                 if ~any(currIx)
0231                     continue;
0232                 end
0233                 
0234                 % adjust radDepth according to range shifter
0235                 currRadDepths = radDepths(currIx) + stf(i).ray(j).rangeShifter(k).eqThickness;
0236 
0237                 % calculate initial focus sigma
0238                 sigmaIni = matRad_interp1(machine.data(energyIx).initFocus.dist (stf(i).ray(j).focusIx(k),:)', ...
0239                                              machine.data(energyIx).initFocus.sigma(stf(i).ray(j).focusIx(k),:)',stf(i).ray(j).SSD);
0240                 sigmaIni_sq = sigmaIni^2;
0241                 
0242                 % consider range shifter for protons if applicable
0243                 if stf(i).ray(j).rangeShifter(k).eqThickness > 0 && strcmp(pln.radiationMode,'protons')
0244                     
0245                     % compute!
0246                     sigmaRashi = matRad_calcSigmaRashi(machine.data(energyIx).energy, ...
0247                                                        stf(i).ray(j).rangeShifter(k), ...
0248                                                        stf(i).ray(j).SSD);
0249                               
0250                     % add to initial sigma in quadrature
0251                     sigmaIni_sq = sigmaIni_sq +  sigmaRashi^2;
0252                     
0253                 end
0254                                 
0255                 % calculate particle dose for bixel k on ray j of beam i
0256                 bixelDose = matRad_calcParticleDoseBixel(...
0257                     currRadDepths, ...
0258                     radialDist_sq(currIx), ...
0259                     sigmaIni_sq, ...
0260                     machine.data(energyIx));                 
0261   
0262                 % dij sampling is exluded for particles until we investigated the influence of voxel sampling for particles
0263                 %relDoseThreshold   =  0.02;   % sample dose values beyond the relative dose
0264                 %Type               = 'dose';
0265                 %[currIx,bixelDose] = matRad_DijSampling(currIx,bixelDose,radDepths(currIx),radialDist_sq(currIx),Type,relDoseThreshold);
0266                 
0267                 % Save dose for every bixel in cell array
0268                 doseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1} = sparse(VdoseGrid(ix(currIx)),1,bixelDose,dij.doseGrid.numOfVoxels,1);
0269 
0270                 if isfield(dij,'mLETDose')
0271                   % calculate particle LET for bixel k on ray j of beam i
0272                   depths = machine.data(energyIx).depths + machine.data(energyIx).offset; 
0273                   bixelLET = matRad_interp1(depths,machine.data(energyIx).LET,currRadDepths); 
0274 
0275                   % Save LET for every bixel in cell array
0276                   letDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1} = sparse(VdoseGrid(ix(currIx)),1,bixelLET.*bixelDose,dij.doseGrid.numOfVoxels,1);
0277                 end
0278                              
0279                 if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
0280                     && strcmp(pln.radiationMode,'carbon')
0281                     % calculate alpha and beta values for bixel k on ray j of
0282                     [bixelAlpha, bixelBeta] = matRad_calcLQParameter(...
0283                         currRadDepths,...
0284                         vTissueIndex_j(currIx,:),...
0285                         machine.data(energyIx));
0286                     
0287                     alphaDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1} = sparse(VdoseGrid(ix(currIx)),1,bixelAlpha.*bixelDose,dij.doseGrid.numOfVoxels,1);
0288                     betaDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1}  = sparse(VdoseGrid(ix(currIx)),1,sqrt(bixelBeta).*bixelDose,dij.doseGrid.numOfVoxels,1);
0289                 end
0290                 
0291                 matRad_calcDoseFillDij;                
0292 
0293             end
0294             
0295         end
0296         
0297     end
0298 end
0299 
0300 %Close Waitbar
0301 if ishandle(figureWait)
0302     delete(figureWait);
0303 end
0304 
0305