Full RO added and it tenatitvely works! Still need to verify a few thigns, but at least the function runs through entire thing.

WiscPlanPhotonkV125/matlab_frontend/FullRO_NLP_optimizer_v3.m

@@ -0,0 +1,580 @@
+
+
+% paths.patient_dir
+% paths.Goal_dir (previously called DP_dir)
+% paths.patient
+% paths.goalsName
+
+% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
+% orthoslice(D_full, [0,70])
+
+function [D_full, w_fin, Geometry, optGoal] = FullRO_NLP_optimizer_v3(varargin)
+% This function performs the beamlet optimization
+% [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
+% 
+% Inputs:
+%       () OR
+%       (Pat_path, path2goal) OR
+%       (Pat_path, path2goal, beamlet_weights)
+%   Pat_path, path2goal = strings to patient folder and optimal goals
+%   beamlet_weights = initial beamlet weights
+%   
+% Outputs:
+%       full dose image dose: [D_full, w_fin, Geometry, optGoal]
+% 
+% Made by Peter Ferjancic 1. May 2018
+% Last updated: 1. April 2019
+
+if nargin<2
+    load('WiscPlan_preferences.mat')
+    [Pat_path] = uigetdir([WiscPlan_preferences.patientDataPath ], 'Select Patient folder');
+    [Goal_file,Goal_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select OptGoal file');
+    
+    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
+    path2goal = [Goal_path, Goal_file];
+    
+    [path2scBeams] = uigetdir([WiscPlan_preferences.patientDataPath ], 'Select folder with scenario beamlets');
+else
+    Pat_path = varargin{1};
+    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
+    path2goal = varargin{2};
+    [Goal_path,Goal_file,ext] = fileparts(path2goal);
+    path2scBeams = varargin{3};
+end
+
+dialogue_box = 'no'
+switch dialogue_box
+    case 'yes'
+        str = inputdlg({'N of iterations for initial calc', 'N of iterations for full calc', ...
+            'Use pre-existing NLP_result to initiate? (y/n)'}, 'input', [1,35], {'100000', '500000', 'n'});
+        N_fcallback1 = str2double(str{1}); % 100000  is a good guesstimate
+        N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
+        pre_beamWeights = str{3};
+    case 'no'
+        disp('dialogue box skipped')
+        N_fcallback1 = 1e5;
+        N_fcallback2 = 2e6; % 500000;
+        pre_beamWeights = 'n';
+end
+
+
+switch pre_beamWeights
+    case 'y'
+        [NLP_file,NLP_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select NLP_result file');
+        load([NLP_path, NLP_file])
+        w_beamlets = NLP_result.weights;
+        
+        load([path2scBeams '\scenario1\all_beams.mat'])
+        if numel(all_beams) ~= numel(w_beamlets)
+            error('Provided weight number does not match beamlet number!')
+        end
+    case 'n'
+        disp('Initial beam weights will be calculated.')
+end
+
+%% PROGRAM STARTS HERE
+% - no tocar lo que hay debajo -
+fprintf('starting NLP optimization process... \n')
+
+% % -- LOAD GEOMETRY, GOALS, BEAMLETS --
+load(path2geometry)
+load(path2goal)
+%%%% FIGURE out this part for the full RO appraoch.
+% [beamlets, numBeamlet] = get_beamlets(Geometry, Pat_path, OptGoals);
+% [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, Pat_path);
+
+
+%% -- OPTIMIZATION TARGETS --
+% -- make the optimization optGoal structure --
+for i_goal = 1:size(OptGoals.goals,1)
+        
+    if isfield(OptGoals.data{i_goal}, 'SupVox_num')
+        SupVox_num = OptGoals.data{i_goal}.SupVox_num;
+    else
+        answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
+        SupVox_num = str2double(answer{1})
+    end
+    optGoal{i_goal} = OptGoals.data{i_goal};
+    optGoal{i_goal}.sss_scene_list = OptGoals.sss_scene_list;
+    optGoal{i_goal}.maxModulation = OptGoals.maxModulation;
+    optGoal{i_goal}.BeamSmoothMax = OptGoals.BeamSmoothMax;
+    
+    optGoal{i_goal}.optFuncNum = OptGoals.optFuncNum;
+    
+    % modulation
+    
+    if strcmp( '-||-' , optGoal{i_goal}.ROI_name)
+%         optGoal{i_goal}.ROI_idx_old = optGoal{i_goal-1}.ROI_idx_old; % copy old index data
+%         optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
+        optGoal{i_goal}.ROI_idx = optGoal{i_goal-1}.ROI_idx
+
+        if isfield(OptGoals.data{i_goal}, 'wgt_map')
+            optGoal{i_goal}.vox_wgt = optGoal{i_goal-1}.vox_wgt;
+        end
+        optGoal{i_goal}.beamlets_pruned = optGoal{i_goal-1}.beamlets_pruned;
+
+    else
+        switch SupVox_num
+            case 0
+                error('You should always specify supervoxels now')
+            % if not supervoxel, just select provided ROI_idx
+%             optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
+
+            otherwise
+            % -- if supervoxel, merge given columns
+            % - make supervoxel map
+            switch OptGoals.goals{i_goal, 3}
+                case 'Fixed dose'
+                    mask = zeros(OptGoals.data{i_goal}.imgDim);
+                    mask(OptGoals.data{i_goal}.ROI_idx) = 1;
+                    
+                case 'Dose map'
+                    mask = zeros(OptGoals.data{i_goal}.imgDim);
+                    mask(OptGoals.data{i_goal}.ROI_idx) = OptGoals.data{i_goal}.D_final;
+                    mask2 = round(mask/3)*3;
+            end
+            % group superpixels
+            superMask = superpix_group(mask, SupVox_num, 'no');
+            optGoal{i_goal}.superMask = superMask;
+%             orthoslice(superMask)
+        end
+    end
+end
+
+% -- make them robust --
+RO_params=0;
+[optGoal, numBeamlet] = make_fullRO_optGoal(optGoal, RO_params, path2scBeams);
+% save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
+
+% -- get beamlet indeces --
+load([path2scBeams '\scenario1\all_beams.mat'])
+Nbeamlets = all_beams{1}.Mxp;   % number of beamlets in a beam - usually 64 or 32
+weightTable = zeros(3,Nbeamlets);
+for ind_bmlt = 1:numel(all_beams)
+    bLet_idx.y(ind_bmlt) = floor(all_beams{1, ind_bmlt}.num/Nbeamlets)+1;
+    bLet_idx.x(ind_bmlt) = rem(all_beams{1, ind_bmlt}.num, Nbeamlets)+1;
+    weightTable(bLet_idx.y(ind_bmlt),bLet_idx.x(ind_bmlt)) = 1;
+end
+bLet_idx.idx = find(weightTable>0);
+bLet_idx.Nbeamlets = Nbeamlets;
+disp('.')
+
+% -- CALLBACK OPTIMIZATION FUNCTION --
+fun1 = @(x) get_penalty(x, optGoal_beam, bLet_idx);
+fun2 = @(x) get_penalty(x, optGoal,     bLet_idx);
+
+% -- OPTIMIZATION PARAMETERS --
+% define optimization parameters
+A = [];
+b = [];
+Aeq = [];
+beq = [];
+lb = zeros(1, numBeamlet);
+% lb_beam = zeros(1, numBeam);
+ub = [];
+nonlcon = [];
+
+% define opt limits, and make it fmincon progress
+options = optimoptions('fmincon');
+options.MaxFunctionEvaluations = N_fcallback1;
+options.Display = 'iter';
+options.PlotFcn = 'optimplotfval';
+% options.UseParallel = true;
+options.UseParallel = false;
+% options.OptimalityTolerance = 1e-9;
+
+%% -- INITIALIZE BEAMLET WEIGHTS --
+switch pre_beamWeights
+    case 'y'
+        % should have been assigned previously.
+        disp('Provided beamlet weights used for initial comparison')
+    case 'n'
+        % if initial beamlet weights are not provided, get quick estimate
+%         fprintf('\n running initial optimizer:')
+        % initialize beamlet weights, OR
+        w0 = ones(numBeamlet,1);
+%         w0 = mean(optGoal{1}.D_final(optGoal{1}.ROI_idx) ./ (optGoal{1}.beamlets_pruned*w0+0.1)) * w0; % old
+        w0 = mean(optGoal{1}.D_final(:)) ./ mean(optGoal{1}.nrs{1}.sss{1}.rgs{1}.beamlets_pruned*w0+0.05) * w0;
+        w_beamlets = double(w0);
+
+        % -- GET BEAM WEIGHTS --
+%         tic
+%         w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb,ub,nonlcon,options);
+%         fprintf('  done!:')
+%         t=toc;
+%         disp(['Optimization time for beams = ',num2str(t)]);
+% 
+%         w_beamlets = ones(numBeamlet,1);
+%         numBeam=numel(unique(beam_i_list));
+%         for beam_i = 1:numBeam % assign weights to beamlets
+%             % beamlets from same beam get same initial weights
+%             w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
+%         end
+end
+
+%% FULL OPTIMIZATION
+
+% -- GET FULL BEAMLET WEIGHTS --
+options.MaxFunctionEvaluations = N_fcallback2;
+tic
+fprintf('\n running full optimizer:')
+w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
+fprintf('  done!:')
+t=toc;
+disp(['Optimization time for beamlets = ',num2str(t)]);
+
+%% evaluate the results
+[beamlets, numBeamlet] = get_beamlets_fullRO(optGoal{1}.imgDim, [path2scBeams '\scenario1'], optGoal);
+D_full = reshape(beamlets * w_fin, size(Geometry.data));
+
+%% save outputs
+NLP_result.dose = D_full;
+NLP_result.weights = w_fin;
+NLP_result.sss_scene_list = optGoal{1}.sss_scene_list;
+NLP_result.maxModulation = OptGoals.maxModulation;
+NLP_result.BeamSmoothMax = OptGoals.BeamSmoothMax;
+save([path2scBeams, '\FRO_result_' Goal_file '.mat'], 'NLP_result');
+
+
+
+% plot_DVH(Geometry, D_full)
+DQM_DVH(Geometry, D_full)
+colorwash(Geometry.data, D_full, [500, 1500], [0, 90]);
+
+
+end
+
+%% support functions
+% ---- PENALTY FUNCTION ----
+function penalty = get_penalty(x, optGoal, bLet_idx)
+    % this function gets called by the optimizer. It checks the penalty for
+    % all the robust implementation and returns the worst result.
+    
+    NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
+    fobj = zeros(NumScenarios,numel(optGoal)+2);  
+    sc_i = 1;
+    
+
+    
+    for nrs_i = 1:optGoal{1}.NbrRandScenarios 
+        for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
+            for rgs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
+%                 fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx);
+                fobj(sc_i,:) = eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx);
+                sc_i = sc_i + 1;
+            end
+        end
+    end
+    
+        
+    switch optGoal{1}.optFuncNum
+        case 1 % "RO max"
+            penalty=max(sum(fobj,2));
+        case 2 % "RO mean"
+            penalty=mean(sum(fobj,2));
+        case 3 % "RO objective-based"
+            penalty=sum(max(fobj,[],1));
+        case 4 % "Classical"
+            penalty=sum(fobj(1,:));
+    end 
+    % take the worst case penalty of evaluated scenarios
+%     penalty=max(fobj);
+    
+    % take the median worst case (stochastic robust)
+%     penalty=median(fobj);
+
+
+    
+end
+% ------ supp: penalty for single scenario ------
+function penArray = eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx)
+%     penalty = 0;
+    penArray = zeros(1,numel(optGoal)+2); % all the optGoals, plus modulation plus smoothing
+    % for each condition
+    for goal_i = 1:numel(optGoal)
+        switch optGoal{goal_i}.function
+            % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
+            case 'min'
+                % penalize if achieved dose is lower than target dose
+                d_penalty = 1.0e0 * sum(max(0, ...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)));
+            case 'max'
+                % penalize if achieved dose is higher than target dose
+                d_penalty = 1.0e0 * sum(max(0, ...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target)));
+            case 'min_sq'
+                % penalize if achieved dose is lower than target dose
+                temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e0 * sum(temp1.*temp1);
+            case 'max_sq'
+                % penalize if achieved dose is higher than target dose
+                temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e0 * sum(temp1.*temp1);
+            case 'min_step'
+                % penalize if achieved dose is lower than target dose
+                temp1=-min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e1 * (sum(temp1.*temp1) + 1.0e3* sum(temp1>0)/numel(temp1));
+            case 'max_step'
+                % penalize if achieved dose is higher than target dose
+                temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e1 * (sum(temp1.*temp1) + 1.0e3* sum(temp1>0)/numel(temp1));
+            case 'LeastSquare'
+                % penalize with sum of squares any deviation from target
+                % dose
+                temp1 = (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) - ...
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target;
+                d_penalty = 1.0e0* sum(temp1.^2);
+            case 'min_perc_Volume'
+                % penalize by amount of volume under threshold
+                perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) > 0)) ...
+                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
+                d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
+                
+            case 'max_perc_Volume'
+                % penalize by amount of volume under threshold
+                perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) < 0)) ...
+                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
+                d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
+            case 'min_sq_voxwgt'
+                % penalize if achieved dose is lower than target dose
+                temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e0 * sum(temp1.*temp1.* optGoal{goal_i}.vox_wgt');
+            case 'max_sq_voxwgt'
+                % penalize if achieved dose is lower than target dose
+                temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
+                d_penalty = 1.0e0 * sum(temp1.*temp1.* optGoal{goal_i}.vox_wgt'); 
+        end
+%         penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
+        penArray(goal_i) = d_penalty * optGoal{goal_i}.opt_weight;
+    end
+    
+    %% add modulation penalty
+    if true
+        if isfield(optGoal{goal_i}, 'maxModulation')
+            max_modulation = optGoal{goal_i}.maxModulation;
+        else
+            error('no Max modulation parameter enterd!')
+            max_modulation = 5;
+        end
+        
+        mod_pen_weight = 1.0e10;
+        % calc the penalty
+        mod_excess = max(0, x-max_modulation*mean(x));
+%         mod_excess = max(0, x-max_modulation*mean(x(x>1)));
+        mod_pen1 = mod_pen_weight*(sum(mod_excess) + numel(mod_excess)* any(mod_excess));
+%         penalty = penalty + mod_pen1;
+        penArray(end-1) = mod_pen1;
+    end
+    
+    %% add penlty for single off beamlets - version 2
+    if false
+        if isfield(optGoal{goal_i}, 'BeamSmoothMax')
+            BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
+        else
+            error('no Max beam smooth parameter enterd!')
+            BeamSmoothMax = 1e6;
+        end
+        mod_pen_weight = BeamSmoothMax; %1.0e6
+        x_down = zeros(size(x));
+        x_down(2:end) = x(1:end-1);
+        x_down2 = zeros(size(x));
+        x_down2(3:end) = x(1:end-2);
+        x_up = zeros(size(x));
+        x_up(1:end-1) = x(2:end);
+        x_up2 = zeros(size(x));
+        x_up2(1:end-2) = x(3:end);
+        mod_pen2 = mod_pen_weight*sum((x_down+x_up-2*x).^2 + 0.5*(x_down2+x_up2-2*x).^2)/(mean(x)^(2)*numel(x));
+%         penalty = penalty + mod_pen2;
+        penArray(end) = mod_pen2;
+    end
+    %% add penlty for single off beamlets - version 3
+    if true
+        if isfield(optGoal{goal_i}, 'BeamSmoothMax')
+            BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
+        else
+            error('no Max beam smooth parameter enterd!')
+            BeamSmoothMax = 1e6;
+        end
+        mod_pen_weight = BeamSmoothMax; %1.0e6
+        
+        % make 2D beamletWeight map
+        maxY = max(bLet_idx.y);
+        tabula = zeros(maxY, bLet_idx.Nbeamlets);
+        tabula(bLet_idx.idx) = x;
+        % for each index calculate laplace
+        myWeights = 1/12*[1,3,1; 1,-12,1; 1,3,1];
+        i=2:maxY-1;
+        j=2:bLet_idx.Nbeamlets-1;
+        tabula2(i,j) = ...
+             myWeights(1)*tabula(i-1,j-1) + myWeights(4)*tabula(i-1,j) + myWeights(7)*tabula(i-1,j+1)...
+            +myWeights(2)*tabula(i,j-1)  + myWeights(8)*tabula(i,j+1) ...
+            +myWeights(3)*tabula(i+1,j-1) + myWeights(6)*tabula(i+1,j) + myWeights(9)*tabula(i+1,j+1)...
+            +myWeights(5)*tabula(i,j);
+        tabula2(1,j) = ...
+             myWeights(2)*tabula(1,j-1) + myWeights(8)*tabula(1,j+1)...
+            +myWeights(3)*tabula(2,j-1) + myWeights(6)*tabula(2,j) + myWeights(9)*tabula(2,j+1)...
+            +myWeights(5)*tabula(1,j);
+        tabula2(maxY,j) =...
+             myWeights(1)*tabula(maxY-1,j-1) + myWeights(4)*tabula(maxY-1,j) + myWeights(7)*tabula(maxY-1,j+1)...
+            +myWeights(2)*tabula(maxY,j-1)  + myWeights(8)*tabula(maxY,j+1) ...
+            +myWeights(5)*tabula(maxY,j);
+        % make sum of squares
+        mod_pen2 = mod_pen_weight*sum((tabula2(:)).^2)/(mean(x.^2)*numel(x));
+%         penalty = penalty + mod_pen2;
+        penArray(end) = mod_pen2;
+    end
+end
+
+
+
+% ---- MAKE ROI ROBUST ----
+function [optGoal, numBeamlet] = make_fullRO_optGoal(optGoal, RO_params, path2scBeams);
+    % take regular optimal goal and translate it into several robust cases
+    
+    % RO_params - should have the information below
+    % nrs - random scenarios
+    % sss - system setup scenarios
+    % rgs - random range scenarios
+    
+    % X - X>0 moves image right
+    % Y - Y>0 moves image down
+    % Z - in/out.
+    
+    shift_X = 2; % vox of shift
+    shift_Y = 2; % vox of shift
+    shift_Z = 1; % vox of shift
+    nrs_scene_list={[0,0,0]};
+
+    
+    % ----====#### CHANGE ROBUSTNESS HERE ####====----
+    if isfield(optGoal{1}, 'sss_scene_list')
+        sss_scene_list = optGoal{1}.sss_scene_list;
+    else
+        error('New OptGoals should no longer reach this code!')
+%         sss_scene_list={[0,0,0], [-shift_Y,0,0], [shift_Y,0,0], [0,-shift_X,0], [0,shift_X,0], [0,0,-shift_Z], [0,0,shift_Z]};
+%         optGoal{1}.sss_scene_list = sss_scene_list;
+    end
+%     sss_scene_list={[0,0,0]};
+
+    % ----====#### CHANGE ROBUSTNESS HERE ####====----
+    rgs_scene_list={[0,0,0]};
+    for i = 1:numel(optGoal)
+        optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
+        optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
+        optGoal{i}.NbrRangeScenarios    =numel(rgs_scene_list);
+    end
+    
+    % for each scenario
+    for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
+    % load new beamlets
+    
+    
+    patient_dir = [path2scBeams '\scenario' num2str(sss_i)];
+    [beamlets, numBeamlet] = get_beamlets_fullRO(optGoal{1}.imgDim, patient_dir, optGoal);
+    
+    disp(['-  scenario ' num2str(sss_i) ' beamlets loaded.'])
+        for goal_i = 1:numel(optGoal)
+            % get target
+            idx=optGoal{goal_i}.ROI_idx;
+            targetImg1=zeros(optGoal{goal_i}.imgDim);
+            targetImg1(idx)=1;
+            % get beamlets
+            for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios          % num. of random scenarios
+                % modify target and beamlets
+%                 targetImg3=targetImg1;
+                % beamlets stay the same
+            
+                for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rgs
+                    % modify target and beamlets
+                    targetImg4=targetImg1;
+                    % beamlets stay the same
+                    
+                    % function to put in superMask and get out updated optGoal
+                    [beamlets_pruned, targetD] = superVoxMerge(optGoal{goal_i}, beamlets);
+                    % needs to return "beamlets_pruned" and "target" arrays
+
+                    % save to optGoal output
+%                     optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx          = ROI_idx; % some old code, not used?
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned    = beamlets_pruned; % beamlets, only within target and supervoxeled
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target             = targetD'; % target does that we try to achieve
+                    
+                end
+            end
+        end
+    end
+end
+
+%% ---- Supervoxel merge optGoals ----
+function [beamlets_pruned, targetD]  = superVoxMerge(optGoal, beamlets)
+    superMask = optGoal.superMask;
+    superVoxList = unique(superMask);
+    superVoxList = superVoxList(superVoxList>0);
+
+    optGoal.ROI_idx_full = optGoal.ROI_idx; % copy old index data
+    % optGoal.ROI_idx = zeros(numel(superVoxList), 1);
+    optGoal.opt_weight = optGoal.opt_weight * numel(optGoal.ROI_idx_full)/numel(optGoal.ROI_idx);
+
+
+    if isfield(optGoal, 'wgt_map')
+        tabula_wgtmap = zeros(size(superMask));
+        tabula_wgtmap(optGoal.optGoal.ROI_idx_full) = optGoal.wgt_map;
+    end
+
+    h_w1 = waitbar(0, 'merging superboxels');
+    
+    dose_tabula = zeros(optGoal.imgDim);
+    dose_tabula(optGoal.ROI_idx_full) = optGoal.D_final;
+    
+    for i_supVox = 1:numel(superVoxList)
+        waitbar(i_supVox/numel(superVoxList), h_w1)
+        supVox_idx = superVoxList(i_supVox);
+        idxList = find(superMask == supVox_idx);
+        % merge beamlet doses
+        beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
+
+        targetD(i_supVox) = mean(dose_tabula(idxList),1);
+        
+        if isfield(optGoal, 'vox_wgt')
+            optGoal.vox_wgt(i_supVox) = sum(tabula_wgtmap(idxList));
+        end
+
+        % -- make new indeces
+    %     optGoal.ROI_idx(i_supVox) = idxList(1);
+    end
+    close(h_w1)
+
+end
+
+%% ------ supp: RO case SSS ------
+function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
+    % translate the target image 
+    
+    targetImg3 = imtranslate(targetImg2,sss_scene_shift);
+    
+    % now we need to figure out if any target voxels fell out during the
+    % shift
+    imgValid = imtranslate(targetImg3,-sss_scene_shift);
+    imgInvalid = (targetImg2-imgValid);
+    
+    idx_1 = find(targetImg2);
+    idx_2 = find(imgInvalid);
+    
+    [idxValid,ia] = setdiff(idx_1,idx_2);
+    
+    [C,ia, ib] = intersect(idx_1,idxValid);
+    
+end
+
+
+
+

WiscPlanPhotonkV125/matlab_frontend/dicomrt/dicomrt2geometry.m

@@ -151,6 +151,7 @@ if downsample_flag == true
     % Geometry.data = resamplegrid(Geometry.data, 0.5);
     % note that ct_xmesh, ct_ymesh, ct_zmesh is already downsampled
     Geometry.start(1:2) = Geometry.start(1:2) + Geometry.voxel_size(1:2) / 2;
+    Geometry.voxel_size = Geometry.voxel_size/0.5;
 end
 
 Geometry.rhomw = CT2dens(Geometry.data, 'pinn');

WiscPlanPhotonkV125/matlab_frontend/fullRO_mergeBeams.m

@@ -2,11 +2,7 @@
 function fullRO_mergeBeams(Geometry)
 % this function loops through a full RO folder and merges all the beamlets
 
-    if isfield(Geometry, 'num_batches')
-        str = num2str(obj.handles.hSVPS.Geometry.num_batches);
-    else
-        str = input('Enter num of beam batches: ','s');
-    end
+    str = input('Enter num of beam batches: ','s');
     
     % get the full RO folder
     beamlet_dir = uigetdir([Geometry.patient_dir ], 'Select full RO folder' );

WiscPlanPhotonkV125/matlab_frontend/get_beamlets_fullRO.m

@@ -0,0 +1,87 @@
+
+
+function [beamlets, numBeamlet] = get_beamlets_fullRO(imgDim, patient_dir, optGoal)
+% this function loads and returns beamlets and joined beams
+
+% ROI_idxList = [];
+% for i = 1:numel(OptGoals.data)
+%     ROI_idxList = [ROI_idxList; OptGoals.data{i}.ROI_idx];
+% end
+% ROI_idxList = unique(ROI_idxList);
+
+tabula = zeros(imgDim);
+for i = 1:numel(optGoal)
+    tabula(optGoal{i}.ROI_idx) = 1;
+end
+ROI_idxList = find(tabula>0);
+
+
+%% add idx lists for each goal!
+
+beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
+% beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
+beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'peter', ROI_idxList);
+
+numVox  = prod(imgDim);
+numBeamlet = size(beamlet_cell_array,2);
+batchSize = 200;
+beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
+% beamlets = beamlet_cell_array;
+
+end
+
+% ---- GET BEAMLETS ----
+function beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
+    wbar1 = waitbar(0, 'Creating beamlet array');
+    numBeam = size(beamlet_cell_array,2);
+    beamlets = sparse(0, 0);
+    for beam_i=1:numBeam
+        % for each beam define how much dose it delivers on each voxel
+        idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
+        [ROI_idxIntersect, ia, ~] = intersect (idx, ROI_idxList);
+%         if ~isequal(ROI_idxIntersect,idx)
+%             warning('this')
+%         end
+        
+%         ROI_idxIntersect = beamlet_cell_array{1, beam_i}.non_zero_indices;
+        
+        if isempty(ROI_idxIntersect)
+            warning(['Beamlet ' num2str(beam_i) ' is empty!'])
+        end
+        
+        % break the beamlets into multiple batches
+        if rem(beam_i, batchSize)==1;
+            beamlet_batch = sparse(numVox, batchSize);
+            beam_i_temp=1;
+        end
+
+%         beamlet_batch(idx, beam_i_temp) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
+        beamlet_batch(ROI_idxIntersect, beam_i_temp) = beamlet_cell_array{1, beam_i}.non_zero_values(ia);
+        waitbar(beam_i/numBeam, wbar1, ['Adding beamlet array: #', num2str(beam_i)])
+
+        % add the batch to full set when filled
+        if rem(beam_i, batchSize)==0;
+            beamlets =[beamlets, beamlet_batch];
+        end
+        % crop and add the batch to full set when completed
+        if beam_i==numBeam;
+            beamlet_batch=beamlet_batch(:, 1:beam_i_temp);
+            beamlets =[beamlets, beamlet_batch];
+        end
+        beam_i_temp=beam_i_temp+1;
+
+    end
+    close(wbar1)
+
+end
+function show_joint_beamlets(beamlets, IMGsize, Beam_list)
+    % this function overlays and plots multiple beamlets. The goal is to
+    % check whether some beamlets are part of the same beam manually.
+    
+    beam=sum(beamlets(:, Beam_list), 2);
+    
+    beamImg=reshape(full(beam), IMGsize);
+        
+    orthoslice(beamImg)
+    
+end

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -70,7 +70,7 @@ end
 for i = 1: numel(Geometry.ROIS)
     handles.uitable1.ColumnFormat{2}{i} = Geometry.ROIS{i}.name;
 end
-handles.uitable1.ColumnFormat{2}{numel(Geometry.ROIS)} = '-||-';
+handles.uitable1.ColumnFormat{2}{numel(Geometry.ROIS)+1} = '-||-';
 % == populate function options ==
 handles.uitable1.ColumnFormat{5}{1} = 'min';
 handles.uitable1.ColumnFormat{5}{2} = 'max';

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -11,6 +11,7 @@ classdef XTPS < handle
 % See also SliceViewerPanel
 %
 % Author Xiaohu Mo
+% Robust optimization by Peter Ferjancic
 
     properties
         handles
@@ -76,9 +77,9 @@ classdef XTPS < handle
             
             % Full RO menu
             obj.handles.hFullROMenu   = uimenu(obj.handles.hMainFigure, 'Label', 'Full RO');     
-            obj.handles.hFullRO_beamCalc = uimenu(obj.handles.hFullROMenu,  'Label', 'Full RO beam calc');
-            obj.handles.hFullRO_mergeBeams = uimenu(obj.handles.hFullROMenu,  'Label', 'mergeBeams');
-            obj.handles.hFullRO_Opt = uimenu(obj.handles.hFullROMenu,  'Label', 'Full RO optimization');
+            obj.handles.hFullRO_beamCalc = uimenu(obj.handles.hFullROMenu,  'Label', 'FRO beam calc');
+            obj.handles.hFullRO_mergeBeams = uimenu(obj.handles.hFullROMenu,  'Label', 'FRO mergeBeams');
+            obj.handles.hFullRO_Opt = uimenu(obj.handles.hFullROMenu,  'Label', 'FRO optimization');
             
             % Tools menu
             obj.handles.hToolsMenu   = uimenu(obj.handles.hMainFigure, 'Label', 'Tools');            
@@ -701,6 +702,8 @@ classdef XTPS < handle
         end
 %-------------------------------------------------------------------------------    
         function hFullRO_mergeBeams_Callback(obj, src, evt)
+            % loop through the scenario folders and merge and shift
+            % beamlets
             disp('Merging full RO beamlets')
             fullRO_mergeBeams(obj.handles.hSVPS.Geometry)
         end
@@ -709,6 +712,12 @@ classdef XTPS < handle
             % call function to do full RO optimization based on multiple
             % beamlets
             disp('Full RO optimization called')
+            Pat_path = [obj.handles.hSVPS.Geometry.patient_dir];
+            [Obj_file,Obj_path,indx] = uigetfile([obj.handles.hSVPS.Geometry.patient_dir '\matlab_files\*.mat'], 'Select OptGoal file' );
+            path2goal = [Obj_path, Obj_file];
+            
+            path2scBeams = uigetdir([obj.handles.hSVPS.Geometry.patient_dir ], 'Select folder with scenario beamlets');
+            [D_full, w_fin, Geometry, optGoal] = FullRO_NLP_optimizer_v3(Pat_path, path2goal, path2scBeams);
         end
 %-------------------------------------------------------------------------------
         function DoseModeDropdown_Callback(obj, src, evt)

view_beamlet.m

@@ -1,15 +1,17 @@
 
 
-% data = read_ryan_beamlets('F:\021_WiscPlan_data\Hippocampus_HD_2\batch_dose.bin')
+blet1 = read_ryan_beamlets('F:\021_WiscPlan_data\DELETEME3\FRO_Beam1\scenario1\batch_dose.bin');
+blet2 = read_ryan_beamlets('F:\021_WiscPlan_data\DELETEME3\FRO_Beam1\scenario2\batch_dose.bin');
 
 
-idx = 3000;
+idx = 110;
 
-tabula = zeros([data{idx}.x_count, data{idx}.y_count, data{idx}.z_count]);
+tabula = zeros([blet1{idx}.x_count, blet1{idx}.y_count, blet1{idx}.z_count]);
+tabula(blet1{idx}.non_zero_indices)   = tabula(blet1{idx}.non_zero_indices)+blet1{idx}.non_zero_values;
+orthoslice(tabula);
 
-tabula(data{idx}.non_zero_indices)   = tabula(data{idx}.non_zero_indices)+data{idx}.non_zero_values;
-tabula(data{idx+1}.non_zero_indices) = tabula(data{idx+1}.non_zero_indices)+data{idx+1}.non_zero_values;
+tabula = zeros([blet2{idx}.x_count, blet2{idx}.y_count, blet2{idx}.z_count]);
+tabula(blet2{idx}.non_zero_indices)   = tabula(blet2{idx}.non_zero_indices)+blet2{idx}.non_zero_values;
+orthoslice(tabula);
 
-tabula(data{idx+2}.non_zero_indices) = tabula(data{idx+2}.non_zero_indices)+data{idx+2}.non_zero_values;
 
-orthoslice(tabula);