Skip to content

Commit 6f7b87e

Browse files
wahlnwahln
committed
* consistent treatment of weights in pencil beam engine for particles
* stf generator does now consider airOffset * updated test dataset
1 parent 1308623 commit 6f7b87e

File tree

8 files changed

+83
-69
lines changed

8 files changed

+83
-69
lines changed

‎matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -318,13 +318,23 @@ function assignBioModelPropertiesFromPln(this, plnModel, warnWhenPropertyChanged
318318
end
319319
end
320320
else
321-
resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));
321+
if this.multScen.totNumScen > 1
322+
resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));
323+
end
322324
end
323325
end
326+
327+
if isfield(dij,'w')
328+
resultGUI.w = dij.w;
329+
else
330+
resultGUI.w = w;
331+
end
324332

333+
if isfield(dij,'MU')
334+
resultGUI.MU = dij.MU;
335+
end
325336

326-
resultGUI.w = w;
327-
337+
resultGUI = orderfields(resultGUI);
328338
end
329339

330340
function dij = calcDoseInfluence(this,ct,cst,stf)

‎matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -50,8 +50,6 @@
5050

5151
this = [email protected]_PencilBeamEngineAbstract(pln);
5252
end
53-
54-
5553
end
5654

5755
% Should be abstract methods but in order to satisfy the compatibility
@@ -164,6 +162,15 @@ function chooseLateralModel(this)
164162
bixel.numParticlesPerMU = currRay.numParticlesPerMU(k);
165163
end
166164

165+
if this.calcDoseDirect
166+
if ~isfield(currRay,'weight') || numel(currRay.weight) < k
167+
matRad_cfg = MatRad_Config.instance();
168+
matRad_cfg.dispError('No weight available for beam %d, ray %d, bixel %d',bixel.beamIndex,bixel.rayIndex,bixel.bixelIndex);
169+
end
170+
bixel.weight = currRay.weight(k);
171+
bixel.MU = (bixel.weight.*1e6) ./ bixel.numParticlesPerMU;
172+
end
173+
167174
% find energy index in base data
168175
energyIx = find(this.round2(currRay.energy(k),4) == this.round2([this.machine.data.energy],4));
169176
bixel.energyIx = energyIx;
@@ -942,7 +949,9 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement)
942949
dij = this.fillDij@DoseEngines.matRad_PencilBeamEngineAbstract(bixel,dij,stf,scenIdx,currBeamIdx,currRayIdx,currBixelIdx,bixelCounter);
943950

944951
% Add MU information
945-
if ~this.calcDoseDirect
952+
if this.calcDoseDirect
953+
dij.MU(bixelCounter,1) = bixel.MU;
954+
else
946955
dij.minMU(bixelCounter,1) = bixel.minMU;
947956
dij.maxMU(bixelCounter,1) = bixel.maxMU;
948957
dij.numParticlesPerMU(bixelCounter,1) = bixel.numParticlesPerMU;

‎matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m‎

Lines changed: 2 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -442,16 +442,9 @@ function setDefaults(this)
442442

443443
dijColIx = (ceil(counter/this.numOfBixelsContainer)-1)*this.numOfBixelsContainer+1:counter;
444444
containerIx = 1:bixelContainerColIx;
445-
weight = 1;
446445
else
447446
dijColIx = currBeamIdx;
448447
containerIx = 1;
449-
if isfield(stf(currBeamIdx).ray(currRayIdx),'weight') && numel(stf(currBeamIdx).ray(currRayIdx).weight) >= currBixelIdx
450-
weight = stf(currBeamIdx).ray(currRayIdx).weight(currBixelIdx);
451-
else
452-
matRad_cfg = MatRad_Config.instance();
453-
matRad_cfg.dispError('No weight available for beam %d, ray %d, bixel %d',currBeamIdx,currRayIdx,currBixelIdx);
454-
end
455448
end
456449

457450
% Iterate through all quantities
@@ -463,7 +456,7 @@ function setDefaults(this)
463456
this.tmpMatrixContainers.(qName)(containerIx,subScenIdx{:}) = cell(numel(containerIx,subScenIdx{:}));
464457
else
465458
%dij.(qName){1}(this.VdoseGrid(bixel.ix),dijColIx) = dij.(qName){1}(this.VdoseGrid(bixel.ix),dijColIx) + weight * this.tmpMatrixContainers.(qName){containerIx,1}(this.VdoseGrid(bixel.ix));
466-
dij.(qName){scenIdx}(bixel.ix,dijColIx) = dij.(qName){scenIdx}(bixel.ix,dijColIx) + weight * bixel.(qName);
459+
dij.(qName){scenIdx}(bixel.ix,dijColIx) = dij.(qName){scenIdx}(bixel.ix,dijColIx) + bixel.weight * bixel.(qName);
467460
end
468461
end
469462
end
@@ -475,35 +468,13 @@ function setDefaults(this)
475468
dij.beamNum(currBeamIdx) = currBeamIdx;
476469
dij.rayNum(currBeamIdx) = currBeamIdx;
477470
dij.bixelNum(currBeamIdx) = currBeamIdx;
471+
dij.w(counter,1) = bixel.weight;
478472
else
479473
dij.beamNum(counter) = currBeamIdx;
480474
dij.rayNum(counter) = currRayIdx;
481475
dij.bixelNum(counter) = currBixelIdx;
482476
end
483477
end
484-
485-
%{
486-
function ray = computeRaySSD(this,ray)
487-
[alpha,~,rho,d12,~] = matRad_siddonRayTracer(ray.isoCenter, ...
488-
ct.resolution, ...
489-
ray.sourcePoint, ...
490-
ray.targetPoint, ...
491-
this.cubeWED(1));
492-
ixSSD = find(rho{1} > this.ssdDensityThreshold,1,'first');
493-
494-
495-
if isempty(ixSSD)
496-
matRad_cfg.dispError('ray does not hit patient. Trying to fix afterwards...');
497-
boolShowWarning = false;
498-
elseif ixSSD(1) == 1
499-
matRad_cfg.dispWarning('Surface for SSD calculation starts directly in first voxel of CT!');
500-
boolShowWarning = false;
501-
end
502-
503-
% calculate SSD
504-
ray.SSD = double(d12* alpha(ixSSD));
505-
end
506-
%}
507478

508479
function dij = finalizeDose(this,dij)
509480
%TODO: We could also do this by default for all engines, but

‎matRad/steering/matRad_StfGeneratorParticleIMPT.m‎

Lines changed: 27 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
name = 'Particle IMPT stf Generator';
2020
shortName = 'ParticleIMPT';
2121
possibleRadiationModes = {'protons','helium','carbon'};
22+
airOffsetCorrection = true;
2223
end
2324

2425
properties
@@ -44,6 +45,7 @@
4445

4546
function beam = setBeamletEnergies(this,beam)
4647
%Assigns the max particle machine energy layers to all rays
48+
matRad_cfg = MatRad_Config.instance();
4749

4850
isoCenterInCubeCoords = matRad_world2cubeCoords(beam.isoCenter,this.ct);
4951

@@ -52,11 +54,24 @@
5254
else
5355
LUTspotSize = this.machine.meta.LUTspotSize;
5456
end
55-
57+
58+
%Air Offset Correction
59+
if this.airOffsetCorrection
60+
if ~isfield(this.machine.meta, 'fitAirOffset')
61+
this.machine.meta.fitAirOffset = 0; %By default we assume that the base data was fitted to a phantom with surface at isocenter
62+
matRad_cfg.dispDebug('Asked for correction of Base Data Air Offset, but no value found. Using default value of %f mm.\n',this.machine.meta.fitAirOffset);
63+
end
64+
else
65+
this.machine.meta.fitAirOffset = 0;
66+
end
5667

5768
beam.numOfBixelsPerRay = zeros(1,beam.numOfRays);
5869

5970
for j = beam.numOfRays:-1:1
71+
72+
ctEntryPoint = zeros(this.multScen.totNumShiftScen,1);
73+
74+
radDepthOffset = zeros(this.multScen.totNumShiftScen,1);
6075

6176
for shiftScen = 1:this.multScen.totNumShiftScen
6277
% ray tracing necessary to determine depth of the target
@@ -67,7 +82,14 @@
6782
[this.ct.cube {this.voiTarget}]);
6883

6984
%Used for generic range-shifter placement
70-
ctEntryPoint = alphas(1) * d12;
85+
ctEntryPoint(shiftScen) = alphas(1) * d12;
86+
87+
if this.airOffsetCorrection
88+
nozzleToSkin = (ctEntryPoint(shiftScen) + this.machine.meta.BAMStoIsoDist) - this.machine.meta.SAD;
89+
radDepthOffset(shiftScen) = 0.0011 * (nozzleToSkin - this.machine.meta.fitAirOffset);
90+
else
91+
radDepthOffset(shiftScen) = 0;
92+
end
7193
end
7294

7395
% target hit
@@ -90,8 +112,8 @@
90112

91113
% compute radiological depths
92114
% http://www.ncbi.nlm.nih.gov/pubmed/4000088, eq 14
93-
rSP = l{shiftScen} .* rho{shiftScen}{ctScen};
94-
radDepths = cumsum(rSP) - 0.5*rSP;
115+
rSP = l{shiftScen} .* rho{shiftScen}{ctScen} ;
116+
radDepths = cumsum(rSP) - 0.5*rSP + radDepthOffset(shiftScen);
95117

96118
if this.multScen.relRangeShift(rangeShiftScen) ~= 0 || this.multScen.absRangeShift(rangeShiftScen) ~= 0
97119
radDepths = radDepths +... % original cube
@@ -154,7 +176,7 @@
154176

155177
raShi.ID = 1;
156178
raShi.eqThickness = rangeShifterEqD;
157-
raShi.sourceRashiDistance = round(ctEntryPoint - 2*rangeShifterEqD,-1); %place a little away from entry, round to cms to reduce number of unique settings
179+
raShi.sourceRashiDistance = round(min(ctEntryPoint) - 2*rangeShifterEqD,-1); %place a little away from entry, round to cms to reduce number of unique settings
158180

159181
beam.ray(j).energy = [beam.ray(j).energy raShiEnergies];
160182
beam.ray(j).rangeShifter = [beam.ray(j).rangeShifter repmat(raShi,1,length(raShiEnergies))];

‎test/testData/carbon_testData.mat‎

-3.15 KB
Binary file not shown.

‎test/testData/helium_testData.mat‎

-1.98 KB
Binary file not shown.

‎test/testData/helper_testDataCreater.m‎

Lines changed: 29 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@
22
% therapy, that can be red in and used in testing
33

44
%% create ct
5-
65
ct = struct();
76
ct.cubeDim = [20,10,10];
87
ct.resolution.x = 10;
@@ -40,29 +39,32 @@
4039

4140
clear VolHelper ixBody ixTarget i
4241
%% create pln, stf
43-
radMode = 'carbon'; %protons,helium,carbon;
44-
45-
pln.radiationMode = radMode;
46-
pln.machine = 'Generic';
47-
pln.numOfFractions = 30;
48-
pln.propStf.gantryAngles = [0,180];
49-
pln.propStf.couchAngles = zeros(size(pln.propStf.gantryAngles));
50-
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
51-
pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
52-
53-
pln.propStf.longitudinalSpotSpacing = 8;
54-
pln.propStf.bixelWidth = 10;
55-
pln.propDoseCalc.doseGrid.resolution = struct('x',10,'y',10,'z',10); %[mm]
56-
57-
%pln.bioModel = matRad_bioModel(pln.radiationMode,'none');
58-
59-
%% Generate Beam Geometry STF
60-
pln.propStf.addMargin = false; %to make smaller stf, les bixel
61-
stf = matRad_generateStf(ct,cst,pln);
62-
ct = matRad_calcWaterEqD(ct, pln);
63-
%% Dose Calculation
64-
%dij = matRad_calcDoseInfluence(ct,cst,stf,pln);
65-
%resultGUI = matRad_calcCubes(ones(dij.totalNumOfBixels,1),dij);
66-
67-
%% save basic data
68-
save([radMode '_testData.mat'],'ct','cst','pln','stf','-v7')
42+
radModes = ["protons","helium","carbon","VHEE"];
43+
for radMode = radModes
44+
%radMode = 'carbon'; %protons,helium,carbon;
45+
46+
pln.radiationMode = char(radMode);
47+
pln.machine = 'Generic';
48+
pln.numOfFractions = 30;
49+
pln.propStf.gantryAngles = [0,180];
50+
pln.propStf.couchAngles = zeros(size(pln.propStf.gantryAngles));
51+
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
52+
pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
53+
54+
pln.propStf.longitudinalSpotSpacing = 8;
55+
pln.propStf.bixelWidth = 10;
56+
pln.propDoseCalc.doseGrid.resolution = struct('x',10,'y',10,'z',10); %[mm]
57+
58+
%pln.bioModel = matRad_bioModel(pln.radiationMode,'none');
59+
60+
%% Generate Beam Geometry STF
61+
pln.propStf.addMargin = false; %to make smaller stf, les bixel
62+
stf = matRad_generateStf(ct,cst,pln);
63+
ct = matRad_calcWaterEqD(ct, pln);
64+
%% Dose Calculation
65+
dij = matRad_calcDoseInfluence(ct,cst,stf,pln);
66+
resultGUI = matRad_calcCubes(ones(dij.totalNumOfBixels,1),dij);
67+
68+
%% save basic data
69+
save([char(radMode) '_testData.mat'],'ct','cst','pln','stf','dij','resultGUI','-v7');
70+
end

‎test/testData/protons_testData.mat‎

12.1 KB
Binary file not shown.

0 commit comments

Comments
 (0)