Skip to content

Commit 96c3c0c

Browse files
wahlnwahln
authored
Merge branch 'dev' into feature/vhee
2 parents 2ab1746 + 1308623 commit 96c3c0c

File tree

11 files changed

+319
-178
lines changed

11 files changed

+319
-178
lines changed

‎matRad/bioModels/LQbasedModels/kernelBasedModels/matRad_LQKernelBasedModel.m‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,4 +80,22 @@
8080

8181
end
8282
end
83+
84+
methods (Static)
85+
86+
87+
function [alphaX, betaX] = getAvailableTissueParameters(pln)
88+
89+
% load machine
90+
machine = matRad_loadMachine(pln);
91+
if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX')
92+
alphaX = machine.data(1).alphaX;
93+
betaX = machine.data(1).betaX;
94+
else
95+
matRad_cfg = MatRad_Config.instance();
96+
matRad_cfg.dispError('The selected biological model requires AlphaX and BetaX to be set in the machine file but none was found.');
97+
end
98+
99+
end
100+
end
83101
end

‎matRad/bioModels/matRad_BiologicalModel.m‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -251,6 +251,13 @@
251251
end
252252
end
253253

254+
function [alphaX, betaX] = getAvailableTissueParameters(pln)
255+
% empty values in standard implementation, needs to be
256+
% overwritten in subclasses
257+
alphaX = [];
258+
betaX = [];
259+
end
260+
254261

255262
end
256263

‎matRad/bioModels/matRad_bioModel.m‎

Lines changed: 10 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,26 @@
1-
function model = matRad_bioModel(sRadiationMode, sModel)
1+
function model = matRad_bioModel(radiationMode, model, providedQuantities)
22
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33
% matRad_bioModel
44
% This is a helper function to instantiate a matRad_BiologicalModel. This
55
% function currently exists for downwards compatability, as the new
66
% Biological Models will follow a polymorphic software architecture
77
%
88
% call
9-
% matRad_bioModel(sRadiationMode, sModel)
9+
% matRad_bioModel(radiationMode, model)
1010
%
1111
% e.g. pln.bioModel = matRad_bioModel('protons','MCN')
1212
%
1313
% input
14-
% sRadiationMode: radiation modality 'photons' 'protons' 'helium' 'carbon' 'brachy'
14+
% radiationMode: radiation modality 'photons' 'protons' 'helium' 'carbon' 'brachy'
1515
%
16-
% sModel: string to denote which biological model is used
16+
% model: string to denote which biological model is used
1717
% 'none': for photons, protons, carbon 'constRBE': constant RBE for photons and protons
1818
% 'MCN': McNamara-variable RBE model for protons 'WED': Wedenberg-variable RBE model for protons
1919
% 'LEM': Local Effect Model for carbon ions
2020
%
21+
% providedQuantities: optional cell string of provided quantities to
22+
% check if the model can be evaluated
23+
%
2124
% output
2225
% model: instance of a biological model
2326
%
@@ -36,51 +39,10 @@
3639
%
3740
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3841

39-
matRad_cfg = MatRad_Config.instance();
40-
41-
% Look for the correct inputs
42-
p = inputParser;
43-
addRequired(p, 'sRadiationMode', @ischar);
44-
addRequired(p, 'sModel',@ischar);
45-
46-
p.KeepUnmatched = true;
47-
48-
%Check for the available models
49-
mainFolder = fullfile(matRad_cfg.matRadSrcRoot,'bioModels');
50-
userDefinedFolder = fullfile(matRad_cfg.primaryUserFolder, 'bioModels');
51-
52-
if ~exist(userDefinedFolder,"dir")
53-
folders = {mainFolder};
42+
if nargin < 3
43+
model = matRad_BiologicalModel.validate(model,radiationMode);
5444
else
55-
folders = {mainFolder,userDefinedFolder};
56-
end
57-
58-
availableBioModelsClassList = matRad_findSubclasses('matRad_BiologicalModel', 'folders', folders , 'includeSubfolders',true);
59-
modelInfos = matRad_identifyClassesByConstantProperties(availableBioModelsClassList,'model');
60-
modelNames = {modelInfos.model};
61-
62-
if numel(unique({modelInfos.model})) ~= numel(modelInfos)
63-
matRad_cfg.dispError('Multiple biological models with the same name available.');
45+
model = matRad_BiologicalModel.validate(model,radiationMode,providedQuantities);
6446
end
65-
66-
selectedModelIdx = find(strcmp(sModel, modelNames));
67-
68-
% Create first instance of the selected model
69-
if ~isempty(selectedModelIdx)
70-
tmpBioParam = modelInfos(selectedModelIdx).handle();
71-
else
72-
matRad_cfg.dispError('Unrecognized biological model: %s', sModel);
73-
end
74-
75-
% For the time being I do not assigne the model specific parameters, they
76-
% can be assigned by the user later
77-
78-
correctRadiationModality = any(strcmp(tmpBioParam.possibleRadiationModes, sRadiationMode));
79-
80-
if ~correctRadiationModality
81-
matRad_cfg.dispError('Incorrect radiation modality for the required biological model');
82-
end
83-
84-
model = tmpBioParam;
8547

8648
end % end function

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

Lines changed: 103 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -72,107 +72,8 @@
7272
end
7373
end
7474

75-
methods (Access = protected)
76-
function dij = initDoseCalc(this,ct,cst,stf)
77-
78-
matRad_cfg = MatRad_Config.instance();
79-
80-
if this.calcLET == true
81-
matRad_cfg.dispWarning('Engine does not support LET calculation! Disabling!');
82-
this.calcLET = false;
83-
end
84-
85-
if this.calcBioDose == true
86-
matRad_cfg.dispWarning('Engine does not support BioDose calculation! Disabling!');
87-
this.calcBioDose = false;
88-
end
89-
90-
dij = this.initDoseCalc@DoseEngines.matRad_ParticlePencilBeamEngineAbstract(ct,cst,stf);
91-
end
92-
93-
function chooseLateralModel(this)
94-
%Now check if we need tho chose the lateral model because it
95-
%was set to auto
96-
matRad_cfg = MatRad_Config.instance();
97-
if strcmp(this.lateralModel,'auto')
98-
this.lateralModel = 'single';
99-
elseif ~strcmp(this.lateralModel,'single')
100-
matRad_cfg.dispWarning('Engine only supports analytically computed singleGaussian lateral Model!');
101-
this.lateralModel = 'single';
102-
end
103-
matRad_cfg.dispInfo('Using an analytically computed %s Gaussian pencil-beam kernel model!\n');
104-
end
105-
106-
function [currBixel] = getBixelIndicesOnRay(this,currBixel,currRay)
107-
108-
% create offset vector to account for additional offsets modelled in the base data and a potential
109-
% range shifter. In the following, we only perform dose calculation for voxels having a radiological depth
110-
% that is within the limits of the base data set (-> machine.data(i).dephts). By this means, we only allow
111-
% interpolations in this.calcParticleDoseBixel() and avoid extrapolations.
112-
%urrBixel.offsetRadDepth = currBixel.baseData.offset + currBixel.radDepthOffset;
113-
tmpOffset = currBixel.baseData.offset - currBixel.radDepthOffset;
114-
115-
maxDepth = 1.15 * currBixel.baseData.range;
116-
117-
% find depth depended lateral cut off
118-
if this.dosimetricLateralCutOff == 1
119-
currIx = currRay.radDepths <= maxDepth + tmpOffset;
120-
elseif this.dosimetricLateralCutOff < 1 && this.dosimetricLateralCutOff > 0
121-
currIx = currRay.radDepths <= maxDepth + tmpOffset;
122-
sigmaSq = this.calcSigmaLatMCS(currRay.radDepths(currIx) - tmpOffset, currBixel.baseData.energy).^2 + currBixel.sigmaIniSq;
123-
currIx(currIx) = currRay.radialDist_sq(currIx) < currBixel.baseData.LatCutOff.numSig.^2*sigmaSq;
124-
else
125-
matRad_cfg = MatRad_Config.instance();
126-
matRad_cfg.dispError('Cutoff must be a value between 0 and 1!')
127-
end
128-
129-
currBixel.subRayIx = currIx;
130-
currBixel.ix = currRay.ix(currIx);
131-
end
132-
133-
function X = interpolateKernelsInDepth(this,bixel)
134-
baseData = bixel.baseData;
135-
136-
% calculate particle dose for bixel k on ray j of beam i
137-
% convert from MeV cm^2/g per primary to Gy mm^2 per 1e6 primaries
138-
conversionFactor = 1.6021766208e-02;
139-
140-
radDepthOffset = bixel.radDepthOffset;
141-
142-
if isfield(baseData,'energySpectrum')
143-
energyMean = baseData.energySpectrum.mean;
144-
energySpread = baseData.energySpectrum.sigma/100 * baseData.energySpectrum.mean;
145-
else
146-
energyMean = baseData.energy;
147-
energySpread = baseData.energy * this.sigmaEnergy;
148-
end
149-
150-
if isfield(baseData,'offset') && ~isfield(baseData,'energySpectrum')
151-
radDepthOffset = radDepthOffset - baseData.offset;
152-
end
153-
154-
X.Z = conversionFactor * this.calcAnalyticalBragg(energyMean, bixel.radDepths + radDepthOffset, energySpread);
155-
X.sigma = this.calcSigmaLatMCS(bixel.radDepths + radDepthOffset, baseData.energy);
156-
end
157-
158-
function bixel = calcParticleBixel(this,bixel)
159-
160-
kernel = this.interpolateKernelsInDepth(bixel);
161-
162-
%compute lateral sigma
163-
sigmaSq = kernel.sigma.^2 + bixel.sigmaIniSq;
164-
165-
% calculate dose
166-
bixel.physicalDose = bixel.baseData.LatCutOff.CompFac * exp( -bixel.radialDist_sq ./ (2*sigmaSq)) .* kernel.Z ./(2*pi*sigmaSq);
167-
168-
% check if we have valid dose values
169-
if any(isnan(bixel.physicalDose)) || any(bixel.physicalDose<0)
170-
matRad_cfg = MatRad_Config.instance();
171-
matRad_cfg.dispError('Error in particle dose calculation.');
172-
end
173-
end
174-
175-
function doseVector = calcAnalyticalBragg(this, primaryEnergy, depthZ, energySpread)
75+
methods (Access = public)
76+
function [doseVector,hatD] = calcAnalyticalBragg(this, primaryEnergy, depthZ, energySpread)
17677
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17778
% call
17879
% this.calcAnalyticalBragg(PrimaryEnergy, depthz, WidthMod)
@@ -319,6 +220,107 @@ function chooseLateralModel(this)
319220

320221
sigmaMCS(depthZ==0) = 0;
321222
end
223+
end
224+
225+
methods (Access = protected)
226+
function dij = initDoseCalc(this,ct,cst,stf)
227+
228+
matRad_cfg = MatRad_Config.instance();
229+
230+
if this.calcLET == true
231+
matRad_cfg.dispWarning('Engine does not support LET calculation! Disabling!');
232+
this.calcLET = false;
233+
end
234+
235+
if this.calcBioDose == true
236+
matRad_cfg.dispWarning('Engine does not support BioDose calculation! Disabling!');
237+
this.calcBioDose = false;
238+
end
239+
240+
dij = this.initDoseCalc@DoseEngines.matRad_ParticlePencilBeamEngineAbstract(ct,cst,stf);
241+
end
242+
243+
function chooseLateralModel(this)
244+
%Now check if we need tho chose the lateral model because it
245+
%was set to auto
246+
matRad_cfg = MatRad_Config.instance();
247+
if strcmp(this.lateralModel,'auto')
248+
this.lateralModel = 'single';
249+
elseif ~strcmp(this.lateralModel,'single')
250+
matRad_cfg.dispWarning('Engine only supports analytically computed singleGaussian lateral Model!');
251+
this.lateralModel = 'single';
252+
end
253+
matRad_cfg.dispInfo('Using an analytically computed %s Gaussian pencil-beam kernel model!\n');
254+
end
255+
256+
function [currBixel] = getBixelIndicesOnRay(this,currBixel,currRay)
257+
258+
% create offset vector to account for additional offsets modelled in the base data and a potential
259+
% range shifter. In the following, we only perform dose calculation for voxels having a radiological depth
260+
% that is within the limits of the base data set (-> machine.data(i).dephts). By this means, we only allow
261+
% interpolations in this.calcParticleDoseBixel() and avoid extrapolations.
262+
%urrBixel.offsetRadDepth = currBixel.baseData.offset + currBixel.radDepthOffset;
263+
tmpOffset = currBixel.baseData.offset - currBixel.radDepthOffset;
264+
265+
maxDepth = 1.15 * currBixel.baseData.range;
266+
267+
% find depth depended lateral cut off
268+
if this.dosimetricLateralCutOff == 1
269+
currIx = currRay.radDepths <= maxDepth + tmpOffset;
270+
elseif this.dosimetricLateralCutOff < 1 && this.dosimetricLateralCutOff > 0
271+
currIx = currRay.radDepths <= maxDepth + tmpOffset;
272+
sigmaSq = this.calcSigmaLatMCS(currRay.radDepths(currIx) - tmpOffset, currBixel.baseData.energy).^2 + currBixel.sigmaIniSq;
273+
currIx(currIx) = currRay.radialDist_sq(currIx) < currBixel.baseData.LatCutOff.numSig.^2*sigmaSq;
274+
else
275+
matRad_cfg = MatRad_Config.instance();
276+
matRad_cfg.dispError('Cutoff must be a value between 0 and 1!')
277+
end
278+
279+
currBixel.subRayIx = currIx;
280+
currBixel.ix = currRay.ix(currIx);
281+
end
282+
283+
function X = interpolateKernelsInDepth(this,bixel)
284+
baseData = bixel.baseData;
285+
286+
% calculate particle dose for bixel k on ray j of beam i
287+
% convert from MeV cm^2/g per primary to Gy mm^2 per 1e6 primaries
288+
conversionFactor = 1.6021766208e-02;
289+
290+
radDepthOffset = bixel.radDepthOffset;
291+
292+
if isfield(baseData,'energySpectrum')
293+
energyMean = baseData.energySpectrum.mean;
294+
energySpread = baseData.energySpectrum.sigma/100 * baseData.energySpectrum.mean;
295+
else
296+
energyMean = baseData.energy;
297+
energySpread = baseData.energy * this.sigmaEnergy;
298+
end
299+
300+
if isfield(baseData,'offset') && ~isfield(baseData,'energySpectrum')
301+
radDepthOffset = radDepthOffset - baseData.offset;
302+
end
303+
304+
X.Z = conversionFactor * this.calcAnalyticalBragg(energyMean, bixel.radDepths + radDepthOffset, energySpread);
305+
X.sigma = this.calcSigmaLatMCS(bixel.radDepths + radDepthOffset, baseData.energy);
306+
end
307+
308+
function bixel = calcParticleBixel(this,bixel)
309+
310+
kernel = this.interpolateKernelsInDepth(bixel);
311+
312+
%compute lateral sigma
313+
sigmaSq = kernel.sigma.^2 + bixel.sigmaIniSq;
314+
315+
% calculate dose
316+
bixel.physicalDose = bixel.baseData.LatCutOff.CompFac * exp( -bixel.radialDist_sq ./ (2*sigmaSq)) .* kernel.Z ./(2*pi*sigmaSq);
317+
318+
% check if we have valid dose values
319+
if any(isnan(bixel.physicalDose)) || any(bixel.physicalDose<0)
320+
matRad_cfg = MatRad_Config.instance();
321+
matRad_cfg.dispError('Error in particle dose calculation.');
322+
end
323+
end
322324

323325
function calcLateralParticleCutOff(this,cutOffLevel,~)
324326
calcRange = false;

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

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,8 @@
4545
end
4646

4747
this = [email protected]_ParticlePencilBeamEngineAbstract(pln);
48+
49+
this.restrictBeamRadDepthsByMaxEnergy = false;
4850
end
4951

5052
function setDefaults(this)
@@ -66,6 +68,14 @@ function setDefaults(this)
6668
ray.rotMat_system_T = beam.rotMat_system_T;
6769
end
6870

71+
function [currBixel] = getBixelIndicesOnRay(this,currBixel,currRay)
72+
% In Finesampling we can not reduce the number of points to be
73+
% computed that easily based on the rad depth, so we overload
74+
% this method to just give all indices stored in the ray
75+
currBixel.subRayIx = true(size(currRay.ix));
76+
currBixel.ix = currRay.ix;
77+
end
78+
6979
% We override this function to get full lateral distances
7080
function ray = getRayGeometryFromBeam(this,ray,currBeam)
7181
lateralRayCutOff = this.getLateralDistanceFromDoseCutOffOnRay(ray);

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

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,8 @@
3838
vBetaX; % Stores Photon Beta
3939

4040
bioKernelQuantities; % Kernel quantites to request from the machine data for biological dose calculation
41+
42+
restrictBeamRadDepthsByMaxEnergy = true;
4143
end
4244

4345
methods
@@ -449,12 +451,12 @@ function chooseLateralModel(this)
449451
radDepthOffset = this.machine.data(maxEnergyIx).offset + minRaShi;
450452

451453
% apply limit in depth
452-
%subSelectIx = currBeam.radDepths{1} < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset);
453-
454-
subSelectIx = cellfun(@(rD) rD < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset),currBeam.radDepths,'UniformOutput',false);
455-
currBeam.validCoords = cellfun(@and,subSelectIx,currBeam.validCoords,'UniformOutput',false);
456-
currBeam.validCoordsAll = any(cell2mat(currBeam.validCoords),2);
457-
454+
if this.restrictBeamRadDepthsByMaxEnergy
455+
subSelectIx = cellfun(@(rD) rD < (this.machine.data(maxEnergyIx).depths(end) - radDepthOffset),currBeam.radDepths,'UniformOutput',false);
456+
currBeam.validCoords = cellfun(@and,subSelectIx,currBeam.validCoords,'UniformOutput',false);
457+
currBeam.validCoordsAll = any(cell2mat(currBeam.validCoords),2);
458+
end
459+
458460
%currBeam.ixRadDepths = currBeam.ixRadDepths(subSelectIx);
459461
%currBeam.subIxVdoseGrid = currBeam.subIxVdoseGrid(subSelectIx);
460462
%currBeam.radDepths = cellfun(@(rd) rd(subSelectIx),currBeam.radDepths,'UniformOutput',false);

0 commit comments

Comments
 (0)