Skip to content

Commit e9411f8

Browse files
wahlnwahln
authored
Merge pull request #834 from JenHardt/dev
TOPAS 4D Monte Carlo
2 parents ab21a6b + 8c08414 commit e9411f8

File tree

7 files changed

+835
-391
lines changed

7 files changed

+835
-391
lines changed

‎examples/matRad_example10_4DphotonRobust.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -192,7 +192,7 @@
192192
stf = matRad_generateStf(ct,cst,pln);
193193

194194
%% Dose Calculation
195-
dij = matRad_calcPhotonDose(ct,stf,pln,cst);
195+
dij = matRad_calcDoseInfluence(ct,cst,stf,pln);
196196

197197
%% Inverse Optimization for IMPT based on RBE-weighted dose
198198
% The goal of the fluence optimization is to find a set of bixel/spot

‎matRad/4D/matRad_addMovement.m‎

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@
4343

4444
expectedDVF = {'pull', 'push'};
4545

46-
p = inputParser;
46+
p = inputParser;
4747
addParameter(p,'dvfType','pull', @(x) any(validatestring(x,expectedDVF)))
4848
addParameter(p,'visBool',false, @islogical);
4949
parse(p,varargin{:});
@@ -64,64 +64,64 @@
6464

6565
% generate scenarios
6666
for i = 1:numOfCtScen
67-
67+
6868
if isfield(ct,'hlut')
6969
padValue = min(ct.hlut(:,2));
7070
else
7171
padValue = -1024;
7272
end
73-
73+
7474
ct.dvf{i} = zeros([ct.cubeDim, 3]);
75-
75+
7676
dVec = arrayfun(@(A) A*sin((i-1)*pi / numOfCtScen)^2, amp);
77-
77+
7878
ct.dvf{i}(:,:,:,1) = dVec(1); % deformation along x direction (i.e. 2nd coordinate in dose/ct)
7979
ct.dvf{i}(:,:,:,2) = dVec(2);
8080
ct.dvf{i}(:,:,:,3) = dVec(3);
81-
81+
8282
matRad_cfg.dispInfo('Deforming ct phase %d with [dx,dy,dz] = [%f,%f,%f] voxels\n',i,dVec(1),dVec(2),dVec(3));
83-
83+
8484
% warp ct
8585
switch env
8686
case 'MATLAB'
8787
ct.cubeHU{i} = imwarp(ct.cubeHU{1}, ct.dvf{i},'FillValues',padValue);
88-
88+
8989
if isfield(ct,'cube')
9090
ct.cube{i} = imwarp(ct.cube{1}, ct.dvf{i},'FillValues',0);
9191
end
92-
92+
9393
% warp cst
9494
for j = 1:size(cst,1)
9595
tmp = zeros(ct.cubeDim);
9696
tmp(cst{j,4}{1}) = 1;
9797
tmpWarp = imwarp(tmp, ct.dvf{i});
98-
98+
9999
cst{j,4}{i} = find(tmpWarp > .5);
100100
end
101101
case 'OCTAVE'
102102
ct.cubeHU{i} = displaceOctave(ct.cubeHU{1}, ct.dvf{i},'linear',padValue);
103-
103+
104104
if isfield(ct,'cube')
105105
ct.cube{i} = displaceOctave(ct.cube{1},ct.dvf{i},'linear',0);
106106
end
107-
107+
108108
% warp cst
109109
for j = 1:size(cst,1)
110110
tmp = zeros(ct.cubeDim);
111111
tmp(cst{j,4}{1}) = 1;
112112
tmpWarp = displaceOctave(tmp, ct.dvf{i},'linear',0);
113-
113+
114114
cst{j,4}{i} = find(tmpWarp > .5);
115115
end
116-
otherwise
116+
otherwise
117117
end
118-
118+
119119
% convert dvfs to [mm]
120120
ct.dvf{i}(:,:,:,1) = ct.dvf{i}(:,:,:,1).* ct.resolution.x;
121121
ct.dvf{i}(:,:,:,2) = ct.dvf{i}(:,:,:,2).*ct.resolution.y;
122122
ct.dvf{i}(:,:,:,3) = ct.dvf{i}(:,:,:,3).* ct.resolution.z;
123-
124-
ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]);
123+
124+
ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]);
125125
end
126126

127127

@@ -139,8 +139,8 @@
139139
end
140140

141141
function newCube = displaceOctave(cube,vectorfield,interpMethod,fillValue)
142-
x = 1:size(cube,1);
143-
y = 1:size(cube,2);
142+
x = 1:size(cube,2);
143+
y = 1:size(cube,1);
144144
z = 1:size(cube,3);
145145

146146
[X,Y,Z] = meshgrid(x,y,z);

‎matRad/4D/matRad_calc4dDose.m‎

Lines changed: 73 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@
3333
% LICENSE file.
3434
%
3535
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36-
36+
matRad_cfg = MatRad_Config.instance();
3737
if ~exist('accType','var')
3838
accType = 'DDM';
3939
end
@@ -74,47 +74,85 @@
7474
tmpResultGUI = matRad_calcCubes(totalPhaseMatrix(:,i),dij,i);
7575

7676
% compute physical dose for physical opt
77-
if isa(pln.bioModel,'matRad_EmptyBiologicalModel')
78-
resultGUI.phaseDose{i} = tmpResultGUI.physicalDose;
79-
% compute RBExDose with const RBE
80-
elseif isa(pln.bioModel,'matRad_ConstantRBE')
81-
resultGUI.phaseRBExDose{i} = tmpResultGUI.RBExDose;
82-
% compute all fields
83-
elseif isa(pln.bioModel,'matRad_LQBasedModel')
84-
resultGUI.phaseAlphaDose{i} = tmpResultGUI.alpha .* tmpResultGUI.physicalDose;
85-
resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose;
86-
ix = ax{i} ~=0;
87-
resultGUI.phaseEffect{i} = resultGUI.phaseAlphaDose{i} + resultGUI.phaseSqrtBetaDose{i}.^2;
88-
resultGUI.phaseRBExDose{i} = zeros(ct.cubeDim);
89-
resultGUI.phaseRBExDose{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
90-
else
91-
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
77+
resultGUI.phaseDose{i} = tmpResultGUI.physicalDose;
78+
if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
79+
% compute RBExDose
80+
if isa(pln.bioModel,'matRad_ConstantRBE')
81+
resultGUI.phaseRBExDose{i} = tmpResultGUI.RBExDose;
82+
elseif isa(pln.bioModel,'matRad_LQBasedModel')
83+
resultGUI.phaseAlphaDose{i} = tmpResultGUI.alpha .* tmpResultGUI.physicalDose;
84+
resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose;
85+
ix = ax{i} ~=0;
86+
resultGUI.phaseEffect{i} = resultGUI.phaseAlphaDose{i} + resultGUI.phaseSqrtBetaDose{i}.^2;
87+
resultGUI.phaseRBExDose{i} = zeros(ct.cubeDim);
88+
resultGUI.phaseRBExDose{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
89+
else
90+
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
91+
end
92+
end
93+
94+
for beamIx = 1:dij.numOfBeams
95+
resultGUI.(['phaseDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
96+
if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
97+
% compute RBExD
98+
if isa(pln.bioModel,'matRad_ConstantRBE')
99+
resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['RBExDose_beam', num2str(beamIx)]);
100+
elseif isa(pln.bioModel,'matRad_LQBasedModel')
101+
resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['alpha_beam', num2str(beamIx)]).*tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
102+
resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i} = sqrt(tmpResultGUI.(['beta_beam', num2str(beamIx)])).*tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
103+
ix = ax{i} ~=0;
104+
resultGUI.(['phaseEffect_beam', num2str(beamIx)]){i} = resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} + resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i}.^2;
105+
resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = zeros(ct.cubeDim);
106+
resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.(['phaseEffect_beam', num2str(beamIx)]){i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
107+
else
108+
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
109+
end
110+
end
92111
end
93112
end
94113

95114
% accumulation
96-
if isa(pln.bioModel,'matRad_EmptyBiologicalModel')
97-
98-
resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
99-
100-
elseif isa(pln.bioModel,'matRad_ConstantRBE')
101-
102-
resultGUI.accRBExDose = matRad_doseAcc(ct,resultGUI.phaseRBExDose, cst, accType);
103-
104-
elseif isa(pln.bioModel,'matRad_LQBasedModel')
115+
resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
105116

106-
resultGUI.accAlphaDose = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
107-
resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
117+
if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
118+
if isa(pln.bioModel,'matRad_ConstantRBE')
119+
120+
resultGUI.accRBExDose = matRad_doseAcc(ct,resultGUI.phaseRBExDose, cst, accType);
121+
122+
elseif isa(pln.bioModel,'matRad_LQBasedModel')
123+
124+
resultGUI.accAlphaDose = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
125+
resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
126+
127+
% only compute where we have biologically defined tissue
128+
ix = (ax{1} ~= 0);
129+
130+
resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
131+
132+
resultGUI.accRBExDose = zeros(ct.cubeDim);
133+
resultGUI.accRBExDose(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
134+
else
135+
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
136+
end
137+
end
108138

109-
% only compute where we have biologically defined tissue
110-
ix = (ax{1} ~= 0);
139+
for beamIx = 1:dij.numOfBeams
140+
resultGUI.(['accPhysicalDose_beam', num2str(beamIx)])= matRad_doseAcc(ct,resultGUI.(['phaseDose_beam', num2str(beamIx)]), cst, accType);
141+
if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
142+
if isa(pln.bioModel,'matRad_ConstantRBE')
143+
resultGUI.(['accRBExDose_beam', num2str(beamIx)]) = matRad_doseAcc(ct,resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]), cst, accType);
144+
elseif isa(pln.bioModel,'matRad_LQBasedModel')
145+
resultGUI.(['accAlphaDose_beam', num2str(beamIx)]) = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
146+
resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)]) = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
147+
resultGUI.(['accEffect_beam', num2str(beamIx)]) = resultGUI.(['accAlphaDose_beam', num2str(beamIx)]) + resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)]).^2;
148+
resultGUI.(['accRBExDose_beam', num2str(beamIx)]){i} = zeros(ct.cubeDim);
149+
resultGUI.(['accRBExDose_beam', num2str(beamIx)]){i} = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.(['accEffect_beam', num2str(beamIx)]){i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
150+
else
151+
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
152+
end
153+
end
154+
end
111155

112-
resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
113156

114-
resultGUI.accRBExDose = zeros(ct.cubeDim);
115-
resultGUI.accRBExDose(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
116-
else
117-
matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
118-
end
119157
end
120158

‎matRad/4D/matRad_makePhaseMatrix.m‎

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

7171
% permuatation of phaseMatrix from SS order to STF order
7272
timeSequence(i).phaseMatrix = timeSequence(i).phaseMatrix(timeSequence(i).orderToSTF,:);
73-
73+
[timeSequence(i).phaseNum,~] = find(timeSequence(i).phaseMatrix');
7474
% inserting the fluence in phaseMatrix
7575
timeSequence(i).phaseMatrix = timeSequence(i).phaseMatrix .* timeSequence(i).w;
7676

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

Lines changed: 26 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -290,17 +290,41 @@ function assignBioModelPropertiesFromPln(this, plnModel, warnWhenPropertyChanged
290290
if ~isa(this.multScen,'matRad_ScenarioModel')
291291
this.multScen = matRad_ScenarioModel.create(this.multScen,struct('numOfCtScen',ct.numOfCtScen));
292292
end
293-
293+
294294
for i = 1:this.multScen.totNumScen
295295
scenSubIx = this.multScen.linearMask(i,:);
296296
resultGUItmp = matRad_calcCubes(ones(dij.numOfBeams,1),dij,this.multScen.sub2scenIx(scenSubIx(1),scenSubIx(2),scenSubIx(3)));
297297
if i == 1
298298
resultGUI = resultGUItmp;
299299
end
300-
resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));
300+
if isvector(this.multScen.scenMask) && this.multScen.numOfCtScen>1%ctScen
301+
resultGUI.phaseDose{i} = resultGUItmp.physicalDose;
302+
for beamIx = 1:dij.numOfBeams
303+
resultGUI.(['phaseDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
304+
end
305+
if isfield(resultGUItmp, 'alphaDoseCube') && isfield(resultGUItmp, 'SqrtBetaDoseCube')
306+
resultGUI.phaseAlphaDose{i} = resultGUItmp.alpha .* resultGUItmp.physicalDose;
307+
resultGUI.phaseSqrtBetaDose{i} = sqrt(resultGUItmp.beta) .* resultGUItmp.physicalDose;
308+
resultGUI.phaseRBExDose{i} = resultGUItmp.RBExDose;
309+
for beamIx = 1:dij.numOfBeams
310+
resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['alpha_beam', num2str(beamIx)]).*resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
311+
resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i} = sqrt(resultGUItmp.(['beta_beam', num2str(beamIx)])).*resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
312+
resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['RBExDose_beam', num2str(beamIx)]);
313+
end
314+
elseif isfield(resultGUItmp,'RBExDose')
315+
resultGUI.phaseRBExDose{i} = resultGUItmp.RBExDose;
316+
for beamIx = 1:dij.numOfBeams
317+
resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['RBExDose_beam', num2str(beamIx)]);
318+
end
319+
end
320+
else
321+
resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));
322+
end
301323
end
302324

325+
303326
resultGUI.w = w;
327+
304328
end
305329

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

0 commit comments

Comments
 (0)