Skip to content

Commit 76cb23a

Browse files
wahlnwahln
committed
Merge branch 'dev' into pr/828
# Conflicts: # test/doseCalc/test_HongPB.m
2 parents 4079c5a + 6e2cd3d commit 76cb23a

File tree

13 files changed

+281
-73
lines changed

13 files changed

+281
-73
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_ParticleFREDEngine/matRad_ParticleFREDEngine.m‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -468,12 +468,20 @@ function writeHlut(this,hLutFile)
468468
matRad_cfg = MatRad_Config.instance();
469469

470470
try
471-
[status, cmdOut] = system([DoseEngines.matRad_ParticleFREDEngine.cmdCall,' -vn']);
471+
[status, cmdOut] = system([DoseEngines.matRad_ParticleFREDEngine.cmdCall, ' -vn']);
472472

473473
if status == 0
474-
version = cmdOut(1:end-1);
474+
% Extract the version number using a regular expression
475+
versionMatch = regexp(cmdOut, '\d+\.\d+\.\d+', 'match', 'once');
476+
if ~isempty(versionMatch)
477+
version = versionMatch;
478+
else
479+
matRad_cfg.dispError('Unable to parse FRED version from output. Please check FRED installation.');
480+
version = [];
481+
end
475482
else
476-
matRad_cfg.dispError('Something wrong occured in checking FRED installation. Please check correct FRED installation');
483+
matRad_cfg.dispError('Error occurred while checking FRED installation. Please check FRED installation.');
484+
version = [];
477485
end
478486
catch
479487
matRad_cfg.dispWarning('Something wrong occured in checking FRED installation. Please check correct FRED installation');

‎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
@@ -176,6 +174,15 @@ function chooseLateralModel(this)
176174
bixel.numParticlesPerMU = currRay.numParticlesPerMU(k);
177175
end
178176

177+
if this.calcDoseDirect
178+
if ~isfield(currRay,'weight') || numel(currRay.weight) < k
179+
matRad_cfg = MatRad_Config.instance();
180+
matRad_cfg.dispError('No weight available for beam %d, ray %d, bixel %d',bixel.beamIndex,bixel.rayIndex,bixel.bixelIndex);
181+
end
182+
bixel.weight = currRay.weight(k);
183+
bixel.MU = (bixel.weight.*1e6) ./ bixel.numParticlesPerMU;
184+
end
185+
179186
% find energy index in base data
180187
energyIx = find(this.round2(currRay.energy(k),4) == this.round2([this.machine.data.energy],4));
181188
bixel.energyIx = energyIx;
@@ -996,7 +1003,9 @@ function calcLateralParticleCutOff(this,cutOffLevel,stfElement)
9961003
dij = this.fillDij@DoseEngines.matRad_PencilBeamEngineAbstract(bixel,dij,stf,scenIdx,currBeamIdx,currRayIdx,currBixelIdx,bixelCounter);
9971004

9981005
% Add MU information
999-
if ~this.calcDoseDirect
1006+
if this.calcDoseDirect
1007+
dij.MU(bixelCounter,1) = bixel.MU;
1008+
else
10001009
dij.minMU(bixelCounter,1) = bixel.minMU;
10011010
dij.maxMU(bixelCounter,1) = bixel.maxMU;
10021011
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
@@ -21,6 +21,7 @@
2121
name = 'Particle IMPT stf Generator';
2222
shortName = 'ParticleIMPT';
2323
possibleRadiationModes = {'protons','helium','carbon'};
24+
airOffsetCorrection = true;
2425
end
2526

2627
properties
@@ -46,6 +47,7 @@
4647

4748
function beam = setBeamletEnergies(this,beam)
4849
%Assigns the max particle machine energy layers to all rays
50+
matRad_cfg = MatRad_Config.instance();
4951

5052
isoCenterInCubeCoords = matRad_world2cubeCoords(beam.isoCenter,this.ct);
5153

@@ -54,11 +56,24 @@
5456
else
5557
LUTspotSize = this.machine.meta.LUTspotSize;
5658
end
57-
59+
60+
%Air Offset Correction
61+
if this.airOffsetCorrection
62+
if ~isfield(this.machine.meta, 'fitAirOffset')
63+
this.machine.meta.fitAirOffset = 0; %By default we assume that the base data was fitted to a phantom with surface at isocenter
64+
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);
65+
end
66+
else
67+
this.machine.meta.fitAirOffset = 0;
68+
end
5869

5970
beam.numOfBixelsPerRay = zeros(1,beam.numOfRays);
6071

6172
for j = beam.numOfRays:-1:1
73+
74+
ctEntryPoint = zeros(this.multScen.totNumShiftScen,1);
75+
76+
radDepthOffset = zeros(this.multScen.totNumShiftScen,1);
6277

6378
for shiftScen = 1:this.multScen.totNumShiftScen
6479
% ray tracing necessary to determine depth of the target
@@ -69,7 +84,14 @@
6984
[this.ct.cube {this.voiTarget}]);
7085

7186
%Used for generic range-shifter placement
72-
ctEntryPoint = alphas(1) * d12;
87+
ctEntryPoint(shiftScen) = alphas(1) * d12;
88+
89+
if this.airOffsetCorrection
90+
nozzleToSkin = (ctEntryPoint(shiftScen) + this.machine.meta.BAMStoIsoDist) - this.machine.meta.SAD;
91+
radDepthOffset(shiftScen) = 0.0011 * (nozzleToSkin - this.machine.meta.fitAirOffset);
92+
else
93+
radDepthOffset(shiftScen) = 0;
94+
end
7395
end
7496

7597
% target hit
@@ -92,8 +114,8 @@
92114

93115
% compute radiological depths
94116
% http://www.ncbi.nlm.nih.gov/pubmed/4000088, eq 14
95-
rSP = l{shiftScen} .* rho{shiftScen}{ctScen};
96-
radDepths = cumsum(rSP) - 0.5*rSP;
117+
rSP = l{shiftScen} .* rho{shiftScen}{ctScen} ;
118+
radDepths = cumsum(rSP) - 0.5*rSP + radDepthOffset(shiftScen);
97119

98120
if this.multScen.relRangeShift(rangeShiftScen) ~= 0 || this.multScen.absRangeShift(rangeShiftScen) ~= 0
99121
radDepths = radDepths +... % original cube
@@ -156,7 +178,7 @@
156178

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

161183
beam.ray(j).energy = [beam.ray(j).energy raShiEnergies];
162184
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+

‎test/doseCalc/test_HongPB.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323

2424
function test_calcDoseHongPBprotons
2525
testData = load('protons_testData.mat');
26-
26+
2727
assertTrue(DoseEngines.matRad_ParticleHongPencilBeamEngine.isAvailable(testData.pln));
2828

2929
testData.pln.propDoseCalc.engine = 'HongPB';

0 commit comments

Comments
 (0)