function [ machine ] = matRad_calcLateralParticleCutOff(machine,cutOffLevel,stf,visBool)
0001 function [ machine ] = matRad_calcLateralParticleCutOff(machine,cutOffLevel,stf,visBool)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
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
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
0049 vX = [0 logspace(-1,3,1200)];
0050
0051
0052 r_mid = 0.5*(vX(1:end-1) + vX(2:end))';
0053 dr = (vX(2:end) - vX(1:end-1))';
0054 radialDist_sq = r_mid.^2;
0055
0056
0057 numDepthVal = 35;
0058
0059
0060 round2 = @(a,b)round(a*10^b)/10^b;
0061
0062
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
0071 [energySigmaLUT,ixUnique] = unique([[stf.ray(:).energy]; [stf.ray(:).focusIx] ; vSSD]','rows');
0072 rangeShifterLUT = [stf.ray(:).rangeShifter];
0073 rangeShifterLUT = rangeShifterLUT(1,ixUnique);
0074
0075
0076 for i = 1:size(energySigmaLUT,1)
0077
0078
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
0088 if strcmp(machine.meta.radiationMode,'protons') && rangeShifterLUT(i).eqThickness > 0 && ~strcmp(machine.meta.machine,'Generic')
0089
0090
0091 sigmaRashi = matRad_calcSigmaRashi(machine.data(energyIx).energy, ...
0092 rangeShifterLUT(i), ...
0093 energySigmaLUT(i,3));
0094
0095
0096 sigmaIni_sq = sigmaIni_sq + sigmaRashi.^2;
0097
0098 end
0099
0100 energySigmaLUT(i,4) = sigmaIni_sq;
0101
0102 end
0103
0104
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
0113 vEnergiesIx = find(ismember([machine.data(:).energy],uniqueEnergies(:,1)));
0114 cnt = 0;
0115
0116
0117 for energyIx = vEnergiesIx
0118
0119
0120 depthDoseCutOff = inf;
0121
0122
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
0134 cumIntEnergy = cumtrapz(machine.data(energyIx).depths,idd_org);
0135
0136 peakTailRelation = 0.5;
0137 numDepthValToPeak = ceil(numDepthVal*peakTailRelation);
0138 numDepthValTail = ceil(numDepthVal*(1-peakTailRelation));
0139 energyStepsToPeak = cumIntEnergy(peakIxOrg)/numDepthValToPeak;
0140 energyStepsTail = (cumIntEnergy(end)-cumIntEnergy(peakIxOrg))/numDepthValTail;
0141
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
0158 baseData = machine.data(energyIx);
0159 baseData.LatCutOff.CompFac = 1;
0160
0161 for j = 1:numel(depthValues)
0162
0163
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
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;
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
0197 if visBool
0198
0199
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
0211 sStep = 0.5;
0212 vLatX = -100 : sStep : 100;
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
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
0231 if rangeShifter.eqThickness > 0 && strcmp(pln.radiationMode,'protons')
0232
0233
0234 sigmaRashi = matRad_calcSigmaRashi(machine.data(energyIx).energy,rangeShifter,maxSSD);
0235
0236
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
0248 r_mid_Cut = (0.5*(vXCut(1:end-1) + vXCut(2:end)))';
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
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
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