Merge pull request #685 from JenHardt/dev_varRBErobOpt_doseAcc

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, varargin)
 % adds artificial sinosodal patient motion by creating a deformation vector
 % field and applying it to the ct.cube by geometric transformation
 %
@@ -11,7 +11,8 @@
 %   motionPeriod:   the length of a whole breathing cycle (in seconds)
 %   numOfCtScen:    number of ct phases
 %   amp:            amplitude of the sinosoidal movement (in pixels)
-%   visBool         boolean flag for visualization
+%   varargin:   dvfType:        push or pull dvf
+%               visBool         boolean flag for visualization
 %
 %   note:           1st dim --> x LPS coordinate system
 %                   2nd dim --> y LPS coordinate system
@@ -40,9 +41,12 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-if ~exist('visBool','var')
-    visBool = false;
-end
+expectedDVF = {'pull', 'push'};
+
+p = inputParser; 
+addParameter(p,'dvfType','pull', @(x) any(validatestring(x,expectedDVF)))
+addParameter(p,'visBool',false, @islogical);
+parse(p,varargin{:});
 
 matRad_cfg = MatRad_Config.instance();
 
@@ -51,7 +55,10 @@
 ct.numOfCtScen = numOfCtScen;
 
 % set type
-ct.dvfType = 'pull'; % push or pull
+ct.dvfMetadata.dvfType = p.Results.dvfType;
+if strcmp(p.Results.dvfType,'push')
+    amp = -amp; %dvf_pull = -dvf_push;
+end
 
 env = matRad_getEnvironment();
 
@@ -110,17 +117,16 @@
     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
 
 
 
-if visBool
+if p.Results.visBool
     slice = round(ct.cubeDim(3)/2);
     figure,
     for i = 1:numOfCtScen

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')

examples/matRad_example10_4DphotonRobust.m

@@ -122,7 +122,7 @@
 numOfCtScen  = 5;
 motionPeriod = 2.5; % [s] 
 
-[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,1);
+[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'visBool',true);
 
 % show the deformation vector field
 slice = 25; % select a specific slice and to plot the vector field

examples/matRad_example9_4DDoseCalcMinimal.m

@@ -29,7 +29,7 @@
 numOfCtScen  = 5;
 motionPeriod = 2.5; % [s] 
 
-[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude);
+[ct,cst] = matRad_addMovement(ct, cst,motionPeriod, numOfCtScen, amplitude,'dvfType','pull');
 % Set up a plan, compute dose influence on all phases, conventional optimization
 % meta information for treatment plan
 pln.numOfFractions  = 30;

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