4D Update to scenario handling with carbon

e0404 · Apr 19, 2024 · 4d7d1ec · 4d7d1ec
1 parent fd85a05
commit 4d7d1ec
Show file tree

Hide file tree

Showing 4 changed files with 362 additions and 8 deletions.
diff --git a/4D/matRad_calc4dDose.m b/4D/matRad_calc4dDose.m
@@ -55,6 +55,11 @@
    timeSequence = [];
 
 end
+
+if any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'}))
+    [ax,bx] = matRad_getPhotonLQMParameters(cst,numel(resultGUI.physicalDose));
+end
+
 % compute all phases
 for i = 1:ct.numOfCtScen
 
@@ -70,6 +75,10 @@
     elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'}))
         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.phaseRBExD{i}     = zeros(ct.cubeDim);
+        resultGUI.phaseRBExD{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));
     end
 
 end
@@ -89,12 +98,12 @@
     resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
 
     % only compute where we have biologically defined tissue
-    ix = any(cell2mat(cellfun(@(ax) ax ~= 0,dij.ax','UniformOutput',false)),2);
+    ix = (ax{1} ~= 0);
 
     resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
 
     resultGUI.accRBExD     = zeros(ct.cubeDim);
-    resultGUI.accRBExD(ix) = ((sqrt(dij.ax(ix).^2 + 4 .* dij.bx(ix) .* resultGUI.accEffect(ix)) - dij.ax(ix))./(2.*dij.bx(ix)));
+    resultGUI.accRBExD(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
 
 end
 
diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m
@@ -48,6 +48,8 @@
 beamInfo(dij.numOfBeams+1).suffix = '';
 beamInfo(dij.numOfBeams+1).logIx  = true(size(resultGUI.w,1),1);
 
+[ctScen,~] = ind2sub(size(dij.physicalDose),scenNum);
+
 
 %% Physical Dose
 doseFields = {'physicalDose','doseToWater'};
@@ -112,15 +114,15 @@
                 wBeam = (resultGUI.w .* beamInfo(i).logIx);
 
                 % consider biological optimization
-                ix = dij.bx{scenNum}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0;
+                ix = dij.bx{ctScen}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0;
 
                 % Calculate effect from alpha- and sqrtBetaDose
                 resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])       = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam + (dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam).^2);
                 resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])       = reshape(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]),dij.doseGrid.dimensions);
 
                 % Calculate RBExD from the effect
                 resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])        = zeros(size(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])));
-                resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix)    = (sqrt(dij.ax{scenNum}(ix).^2 + 4 .* dij.bx{scenNum}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{scenNum}(ix))./(2.*dij.bx{scenNum}(ix));
+                resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix)    = (sqrt(dij.ax{ctScen}(ix).^2 + 4 .* dij.bx{ctScen}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{ctScen}(ix))./(2.*dij.bx{ctScen}(ix));
 
                 % Divide RBExD with the physicalDose to get the plain RBE cube
                 resultGUI.(['RBE', RBE_model{j}, beamInfo(i).suffix])          = resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])./resultGUI.(['physicalDose', beamInfo(i).suffix]);