Merge branch 'dev_varRBErobOpt' of https://github.com/e0404/matRad in…

…to dev_MixedModality_mergeVarRBE

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

dicom/matRad_convRtssContours2Indices.m

@@ -29,7 +29,7 @@
 % LICENSE file.
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+matRad_cfg = MatRad_Config.instance();
 voiCube = zeros(ct.cubeDim);
 
 % loop over all closed contour items
@@ -40,7 +40,7 @@
         dicomCtSlicePos = unique(structure.item(i).points(:,3));
 
         if numel(dicomCtSlicePos) > 1 || isempty(dicomCtSlicePos)
-            error('Contour defined over multiple planes!');
+            matRad_cfg.dispError('Contour defined over multiple planes!');
         end
 
         round2 = @(a,b) round(a*10^b)/10^b;
@@ -49,7 +49,7 @@
         %Sanity check
         msg = checkSliceThickness(dicomCtSliceThickness);
         if ~isempty(msg)
-            error('Slice Thickness of slice at %f could not be identified: %s',dicomCtSlicePos,msg);
+            matRad_cfg.dispError('Slice Thickness of slice at %f could not be identified: %s',dicomCtSlicePos,msg);
         end
 
         slicesInMatradCt = find(dicomCtSlicePos+dicomCtSliceThickness/2 > ct.z & dicomCtSlicePos-dicomCtSliceThickness/2 <= ct.z);
@@ -67,6 +67,33 @@
     end
 
 end
+% The x- & y-direction in lps-coordinates are specified in:
+% ImageOrientationPatient
+
+xDir = ct.dicomInfo.ImageOrientationPatient(1:3); % lps: [1;0;0]
+yDir = ct.dicomInfo.ImageOrientationPatient(4:6); % lps: [0;1;0]
+nonStandardDirection = false;
+
+if xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0
+%     matRad_cfg.dispInfo('x-direction OK\n')
+elseif xDir(1) == -1 && xDir(2) == 0 && xDir(3) == 0
+    voiCube = flip(voiCube,1);
+else
+    nonStandardDirection = true;
+end
+
+if yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0
+%     matRad_cfg.dispInfo('y-direction OK\n')
+elseif yDir(1) == 0 && yDir(2) == -1 && yDir(3) == 0
+    voiCube = flip(voiCube,2);
+else
+    nonStandardDirection = true;
+end
+
+if nonStandardDirection
+    matRad_cfg.dispWarning(['Non-standard patient orientation.\n'...
+        'CT might not fit to contoured structures\n'])
+end
 
 indices = find(voiCube(:));
 

dicom/matRad_importDicom.m

@@ -82,6 +82,22 @@
     h = waitbar(0,'Please wait...');
     %h.WindowStyle = 'Modal';
     steps = numel(structures);
+
+
+    % The x- & y-direction in lps-coordinates are specified in:
+    % ImageOrientationPatient
+
+xDir = ct.dicomInfo.ImageOrientationPatient(1:3); % lps: [1;0;0]
+yDir = ct.dicomInfo.ImageOrientationPatient(4:6); % lps: [0;1;0]
+
+
+if ~(xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0)
+     matRad_cfg.dispInfo('\nNonstandard image orientation: tring to Mirror RTSS x-direction...')
+end
+
+if ~(yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0)
+    matRad_cfg.dispInfo('\nNonstandard image orientation: trying to Mirror RTSS y direction...')
+end
     for i = 1:numel(structures)
         % computations take place here
         waitbar(i / steps)
@@ -162,15 +178,30 @@
 end
 
 
+
 %% save ct, cst, pln, dose
 matRadFileName = [files.ct{1,3} '.mat']; % use default from dicom
 [FileName,PathName] = uiputfile('*','Save as...',matRadFileName);
 if ischar(FileName)
     % delete unnecessary variables
     if matRad_cfg.isMatlab
+        varNames=who;
+        for var=1:length(varNames)
+            cleanVar=eval([varNames{var}]);
+            if isempty(cleanVar)
+                eval(['clear ' varNames{var} ';'])
+            end
+        end
         clearvars -except ct cst pln stf resultGUI FileName PathName;
         save([PathName, FileName], '-regexp', '^(?!(FileName|PathName)$).','-v7');
     elseif matRad_cfg.isOctave
+        varNames=who;
+        for var=1:length(varNames)
+            cleanVar=eval([varNames{var}]);
+            if isempty(cleanVar)
+                eval(['clear ' varNames{var} ';'])
+            end
+        end
         clear -x ct cst pln stf resultGUI FileName PathName;
         save([PathName, FileName],'-v6');
     else

dicom/matRad_importDicomCt.m

@@ -123,13 +123,13 @@
 % patients in a supine position. Other orientations (e.g. prone, decubitus
 % left/right) are not supported.
 % Defined Terms:
-% HFP     Head First-Prone                  (not supported)
+% HFP     Head First-Prone                  (supported)
 % HFS     Head First-Supine                 (supported)
 % HFDR    Head First-Decubitus Right        (not supported)
 % HFDL    Head First-Decubitus Left         (not supported)
 % FFDR    Feet First-Decubitus Right        (not supported)
 % FFDL    Feet First-Decubitus Left         (not supported)
-% FFP     Feet First-Prone                  (not supported)
+% FFP     Feet First-Prone                  (supported)
 % FFS     Feet First-Supine                 (supported)
 
 if isempty(regexp(ctInfo(1).PatientPosition,{'S','P'}, 'once'))
@@ -185,29 +185,29 @@
 nonStandardDirection = false;
 
 % correct x- & y-direction
-% 
-% if xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0
-%     matRad_cfg.dispInfo('x-direction OK\n')
-% elseif xDir(1) == -1 && xDir(2) == 0 && xDir(3) == 0
-%     matRad_cfg.dispInfo('\nMirroring x-direction...')
-%     origCt = flip(origCt,1);
-%     matRad_cfg.dispInfo('finished!\n')
-% else
-%     nonStandardDirection = true;
-% end
-%     
-% if yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0
-%     matRad_cfg.dispInfo('y-direction OK\n')
-% elseif yDir(1) == 0 && yDir(2) == -1 && yDir(3) == 0
-%     matRad_cfg.dispInfo('\nMirroring y-direction...')
-%     origCt = flip(origCt,2);
-%     matRad_cfg.dispInfo('finished!\n')
-% else
-%     nonStandardDirection = true;
-% end
+
+if xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0
+    matRad_cfg.dispInfo('x-direction OK\n')
+elseif xDir(1) == -1 && xDir(2) == 0 && xDir(3) == 0
+    matRad_cfg.dispInfo('\nMirroring x-direction...')
+    origCt = flip(origCt,1);
+    matRad_cfg.dispInfo('finished!\n')
+else
+    nonStandardDirection = true;
+end
+
+if yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0
+    matRad_cfg.dispInfo('y-direction OK\n')
+elseif yDir(1) == 0 && yDir(2) == -1 && yDir(3) == 0
+    matRad_cfg.dispInfo('\nMirroring y-direction...')
+    origCt = flip(origCt,2);
+    matRad_cfg.dispInfo('finished!\n')
+else
+    nonStandardDirection = true;
+end
 
 if nonStandardDirection
-    matRad_cfg.dispInfo(['Non-standard patient orientation.\n'...
+    matRad_cfg.dispWarning(['Non-standard patient orientation.\n'...
         'CT might not fit to contoured structures\n'])
 end
 

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;