Merge branch 'dev' of https://github.com/e0404/matRad into dev

.github/workflows/tests.yml

@@ -20,7 +20,9 @@ jobs:
     steps:
       - uses: actions/checkout@v2 # Checks-out repository under $GITHUB_WORKSPACE
       - name: Install OCTAVE
-        run: sudo apt-get install -y gdb gfortran fonts-freefont-otf gnuplot-x11 libgdcm-dev octave liboctave-dev
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y gdb gfortran fonts-freefont-otf gnuplot-x11 libgdcm-dev octave liboctave-dev
       - name: Prepare Test Environment
         run: |
           sudo chmod +x .travis/before_install_linux.sh

dicom/matRad_convRtssContours2Indices.m

@@ -39,20 +39,26 @@
 
         dicomCtSlicePos = unique(structure.item(i).points(:,3));
 
-        if numel(dicomCtSlicePos) > 1
-            error('Contour defined over multiple planes\n');
+        if numel(dicomCtSlicePos) > 1 || isempty(dicomCtSlicePos)
+            error('Contour defined over multiple planes!');
         end
 
         round2 = @(a,b) round(a*10^b)/10^b;
         dicomCtSliceThickness = ct.dicomInfo.SliceThickness(round2(ct.dicomInfo.SlicePositions,1)==round2(dicomCtSlicePos,1));
 
+        %Sanity check
+        msg = checkSliceThickness(dicomCtSliceThickness);
+        if ~isempty(msg)
+            error('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);
+
         coords1 = interp1(ct.x,1:ct.cubeDim(2),structure.item(i).points(:,1),'linear','extrap');
         coords2 = interp1(ct.y,1:ct.cubeDim(1),structure.item(i).points(:,2),'linear','extrap');
 
         binIn = poly2mask(coords1,coords2,ct.cubeDim(1),ct.cubeDim(2));
 
-        slicesInMatradCt = find(dicomCtSlicePos+dicomCtSliceThickness/2 > ct.z & dicomCtSlicePos-dicomCtSliceThickness/2 <= ct.z);
-
         % loop over all slices in matRad ct
         for j = 1:numel(slicesInMatradCt)
             voiCube(:,:,slicesInMatradCt(j)) = voiCube(:,:,slicesInMatradCt(j)) | binIn;
@@ -63,3 +69,17 @@
 end
 
 indices = find(voiCube(:));
+
+end
+
+function msg = checkSliceThickness(dicomCtSliceThickness)
+    if isempty(dicomCtSliceThickness)
+        msg = 'Slice could not be identified (empty)';
+    elseif ~isscalar(dicomCtSliceThickness)
+        msg = 'Slice thickness not unique';
+    elseif ~isnumeric(dicomCtSliceThickness)
+        msg = 'unexpected value';
+    else
+        msg = '';
+    end
+end

dicom/matRad_importDicom.m

@@ -1,10 +1,10 @@
-function [ct, cst, pln, resultGUI] = matRad_importDicom( files, dicomMetaBool )
+function [ct,cst,pln,stf,resultGUI] = matRad_importDicom( files, dicomMetaBool )
 % matRad wrapper function to import a predefined set of dicom files 
 % into matRad's native data formats
 % 
 % call
-% [ct, cst, pln, resultGUI] = matRad_importDicom( files )  
-% [ct, cst, pln, resultGUI] = matRad_importDicom( files, dicomMetaBool )
+% [ct, cst, pln, stf, resultGUI] = matRad_importDicom( files )  
+% [ct, cst, pln, stf, resultGUI] = matRad_importDicom( files, dicomMetaBool )
 %
 % input
 %   files:          list of files to be imported (will contain cts and rt
@@ -16,6 +16,7 @@
 %   ct:        matRad ct struct
 %   cst:       matRad cst struct
 %   pln:       matRad plan struct
+%   stf:       matRad stf struct
 %   resultGUI: matRad result struct holding data for visualization in GUI
 %
 % References
@@ -34,7 +35,7 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-[env, ~] = matRad_getEnvironment();
+%[env, ~] = matRad_getEnvironment();
 
 %%
 if ~exist('dicomMetaBool','var')
@@ -85,8 +86,14 @@
     for i = 1:numel(structures)
         % computations take place here
         waitbar(i / steps)
-        fprintf('creating cube for %s volume...\n', structures(i).structName);
-        structures(i).indices = matRad_convRtssContours2Indices(structures(i),ct);
+        fprintf('creating cube for %s volume... ', structures(i).structName);
+        try
+            structures(i).indices = matRad_convRtssContours2Indices(structures(i),ct);
+            fprintf('\n');
+        catch ME
+            warning('matRad:dicomImport','could not be imported: %s',ME.message);
+            structures(i).indices = [];
+        end      
     end
     fprintf('finished!\n');
     close(h)
@@ -105,6 +112,8 @@
     if ~(cellfun(@isempty,files.rtplan(1,:)))
         pln = matRad_importDicomRTPlan(ct, files.rtplan, dicomMetaBool);
     end
+else
+    pln = struct([]);
 end
 
 %% import stf
@@ -121,6 +130,8 @@
             warning('No support for DICOM import of steering information for this modality.');
         end
     end
+else
+    stf = struct([]);
 end
 
 %% import dose cube
@@ -138,30 +149,17 @@
             resultGUI = matRad_importDicomRTDose(ct, files.rtdose);
         end
         if size(resultGUI) == 0
-           clear resultGUI;
+           resultGUI = struct([]);
         end
     end
+else
+    resultGUI = struct([]);
 end
 
 %% put weight also into resultGUI
-if exist('stf','var') && exist('resultGUI','var')
+if ~isempty(stf) && ~isempty(resultGUI)
     resultGUI.w = [];
     for i = 1:size(stf,2)
         resultGUI.w = [resultGUI.w; [stf(i).ray.weight]'];
     end
 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
-    switch env
-    case 'MATLAB'
-        clearvars -except ct cst pln stf resultGUI FileName PathName;
-    case 'OCTAVE' 
-        clear -x ct cst pln stf resultGUI FileName PathName;
-    end
-    % save all except FileName and PathName
-    save([PathName, FileName], '-regexp', '^(?!(FileName|PathName)$).','-v7.3');
-end

dicom/matRad_importDicomGUI.m

@@ -327,7 +327,23 @@ function import_button_Callback(hObject, eventdata, handles)
 
 % dicomMetaBool: store complete DICOM information and patientName or not
 dicomMetaBool = logical(get(handles.checkPatientName,'Value'));
-matRad_importDicom(files, dicomMetaBool);
+[ct,cst,pln,stf,resultGUI] = matRad_importDicom(files, dicomMetaBool);
+
+
+% Save file
+matRadFileName = [files.ct{1,3} '.mat']; % use default from dicom
+[FileName,PathName] = uiputfile('*','Save as...',matRadFileName);
+if ischar(FileName)
+    %Hack to get non-empty variables
+    varNames = {'ct','cst','pln','stf','resultGUI'};
+    ix = cellfun(@(s) ~isempty(evalin('caller',s)),varNames);
+    varNames = varNames(ix);
+
+    FileName = fullfile(PathName,FileName);
+
+    % save all non-empty variables imported
+    save(FileName,'-v7.3',varNames{:});
+end
 
 
 % --- Executes on button press in cancel_button.

matRadGUI.m

@@ -3972,7 +3972,11 @@ function checkbox_lockColormap_Callback(hObject, eventdata, handles)
     classNames(:,clIx) = {cl.Name; pName}; %Store class name and display name
 end
 
-numOfObjectives = sum(cellfun(@numel,cst(:,6)));
+if size(cst,2) < 6
+    numOfObjectives = 0;
+else 
+    numOfObjectives = sum(cellfun(@numel,cst(:,6)));
+end
 
 cnt = 0;
 
@@ -4005,7 +4009,12 @@ function checkbox_lockColormap_Callback(hObject, eventdata, handles)
 cnt = cnt + 1;
 
 %Create Objectives / Constraints controls
-for i = 1:size(cst,1)   
+for i = 1:size(cst,1)
+    %Safety break in case the 6th column is empty, because it allows us to
+    %run less checks afterwards
+    if numOfObjectives == 0
+        break;
+    end
     if strcmp(cst(i,3),'IGNORED')~=1
         %Compatibility Layer for old objective format
         if isstruct(cst{i,6})

matRad_calcParticleDose.m

@@ -278,7 +278,7 @@
 
                     % empty bixels may happen during recalculation of error
                     % scenarios -> skip to next bixel
-                    if ~any(currIx)
+                    if ~(any(any(currIx)))
                         continue;
                     end
 
@@ -312,7 +312,7 @@
                     end
 
                     % run over components
-                    for c = 1:numOfSub
+                    for c = 1:numOfSub(k)
                         tmpDose = zeros(size(currIx,1),1);
                         bixelDose = finalWeight(c,k).*matRad_calcParticleDoseBixel(...
                                 radDepths(currIx(:,:,c),1,c), ...

matRad_calcPhotonDose.m

@@ -204,9 +204,9 @@
         else
             %For some reason the use of interpn here is much faster
             %than using interp2 in Octave
-            Interp_kernel1 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx1,x,y,'linear',NaN);
-            Interp_kernel2 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx2,x,y,'linear',NaN);
-            Interp_kernel3 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx3,x,y,'linear',NaN);
+            Interp_kernel1 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx1',x,y,'linear',NaN);
+            Interp_kernel2 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx2',x,y,'linear',NaN);
+            Interp_kernel3 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx3',x,y,'linear',NaN);
         end
     end
 
@@ -247,9 +247,9 @@
             else
                 %For some reason the use of interpn here is much faster
                 %than using interp2 in Octave
-                Interp_kernel1 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx1,x,y,'linear',NaN);
-                Interp_kernel2 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx2,x,y,'linear',NaN);
-                Interp_kernel3 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx3,x,y,'linear',NaN);
+                Interp_kernel1 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx1',x,y,'linear',NaN);
+                Interp_kernel2 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx2',x,y,'linear',NaN);
+                Interp_kernel3 = @(x,y)interpn(convMx_X(1,:),convMx_Z(:,1),convMx3',x,y,'linear',NaN);
             end
 
         end

matRad_fluenceOptimization.m

@@ -1,13 +1,15 @@
-function [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln)
+function [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln,wInit)
 % matRad inverse planning wrapper function
 % 
 % call
 %   [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln)
+%   [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln,wInit)
 %
 % input
 %   dij:        matRad dij struct
 %   cst:        matRad cst struct
 %   pln:        matRad pln struct
+%   wInit:      (optional) custom weights to initialize problems
 %
 % output
 %   resultGUI:  struct containing optimized fluence vector, dose, and (for
@@ -109,13 +111,16 @@
 
 end
 % calculate initial beam intensities wInit
-if  strcmp(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons')
+if exist('wInit','var')
+    %do Nothing
+elseif  strcmp(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons')
 
     % check if a constant RBE is defined - if not use 1.1
     if ~isfield(dij,'RBE')
         dij.RBE = 1.1;
     end
-    bixelWeight =  (doseTarget)/(dij.RBE * mean(dij.physicalDose{1}(V,:)*wOnes)); 
+    doseTmp = dij.physicalDose{1}*wOnes;
+    bixelWeight =  (doseTarget)/(dij.RBE * mean(doseTmp(V)));     
     wInit       = wOnes * bixelWeight;
 
 elseif (strcmp(pln.propOpt.bioOptimization,'LEMIV_effect') || strcmp(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
@@ -145,8 +150,11 @@
     if isequal(pln.propOpt.bioOptimization,'LEMIV_effect')
 
            effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2;
-           p            = (sum(dij.mAlphaDose{1}(V,:)*wOnes)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
-           q            = -(effectTarget * length(V)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
+           aTmp = dij.mAlphaDose{1}*wOnes;
+           bTmp = dij.mSqrtBetaDose{1} * wOnes;
+           p = sum(aTmp(V)) / sum(bTmp(V).^2);
+           q = -(effectTarget * length(V)) / sum(bTmp(V).^2);
+
            wInit        = -(p/2) + sqrt((p^2)/4 -q) * wOnes;
 
     elseif isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')
@@ -155,14 +163,18 @@
            dij.gamma              = zeros(dij.doseGrid.numOfVoxels,1);   
            dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); 
 
-           % calculate current in target
-           CurrEffectTarget = (dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2);
+           % calculate current effect in target
+           aTmp = dij.mAlphaDose{1}*wOnes;
+           bTmp = dij.mSqrtBetaDose{1} * wOnes;
+           doseTmp = dij.physicalDose{1}*wOnes;
+
+           CurrEffectTarget = aTmp(V) + bTmp(V).^2;
            % ensure a underestimated biological effective dose 
            TolEstBio        = 1.2;
            % calculate maximal RBE in target
            maxCurrRBE = max(-cst{ixTarget,5}.alphaX + sqrt(cst{ixTarget,5}.alphaX^2 + ...
-                        4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*(dij.physicalDose{1}(V,:)*wOnes)));
-           wInit    =  ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes;
+                        4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*doseTmp(V)));
+           wInit    =  ((doseTarget)/(TolEstBio*maxCurrRBE*max(doseTmp(V))))* wOnes;
     end
 
 else 

matRad_setOverlapPriorities.m

@@ -32,6 +32,7 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+matRad_cfg = MatRad_Config.instance();
 numOfCtScenarios = unique(cellfun(@(x)numel(x),cst(:,4)));
 
 if numel(numOfCtScenarios) > 1
@@ -56,7 +57,7 @@
         cst{j,4}{i} = idx;
 
         if isempty(cst{j,4}{i}) && ~isempty(cst{j,6})
-            error([cst{j,2} ': Objective(s) and/or constraints for inverse planning defined ' ...
+            matRad_cfg.dispWarning([cst{j,2} ': Objective(s) and/or constraints for inverse planning defined ' ...
                  'but structure overlapped by structure with higher overlap priority.' ...
                  'Objective(s) will not be considered during optimization']); 
         end

optimization/+DoseObjectives/matRad_MeanDose.m

@@ -68,6 +68,16 @@
                     obj.penalty = penalty;
                 end
             end
+
+            %% Downwards compatability / set default values
+            %TODO: maybe move into set method for parameters
+            if numel(obj.parameters) < 1
+                obj.parameters{1} = 0;
+            end
+
+            if numel(obj.parameters) < 2
+                obj.parameters{2} = 1;
+            end
 
         end       
 

optimization/optimizer/matRad_OptimizerIPOPT.m

@@ -64,17 +64,17 @@
             obj.options.print_options_documentation   = 'no';
 
             % Termination (C.2)
-            obj.options.tol                           = 1e-8; % (Opt1)
-            obj.options.dual_inf_tol                  = 1;    % (Opt2)
+            obj.options.tol                           = 1e-10; % (Opt1)
+            obj.options.dual_inf_tol                  = 1e-4;    % (Opt2)
             obj.options.constr_viol_tol               = 1e-4; % (Opt3)
             obj.options.compl_inf_tol                 = 1e-4; % (Opt4), Optimal Solution Found if (Opt1),...,(Opt4) fullfiled
 
-            obj.options.acceptable_iter               = 3;    % (Acc1)
+            obj.options.acceptable_iter               = 5;    % (Acc1)
             obj.options.acceptable_tol                = 1e10; % (Acc2)
-            obj.options.acceptable_constr_viol_tol    = 1e10; % (Acc3)
+            obj.options.acceptable_constr_viol_tol    = 1e-2; % (Acc3)
             obj.options.acceptable_dual_inf_tol       = 1e10; % (Acc4)
             obj.options.acceptable_compl_inf_tol      = 1e10; % (Acc5)
-            obj.options.acceptable_obj_change_tol     = 1e-3; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled
+            obj.options.acceptable_obj_change_tol     = 1e-4; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled
 
             obj.options.max_iter                      = matRad_cfg.propOpt.defaultMaxIter;
             obj.options.max_cpu_time                  = 3000;