Added resample Geometry functionality

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v3.m

@@ -1,543 +1,543 @@
-
-
-% 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] = 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];
-else
-    Pat_path = varargin{1};
-    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
-    path2goal = varargin{2};
-    [Goal_path,Goal_file,ext] = fileparts(path2goal);
-end
-
-dialogue_box = 'yes'
-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([Pat_path, '\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)
-[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
-    switch SupVox_num
-        case 0
-        % 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;
-                % group superpixels
-                superMask = superpix_group(mask, SupVox_num, 'no'); 
-            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;
-                superMask = superpix_group(mask, SupVox_num, 'no'); 
-                orthoslice(superMask)
-        end
-        
-        
-        
-        superVoxList = unique(superMask);
-        superVoxList = superVoxList(superVoxList>0);
-        
-        optGoal{i_goal}.ROI_idx_old = optGoal{i_goal}.ROI_idx; % copy old index data
-        optGoal{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
-        optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
-        
-        if isfield(OptGoals.data{i_goal}, 'wgt_map')
-            tabula_wgtmap = zeros(size(superMask));
-            tabula_wgtmap(OptGoals.data{i_goal}.ROI_idx) = OptGoals.data{i_goal}.wgt_map;
-        end
-        
-        
-        
-        h_w1 = waitbar(0, 'merging superboxels');
-        for i_supVox = 1:numel(superVoxList)
-            waitbar(i_supVox/numel(superVoxList), h_w1)
-            supVox_idx = superVoxList(i_supVox);
-            idxList = find(superMask == supVox_idx);
-            optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
-            
-            if isfield(OptGoals.data{i_goal}, 'wgt_map')
-                optGoal{i_goal}.vox_wgt(i_supVox) = sum(tabula_wgtmap(idxList));
-            end
-            
-            % -- make new indeces
-            optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
-        end
-        close(h_w1)
-    end
-end
-
-% -- make them robust --
-RO_params=0;
-optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
-% save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
-
-% -- get beamlet indeces --
-load([Pat_path, '\all_beams.mat'])
-Nbeamlets = all_beams{1}.Mxp;   % number of beamlets in a beam - usually 64 or 32
-weightTable = zeros(100,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}.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
-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([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
-
-
-plot_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 = make_robust_optGoal(optGoal, RO_params, beamlets);
-    % 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 ####====----
-%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
-%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
-%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
-    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 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
-            targetImg2=targetImg1;
-            % beamlets stay the same
-            for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
-                % modify target and beamlets
-                [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
-                % beamlets stay the same
-                for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rgs
-                    % modify target and beamlets
-                    targetImg4=targetImg3;
-                    % beamlets stay the same
-                    %% make new target and beamlets
-                    ROI_idx=[];
-                    ROI_idx=find(targetImg4>0);
-                    target = optGoal{goal_i}.D_final(idxValid);
-                    beamlets_pruned = beamlets(ROI_idx, :);
-                    % save to optGoal output
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx            = ROI_idx;
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned    = beamlets_pruned;
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target             = target;
-                end
-            end
-        end
-    end
-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
-
-
-
-
+
+
+% 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] = 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];
+else
+    Pat_path = varargin{1};
+    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
+    path2goal = varargin{2};
+    [Goal_path,Goal_file,ext] = fileparts(path2goal);
+end
+
+dialogue_box = 'yes'
+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([Pat_path, '\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)
+[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
+    switch SupVox_num
+        case 0
+        % 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;
+                % group superpixels
+                superMask = superpix_group(mask, SupVox_num, 'no'); 
+            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;
+                superMask = superpix_group(mask, SupVox_num, 'no'); 
+                orthoslice(superMask)
+        end
+        
+        
+        
+        superVoxList = unique(superMask);
+        superVoxList = superVoxList(superVoxList>0);
+        
+        optGoal{i_goal}.ROI_idx_old = optGoal{i_goal}.ROI_idx; % copy old index data
+        optGoal{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
+        optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
+        
+        if isfield(OptGoals.data{i_goal}, 'wgt_map')
+            tabula_wgtmap = zeros(size(superMask));
+            tabula_wgtmap(OptGoals.data{i_goal}.ROI_idx) = OptGoals.data{i_goal}.wgt_map;
+        end
+        
+        
+        
+        h_w1 = waitbar(0, 'merging superboxels');
+        for i_supVox = 1:numel(superVoxList)
+            waitbar(i_supVox/numel(superVoxList), h_w1)
+            supVox_idx = superVoxList(i_supVox);
+            idxList = find(superMask == supVox_idx);
+            optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
+            
+            if isfield(OptGoals.data{i_goal}, 'wgt_map')
+                optGoal{i_goal}.vox_wgt(i_supVox) = sum(tabula_wgtmap(idxList));
+            end
+            
+            % -- make new indeces
+            optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
+        end
+        close(h_w1)
+    end
+end
+
+% -- make them robust --
+RO_params=0;
+optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
+% save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
+
+% -- get beamlet indeces --
+load([Pat_path, '\all_beams.mat'])
+Nbeamlets = all_beams{1}.Mxp;   % number of beamlets in a beam - usually 64 or 32
+weightTable = zeros(100,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}.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
+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([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
+
+
+plot_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 = make_robust_optGoal(optGoal, RO_params, beamlets);
+    % 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 ####====----
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
+    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 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
+            targetImg2=targetImg1;
+            % beamlets stay the same
+            for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
+                % modify target and beamlets
+                [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
+                % beamlets stay the same
+                for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rgs
+                    % modify target and beamlets
+                    targetImg4=targetImg3;
+                    % beamlets stay the same
+                    %% make new target and beamlets
+                    ROI_idx=[];
+                    ROI_idx=find(targetImg4>0);
+                    target = optGoal{goal_i}.D_final(idxValid);
+                    beamlets_pruned = beamlets(ROI_idx, :);
+                    % save to optGoal output
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx            = ROI_idx;
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned    = beamlets_pruned;
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target             = target;
+                end
+            end
+        end
+    end
+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/dicomrt_loadctlist.m

@@ -1,196 +1,196 @@
-function [filelistAXIAL,xlocationAXIAL,ylocationAXIAL,zlocationAXIAL,change,user_option] = dicomrt_loadctlist(filename)
-% dicomrt_loadctlist(filename)
-%
-% Parse CT data set specified in filename. If more than one study are found within the CT data set
-% the user can select a study and dicomrt_loadctlist return a filelist and x,y,z coordinates of the
-% selected CT subset
-%
-% filename contains a list of CT slices to import
-%
-% change is 1 if any change to the filelist have been done, 0 otherwise
-% user_option is 1 if the user select not to continue with the program, 0 otherwise
-%
-% Example:
-% [list,x,y,z,change,user_option]=dicomrt_loadctlist(filename)
-%
-% with filename containint the following:
-% ct1 (group1)
-% ct2 (group2)
-% ct3 (group1)
-% ct4 (group2)
-%
-% if the user select one of them (e.g. group1):
-% list= contain only ct1 and ct3,
-% x= xlocation of ct1 and ct3 ,
-% y= ylocation of ct1 and ct3 ,
-% z= zlocation of ct1 and ct3 ,
-% change= 1,
-% user_option= 0
-%
-% See also dicomrt_loaddose dicomrt_loadct dicomrt_sortct
-%
-% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)
-
-% Retrieve the number of images
-nlines=dicomrt_nASCIIlines(filename);
-
-% Get CT images and create 3D Volume
-fid=fopen(filename);
-nct=0;
-counter=0;
-
-% Initialize variable
-filelist=' ';
-studyUIDlist=' ';
-listUID=' ';
-imagetype=' ';
-user_option=0;
-
-% Progress bar
-h = waitbar(0,['Loading progress:']);
-set(h,'Name','dicomrt_loadctlist: loading CT objects tags');
-
-%loop until the end-of-file is reached and build 3D CT matrix
-while (feof(fid)~=1);
-    nct=nct+1; % counting
-
-    nctcheck=nct; % check for eof
-
-    ct_file_location{1,nct}=fgetl(fid);
-
-    if isnumeric(ct_file_location{1,nct}), nct=nct-1, break, end %end of line reached
-
-    dictFlg = checkDictUse;
-    if dictFlg
-        info_temp=dicominfo(ct_file_location{1,nct}, 'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
-    else
-        info_temp=dicominfo(ct_file_location{1,nct});
-    end
-
-
-
-    xlocation(nct)=info_temp.ImagePositionPatient(1);
-    ylocation(nct)=info_temp.ImagePositionPatient(2);
-    zlocation(nct)=info_temp.ImagePositionPatient(3);
-
-    xlocation=xlocation';
-    ylocation=ylocation';
-    zlocation=zlocation';
-
-    filelist=char(filelist,info_temp.Filename);
-
-    if isfield(info_temp,'ImageType')~=1
-        warning('dicomrt_loadctlist: no DICOM ImageType was found. Assuming AXIAL CT Images');
-        imagetype='AXIAL';
-    else
-        imagetype=char(imagetype,info_temp.ImageType);
-    end
-
-    listUID=char(listUID,info_temp.StudyInstanceUID);
-
-    studyUID=info_temp.StudyInstanceUID;
-    if isequal(studyUID, studyUIDlist(size(studyUIDlist,1),:))==0
-        studyUIDlist=char(studyUIDlist,studyUID);
-    end
-    waitbar(nct/nlines,h);
-end
-filelist(1,:)=[];
-studyUIDlist(1,:)=[];
-listUID(1,:)=[];
-imagetype(1,:)=[];
-
-if size(studyUIDlist,1)>=2
-    change=1;
-    disp(' ');
-    warning([int2str(size(studyUIDlist,1)),' studies was found among the ct slices you want to import']);
-    disp(' ');
-    leave = input('Do you want to leave (Y/N) ? [N] ','s');
-    if leave == 'Y' | leave == 'y';
-        user_option=1;
-        return
-    else
-        disp('Available studies:');
-        for j=1:size(studyUIDlist,1)
-            disp([int2str(j), ' - ', studyUIDlist(j,:)]);
-        end
-        chooseUID = input(['Select a study to be imported from 1 to ',int2str(size(studyUIDlist,1)),': ']);
-        if isempty(chooseUID)==1 | isnumeric(chooseUID)~=1 | chooseUID>size(studyUIDlist,1)
-            error('dicomrt_loadctlist: There is no default to this answer or the number to entered is invalid. Exit now !');
-            user_option=1;
-        else
-            filelistUID=' ';
-            imagetypeUID=' ';
-            for k=1:size(filelist,1)
-                if listUID(k,:)==studyUIDlist(chooseUID,:);
-                    counter=counter+1;
-                    filelistUID=char(filelistUID,filelist(k,:));
-                    imagetypeUID=char(imagetypeUID,imagetype(k,:));
-                    xlocationUID(counter)=xlocation(k);
-                    ylocationUID(counter)=ylocation(k);
-                    zlocationUID(counter)=zlocation(k);
-                end
-            end
-            filelistUID(1,:)=[];
-            imagetypeUID(1,:)=[];
-            counter=0; % reset counter
-        end
-    end
-else
-    filelistUID=filelist;
-    imagetypeUID=imagetype;
-    xlocationUID=xlocation;
-    ylocationUID=ylocation;
-    zlocationUID=zlocation;
-end
-
-% Check for scout images (not AXIAL)
-imagetypeAXIAL=' ';
-filelistAXIAL=' ';
-
-for i=1:size(filelistUID,1)
-    if isempty(findstr('AXIAL',imagetypeUID(i,:)))==1 % Scout image found
-        disp(['The following image :',filelistUID(i,:),' is not AXIAL. Skipped ...']);
-        change=1; % just make sure we return alterations to the filelist
-        
-        % added to try read uulke prosate data
-        counter=counter+1;
-        imagetypeAXIAL=char(imagetypeAXIAL,imagetypeUID(i,:));
-        filelistAXIAL=char(filelistAXIAL,filelistUID(i,:));
-        xlocationAXIAL(counter)=xlocationUID(i);
-        ylocationAXIAL(counter)=ylocationUID(i);
-        zlocationAXIAL(counter)=zlocationUID(i);
-        % added to here
-    else
-        counter=counter+1;
-        imagetypeAXIAL=char(imagetypeAXIAL,imagetypeUID(i,:));
-        filelistAXIAL=char(filelistAXIAL,filelistUID(i,:));
-        xlocationAXIAL(counter)=xlocationUID(i);
-        ylocationAXIAL(counter)=ylocationUID(i);
-        zlocationAXIAL(counter)=zlocationUID(i);
-    end
-end
-
-if exist('change')~=1
-    change=0;
-end
-
-filelistAXIAL(1,:)=[];
-imagetypeAXIAL(1,:)=[];
-
-if change ==1 % export changes to file
-    newfilename=[filename,'.sort.txt'];
-    newfile=fopen(newfilename,'w');
-
-    for i=1:size(filelistAXIAL,1)
-        fprintf(newfile,'%c',deblank(filelistAXIAL(i,:))); fprintf(newfile,'\n');
-    end
-
-    disp(['A new file list has been written by dicomrt_loadctlist with name: ',newfilename]);
-    disp('This file will be used to import ct data instead');
-
-    fclose(newfile);
-end
-
-% Close progress bar
-close(h);
+function [filelistAXIAL,xlocationAXIAL,ylocationAXIAL,zlocationAXIAL,change,user_option] = dicomrt_loadctlist(filename)
+% dicomrt_loadctlist(filename)
+%
+% Parse CT data set specified in filename. If more than one study are found within the CT data set
+% the user can select a study and dicomrt_loadctlist return a filelist and x,y,z coordinates of the
+% selected CT subset
+%
+% filename contains a list of CT slices to import
+%
+% change is 1 if any change to the filelist have been done, 0 otherwise
+% user_option is 1 if the user select not to continue with the program, 0 otherwise
+%
+% Example:
+% [list,x,y,z,change,user_option]=dicomrt_loadctlist(filename)
+%
+% with filename containint the following:
+% ct1 (group1)
+% ct2 (group2)
+% ct3 (group1)
+% ct4 (group2)
+%
+% if the user select one of them (e.g. group1):
+% list= contain only ct1 and ct3,
+% x= xlocation of ct1 and ct3 ,
+% y= ylocation of ct1 and ct3 ,
+% z= zlocation of ct1 and ct3 ,
+% change= 1,
+% user_option= 0
+%
+% See also dicomrt_loaddose dicomrt_loadct dicomrt_sortct
+%
+% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)
+
+% Retrieve the number of images
+nlines=dicomrt_nASCIIlines(filename);
+
+% Get CT images and create 3D Volume
+fid=fopen(filename);
+nct=0;
+counter=0;
+
+% Initialize variable
+filelist=' ';
+studyUIDlist=' ';
+listUID=' ';
+imagetype=' ';
+user_option=0;
+
+% Progress bar
+h = waitbar(0,['Loading progress:']);
+set(h,'Name','dicomrt_loadctlist: loading CT objects tags');
+
+%loop until the end-of-file is reached and build 3D CT matrix
+while (feof(fid)~=1);
+    nct=nct+1; % counting
+
+    nctcheck=nct; % check for eof
+
+    ct_file_location{1,nct}=fgetl(fid);
+
+    if isnumeric(ct_file_location{1,nct}), nct=nct-1, break, end %end of line reached
+
+    dictFlg = checkDictUse;
+    if dictFlg
+        info_temp=dicominfo(ct_file_location{1,nct}, 'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
+    else
+        info_temp=dicominfo(ct_file_location{1,nct});
+    end
+
+
+
+    xlocation(nct)=info_temp.ImagePositionPatient(1);
+    ylocation(nct)=info_temp.ImagePositionPatient(2);
+    zlocation(nct)=info_temp.ImagePositionPatient(3);
+
+    xlocation=xlocation';
+    ylocation=ylocation';
+    zlocation=zlocation';
+
+    filelist=char(filelist,info_temp.Filename);
+
+    if isfield(info_temp,'ImageType')~=1
+        warning('dicomrt_loadctlist: no DICOM ImageType was found. Assuming AXIAL CT Images');
+        imagetype='AXIAL';
+    else
+        imagetype=char(imagetype,info_temp.ImageType);
+    end
+
+    listUID=char(listUID,info_temp.StudyInstanceUID);
+
+    studyUID=info_temp.StudyInstanceUID;
+    if isequal(studyUID, studyUIDlist(size(studyUIDlist,1),:))==0
+        studyUIDlist=char(studyUIDlist,studyUID);
+    end
+    waitbar(nct/nlines,h);
+end
+filelist(1,:)=[];
+studyUIDlist(1,:)=[];
+listUID(1,:)=[];
+imagetype(1,:)=[];
+
+if size(studyUIDlist,1)>=2
+    change=1;
+    disp(' ');
+    warning([int2str(size(studyUIDlist,1)),' studies was found among the ct slices you want to import']);
+    disp(' ');
+    leave = input('Do you want to leave (Y/N) ? [N] ','s');
+    if leave == 'Y' | leave == 'y';
+        user_option=1;
+        return
+    else
+        disp('Available studies:');
+        for j=1:size(studyUIDlist,1)
+            disp([int2str(j), ' - ', studyUIDlist(j,:)]);
+        end
+        chooseUID = input(['Select a study to be imported from 1 to ',int2str(size(studyUIDlist,1)),': ']);
+        if isempty(chooseUID)==1 | isnumeric(chooseUID)~=1 | chooseUID>size(studyUIDlist,1)
+            error('dicomrt_loadctlist: There is no default to this answer or the number to entered is invalid. Exit now !');
+            user_option=1;
+        else
+            filelistUID=' ';
+            imagetypeUID=' ';
+            for k=1:size(filelist,1)
+                if listUID(k,:)==studyUIDlist(chooseUID,:);
+                    counter=counter+1;
+                    filelistUID=char(filelistUID,filelist(k,:));
+                    imagetypeUID=char(imagetypeUID,imagetype(k,:));
+                    xlocationUID(counter)=xlocation(k);
+                    ylocationUID(counter)=ylocation(k);
+                    zlocationUID(counter)=zlocation(k);
+                end
+            end
+            filelistUID(1,:)=[];
+            imagetypeUID(1,:)=[];
+            counter=0; % reset counter
+        end
+    end
+else
+    filelistUID=filelist;
+    imagetypeUID=imagetype;
+    xlocationUID=xlocation;
+    ylocationUID=ylocation;
+    zlocationUID=zlocation;
+end
+
+% Check for scout images (not AXIAL)
+imagetypeAXIAL=' ';
+filelistAXIAL=' ';
+
+for i=1:size(filelistUID,1)
+    if isempty(findstr('AXIAL',imagetypeUID(i,:)))==1 % Scout image found
+        disp(['The following image :',filelistUID(i,:),' is not AXIAL. Skipped ...']);
+        change=1; % just make sure we return alterations to the filelist
+        
+        % added to try read uulke prosate data
+        counter=counter+1;
+        imagetypeAXIAL=char(imagetypeAXIAL,imagetypeUID(i,:));
+        filelistAXIAL=char(filelistAXIAL,filelistUID(i,:));
+        xlocationAXIAL(counter)=xlocationUID(i);
+        ylocationAXIAL(counter)=ylocationUID(i);
+        zlocationAXIAL(counter)=zlocationUID(i);
+        % added to here
+    else
+        counter=counter+1;
+        imagetypeAXIAL=char(imagetypeAXIAL,imagetypeUID(i,:));
+        filelistAXIAL=char(filelistAXIAL,filelistUID(i,:));
+        xlocationAXIAL(counter)=xlocationUID(i);
+        ylocationAXIAL(counter)=ylocationUID(i);
+        zlocationAXIAL(counter)=zlocationUID(i);
+    end
+end
+
+if exist('change')~=1
+    change=0;
+end
+
+filelistAXIAL(1,:)=[];
+imagetypeAXIAL(1,:)=[];
+
+if change ==1 % export changes to file
+    newfilename=[filename,'.sort.txt'];
+    newfile=fopen(newfilename,'w');
+
+    for i=1:size(filelistAXIAL,1)
+        fprintf(newfile,'%c',deblank(filelistAXIAL(i,:))); fprintf(newfile,'\n');
+    end
+
+    disp(['A new file list has been written by dicomrt_loadctlist with name: ',newfilename]);
+    disp('This file will be used to import ct data instead');
+
+    fclose(newfile);
+end
+
+% Close progress bar
+close(h);
 clear info_temp

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -77,6 +77,7 @@ classdef XTPS < handle
             % Tools menu
             obj.handles.hToolsMenu   = uimenu(obj.handles.hMainFigure, 'Label', 'Tools');            
             obj.handles.hPlotDVHMenu = uimenu(obj.handles.hToolsMenu,  'Label', 'Plot DVH');
+            obj.handles.hResampleGeometry = uimenu(obj.handles.hToolsMenu,  'Label', 'Resample Geometry');
             
 
             % Associate callbacks
@@ -108,6 +109,7 @@ classdef XTPS < handle
             set(obj.handles.hSaveAsOptResultsMenu, 'Callback',      @obj.SaveAsOptResultsMenu_Callback);
             set(obj.handles.hExportPinnacleMenu, 'Callback',        @obj.ExportPinnacleMenu_Callback);
             set(obj.handles.hPlotDVHMenu, 'Callback',               @obj.PlotDVHMenu_Callback);
+            set(obj.handles.hResampleGeometry, 'Callback',          @obj.ResampleGeometry_Callback);
             set(obj.handles.hDoseModeDropdown, 'Callback',          @obj.DoseModeDropdown_Callback);
             set(obj.handles.hWWSlider, 'Callback',                  @obj.WWSlider_Callback);
             set(obj.handles.hWWEdit, 'Callback',                    @obj.WWEdit_Callback);
@@ -671,6 +673,17 @@ classdef XTPS < handle
             ylabel('% volume');
             hold off;
         end
+        
+%-------------------------------------------------------------------------------
+        function ResampleGeometry_Callback(obj, src, evt)
+            %% --- make the figure prompt for number of angles and beamlets
+            
+            
+            % call the resample function
+            resample_geometry(obj.handles.hSVPS.Geometry)
+            disp('Resampled geometry created!')
+            
+        end
 %-------------------------------------------------------------------------------
         function DoseModeDropdown_Callback(obj, src, evt)
             str = get(obj.handles.hDoseModeDropdown, 'String');

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -1,35 +1,35 @@
-function B = merge_beamlets(num_batches, pat_dir)
-
-% num_batches = 4;
-
-
-B = [];
-for k = 1:num_batches
-    batch_doses = read_ryan_beamlets([pat_dir, '\batch_dose' num2str(k-1) '.bin']);
-    B = [B batch_doses];
-end
-for k = 1:numel(B)
-    B{k}.num = k-1;
-end
-write_ryan_beamlets([pat_dir '\batch_dose_unshifted.bin'],B);
-
-% apply the shift (0.5 vox up and left) to make dose match geometry
-if true 
-    f = waitbar(0, 'Correcting beamlet shift');
-    for beam_i = 1:numel(B)
-        new_dose_list = B{beam_i};
-        waitbar(beam_i/numel(B),f, ['Correcting beamlet shift: ' num2str(beam_i)])
-
-        tabula = zeros(B{1, beam_i}.y_count, B{1, beam_i}.x_count, B{1, beam_i}.z_count);
-        tabula(B{1, beam_i}.non_zero_indices) = B{1, beam_i}.non_zero_values;
-
-        shifted_img = imtranslate(tabula, [-0.5,-0.5,-0.5]);
-        nonzero_idx = find(shifted_img>0.005*max(shifted_img(:)));
-
-        B{beam_i}.non_zero_indices = nonzero_idx;
-        B{beam_i}.non_zero_values = shifted_img(nonzero_idx);
-    end
-    close(f)
-end
-write_ryan_beamlets([pat_dir '\batch_dose.bin'],B);
-write_ryan_beamlets([pat_dir '\batch_dose_backup.bin'],B);
+function B = merge_beamlets(num_batches, pat_dir)
+
+% num_batches = 4;
+
+
+B = [];
+for k = 1:num_batches
+    batch_doses = read_ryan_beamlets([pat_dir, '\batch_dose' num2str(k-1) '.bin']);
+    B = [B batch_doses];
+end
+for k = 1:numel(B)
+    B{k}.num = k-1;
+end
+write_ryan_beamlets([pat_dir '\batch_dose_unshifted.bin'],B);
+
+% apply the shift (0.5 vox up and left) to make dose match geometry
+if true 
+    f = waitbar(0, 'Correcting beamlet shift');
+    for beam_i = 1:numel(B)
+        new_dose_list = B{beam_i};
+        waitbar(beam_i/numel(B),f, ['Correcting beamlet shift: ' num2str(beam_i)])
+
+        tabula = zeros(B{1, beam_i}.y_count, B{1, beam_i}.x_count, B{1, beam_i}.z_count);
+        tabula(B{1, beam_i}.non_zero_indices) = B{1, beam_i}.non_zero_values;
+
+        shifted_img = imtranslate(tabula, [-0.5,-0.5,-0.5]);
+        nonzero_idx = find(shifted_img>0.005*max(shifted_img(:)));
+
+        B{beam_i}.non_zero_indices = nonzero_idx;
+        B{beam_i}.non_zero_values = shifted_img(nonzero_idx);
+    end
+    close(f)
+end
+write_ryan_beamlets([pat_dir '\batch_dose.bin'],B);
+write_ryan_beamlets([pat_dir '\batch_dose_backup.bin'],B);

WiscPlanPhotonkV125/matlab_frontend/resample_geometry.m

@@ -0,0 +1,99 @@
+
+function resample_geometry(Geometry, Nx, Ny, Nz)
+    Geometry_old = Geometry;
+    str = inputdlg({'Enter X voxel dimension (cm):', 'Enter Y voxel dimension (cm):', ...
+        'Enter Z voxel dimension (cm):'}, 'input', [1,8], {'0.2', '0.2', '0.2'});
+    % dimensions in mm
+    X_dim = str2double(str{1});
+    Y_dim = str2double(str{2});
+    Z_dim = str2double(str{3});
+    disp(['Desired voxel dimensions: ' num2str(X_dim) ' ' num2str(Y_dim) ' ' num2str(Z_dim) ' mm'])
+    
+    % round to appropriate number of voxels
+    Nx = round(size(Geometry_old.data,1)*Geometry_old.voxel_size(1)/X_dim);
+    Ny = round(size(Geometry_old.data,2)*Geometry_old.voxel_size(2)/Y_dim);
+    Nz = round(size(Geometry_old.data,3)*Geometry_old.voxel_size(3)/Z_dim);
+    
+    % find the correct scale factor
+    fx = Nx/size(Geometry_old.data,1);
+    fy = Ny/size(Geometry_old.data,2);
+    fz = Nz/size(Geometry_old.data,3);
+    
+    %% start converting the Geometry file
+    Geometry.voxel_size = Geometry_old.voxel_size./[fx, fy, fz];
+    disp(['New voxel sizes: ' num2str(Geometry.voxel_size(1)) ' ' num2str(Geometry.voxel_size(2)) ...
+        ' ' num2str(Geometry.voxel_size(3)) ' mm'])
+    
+    load('WiscPlan_preferences.mat')
+    if ~isfield(WiscPlan_preferences,'inDataPath')
+        Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
+    elseif WiscPlan_preferences.inDataPath ~= 0
+        Geometry.patient_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the new patient directory');
+    else
+        Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
+    end
+    
+    Geometry.data = imresize3(Geometry_old.data, [Nx, Ny, Nz]);
+    Geometry.rhomw = imresize3(Geometry_old.rhomw, [Nx, Ny, Nz]);
+    Geometry.Smw = imresize3(Geometry_old.Smw, [Nx, Ny, Nz]);
+    Geometry.Fmw2 = imresize3(Geometry_old.Fmw2, [Nx, Ny, Nz]);
+    Geometry.BTV = logical(round(imresize3(double(Geometry_old.BTV), [Nx, Ny, Nz])));
+    Geometry.Ring = logical(round(imresize3(double(Geometry_old.Ring), [Nx, Ny, Nz])));
+
+    Geometry.x = Geometry.start(1) +linspace(0, Nx-1, Nx)*Geometry.voxel_size(1);
+    Geometry.y = Geometry.start(2) +linspace(0, Ny-1, Ny)*Geometry.voxel_size(2);
+    Geometry.z = Geometry.start(3) +linspace(0, Nz-1, Nz)*Geometry.voxel_size(3);
+    
+    Geometry.ROIS = {};
+    for ROI_i = 1:numel(Geometry_old.ROIS)
+       Geometry.ROIS{ROI_i}.name = Geometry_old.ROIS{ROI_i}.name;
+       tabula = zeros(size(Geometry_old.data));
+       tabula(Geometry_old.ROIS{ROI_i}.ind) = 1;
+       tabula = imresize3(tabula, [Nx, Ny, Nz]);
+       tabula = imgaussfilt3(tabula, 0.5);
+       Geometry.ROIS{ROI_i}.ind = find(tabula>0.3);
+    end
+    
+    %% create export folder and start populating with files
+    mkdir([Geometry.patient_dir '\matlab_files']);
+    save([Geometry.patient_dir '\matlab_files\Geometry.mat'], 'Geometry');
+    
+    
+    if isfile([Geometry_old.patient_dir '\batch_dose_backup.bin'])
+        disp('calculated beamlet file found - resampling')
+        
+        % resample the beamlets file
+        f = waitbar(0, 'Resampling beamlets');
+        
+        beamlets_in = read_ryan_beamlets([Geometry_old.patient_dir '\batch_dose_backup.bin'], 'ryan');
+        for beam_i = 1:numel(beamlets_in)
+            
+            waitbar(beam_i/numel(beamlets_in),f, ['Correcting beamlet shift: ' num2str(beam_i)])
+            tabula = zeros(beamlets_in{beam_i}.x_count, beamlets_in{beam_i}.y_count, beamlets_in{beam_i}.z_count);
+            
+            beamlets_out{beam_i}.num = beamlets_in{beam_i}.num-1; %reading beamlets adds +1 to number, this fixes it
+            beamlets_out{beam_i}.x_count = Nx;
+            beamlets_out{beam_i}.y_count = Ny;
+            beamlets_out{beam_i}.z_count = Nz;
+            
+            tabula(beamlets_in{beam_i}.non_zero_indices) = beamlets_in{beam_i}.non_zero_values;
+            tabula = imresize3(tabula, [Nx, Ny, Nz]);
+            
+            beamlets_out{beam_i}.non_zero_indices = find(tabula>0.005*max(tabula(:)));
+            beamlets_out{beam_i}.non_zero_values = tabula( beamlets_out{beam_i}.non_zero_indices);
+
+        end
+        
+        write_ryan_beamlets([Geometry.patient_dir '\batch_dose.bin'],beamlets_out);
+        write_ryan_beamlets([Geometry.patient_dir '\batch_dose_backup.bin'],beamlets_out);
+        
+        close(f)
+        
+    else
+        disp('no beamlet file found')
+    end
+    disp('done!')
+
+end
+
+

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -1,82 +1,82 @@
-
-
-function superMask = superpix_group(mask, N_svox_in, orthoPlot)
-% this function contains supervoxel grouping
-fprintf(['\n' 'Creating ' num2str(N_svox_in) ' SupVox:\n' ])
-canvas2 = mask;
-% canvas2 = sqrt(bwdist(1-mask));
-% orthoslice (canvas2)
-
-% find box crop of ROI
-ymin = 0;
-ymax = 0;
-for yi = 1:size(mask,1)
-    data = mask(yi, :,:);
-    if max(data(:)) >0
-        if ymin == 0
-            ymin = yi;
-        end
-        ymax = yi;
-    end
-end
-xmin = 0;
-xmax = 0;
-for xi = 1:size(mask,2)
-    data = mask(:,xi,:);
-    if max(data(:)) >0
-        if xmin == 0
-            xmin = xi;
-        end
-        xmax = xi;
-    end
-end
-zmin = 0;
-zmax = 0;
-for zi = 1:size(mask,3)
-    data = mask(:,:,zi);
-    if max(data(:)) >0
-        if zmin == 0
-            zmin = zi;
-        end
-        zmax = zi;
-    end
-end
-
-canvas3 = canvas2(ymin:ymax, xmin:xmax, zmin:zmax);
-
-N_svox = N_svox_in; % number of supervoxels to give as param to parse
-converge_factor = 1;
-
-while true
-    [L,NumLabels] = superpixels3(canvas3,N_svox);
-
-    superMask = zeros(size(mask));
-    superMask(ymin:ymax, xmin:xmax, zmin:zmax) = L;
-    superMask(mask==0) = 0;
-
-    numSupVox = numel(unique(superMask))-1; % number of created supervoxels
-    fprintf([num2str(numSupVox) ' areas created' ])
-    
-    if abs(numSupVox-N_svox_in)/N_svox_in < 0.2/sqrt(converge_factor)
-        
-        switch orthoPlot
-            case 'yes'
-                orthoslice(superMask)
-            case 'no'
-                fprintf('\n supervoxel volume created! (orthoslice skipped)')
-        end
-        break
-    end
-    N_svox2 = N_svox * N_svox_in/numSupVox;
-    N_svox = round(converge_factor * N_svox2 + (1-converge_factor)*N_svox);
-    converge_factor = converge_factor * 0.95;
-    fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target. New start size: ' num2str(N_svox) '\n'])
-    
-    if converge_factor < 0.01
-%         orthoslice(superMask)
-        break
-    end
-end 
-fprintf([' -  ok.\n'])
-
+
+
+function superMask = superpix_group(mask, N_svox_in, orthoPlot)
+% this function contains supervoxel grouping
+fprintf(['\n' 'Creating ' num2str(N_svox_in) ' SupVox:\n' ])
+canvas2 = mask;
+% canvas2 = sqrt(bwdist(1-mask));
+% orthoslice (canvas2)
+
+% find box crop of ROI
+ymin = 0;
+ymax = 0;
+for yi = 1:size(mask,1)
+    data = mask(yi, :,:);
+    if max(data(:)) >0
+        if ymin == 0
+            ymin = yi;
+        end
+        ymax = yi;
+    end
+end
+xmin = 0;
+xmax = 0;
+for xi = 1:size(mask,2)
+    data = mask(:,xi,:);
+    if max(data(:)) >0
+        if xmin == 0
+            xmin = xi;
+        end
+        xmax = xi;
+    end
+end
+zmin = 0;
+zmax = 0;
+for zi = 1:size(mask,3)
+    data = mask(:,:,zi);
+    if max(data(:)) >0
+        if zmin == 0
+            zmin = zi;
+        end
+        zmax = zi;
+    end
+end
+
+canvas3 = canvas2(ymin:ymax, xmin:xmax, zmin:zmax);
+
+N_svox = N_svox_in; % number of supervoxels to give as param to parse
+converge_factor = 1;
+
+while true
+    [L,NumLabels] = superpixels3(canvas3,N_svox);
+
+    superMask = zeros(size(mask));
+    superMask(ymin:ymax, xmin:xmax, zmin:zmax) = L;
+    superMask(mask==0) = 0;
+
+    numSupVox = numel(unique(superMask))-1; % number of created supervoxels
+    fprintf([num2str(numSupVox) ' areas created' ])
+    
+    if abs(numSupVox-N_svox_in)/N_svox_in < 0.2/sqrt(converge_factor)
+        
+        switch orthoPlot
+            case 'yes'
+                orthoslice(superMask)
+            case 'no'
+                fprintf('\n supervoxel volume created! (orthoslice skipped)')
+        end
+        break
+    end
+    N_svox2 = N_svox * N_svox_in/numSupVox;
+    N_svox = round(converge_factor * N_svox2 + (1-converge_factor)*N_svox);
+    converge_factor = converge_factor * 0.95;
+    fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target. New start size: ' num2str(N_svox) '\n'])
+    
+    if converge_factor < 0.01
+%         orthoslice(superMask)
+        break
+    end
+end 
+fprintf([' -  ok.\n'])
+
 end

WiscPlanPhotonkV125/matlab_frontend/update_geometry_path.m

@@ -1,10 +1,10 @@
-
-
-path = 'C:\021_WiscPlan_data\patient013\WiscPlan';
-
-load([path '\matlab_files\Geometry.mat'])
-Geometry.patient_dir = path;
-Geometry.data_dir = path;
-save([path '\matlab_files\Geometry.mat'], 'Geometry')
-
-disp('Done!')
+
+
+path = 'C:\021_WiscPlan_data\patient013\WiscPlan';
+
+load([path '\matlab_files\Geometry.mat'])
+Geometry.patient_dir = path;
+Geometry.data_dir = path;
+save([path '\matlab_files\Geometry.mat'], 'Geometry')
+
+disp('Done!')