TOPAS 4D Monte Carlo

examples/matRad_example10_4DphotonRobust.m

@@ -192,7 +192,7 @@
 stf = matRad_generateStf(ct,cst,pln);
 
 %% Dose Calculation
-dij = matRad_calcPhotonDose(ct,stf,pln,cst);
+dij =  matRad_calcDoseInfluence(ct,cst,stf,pln);
 
 %% Inverse Optimization  for IMPT based on RBE-weighted dose
 % The goal of the fluence optimization is to find a set of bixel/spot 

matRad/4D/matRad_addMovement.m

@@ -43,7 +43,7 @@
 
 expectedDVF = {'pull', 'push'};
 
-p = inputParser; 
+p = inputParser;
 addParameter(p,'dvfType','pull', @(x) any(validatestring(x,expectedDVF)))
 addParameter(p,'visBool',false, @islogical);
 parse(p,varargin{:});
@@ -64,64 +64,64 @@
 
 % generate scenarios
 for i = 1:numOfCtScen
-    
+
     if isfield(ct,'hlut')
         padValue = min(ct.hlut(:,2));
     else
         padValue = -1024;
     end
-    
+
     ct.dvf{i} = zeros([ct.cubeDim, 3]);
-    
+
     dVec = arrayfun(@(A)  A*sin((i-1)*pi / numOfCtScen)^2, amp);
-    
+
     ct.dvf{i}(:,:,:,1) = dVec(1); % deformation along x direction (i.e. 2nd coordinate in dose/ct)
     ct.dvf{i}(:,:,:,2) = dVec(2);
     ct.dvf{i}(:,:,:,3) = dVec(3);
-    
+
     matRad_cfg.dispInfo('Deforming ct phase %d with [dx,dy,dz] = [%f,%f,%f] voxels\n',i,dVec(1),dVec(2),dVec(3));
-    
+
     % warp ct
     switch env
         case 'MATLAB'
             ct.cubeHU{i} = imwarp(ct.cubeHU{1}, ct.dvf{i},'FillValues',padValue);
-            
+
             if isfield(ct,'cube')
                 ct.cube{i}   = imwarp(ct.cube{1},   ct.dvf{i},'FillValues',0);
             end
-            
+
             % warp cst
             for j = 1:size(cst,1)
                 tmp = zeros(ct.cubeDim);
                 tmp(cst{j,4}{1}) = 1;
                 tmpWarp     = imwarp(tmp, ct.dvf{i});
-                
+
                 cst{j,4}{i} = find(tmpWarp > .5);
             end
         case 'OCTAVE'
             ct.cubeHU{i} = displaceOctave(ct.cubeHU{1}, ct.dvf{i},'linear',padValue);
-            
+
             if isfield(ct,'cube')
                 ct.cube{i}   = displaceOctave(ct.cube{1},ct.dvf{i},'linear',0);
             end
-            
+
             % warp cst
             for j = 1:size(cst,1)
                 tmp = zeros(ct.cubeDim);
                 tmp(cst{j,4}{1}) = 1;
                 tmpWarp     = displaceOctave(tmp, ct.dvf{i},'linear',0);
-                
+
                 cst{j,4}{i} = find(tmpWarp > .5);
             end
-        otherwise   
+        otherwise
     end
-    
+
     % convert dvfs to [mm]
     ct.dvf{i}(:,:,:,1) = ct.dvf{i}(:,:,:,1).* ct.resolution.x;
     ct.dvf{i}(:,:,:,2) = ct.dvf{i}(:,:,:,2).*ct.resolution.y;
     ct.dvf{i}(:,:,:,3) = ct.dvf{i}(:,:,:,3).* ct.resolution.z;
-    
-    ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]);    
+
+    ct.dvf{i} = permute(ct.dvf{i}, [4,1,2,3]);
 end
 
 
@@ -139,8 +139,8 @@
 end
 
 function newCube = displaceOctave(cube,vectorfield,interpMethod,fillValue)
-x = 1:size(cube,1);
-y = 1:size(cube,2);
+x = 1:size(cube,2);
+y = 1:size(cube,1);
 z = 1:size(cube,3);
 
 [X,Y,Z] = meshgrid(x,y,z);

matRad/4D/matRad_calc4dDose.m

@@ -33,7 +33,7 @@
 % LICENSE file.
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+matRad_cfg = MatRad_Config.instance();
 if ~exist('accType','var')
     accType = 'DDM';
 end
@@ -74,47 +74,85 @@
     tmpResultGUI = matRad_calcCubes(totalPhaseMatrix(:,i),dij,i);
 
     % compute physical dose for physical opt
-    if isa(pln.bioModel,'matRad_EmptyBiologicalModel')
-        resultGUI.phaseDose{i} = tmpResultGUI.physicalDose;
-        % compute RBExDose with const RBE
-    elseif isa(pln.bioModel,'matRad_ConstantRBE')
-        resultGUI.phaseRBExDose{i} = tmpResultGUI.RBExDose;
-        % compute all fields
-    elseif isa(pln.bioModel,'matRad_LQBasedModel')
-        resultGUI.phaseAlphaDose{i}    = tmpResultGUI.alpha .* tmpResultGUI.physicalDose;
-        resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose;
-        ix = ax{i} ~=0;
-        resultGUI.phaseEffect{i}    = resultGUI.phaseAlphaDose{i} + resultGUI.phaseSqrtBetaDose{i}.^2;
-        resultGUI.phaseRBExDose{i}     = zeros(ct.cubeDim);
-        resultGUI.phaseRBExDose{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
-    else
-        matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
+    resultGUI.phaseDose{i} = tmpResultGUI.physicalDose;
+    if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
+        % compute RBExDose
+        if isa(pln.bioModel,'matRad_ConstantRBE')
+            resultGUI.phaseRBExDose{i} = tmpResultGUI.RBExDose;
+        elseif isa(pln.bioModel,'matRad_LQBasedModel')
+            resultGUI.phaseAlphaDose{i}    = tmpResultGUI.alpha .* tmpResultGUI.physicalDose;
+            resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose;
+            ix = ax{i} ~=0;
+            resultGUI.phaseEffect{i}    = resultGUI.phaseAlphaDose{i} + resultGUI.phaseSqrtBetaDose{i}.^2;
+            resultGUI.phaseRBExDose{i}     = zeros(ct.cubeDim);
+            resultGUI.phaseRBExDose{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
+        else
+            matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
+        end
+    end
+
+     for beamIx = 1:dij.numOfBeams
+        resultGUI.(['phaseDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
+        if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
+            % compute RBExD 
+            if isa(pln.bioModel,'matRad_ConstantRBE')
+                resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['RBExDose_beam', num2str(beamIx)]);
+            elseif isa(pln.bioModel,'matRad_LQBasedModel')
+                resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} = tmpResultGUI.(['alpha_beam', num2str(beamIx)]).*tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
+                resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i} = sqrt(tmpResultGUI.(['beta_beam', num2str(beamIx)])).*tmpResultGUI.(['physicalDose_beam', num2str(beamIx)]);
+                 ix = ax{i} ~=0;
+                 resultGUI.(['phaseEffect_beam', num2str(beamIx)]){i}    = resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} + resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i}.^2;  
+                 resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i}    = zeros(ct.cubeDim);
+                 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)));
+            else
+                matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
+            end
+        end
     end
 end
 
 % accumulation
-if isa(pln.bioModel,'matRad_EmptyBiologicalModel')
-
-    resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
-
-elseif isa(pln.bioModel,'matRad_ConstantRBE')
-
-    resultGUI.accRBExDose = matRad_doseAcc(ct,resultGUI.phaseRBExDose, cst, accType);
-
-elseif isa(pln.bioModel,'matRad_LQBasedModel')
+resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
 
-    resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
-    resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
+if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
+    if isa(pln.bioModel,'matRad_ConstantRBE')
+
+        resultGUI.accRBExDose = matRad_doseAcc(ct,resultGUI.phaseRBExDose, cst, accType);
+
+    elseif isa(pln.bioModel,'matRad_LQBasedModel')
+
+        resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
+        resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
+
+        % only compute where we have biologically defined tissue
+        ix = (ax{1} ~= 0);
+
+        resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
+
+        resultGUI.accRBExDose     = zeros(ct.cubeDim);
+        resultGUI.accRBExDose(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
+    else
+        matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
+    end
+end
 
-    % only compute where we have biologically defined tissue
-    ix = (ax{1} ~= 0);
+for beamIx = 1:dij.numOfBeams
+    resultGUI.(['accPhysicalDose_beam', num2str(beamIx)])= matRad_doseAcc(ct,resultGUI.(['phaseDose_beam', num2str(beamIx)]), cst, accType);
+    if ~isa(pln.bioModel,'matRad_EmptyBiologicalModel')
+        if isa(pln.bioModel,'matRad_ConstantRBE')
+            resultGUI.(['accRBExDose_beam', num2str(beamIx)]) = matRad_doseAcc(ct,resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]), cst, accType);
+       elseif isa(pln.bioModel,'matRad_LQBasedModel')
+            resultGUI.(['accAlphaDose_beam', num2str(beamIx)])      = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
+            resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)])   = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
+            resultGUI.(['accEffect_beam', num2str(beamIx)])         = resultGUI.(['accAlphaDose_beam', num2str(beamIx)]) + resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)]).^2;
+            resultGUI.(['accRBExDose_beam', num2str(beamIx)]){i}       = zeros(ct.cubeDim);
+            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)));    
+       else
+        matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
+        end
+    end
+end
 
-    resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
 
-    resultGUI.accRBExDose     = zeros(ct.cubeDim);
-    resultGUI.accRBExDose(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
-else
-    matRad_cfg.dispError('Unsupported biological model %s!',pln.bioModel.model);
-end
 end
 

matRad/4D/matRad_makePhaseMatrix.m

@@ -70,7 +70,7 @@
 
     % permuatation of phaseMatrix from SS order to STF order
     timeSequence(i).phaseMatrix = timeSequence(i).phaseMatrix(timeSequence(i).orderToSTF,:);
-
+    [timeSequence(i).phaseNum,~] = find(timeSequence(i).phaseMatrix');
     % inserting the fluence in phaseMatrix
     timeSequence(i).phaseMatrix = timeSequence(i).phaseMatrix .* timeSequence(i).w;
 

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

@@ -290,17 +290,41 @@
             if ~isa(this.multScen,'matRad_ScenarioModel')
                 this.multScen = matRad_ScenarioModel.create(this.multScen,struct('numOfCtScen',ct.numOfCtScen));
             end
-            
+
             for i = 1:this.multScen.totNumScen
                 scenSubIx = this.multScen.linearMask(i,:);
                 resultGUItmp = matRad_calcCubes(ones(dij.numOfBeams,1),dij,this.multScen.sub2scenIx(scenSubIx(1),scenSubIx(2),scenSubIx(3)));
                 if i == 1
                     resultGUI = resultGUItmp;
                 end
-                resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));                
+                if isvector(this.multScen.scenMask) &&  this.multScen.numOfCtScen>1%ctScen
+                    resultGUI.phaseDose{i} = resultGUItmp.physicalDose;
+                    for beamIx = 1:dij.numOfBeams
+                        resultGUI.(['phaseDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
+                    end
+                    if isfield(resultGUItmp, 'alphaDoseCube') && isfield(resultGUItmp, 'SqrtBetaDoseCube')
+                        resultGUI.phaseAlphaDose{i}    = resultGUItmp.alpha .* resultGUItmp.physicalDose;
+                        resultGUI.phaseSqrtBetaDose{i} = sqrt(resultGUItmp.beta) .* resultGUItmp.physicalDose;
+                        resultGUI.phaseRBExDose{i} = resultGUItmp.RBExDose;
+                        for beamIx = 1:dij.numOfBeams
+                            resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['alpha_beam', num2str(beamIx)]).*resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
+                            resultGUI.(['phaseSqrtBetaDose_beam', num2str(beamIx)]){i} = sqrt(resultGUItmp.(['beta_beam', num2str(beamIx)])).*resultGUItmp.(['physicalDose_beam', num2str(beamIx)]);
+                            resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['RBExDose_beam', num2str(beamIx)]);
+                        end
+                    elseif isfield(resultGUItmp,'RBExDose')
+                        resultGUI.phaseRBExDose{i} = resultGUItmp.RBExDose;
+                        for beamIx = 1:dij.numOfBeams
+                          resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]){i} = resultGUItmp.(['RBExDose_beam', num2str(beamIx)]);
+                        end
+                    end
+                else
+                    resultGUI = matRad_appendResultGUI(resultGUI,resultGUItmp,false,sprintf('scen%d',i));
+                end
             end
 
+
             resultGUI.w  = w; 
+
         end
 
         function dij = calcDoseInfluence(this,ct,cst,stf)