Skip to content

Commit 6a7a48c

Browse files
wahlnwahln
authored
More consistent forward dose calculation and new dose engine tests (#869)
* consistent treatment of weights in pencil beam engine for particles * stf generator does now consider airOffset * updated test dataset * add HongPB dose calculation tests * add photon dose calc test and adapt dose engine for tracking weights in bixel struct * add tests for fine sampling dose engine * disable dij sampling for octave in photon engine tests
1 parent 1308623 commit 6a7a48c

File tree

12 files changed

+348
-69
lines changed

12 files changed

+348
-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/doseCalc/+DoseEngines/matRad_PhotonPencilBeamSVDEngine.m‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -295,6 +295,10 @@ function setDefaults(this)
295295

296296
bixel = struct();
297297

298+
if this.calcDoseDirect
299+
bixel.weight = currRay.weight;
300+
end
301+
298302
if isfield(this.tmpMatrixContainers,'physicalDose')
299303
bixel.physicalDose = this.calcSingleBixel(currRay.SAD,...
300304
this.machine.data.m,...

‎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/doseCalc/test_FSPB.m‎

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
function test_suite = test_FSPB
2+
3+
test_functions=localfunctions();
4+
5+
initTestSuite;
6+
7+
function test_getSubsamplingPBEngineFromPln
8+
% Single gaussian lateral model
9+
testData.pln = struct('radiationMode','protons','machine','Generic');
10+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
11+
engine = DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.getEngineFromPln(testData.pln);
12+
assertTrue(isa(engine,'DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine'));
13+
14+
function test_loadMachineForSubsamplingPB
15+
possibleRadModes = DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.possibleRadiationModes;
16+
for i = 1:numel(possibleRadModes)
17+
machine = DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.loadMachine(possibleRadModes{i},'Generic');
18+
assertTrue(isstruct(machine));
19+
assertTrue(isfield(machine, 'meta'));
20+
assertTrue(isfield(machine.meta, 'radiationMode'));
21+
assertTrue(strcmp(machine.meta.radiationMode, possibleRadModes{i}));
22+
end
23+
24+
25+
function test_calcDoseSubsamplingPBprotonsFitCircle
26+
testData = load('protons_testData.mat');
27+
28+
assertTrue(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
29+
30+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
31+
testData.pln.propDoseCalc.dosimetricLateralCutOff = 0.995;
32+
testData.pln.propDoseCalc.geometricLateralCutOff = 50;
33+
testData.pln.propDoseCalc.fineSampling.method = 'fitCircle';
34+
35+
for N = [2,3,8]
36+
testData.pln.propDoseCalc.fineSampling.N = N;
37+
38+
resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1));
39+
40+
assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
41+
assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose)));
42+
assertElementsAlmostEqual(resultGUI.physicalDose,testData.resultGUI.physicalDose,'relative',1e-2,1e-2);
43+
end
44+
45+
function test_calcDoseSubsamplingPBprotonsFitSquare
46+
testData = load('protons_testData.mat');
47+
48+
assertTrue(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
49+
50+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
51+
testData.pln.propDoseCalc.dosimetricLateralCutOff = 0.995;
52+
testData.pln.propDoseCalc.geometricLateralCutOff = 50;
53+
testData.pln.propDoseCalc.fineSampling.method = 'fitSquare';
54+
55+
for N = [2,3]
56+
testData.pln.propDoseCalc.fineSampling.N = N;
57+
resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1));
58+
59+
assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
60+
assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose)));
61+
assertElementsAlmostEqual(resultGUI.physicalDose,testData.resultGUI.physicalDose,'relative',1e-2,1e-2);
62+
end
63+
64+
function test_calcDoseSubsamplingPBprotonsRusso
65+
testData = load('protons_testData.mat');
66+
67+
testData.pln.propDoseCalc.fineSampling.N = 10;
68+
testData.pln.propDoseCalc.fineSampling.sigmaSub = 2;
69+
testData.pln.propDoseCalc.fineSampling.method = 'russo';
70+
71+
resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1));
72+
73+
assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
74+
assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose)));
75+
assertElementsAlmostEqual(resultGUI.physicalDose,testData.resultGUI.physicalDose,'relative',1e-2,1e-2);
76+
77+
function test_calcDoseSubsamplingPBhelium
78+
testData = load('helium_testData.mat');
79+
assertTrue(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
80+
81+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
82+
testData.pln.propDoseCalc.dosimetricLateralCutOff = 0.995;
83+
testData.pln.propDoseCalc.geometricLateralCutOff = 50;
84+
testData.pln.propDoseCalc.fineSampling.N = 10;
85+
testData.pln.propDoseCalc.fineSampling.sigmaSub = 2;
86+
testData.pln.propDoseCalc.fineSampling.method = 'russo';
87+
resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1));
88+
89+
assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
90+
assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose)));
91+
assertElementsAlmostEqual(resultGUI.physicalDose,testData.resultGUI.physicalDose,'relative',1e-2,1e-2);
92+
93+
function test_calcDoseSubsamplingPBcarbon
94+
testData = load('carbon_testData.mat');
95+
assertTrue(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
96+
97+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
98+
testData.pln.propDoseCalc.dosimetricLateralCutOff = 0.995;
99+
testData.pln.propDoseCalc.geometricLateralCutOff = 50;
100+
testData.pln.propDoseCalc.fineSampling.N = 10;
101+
testData.pln.propDoseCalc.fineSampling.sigmaSub = 2;
102+
testData.pln.propDoseCalc.fineSampling.method = 'russo';
103+
resultGUI = matRad_calcDoseForward(testData.ct, testData.cst, testData.stf, testData.pln, ones(sum([testData.stf(:).totalNumOfBixels]),1));
104+
105+
assertTrue(isequal(fieldnames(resultGUI),fieldnames(testData.resultGUI)));
106+
assertTrue(isequal(testData.ct.cubeDim, size(resultGUI.physicalDose)));
107+
assertElementsAlmostEqual(resultGUI.physicalDose,testData.resultGUI.physicalDose,'relative',1e-2,1e-2);
108+
109+
110+
function test_nonSupportedSettings
111+
% Radiation mode other than protons not implemented
112+
testData = load('photons_testData.mat');
113+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
114+
assertFalse(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
115+
116+
% Invalid machine without radiation mode field
117+
testData.pln.machine = 'Empty';
118+
testData.pln.propDoseCalc.engine = 'SubsamplingPB';
119+
assertExceptionThrown(@() DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln));
120+
assertFalse(DoseEngines.matRad_ParticleFineSamplingPencilBeamEngine.isAvailable(testData.pln,[]));
121+
122+
123+
124+

0 commit comments

Comments
 (0)