Update: dose Accumulation

4D/matRad_addMovement.m

@@ -1,4 +1,4 @@
-function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp,visBool)
+function [ct, cst] = matRad_addMovement(ct, cst, motionPeriod, numOfCtScen, amp, dvfType, visBool)
 % adds artificial sinosodal patient motion by creating a deformation vector
 % field and applying it to the ct.cube by geometric transformation
 %
@@ -11,6 +11,7 @@
 %   motionPeriod:   the length of a whole breathing cycle (in seconds)
 %   numOfCtScen:    number of ct phases
 %   amp:            amplitude of the sinosoidal movement (in pixels)
+%   dvfType:        push or pull dvf
 %   visBool         boolean flag for visualization
 %
 %   note:           1st dim --> x LPS coordinate system
@@ -43,15 +44,25 @@
 if ~exist('visBool','var')
     visBool = false;
 end
+if ~exist('dvfType','var')
+   dvfType = 'pull';
+end
 
 matRad_cfg = MatRad_Config.instance();
 
 % book keeping
 ct.motionPeriod = motionPeriod;
 ct.numOfCtScen = numOfCtScen;
 
+if ~(strcmp(dvfType, 'pull') || strcmp(dvfType,'push'))
+    matRad_cfg.disperror('Chosse Push or Pull DVF type')
+end
+
 % set type
-ct.dvfType = 'pull'; % push or pull
+ct.dvfMetadata.dvfType = dvfType;
+if strcmp(dvfType,'push')
+    amp = -amp; %dvf_pull = -dvf_push;
+end
 
 env = matRad_getEnvironment();
 
@@ -110,10 +121,9 @@
     end
 
     % convert dvfs to [mm]
-    tmp = ct.dvf{i}(:,:,:,1);
-    ct.dvf{i}(:,:,:,1) = -ct.dvf{i}(:,:,:,2) * ct.resolution.x;
-    ct.dvf{i}(:,:,:,2) = -tmp * ct.resolution.y;
-    ct.dvf{i}(:,:,:,3) = -ct.dvf{i}(:,:,:,3) * ct.resolution.z;
+    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]);    
 end

4D/matRad_calc4dDose.m

@@ -1,4 +1,4 @@
-function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix)
+function [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI, totalPhaseMatrix,accType)
 % wrapper for the whole 4D dose calculation pipeline and calculated dose
 % accumulation
 %
@@ -13,6 +13,7 @@
 %   cst:             matRad cst struct
 %   resultGUI:       struct containing optimized fluence vector
 %   totalPhaseMatrix optional intput for totalPhaseMatrix
+%   accType:         witch algorithim for dose accumulation
 % output
 %   resultGUI:      structure containing phase dose, RBE weighted dose, etc
 %   timeSequence:   timing information about the irradiation
@@ -33,6 +34,10 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+if ~exist('accType','var')
+    accType = 'DDM';
+end
+
 if ~exist('totalPhaseMatrix','var')
    % make a time sequence for when each bixel is irradiated, the sequence
    % follows the backforth spot scanning
@@ -62,7 +67,7 @@
     elseif strcmp(pln.bioParam.model,'constRBE')
         resultGUI.phaseRBExD{i} = tmpResultGUI.RBExD;
     % compute all fields
-    elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'}))
+    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;
     end
@@ -72,16 +77,16 @@
 % accumulation
 if strcmp(pln.bioParam.model,'none')       
 
-    resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, 'DDM');
+    resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
 
 elseif strcmp(pln.bioParam.model,'constRBE')
 
-    resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, 'DDM');
+    resultGUI.accRBExD = matRad_doseAcc(ct,resultGUI.phaseRBExD, cst, accType);
 
-elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED'}))
+elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'}))
 
-    resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst, 'DDM');
-    resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, 'DDM');
+    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 = dij.ax(:,1)~=0; 

4D/matRad_doseAcc.m

@@ -57,7 +57,7 @@
 
 if strcmp(accMethod,'DDM')
 
-    if ~strcmp(ct.dvfType,'pull')
+    if ~strcmp(ct.dvfMetadata.dvfType,'pull')
         error('dose accumulation via direct dose mapping (DDM) requires pull dvfs');
     end    
 
@@ -75,19 +75,16 @@
         dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y;
         dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z;
 
-        d_ref = matRad_interp3(yGridVec,xGridVec,zGridVec,permute(phaseCubes{i},[2 1 3]), ...
-                         Y(ix) + dvf_y_i(ix), ...     
-                         X(ix) + dvf_x_i(ix), ... 
-                         Z(ix) + dvf_z_i(ix), ...
-                         'linear',0);  
+        d_ref = matRad_interp3(X, Y, Z,...
+                            phaseCubes{i}, X-dvf_x_i,Y-dvf_y_i,Z-dvf_z_i,'linear',0);   
 
-        dAcc(ix) = dAcc(ix) + d_ref;
+        dAcc(ix) = dAcc(ix) + d_ref(ix);
 
     end
 
 elseif strcmp(accMethod,'EMT')   % funktioniert nicht wenn Dosis in einer Phase = 0 ist...
 
-    if ~strcmp(ct.dvfType,'push')
+    if ~strcmp(ct.dvfMetadata.dvfType,'push')
         error('dose accumulation via interpolation requires push dvfs');
     end
 
@@ -99,8 +96,10 @@
         dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y;
         dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z;
 
-        m_i     = ct.cube{i};
-        e_i     = phaseCubes{i}.*m_i;
+        m_i = ct.cube{i};
+        m_i = permute(m_i,[2 1 3]);
+        e_i = phaseCubes{i}.*ct.cube{i};
+        e_i     = permute(e_i,[2 1 3]);
 
         ix = e_i>0;
 
@@ -163,6 +162,7 @@
          k = find(m_ref);
          dAcc(k) = dAcc(k) + e_ref(k)./m_ref(k);
 
+          dAcc = permute(dAcc,[2 1 3]);
     end
 
 % elseif strcmp(accMethod,'DDMM')

tools/matRad_compareDose.m

@@ -142,7 +142,7 @@
 
     % Calculate absolute difference cube and dose windows for plots
     differenceCube  = cube1-cube2;
-    doseDiffWindow  = [-max(differenceCube(:)) max(differenceCube(:))];
+    doseDiffWindow  = [-max(abs(differenceCube(:))) max(abs(differenceCube(:)))];
     %doseGammaWindow = [0 max(gammaCube(:))];
     doseGammaWindow = [0 2]; %We choose 2 as maximum value since the gamma colormap has a sharp cut in the middle