matRad_calcLateralParticleCutOff

Purpose ^

matRad function to calculate a depth dependend lateral cutoff

Synopsis ^

function [ machine ] = matRad_calcLateralParticleCutOff(machine,cutOffLevel,stf,visBool)

Description ^

 matRad function to calculate a depth dependend lateral cutoff 
 for each pristine particle beam
 
 call
   [ machine ] = matRad_calcLateralParticleCutOff( machine,cutOffLevel,stf,visBool )

 input
   machine:        machine base data file
   cutOffLevel:    cut off level - number between 0 and 1
   stf:              matRad steering information struct
   visBool:         toggle visualization (optional)

 output
   machine:        machine base data file including an additional field representing the lateral
                    cutoff
   
 References
   -

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

 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.

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

Cross-reference information ^

This function calls: This function is called by:

Source code ^

0001 function [ machine ] = matRad_calcLateralParticleCutOff(machine,cutOffLevel,stf,visBool)
0002 % matRad function to calculate a depth dependend lateral cutoff
0003 % for each pristine particle beam
0004 %
0005 % call
0006 %   [ machine ] = matRad_calcLateralParticleCutOff( machine,cutOffLevel,stf,visBool )
0007 %
0008 % input
0009 %   machine:        machine base data file
0010 %   cutOffLevel:    cut off level - number between 0 and 1
0011 %   stf:              matRad steering information struct
0012 %   visBool:         toggle visualization (optional)
0013 %
0014 % output
0015 %   machine:        machine base data file including an additional field representing the lateral
0016 %                    cutoff
0017 %
0018 % References
0019 %   -
0020 %
0021 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 %
0023 % Copyright 2015 the matRad development team.
0024 %
0025 % This file is part of the matRad project. It is subject to the license
0026 % terms in the LICENSE file found in the top-level directory of this
0027 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0028 % of the matRad project, including this file, may be copied, modified,
0029 % propagated, or distributed except according to the terms contained in the
0030 % LICENSE file.
0031 %
0032 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 
0034 if cutOffLevel <= 0.98
0035    warning('a lateral cut off below 0.98 may result in an inaccurate dose calculation') 
0036 end
0037 
0038 conversionFactor = 1.6021766208e-02;
0039 
0040 % function handle for calculating depth dose for APM
0041 sumGauss = @(x,mu,SqSigma,w) ((1./sqrt(2*pi*ones(numel(x),1) * SqSigma') .* ...
0042                               exp(-bsxfun(@minus,x,mu').^2 ./ (2* ones(numel(x),1) * SqSigma' ))) * w);
0043                           
0044 if (cutOffLevel < 0 || cutOffLevel > 1)
0045    warning('lateral cutoff is out of range - using default cut off of 0.99') 
0046    cutOffLevel = 0.99;
0047 end
0048 % define some variables needed for the cutoff calculation
0049 vX = [0 logspace(-1,3,1200)]; % [mm]
0050 
0051 % integration steps
0052 r_mid          = 0.5*(vX(1:end-1) +  vX(2:end))'; % [mm]
0053 dr             = (vX(2:end) - vX(1:end-1))';
0054 radialDist_sq  = r_mid.^2;
0055 
0056 % number of depth points for which a lateral cutoff is determined
0057 numDepthVal    = 35; 
0058 
0059 % helper function for energy selection
0060 round2 = @(a,b)round(a*10^b)/10^b;
0061 
0062 % extract SSD for each bixel
0063 vSSD = ones(1,length([stf.ray(:).energy]));
0064 cnt = 1;
0065 for i  = 1:length(stf.ray)
0066     vSSD(cnt:cnt+numel([stf.ray(i).energy])-1) = stf.ray(i).SSD;
0067     cnt = cnt + numel(stf.ray(i).energy);
0068 end
0069 
0070 % setup energy, focus index, sigma look up table - only consider unique rows
0071 [energySigmaLUT,ixUnique]  = unique([[stf.ray(:).energy]; [stf.ray(:).focusIx] ; vSSD]','rows');
0072 rangeShifterLUT = [stf.ray(:).rangeShifter];
0073 rangeShifterLUT = rangeShifterLUT(1,ixUnique);
0074 
0075 % find the largest inital beam width considering focus index, SSD and range shifter for each individual energy
0076 for i = 1:size(energySigmaLUT,1)
0077    
0078     % find index of maximum used energy (round to keV for numerical reasons
0079     energyIx = max(round2(energySigmaLUT(i,1),4)) == round2([machine.data.energy],4);
0080     
0081     currFoci = energySigmaLUT(i,2);
0082     sigmaIni = matRad_interp1(machine.data(energyIx).initFocus.dist(currFoci,:)',...
0083                               machine.data(energyIx).initFocus.sigma(currFoci,:)',...
0084                               energySigmaLUT(i,3));
0085     sigmaIni_sq = sigmaIni^2;
0086     
0087     % consider range shifter for protons if applicable
0088     if  strcmp(machine.meta.radiationMode,'protons') && rangeShifterLUT(i).eqThickness > 0  && ~strcmp(machine.meta.machine,'Generic')
0089 
0090         %get max range shift
0091         sigmaRashi = matRad_calcSigmaRashi(machine.data(energyIx).energy, ...
0092                                            rangeShifterLUT(i), ...
0093                                            energySigmaLUT(i,3));
0094 
0095         % add to initial sigma in quadrature
0096         sigmaIni_sq = sigmaIni_sq +  sigmaRashi.^2;
0097 
0098     end                          
0099                                                          
0100     energySigmaLUT(i,4) = sigmaIni_sq;
0101     
0102 end
0103 
0104 % find for each individual energy the broadest inital beam width
0105 uniqueEnergies                = unique(energySigmaLUT(:,1));
0106 largestSigmaSq4uniqueEnergies = NaN * ones(numel(uniqueEnergies),1);
0107 ix_Max                        = NaN * ones(numel(uniqueEnergies),1);
0108 for i = 1:numel(uniqueEnergies)
0109     [largestSigmaSq4uniqueEnergies(i), ix_Max(i)] = max(energySigmaLUT(uniqueEnergies(i) == energySigmaLUT(:,1),4));
0110 end
0111 
0112 % get energy indices for looping
0113 vEnergiesIx = find(ismember([machine.data(:).energy],uniqueEnergies(:,1)));
0114 cnt         = 0;    
0115 
0116 % loop over all entries in the machine.data struct
0117 for energyIx = vEnergiesIx
0118    
0119     % set default depth cut off - finite value will be set during first iteration
0120     depthDoseCutOff = inf;
0121 
0122     % get the current integrated depth dose profile
0123     if isstruct(machine.data(energyIx).Z)
0124         idd_org = sumGauss(machine.data(energyIx).depths,machine.data(energyIx).Z.mean,...
0125                                    machine.data(energyIx).Z.width.^2,...
0126                                    machine.data(energyIx).Z.weight) * conversionFactor;
0127     else
0128         idd_org = machine.data(energyIx).Z * conversionFactor;
0129     end
0130     
0131     [~,peakIxOrg] = max(idd_org); 
0132     
0133     % get indices for which a lateral cutoff should be calculated
0134     cumIntEnergy = cumtrapz(machine.data(energyIx).depths,idd_org);
0135     
0136     peakTailRelation   = 0.5;
0137     numDepthValToPeak  = ceil(numDepthVal*peakTailRelation);                                                                          % number of depth values from 0 to peak position
0138     numDepthValTail    = ceil(numDepthVal*(1-peakTailRelation));                                                                      % number of depth values behind peak position
0139     energyStepsToPeak  = cumIntEnergy(peakIxOrg)/numDepthValToPeak;
0140     energyStepsTail    = (cumIntEnergy(end)-cumIntEnergy(peakIxOrg))/numDepthValTail;
0141     % make sure to include 0, peak position and end position
0142     vEnergySteps       = unique([0:energyStepsToPeak:cumIntEnergy(peakIxOrg) cumIntEnergy(peakIxOrg) ...
0143                                  cumIntEnergy(peakIxOrg+1):energyStepsTail:cumIntEnergy(end) cumIntEnergy(end)]);
0144 
0145     [cumIntEnergy,ix] = unique(cumIntEnergy);
0146     depthValues       = matRad_interp1(cumIntEnergy,machine.data(energyIx).depths(ix),vEnergySteps);
0147 
0148     if isstruct(machine.data(energyIx).Z)
0149         idd = sumGauss(depthValues,machine.data(energyIx).Z.mean,...
0150                                    machine.data(energyIx).Z.width.^2,...
0151                                    machine.data(energyIx).Z.weight) * conversionFactor;
0152     else
0153         idd  = matRad_interp1(machine.data(energyIx).depths,machine.data(energyIx).Z,depthValues) * conversionFactor; 
0154     end
0155      
0156     cnt = cnt +1 ;
0157     % % calculate dose in spot
0158     baseData                   = machine.data(energyIx);
0159     baseData.LatCutOff.CompFac = 1;   
0160   
0161     for j = 1:numel(depthValues)
0162         
0163         % save depth value
0164         machine.data(energyIx).LatCutOff.depths(j) = depthValues(j);
0165         
0166         if cutOffLevel == 1
0167             machine.data(energyIx).LatCutOff.CompFac   = 1;
0168             machine.data(energyIx).LatCutOff.CutOff(j) = Inf;
0169         else
0170         
0171             % calculate dose
0172             dose_r = matRad_calcParticleDoseBixel(depthValues(j) + baseData.offset, radialDist_sq, largestSigmaSq4uniqueEnergies(cnt), baseData);
0173 
0174             cumArea = cumsum(2*pi.*r_mid.*dose_r.*dr);
0175             relativeTolerance = 0.5; %in [%]
0176             if abs((cumArea(end)./(idd(j)))-1)*100 > relativeTolerance
0177                 warning('LateralParticleCutOff: shell integration is wrong !')
0178             end
0179 
0180             IX = find(cumArea >= idd(j) * cutOffLevel,1, 'first'); 
0181             machine.data(energyIx).LatCutOff.CompFac = cutOffLevel^-1;
0182 
0183             if isempty(IX)
0184                 depthDoseCutOff = Inf;
0185                 warning('LateralParticleCutOff: Couldnt find lateral cut off !')
0186             elseif isnumeric(IX)
0187                 depthDoseCutOff = r_mid(IX);
0188             end
0189 
0190             machine.data(energyIx).LatCutOff.CutOff(j) = depthDoseCutOff;
0191 
0192         end
0193     end    
0194 end    
0195           
0196 %% visualization
0197 if visBool
0198     
0199     % determine which pencil beam should be plotted
0200     subIx    = ceil(numel(vEnergiesIx)/2);
0201     energyIx = vEnergiesIx(subIx);
0202     
0203     baseData       = machine.data(energyIx);
0204     focusIx        = energySigmaLUT(ix_Max(subIx),2);
0205     maxSSD         = energySigmaLUT(ix_Max(subIx),3);
0206     rangeShifter   = rangeShifterLUT(ix_Max(subIx));
0207     TmpCompFac     = baseData.LatCutOff.CompFac;
0208     baseData.LatCutOff.CompFac = 1;
0209     
0210     % plot 3D cutoff at one specific depth on a rather sparse grid
0211     sStep         = 0.5;
0212     vLatX         = -100 : sStep : 100; % [mm]
0213     dimX          = numel(vLatX);
0214     midPos        = round(length(vLatX)/2);
0215     [X,Y]         = meshgrid(vLatX,vLatX);
0216     
0217     radDepths     = [0:sStep:machine.data(energyIx).depths(end)] + machine.data(energyIx).offset;
0218     radialDist_sq = (X.^2 + Y.^2);
0219     radialDist_sq = radialDist_sq(:);
0220     mDose         = zeros(dimX,dimX,numel(radDepths));
0221     vDoseInt      = zeros(numel(radDepths),1);
0222     
0223     for kk = 1:numel(radDepths)    
0224           
0225          % calculate initial focus sigma
0226          sigmaIni = matRad_interp1(machine.data(energyIx).initFocus.dist(focusIx,:)', ...
0227                                    machine.data(energyIx).initFocus.sigma(focusIx,:)',maxSSD);
0228          sigmaIni_sq = sigmaIni^2;
0229 
0230          % consider range shifter for protons if applicable
0231          if rangeShifter.eqThickness > 0 && strcmp(pln.radiationMode,'protons')
0232 
0233               % compute!
0234               sigmaRashi = matRad_calcSigmaRashi(machine.data(energyIx).energy,rangeShifter,maxSSD);
0235 
0236               % add to initial sigma in quadrature
0237               sigmaIni_sq = sigmaIni_sq +  sigmaRashi^2;
0238 
0239          end
0240 
0241          mDose(:,:,kk) = reshape(matRad_calcParticleDoseBixel(radDepths(kk), radialDist_sq, sigmaIni_sq,baseData),[dimX dimX]);
0242           
0243          [~,IX]           = min(abs((machine.data(energyIx).LatCutOff.depths + machine.data(energyIx).offset) - radDepths(kk)));
0244          TmpCutOff        = machine.data(energyIx).LatCutOff.CutOff(IX);    
0245          vXCut            = vX(vX<=TmpCutOff);
0246          
0247          % integration steps
0248          r_mid_Cut        = (0.5*(vXCut(1:end-1) +  vXCut(2:end)))'; % [mm]
0249          dr_Cut           = (vXCut(2:end) - vXCut(1:end-1))';
0250          radialDist_sqCut = r_mid_Cut.^2;    
0251          
0252          dose_r_Cut       = matRad_calcParticleDoseBixel(radDepths(kk), radialDist_sqCut(:), sigmaIni_sq,baseData);
0253          
0254          cumAreaCut = cumsum(2*pi.*r_mid_Cut.*dose_r_Cut.*dr_Cut);  
0255          
0256          if ~isempty(cumAreaCut)
0257              vDoseInt(kk) = cumAreaCut(end);
0258          end
0259     end
0260     
0261     % obtain maximum dose
0262     if isstruct(machine.data(energyIx).Z)
0263         idd = sumGauss(depthValues,machine.data(energyIx).Z.mean,...
0264                                    machine.data(energyIx).Z.width.^2,...
0265                                    machine.data(energyIx).Z.weight) * conversionFactor;
0266     else
0267         idd  = matRad_interp1(machine.data(energyIx).depths,machine.data(energyIx).Z,depthValues) * conversionFactor; 
0268     end
0269     
0270     [~,peakixDepth] = max(idd); 
0271     dosePeakPos = matRad_calcParticleDoseBixel(machine.data(energyIx).depths(peakixDepth), 0, sigmaIni_sq, baseData);   
0272     
0273     vLevelsDose = dosePeakPos.*[0.01 0.05 0.1 0.9];
0274     doseSlice   = squeeze(mDose(midPos,:,:));
0275     figure,set(gcf,'Color',[1 1 1]);
0276     subplot(311),h=imagesc(squeeze(mDose(midPos,:,:)));hold on;
0277     set(h,'AlphaData', .8*double(doseSlice>0));
0278     contour(doseSlice,vLevelsDose,'LevelListMode','manual','LineWidth',2);hold on
0279     
0280     ax = gca;
0281     ax.XTickLabelMode = 'manual';
0282     ax.XTickLabel     = strsplit(num2str(ax.XTick*sStep + machine.data(energyIx).offset),' ')';
0283     ax.YTickLabelMode = 'manual';
0284     ax.YTickLabel     = strsplit(num2str(ax.YTick*sStep + machine.data(energyIx).offset),' ')';
0285   
0286     plot(1+(machine.data(energyIx).LatCutOff.depths)*sStep^-1,...
0287           machine.data(energyIx).LatCutOff.CutOff * sStep^-1 + midPos,'rx');
0288 
0289     legend({'isodose 1%,5%,10% 90%','calculated cutoff'}) ,colorbar,set(gca,'FontSize',12),xlabel('z [mm]'),ylabel('x [mm]');
0290        
0291     entry = machine.data(energyIx);
0292     if isstruct(entry.Z)
0293        idd = sumGauss(entry.depths,entry.Z.mean,entry.Z.width.^2,entry.Z.weight);
0294     else
0295        idd = machine.data(energyIx).Z;
0296     end
0297     subplot(312),plot(machine.data(energyIx).depths,idd*conversionFactor,'k','LineWidth',2),grid on,hold on
0298                  plot(radDepths - machine.data(energyIx).offset,vDoseInt,'r--','LineWidth',2),hold on,
0299                  plot(radDepths - machine.data(energyIx).offset,vDoseInt * TmpCompFac,'bx','LineWidth',1),hold on,
0300     legend({'original IDD',['cut off IDD at ' num2str(cutOffLevel) '%'],'cut off IDD with compensation'},'Location','northwest'),
0301     xlabel('z [mm]'),ylabel('[MeV cm^2 /(g * primary)]'),set(gca,'FontSize',12)     
0302            
0303     totEnergy        = trapz(machine.data(energyIx).depths,idd*conversionFactor) ;
0304     totEnergyCutOff  = trapz(radDepths,vDoseInt * TmpCompFac) ;
0305     relDiff          =  ((totEnergy/totEnergyCutOff)-1)*100;   
0306     title(['rel diff of integral dose ' num2str(relDiff) '%']);
0307     baseData.LatCutOff.CompFac = TmpCompFac;
0308     
0309     subplot(313),
0310     if isfield(machine.data(energyIx),'sigma1')
0311         yyaxis left;
0312         plot(machine.data(energyIx).LatCutOff.depths,machine.data(energyIx).LatCutOff.CutOff,'LineWidth',2),hold on
0313         plot(machine.data(energyIx).depths,(machine.data(energyIx).sigma1),':','LineWidth',2),grid on,hold on,ylabel('mm')
0314         yyaxis right; 
0315         plot(machine.data(energyIx).depths,(machine.data(energyIx).sigma2),'-.','LineWidth',2),grid on,hold on,ylabel('mm')
0316         legend({'Cutoff','sigma1','sigma2'});
0317     else
0318         yyaxis left;plot(machine.data(energyIx).LatCutOff.depths,machine.data(energyIx).LatCutOff.CutOff,'LineWidth',2),hold on,ylabel('mm')
0319         yyaxis right;subplot(313),plot(machine.data(energyIx).depths,machine.data(energyIx).sigma,'LineWidth',2),grid on,hold on
0320         legend({'Cutoff','sigma'});ylabel('mm')
0321     end
0322 
0323     set(gca,'FontSize',12),xlabel('z [mm]'),  ylabel('mm')
0324 
0325     % plot cutoff of different energies
0326     figure,set(gcf,'Color',[1 1 1]);
0327     cnt = 1;
0328     for i = vEnergiesIx
0329         plot(machine.data(i).LatCutOff.depths,machine.data(i).LatCutOff.CutOff,'LineWidth',1.5),hold on
0330         cellLegend{cnt} = [num2str(machine.data(i).energy) ' MeV'];
0331         cnt = cnt + 1;
0332     end
0333     grid on, grid minor,xlabel('depth in [mm]'),ylabel('lateral cutoff in [mm]')
0334     title(['cutoff level = ' num2str(cutOffLevel)]),
0335     ylim = get(gca,'Ylim');    set(gca,'Ylim',[0 ylim(2)+3]),    legend(cellLegend)
0336 end
0337 
0338 
0339 
0340 
0341 end
0342

| Generated by m2html © 2005