matRad_calcPhotonDose

 matRad photon dose calculation wrapper
 
 call
   dij = matRad_calcPhotonDose(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://www.ncbi.nlm.nih.gov/pubmed/8497215

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

 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_calcPhotonDose(ct,stf,pln,cst,calcDoseDirect)
0002 % matRad photon dose calculation wrapper
0003 %
0004 % call
0005 %   dij = matRad_calcPhotonDose(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://www.ncbi.nlm.nih.gov/pubmed/8497215
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 % initialize
0039 matRad_calcDoseInit;
0040 
0041 [env, ~] = matRad_getEnvironment();
0042 
0043 switch env
0044      case 'MATLAB'
0045           rng('default');
0046      case 'OCTAVE'
0047           rand('seed',0)
0048 end
0049 
0050 % issue warning if biological optimization not possible
0051 if sum(strcmp(pln.propOpt.bioOptimization,{'effect','RBExD'}))>0
0052     warndlg('Effect based and RBE optimization not available for photons - physical optimization is carried out instead.');
0053     pln.bioOptimization = 'none';
0054 end
0055 
0056 % initialize waitbar
0057 figureWait = waitbar(0,'calculate dose influence matrix for photons...');
0058 % show busy state
0059 set(figureWait,'pointer','watch');
0060 
0061 % set lateral cutoff value
0062 lateralCutoff = matRad_cfg.propDoseCalc.defaultGeometricCutOff; % [mm]
0063 
0064 % toggle custom primary fluence on/off. if 0 we assume a homogeneous
0065 % primary fluence, if 1 we use measured radially symmetric data
0066 if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'useCustomPrimaryPhotonFluence')
0067     useCustomPrimFluenceBool = matRad_cfg.propDoseCalc.defaultUseCustomPrimaryPhotonFluence;
0068 else
0069     useCustomPrimFluenceBool = pln.propDoseCalc.useCustomPrimaryPhotonFluence;
0070 end
0071 
0072 
0073 % 0 if field calc is bixel based, 1 if dose calc is field based
0074 isFieldBasedDoseCalc = strcmp(num2str(pln.propStf.bixelWidth),'field');
0075 
0076 %% kernel convolution
0077 % set up convolution grid
0078 if isFieldBasedDoseCalc
0079     % get data from DICOM import
0080     intConvResolution = pln.propStf.collimation.convResolution; 
0081     fieldWidth = pln.propStf.collimation.fieldWidth;
0082 else
0083     intConvResolution = .5; % [mm]
0084     fieldWidth = pln.propStf.bixelWidth;
0085 end
0086 
0087 % calculate field size and distances
0088 fieldLimit = ceil(fieldWidth/(2*intConvResolution));
0089 [F_X,F_Z] = meshgrid(-fieldLimit*intConvResolution: ...
0090                   intConvResolution: ...
0091                   (fieldLimit-1)*intConvResolution);    
0092 
0093 % gaussian filter to model penumbra
0094 sigmaGauss = 2.123; % [mm] / see diploma thesis siggel 4.1.2
0095 % use 5 times sigma as the limits for the gaussian convolution
0096 gaussLimit = ceil(5*sigmaGauss/intConvResolution);
0097 [gaussFilterX,gaussFilterZ] = meshgrid(-gaussLimit*intConvResolution: ...
0098                                     intConvResolution: ...
0099                                     (gaussLimit-1)*intConvResolution);   
0100 gaussFilter =  1/(2*pi*sigmaGauss^2/intConvResolution^2) * exp(-(gaussFilterX.^2+gaussFilterZ.^2)/(2*sigmaGauss^2) );
0101 gaussConvSize = 2*(fieldLimit + gaussLimit);
0102 
0103 if ~isFieldBasedDoseCalc   
0104     % Create fluence matrix
0105     F = ones(floor(fieldWidth/intConvResolution));
0106     
0107     if ~useCustomPrimFluenceBool
0108     % gaussian convolution of field to model penumbra
0109         F = real(ifft2(fft2(F,gaussConvSize,gaussConvSize).*fft2(gaussFilter,gaussConvSize,gaussConvSize)));     
0110     end
0111 end
0112 
0113 % get kernel size and distances
0114 kernelLimit = ceil(lateralCutoff/intConvResolution);
0115 [kernelX, kernelZ] = meshgrid(-kernelLimit*intConvResolution: ...
0116                             intConvResolution: ...
0117                             (kernelLimit-1)*intConvResolution);
0118 
0119 % precalculate convoluted kernel size and distances
0120 kernelConvLimit = fieldLimit + gaussLimit + kernelLimit;
0121 [convMx_X, convMx_Z] = meshgrid(-kernelConvLimit*intConvResolution: ...
0122                                 intConvResolution: ...
0123                                 (kernelConvLimit-1)*intConvResolution);
0124 % calculate also the total size and distance as we need this during convolution extensively
0125 kernelConvSize = 2*kernelConvLimit;
0126 
0127 % define an effective lateral cutoff where dose will be calculated. note
0128 % that storage within the influence matrix may be subject to sampling
0129 effectiveLateralCutoff = lateralCutoff + fieldWidth/2;
0130 
0131 counter = 0;
0132 matRad_cfg.dispInfo('matRad: Photon dose calculation...\n');
0133 
0134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0135 for i = 1:dij.numOfBeams % loop over all beams
0136     
0137     matRad_calcDoseInitBeam;
0138     
0139     % get index of central ray or closest to the central ray
0140     [~,center] = min(sum(reshape([stf(i).ray.rayPos_bev],3,[]).^2));
0141     
0142     % get correct kernel for given SSD at central ray (nearest neighbor approximation)
0143     [~,currSSDIx] = min(abs([machine.data.kernel.SSD]-stf(i).ray(center).SSD));
0144     % Display console message.
0145     matRad_cfg.dispInfo('\tSSD = %g mm ...\n',machine.data.kernel(currSSDIx).SSD);
0146     
0147     kernelPos = machine.data.kernelPos;
0148     kernel1 = machine.data.kernel(currSSDIx).kernel1;
0149     kernel2 = machine.data.kernel(currSSDIx).kernel2;
0150     kernel3 = machine.data.kernel(currSSDIx).kernel3;
0151 
0152     % Evaluate kernels for all distances, interpolate between values
0153     kernel1Mx = interp1(kernelPos,kernel1,sqrt(kernelX.^2+kernelZ.^2),'linear',0);
0154     kernel2Mx = interp1(kernelPos,kernel2,sqrt(kernelX.^2+kernelZ.^2),'linear',0);
0155     kernel3Mx = interp1(kernelPos,kernel3,sqrt(kernelX.^2+kernelZ.^2),'linear',0);
0156     
0157     % convolution here if no custom primary fluence and no field based dose calc
0158     if ~useCustomPrimFluenceBool && ~isFieldBasedDoseCalc
0159         
0160         % Display console message.
0161         matRad_cfg.dispInfo('\tUniform primary photon fluence -> pre-compute kernel convolution...\n');   
0162 
0163         % 2D convolution of Fluence and Kernels in fourier domain
0164         convMx1 = real(ifft2(fft2(F,kernelConvSize,kernelConvSize).* fft2(kernel1Mx,kernelConvSize,kernelConvSize)));
0165         convMx2 = real(ifft2(fft2(F,kernelConvSize,kernelConvSize).* fft2(kernel2Mx,kernelConvSize,kernelConvSize)));
0166         convMx3 = real(ifft2(fft2(F,kernelConvSize,kernelConvSize).* fft2(kernel3Mx,kernelConvSize,kernelConvSize)));
0167 
0168         % Creates an interpolant for kernes from vectors position X and Z
0169         if strcmp(env,'MATLAB')
0170             Interp_kernel1 = griddedInterpolant(convMx_X',convMx_Z',convMx1','linear','none');
0171             Interp_kernel2 = griddedInterpolant(convMx_X',convMx_Z',convMx2','linear','none');
0172             Interp_kernel3 = griddedInterpolant(convMx_X',convMx_Z',convMx3','linear','none');
0173         else
0174             Interp_kernel1 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx1,x,y,'linear',NaN);
0175             Interp_kernel2 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx2,x,y,'linear',NaN);
0176             Interp_kernel3 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx3,x,y,'linear',NaN);
0177         end
0178     end
0179     
0180     for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray!
0181 
0182         counter = counter + 1;
0183         bixelsPerBeam = bixelsPerBeam + 1;
0184     
0185         % convolution here if custom primary fluence OR field based dose calc
0186         if useCustomPrimFluenceBool || isFieldBasedDoseCalc
0187             
0188             % overwrite field opening if necessary
0189             if isFieldBasedDoseCalc
0190                 F = stf(i).ray(j).shape;
0191             end
0192             
0193             % prepare primary fluence array
0194             primaryFluence = machine.data.primaryFluence;
0195             r     = sqrt( (F_X-stf(i).ray(j).rayPos(1)).^2 + (F_Z-stf(i).ray(j).rayPos(3)).^2 );
0196             Psi   = interp1(primaryFluence(:,1)',primaryFluence(:,2)',r,'linear',0);
0197                 
0198             % apply the primary fluence to the field
0199             Fx = F .* Psi;
0200             
0201             % convolute with the gaussian
0202             Fx = real( ifft2(fft2(Fx,gaussConvSize,gaussConvSize).* fft2(gaussFilter,gaussConvSize,gaussConvSize)) );
0203 
0204             % 2D convolution of Fluence and Kernels in fourier domain
0205             convMx1 = real( ifft2(fft2(Fx,kernelConvSize,kernelConvSize).* fft2(kernel1Mx,kernelConvSize,kernelConvSize)) );
0206             convMx2 = real( ifft2(fft2(Fx,kernelConvSize,kernelConvSize).* fft2(kernel2Mx,kernelConvSize,kernelConvSize)) );
0207             convMx3 = real( ifft2(fft2(Fx,kernelConvSize,kernelConvSize).* fft2(kernel3Mx,kernelConvSize,kernelConvSize)) );
0208             
0209             % Creates an interpolant for kernes from vectors position X and Z
0210             if strcmp(env,'MATLAB')
0211                 Interp_kernel1 = griddedInterpolant(convMx_X',convMx_Z',convMx1','linear','none');
0212                 Interp_kernel2 = griddedInterpolant(convMx_X',convMx_Z',convMx2','linear','none');
0213                 Interp_kernel3 = griddedInterpolant(convMx_X',convMx_Z',convMx3','linear','none');
0214             else
0215                 Interp_kernel1 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx1,x,y,'linear',NaN);
0216                 Interp_kernel2 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx2,x,y,'linear',NaN);
0217                 Interp_kernel3 = @(x,y)interp2(convMx_X(1,:),convMx_Z(:,1),convMx3,x,y,'linear',NaN);
0218             end
0219 
0220         end
0221 
0222         % Display progress and update text only 200 times
0223         if mod(bixelsPerBeam,max(1,round(stf(i).totalNumOfBixels/200))) == 0
0224             matRad_progress(bixelsPerBeam/max(1,round(stf(i).totalNumOfBixels/200)),...
0225                             floor(stf(i).totalNumOfBixels/max(1,round(stf(i).totalNumOfBixels/200))));
0226         end
0227         % update waitbar only 100 times
0228         if mod(counter,round(dij.totalNumOfBixels/100)) == 0 && ishandle(figureWait)
0229             waitbar(counter/dij.totalNumOfBixels);
0230         end
0231         
0232         % remember beam and bixel number
0233         if ~calcDoseDirect
0234            dij.beamNum(counter)  = i;
0235            dij.rayNum(counter)   = j;
0236            dij.bixelNum(counter) = 1;
0237         else
0238             k = 1;
0239         end
0240         
0241         % Ray tracing for beam i and bixel j
0242         [ix,rad_distancesSq,isoLatDistsX,isoLatDistsZ] = matRad_calcGeoDists(rot_coordsVdoseGrid, ...
0243                                                                stf(i).sourcePoint_bev, ...
0244                                                                stf(i).ray(j).targetPoint_bev, ...
0245                                                                machine.meta.SAD, ...
0246                                                                find(~isnan(radDepthVdoseGrid{1})), ...
0247                                                                effectiveLateralCutoff);
0248 
0249         % empty bixels may happen during recalculation of error
0250         % scenarios -> skip to next bixel
0251         if isempty(ix)
0252             continue;
0253         end
0254 
0255                 
0256         % calculate photon dose for beam i and bixel j
0257         bixelDose = matRad_calcPhotonDoseBixel(machine.meta.SAD,machine.data.m,...
0258                                                    machine.data.betas, ...
0259                                                    Interp_kernel1,...
0260                                                    Interp_kernel2,...
0261                                                    Interp_kernel3,...
0262                                                    radDepthVdoseGrid{1}(ix),...
0263                                                    geoDistVdoseGrid{1}(ix),...
0264                                                    isoLatDistsX,...
0265                                                    isoLatDistsZ);
0266                                                
0267         % sample dose only for bixel based dose calculation
0268         if ~isFieldBasedDoseCalc
0269             r0   = 20 + stf(i).bixelWidth;   % [mm] sample beyond the inner core
0270             Type = 'radius';
0271             [ix,bixelDose] = matRad_DijSampling(ix,bixelDose,radDepthVdoseGrid{1}(ix),rad_distancesSq,Type,r0);
0272         end
0273         
0274         % Save dose for every bixel in cell array
0275         doseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1} = sparse(VdoseGrid(ix),1,bixelDose,dij.doseGrid.numOfVoxels,1);
0276 
0277         matRad_calcDoseFillDij;
0278                
0279     end
0280 end
0281 
0282 %Close Waitbar
0283 if ishandle(figureWait)
0284     delete(figureWait);
0285 end
0286