Dev mix mod effect prescription bug fix

IO/matRad_readNRRD.m

@@ -27,23 +27,25 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+matRad_cfg = MatRad_Config.instance();
+
 % Open file.
 hFile = fopen(filename, 'r');
 if hFile <= 0
-    error('Could not open NRRD file!');
+    matRad_cfg.dispError('Could not open NRRD file!');
 end
 cleaner = onCleanup(@() fclose(hFile));
 
 %% Determine NRRD Version
 nrrdLine = fgetl(hFile);
 regTokens = regexp(nrrdLine,'NRRD00(?:\.)?0([1-5])','tokens');
 if isempty(regTokens)
-    error('Invalid Header line! Could not identify NRRD version!');
+    matRad_cfg.dispError('Invalid Header line! Could not identify NRRD version!');
 end
 nrrdVersion = str2num(regTokens{1}{1});
 
 if nrrdVersion > 5
-    error('NRRD version > 5 not supported!');
+    matRad_cfg.dispError('NRRD version > 5 not supported!');
 end
 
 %% Read header
@@ -60,15 +62,15 @@
         %Parse the line
         lineContent = regexp(currentLine, '(.+):(=|\s)(.+)', 'tokens');
         if isempty(lineContent)
-            warning(['Could not parse line: "' lineContent '"']);
+            matRad_cfg.dispWarning(['Could not parse line: "' currentLine '"']);
         elseif isequal(lineContent{1}{2},' ') %space after colon refers to "field"
             nrrdMetaData.fields{end+1,1} = lineContent{1}{1}; %Fieldname
             nrrdMetaData.fields{end,2} = lineContent{1}{3}; %Information
         elseif isequal(lineContent{1}{2},'=') %= after colon refers to key-value pair
             nrrdMetaData.keys{end+1,1} = lineContent{1}{1}; %Key
             nrrdMetaData.keys{end,2} = lineContent{1}{3}; %Value
         else
-            warning(['Could not parse line: "' lineContent '"']);
+            matRad_cfg.dispWarning(['Could not parse line: "' currentLine '"']);
         end        
     end
     currentLine = fgetl(hFile);
@@ -85,7 +87,7 @@
         endianFieldIx = find(ismember(nrrdMetaData.fields(:,1),'endian'));
         if ~isempty(endianFieldIx)
             if ~isequal(nrrdMetaData.fields{endianFieldIx,2},'little') && ~isequal(nrrdMetaData.fields{endianFieldIx,2},'big')
-                error(['Datatype is ' datatype ', thus endian information is required but could not be interpreted!']);
+                matRad_cfg.dispError(['Datatype is ' datatype ', thus endian information is required but could not be interpreted!']);
             end;
             %Now we compare the file endian to the system endian
             %First acquire system endian
@@ -101,36 +103,36 @@
                 doSwapBytes = true;
             end
         else
-            error(['Datatype is ' datatype ', thus endian information is required but could not be found!']);
+            matRad_cfg.dispError(['Datatype is ' datatype ', thus endian information is required but could not be found!']);
         end
     end
     metadata.datatype = datatype;
 
 else
-    error('Could not find required "type" field!');
+    matRad_cfg.dispError('Could not find required "type" field!');
 end
 
 %Check for the always required image dimension
 dimFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'dimension')); 
 if ~isempty(dimFieldIx)
     [metadata.dimension,success] = str2num(nrrdMetaData.fields{dimFieldIx,2});
     if ~success
-        error('Could not read required dimension field');
+        matRad_cfg.dispError('Could not read required dimension field');
     end
 else
-    error('Could not find required "dimension" field!');
+    matRad_cfg.dispError('Could not find required "dimension" field!');
 end
 
 %Check for size / dim length
 sizeFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'sizes')); 
 if ~isempty(sizeFieldIx)
     sizes = textscan(nrrdMetaData.fields{sizeFieldIx,2},'%d');
     if numel(sizes{1}) ~= metadata.dimension || ~all(sizes{1} > 0) 
-        error('Incorrect size definition!');
+        matRad_cfg.dispError('Incorrect size definition!');
     end
     metadata.cubeDim = sizes{1}';
 else
-    error('Could not find required "dimension" field!');
+    matRad_cfg.dispError('Could not find required "dimension" field!');
 end
 
 %Check for resolution
@@ -140,7 +142,7 @@
 if ~isempty(spacingFieldIx)
     resolutions = textscan(nrrdMetaData.fields{spacingFieldIx,2},'%f');
     if numel(resolutions{1}) ~= metadata.dimension
-        error('Incorrect spacings definition');
+        matRad_cfg.dispError('Incorrect spacings definition');
     end
     metadata.resolution = resolutions{1}';    
 
@@ -167,14 +169,14 @@
         currentAxis = find(vectors{c});
 
         if numel(find(vectors{c})) ~= 1
-            error('Sorry! We currently only support spaces with cartesian basis!'); 
+            matRad_cfg.dispError('Sorry! We currently only support spaces with cartesian basis!'); 
         end        
         metadata.axisPermutation(c) = currentAxis*sign(vectors{c}(currentAxis));
         metadata.resolution(c) = vectors{c}(currentAxis);       
     end
 
 else
-    warning('No Resolution Information available');
+    matRad_cfg.dispWarning('No Resolution Information available');
 end
 
 %find the origin if we have one
@@ -202,7 +204,7 @@
 data_fileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'data file'));
 datafileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'datafile'));
 if ~isempty(data_fileFieldIx) || ~isempty(datafileFieldIx)
-    error('Sorry! We currently do not support detached data files!');
+    matRad_cfg.dispError('Sorry! We currently do not support detached data files!');
     %Proposed workflow:
     %check for data file    
     %close file
@@ -215,15 +217,15 @@
 %Check for encoding
 encodingFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'encoding'));
 if isempty(encodingFieldIx)
-    error('Could not find required "encoding" field!');
+    matRad_cfg.dispError('Could not find required "encoding" field!');
 end
 switch nrrdMetaData.fields{encodingFieldIx,2}
     case 'raw'
         cube = fread(hFile,prod(metadata.cubeDim), metadata.datatype);
     case {'txt','text','ascii'}
         cube = cast(fscanf(hFile,'%f'),metadata.datatype);
     case 'hex'
-        error('Sorry: NRRD hex file not yet supported!');
+        matRad_cfg.dispError('Sorry: NRRD hex file not yet supported!');
     case {'gz','gzip'}
         compressedByteArray = fread(hFile, inf, 'uint8');
 
@@ -242,7 +244,7 @@
                 javaByteOutputStream.close();
                 javaUnpackSuccessful = true;
             catch
-                warning('Java unpacking failed... using temporary files!');
+                matRad_cfg.dispWarning('Java unpacking failed... using temporary files!');
             end
         end
         if ~javaUnpackSuccessful
@@ -252,7 +254,7 @@
             hFileTmp = fopen(tmpFile, 'wb');
 
             if hFileTmp <= 0
-                error('Could not open temporary file for GZIP!');
+                matRad_cfg.dispError('Could not open temporary file for GZIP!');
             end
 
             fwrite(hFileTmp, compressedByteArray, 'uint8');
@@ -264,16 +266,16 @@
             %Read the uncompressed file
             hFileTmp = fopen(tmpName, 'rb');
             if hFileTmp <= 0
-                error('Could not open unpacked file!');
+                matRad_cfg.dispError('Could not open unpacked file!');
             end
             cleanTmpFile = onCleanup(@() fclose(hFileTmp));
 
             cube = fread(hFileTmp, prod(metadata.cubeDim), metadata.datatype);
         end
     case {'bz2','bzip2'}
-        error('Sorry: bzip compression not yet supported!');
+        matRad_cfg.dispError('Sorry: bzip compression not yet supported!');
     otherwise 
-        error(['Undefined NRRD encoding scheme: ' nrrdMetaData.encoding]);
+        matRad_cfg.dispError(['Undefined NRRD encoding scheme: ' nrrdMetaData.encoding]);
 end
 
 %maybe we need to correct the byte ordering (endian)
@@ -327,7 +329,7 @@
     case 'double'
         datatype = 'double';        
     otherwise
-        error('Could not identify datatype for NRRD data');
+        matRad_cfg.dispError('Could not identify datatype for NRRD data');
 end
 
 end

dicom/matRad_scanDicomImportFolder.m

@@ -40,8 +40,11 @@
 %% get all files in search directory
 
 % dicom import needs image processing toolbox -> check if available
+v = ver;
 if ~license('checkout','image_toolbox')
-    matRad_cfg.dispError('image processing toolbox and/or corresponding licence not available');
+    matRad_cfg.dispError('Image Processing Toolbox and/or corresponding license not available');
+elseif ~any(strcmp('Image Processing Toolbox', {v.Name}))
+    matRad_cfg.dispError('Image Processing Toolbox not installed');
 end
 
 fileList = matRad_listAllFiles(patDir);

matRad_fluenceOptimization.m

@@ -286,7 +286,7 @@
         optimizer = matRad_OptimizerIPOPT;
 end
 
-
+optimizer.options.max_iter = 10000;
 optimizer = optimizer.optimize(wInit,optiProb,dij,cst);
 
 wOpt = optimizer.wResult;

matRad_fluenceOptimizationJO.m

@@ -1,4 +1,4 @@
-function [resultGUI,optimizer] = matRad_fluenceOptimizationJO(dij,cst,pln,wInit)
+function [resultGUI,optimizer] = matRad_fluenceOptimizationJO(dij,cst,pln,wInit,loadwInit)
 % matRad inverse planning wrapper function
 % 
 % call
@@ -106,7 +106,11 @@
 else
     numOfModalities = 1;
 end
-if exist('wInit','var')
+
+if ~exist('wInit', 'var')
+    wInit = [];
+end
+if ~isempty(wInit) %can never happen now
     %do nothing as wInit was passed to the function
     matRad_cfg.dispInfo('chosen provided wInit!\n');   
 
@@ -364,7 +368,13 @@
         warning(['Optimizer ''' pln.propOpt.optimizer ''' not known! Fallback to IPOPT!']);
         optimizer = matRad_OptimizerIPOPT;
 end
-
+
+if exist('loadwInit', 'var')
+    if loadwInit ==1
+        load('wInit.mat');
+    end
+end
+optimizer.options.max_iter = 10000;
 optimizer = optimizer.optimize(wInit,optiProb,dij,cst);
 
 wOpt = optimizer.wResult;

optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m

@@ -313,3 +313,21 @@
     %Only implemented for first scenario now
     weightGradient = weightGradient + gProb{1};
 end
+    gradientChecker = 1;
+if gradientChecker == 1
+    f =  matRad_objectiveFunction(optiProb,w,dij,cst);
+    epsilon = 1e-8;
+    ix = unique(randi([1 numel(w)],1,5));
+
+    for i=ix
+        wInit = w;
+        wInit(i) = wInit(i) + epsilon;
+        fDel= matRad_objectiveFunction(optiProb,wInit,dij,cst);
+        numGrad = (fDel - f)/epsilon;
+        diff = (numGrad/weightGradient(i) - 1)*100;
+        fprintf(['grad val #' num2str(i) ' - rel diff numerical and analytical gradient = ' num2str(diff) '\n']);
+        %fprintf([' any nan or zero for photons' num2str(sum(isnan(glog{1}))) ',' num2str(sum(~logical(glog{1}))) ' for protons: ' num2str(sum(isnan(glog{2}))) ',' num2str(sum(~logical(glog{2}))) '\n']);
+    end
+
+end
+end

optimization/@matRad_OptimizationProblemMixMod/matRad_constraintFunctions.m

@@ -32,19 +32,43 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+matRad_cfg = MatRad_Config.instance();
+% DOSE PROJECTION
+bxidx = 1; %modality  bixel index
 
-% get current dose / effect / RBExDose vector
-optiProb.BP.compute(dij,w);
-d = optiProb.BP.GetResult();
+% Obtain cumulative dose cube 
+[d{1}]    = deal(zeros(dij.doseGrid.numOfVoxels,1));
+
+for mod = 1: length(dij.original_Dijs)
+    wt = [];
+    % split the w for current modality
+    STrepmat = (~dij.spatioTemp(mod) + dij.spatioTemp(mod)*dij.numOfSTscen(mod));
+    wt = reshape(w(bxidx: bxidx+STrepmat*dij.original_Dijs{mod}.totalNumOfBixels-1),[dij.original_Dijs{mod}.totalNumOfBixels,STrepmat]);
+    % get current dose / effect / RBExDose vector
+    optiProb.BP.compute(dij.original_Dijs{mod},wt);
+    dt = optiProb.BP.GetResult();
+
+    % also get probabilistic quantities (nearly no overhead if empty)
+    [dExp,dOmega] = optiProb.BP.GetResultProb();                        % NOTE: not sure what exactly to do for the dOmegas 
+
+    %Index of nxt modality
+    bxidx = bxidx + STrepmat*dij.original_Dijs{mod}.totalNumOfBixels;
+    % Accumulat Dose for all scenarios FOR FUTURE REVIEW ON HOW TO COMBINE
+    % DIFFFERENT UNCERTAINTY SCENARIOS FOR DIFFERENT MODALITIES
+    % currently for ST optimization 
+    for scen = 1:numel(dt)
+         d{scen} = d{scen} + sum(dt{scen}.*dij.STfractions{mod}',2);
+    end
+end
 
 % get the used scenarios
 useScen  = optiProb.BP.scenarios;
 scenProb = optiProb.BP.scenarioProb;
 
 % retrieve matching 4D scenarios
-fullScen = cell(ndims(d),1);
-[fullScen{:}] = ind2sub(size(d),useScen);
-contourScen = fullScen{1};
+fullScen      = cell(ndims(dt),1);
+[fullScen{:}] = ind2sub(size(dt),useScen);
+contourScen   = fullScen{1};
 
 % Initializes constraints
 c = [];
@@ -64,21 +88,25 @@
          if isa(constraint,'DoseConstraints.matRad_DoseConstraint')
 
             % rescale dose parameters to biological optimization quantity if required
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter./sum([dij.STfractions{:}]));
+
             constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX);
 
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter.*sum([dij.STfractions{:}]));            
+
             % retrieve the robustness type
             robustness = constraint.robustness;
 
             switch robustness
                case 'none' % if conventional opt: just sum objectives of nominal dose
                    d_i = d{1}(cst{i,4}{1});
                    c = [c; constraint.computeDoseConstraintFunction(d_i)];
-
                case 'PROB' % if prob opt: sum up expectation value of objectives
 
                   d_i = dExp{1}(cst{i,4}{1});
                   c = [c; constraint.computeDoseConstraintFunction(d_i)];
-
                case 'VWWC'  % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR
                   contourIx = unique(contourScen);
                   if ~isscalar(contourIx)