Save before goal maker updates.

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -1,340 +0,0 @@
-
-
-% 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_beamlet_optimizer(paths)
-% This function performs the beamlet optimization
-% Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
-% Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
-% 
-% [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
-% 
-% Made by Peter Ferjancic 1. May 2018
-% Last updated: 14. August 2018
-% Inspired by Ana Barrigan's REGGUI optimization procedures
-%
-% paths.patient_dir
-% paths.Goal_dir (previously called DP_dir)
-% paths.patient
-% paths.goalsName
-
-N_fcallback1 = 5000;
-N_fcallback2 = 200000;
-% paths.patient = 'medivation_01';
-% % paths.goalsName = 'ROI_goals_simple';
-
-if false
-     'medivation_01'
-        paths.patient_dir = 'F:\021_WiscPlan_data\Medivation_0391811_crop1';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2';
-        paths.Goal_dir = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2';
-        paths.goalsName = 'ROI_goals_simple';
-     'temp'
-        paths.patient_dir = 'C:\000-temp\DELETEME';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis'
-     'avastin_009'
-        paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
-     'gbm_005'
-%         paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
-        paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
-     'gbm_015'
-        paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
-     'gbm_022'
-        paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
-        paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
-     'avastin_009'
-    
-        error('invalid case')
-end
-
-% merge_beamlets(4, patient_dir);
-
-
-%% PROGRAM STARTS HERE
-% - no tocar lo que hay debajo -
-fprintf('starting NLP optimization process... \n')
-
-% % -- LOAD GEOMETRY AND BEAMLETS --
-load([paths.patient_dir '\matlab_files\Geometry.mat'])
-
-
-%% -- OPTIMIZATION TARGETS --
-load([paths.Goal_dir '\RODP_files\' paths.goalsName '.mat']);
-
-
-optGoal = ROI_goals.optGoal;
-optGoal_beam = ROI_goals.optGoal_beam;
-optGoal_idx = ROI_goals.optGoal_idx;
-targetMinMax_idx = ROI_goals.targetMinMax_idx;
-
-[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, paths.patient_dir);
-
-
-% -- make them robust --
-RO_params=0;
-optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
-optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
-
-%% -- INITIALIZE BEAMLET WEIGHTS --
-w0_beams = ones(numBeam,1);
-w0_beams = mean(optGoal_beam{1}.target ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
-w0_beams = w0_beams + (2*rand(size(w0_beams))-1) *0.1 .*w0_beams; % random perturbation
-w0_beams = double(w0_beams);
-
-% -- CALLBACK OPTIMIZATION FUNCTION --
-fun1 = @(x) get_penalty(x, optGoal_beam);
-fun2 = @(x) get_penalty(x, optGoal);
-
-% -- 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.OptimalityTolerance = 1e-9;
-
-fprintf('\n running initial optimizer:')
-
-%% Run the optimization
-% -- GET FULL BEAM WEIGHTS --
-tic
-w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,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
-w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % small random perturbation
-w_beamlets = 1.1* w_beamlets; % this just kicks the beamlets up a bit to make it easier for the optimizer to start
-
-% -- 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;
-save([paths.patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
-
-plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
-colorwash(Geometry.data, D_full, [-500, 500], [0, 110]);
-% plot_DVH_robust(D_full, optGoal, optGoal_idx)
-end
-
-%% support functions
-% ---- PENALTY FUNCTION ----
-function penalty = get_penalty(x, optGoal)
-    % 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,1);  
-    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);
-                sc_i = sc_i + 1;
-            end
-        end
-    end
-    % take the worst case penalty of evaluated scenarios
-    penalty=max(fobj);
-end
-% ------ supp: penalty for single scenario ------
-function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
-    penalty = 0;
-    % 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 higher 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.^2);
-            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.^2);
-            case 'LeastSquare'
-                % penalize with sum of squares any deviation from target
-                % dose
-                d_penalty = 1.0e-1* sum(((...
-                    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).^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.0e5 * 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)
-                    
-        end
-        penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
-    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_mag = 2; % vox of shift
-    nrs_scene_list={[0,0,0]};
-
-    
-    % ----====#### CHANGE ROBUSTNESS HERE ####====----
-    
-%     sss_scene_list={[0,0,0]};
-    sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
-%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
-%         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
-
-    % ----====#### CHANGE ROBUSTNESS HERE ####====----
-
-    rgs_scene_list={[0,0,0]};
-
-%     [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');
-    
-    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);
-                    
-                    if isfield(optGoal{goal_i}, 'target_alpha')
-                        target = optGoal{goal_i}.target;
-                        target = target(idxValid);
-%                         disp('exists')
-                    else
-                        target = ones(numel(ROI_idx), 1) .* optGoal{goal_i}.D_final;
-%                         disp('not exists')
-                    end
-                    
-                    
-                    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/NLP_getFullDose.m

@@ -0,0 +1,18 @@
+
+
+function NLP_getFullDose(NLP_path, NLP_file)
+% this function calculates full dose matrix for Robust optimized plans
+
+% get the NLP result
+load([NLP_path '\' NLP_file]);
+
+path2 = NLP_path(1:end-14);
+D = read_ryan_beamlets([path2 '\beamlet_batch_files\beamletbatch0.bin'],'ryan sum', NLP_result.weights);
+
+% orthoslice(D-NLP_result.dose, [0,90])
+
+NLP_result.dose = double(D);
+% NLP_result.weights = double(optResults.weights{end});
+save([NLP_path, NLP_file], 'NLP_result');
+
+end

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v2.m

@@ -1,403 +0,0 @@
-
-
-% 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_v2(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
-
-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], {'10000', '500000', 'n'});
-N_fcallback1 = str2double(str{1}); % 10000  is a good guesstimate
-N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
-pre_beamWeights = str{3};
-
-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, 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
-    
-    
-    switch SupVox_num
-        case 0
-        % if not supervoxel, just select provided ROI_idx
-        optGoal{i_goal} = OptGoals.data{i_goal};
-        optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
-        optGoal_beam{i_goal} = OptGoals.data{i_goal};
-        optGoal_beam{i_goal}.beamlets_pruned = sparse(beamlets_joined(optGoal{i_goal}.ROI_idx, :));
-    
-        otherwise
-        % -- if supervoxel, merge given columns
-        % - make supervoxel map
-        
-        mask = zeros(OptGoals.data{i_goal}.imgDim);
-        mask(OptGoals.data{i_goal}.ROI_idx) = 1;
-        superMask = superpix_group(mask, SupVox_num); 
-        superVoxList = unique(superMask);
-        superVoxList = superVoxList(superVoxList>0);
-        
-        optGoal{i_goal} = OptGoals.data{i_goal};
-        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);
-        
-        optGoal_beam{i_goal} = OptGoals.data{i_goal};
-        optGoal_beam{i_goal}.ROI_idx_old = optGoal_beam{i_goal}.ROI_idx; 
-        optGoal_beam{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
-        optGoal_beam{i_goal}.opt_weight = optGoal_beam{i_goal}.opt_weight * numel(optGoal_beam{i_goal}.ROI_idx_old)/numel(optGoal_beam{i_goal}.ROI_idx);
-        
-        h_w1 = waitbar(0, 'merging superboxels');
-        for i_supVox = 1:numel(superVoxList)
-            supVox_idx = superVoxList(i_supVox);
-            
-            waitbar(i_supVox/numel(superVoxList), h_w1)
-            
-            idxList = find(superMask == supVox_idx);
-            
-            optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
-            optGoal_beam{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets_joined(idxList, :),1));
-            % -- make new indeces
-            optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
-            optGoal_beam{i_goal}.ROI_idx(i_supVox) = idxList(1);
-        end
-        close(h_w1)
-    end
-    
-    
-end
-
-
-
-% -- make them robust --
-RO_params=0;
-optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
-optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
-
-% -- CALLBACK OPTIMIZATION FUNCTION --
-fun1 = @(x) get_penalty(x, optGoal_beam);
-fun2 = @(x) get_penalty(x, optGoal);
-
-% -- 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_beams = ones(numBeam,1);
-        w0_beams = mean(optGoal_beam{1}.D_final(optGoal{1}.ROI_idx) ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
-        w0_beams = double(w0_beams);
-
-        % -- GET BEAM WEIGHTS --
-        tic
-        w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,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;
-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, 36]);
-
-
-end
-
-%% support functions
-% ---- PENALTY FUNCTION ----
-function penalty = get_penalty(x, optGoal)
-    % 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,1);  
-    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);
-                sc_i = sc_i + 1;
-            end
-        end
-    end
-    % take the worst case penalty of evaluated scenarios
-    penalty=max(fobj);
-end
-% ------ supp: penalty for single scenario ------
-function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
-    penalty = 0;
-    % 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_exp'
-                % 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(exp(temp1));
-            case 'max_exp'
-                % 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(exp(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)
-                    
-        end
-        penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
-    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_mag = 1; % vox of shift
-    nrs_scene_list={[0,0,0]};
-
-    
-    % ----====#### CHANGE ROBUSTNESS HERE ####====----
-    
-    sss_scene_list={[0,0,0]};
-%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
-%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
-%         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
-
-    % ----====#### CHANGE ROBUSTNESS HERE ####====----
-
-    rgs_scene_list={[0,0,0]};
-
-%     [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');
-    
-    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/NLP_optimizer_v3.m

@@ -8,7 +8,7 @@
 % colorwash(Geometry.data, D_full, [500, 1500], [0,70])
 % orthoslice(D_full, [0,70])
 
-function [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v2(varargin)
+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;
 % 
@@ -66,7 +66,7 @@ fprintf('starting NLP optimization process... \n')
 % % -- LOAD GEOMETRY, GOALS, BEAMLETS --
 load(path2geometry)
 load(path2goal)
-[beamlets, numBeamlet] = get_beamlets(Geometry, Pat_path);
+[beamlets, numBeamlet] = get_beamlets(Geometry, Pat_path, OptGoals);
 % [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, Pat_path);
 
 
@@ -99,12 +99,24 @@ for i_goal = 1:size(OptGoals.goals,1)
         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
@@ -118,6 +130,8 @@ end
 RO_params=0;
 optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 
+save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
+
 % -- CALLBACK OPTIMIZATION FUNCTION --
 fun1 = @(x) get_penalty(x, optGoal_beam);
 fun2 = @(x) get_penalty(x, optGoal);
@@ -152,7 +166,8 @@ switch pre_beamWeights
 %         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;
+%         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 --
@@ -192,7 +207,7 @@ 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, 36]);
+colorwash(Geometry.data, D_full, [500, 1500], [0, 90]);
 
 
 end
@@ -217,6 +232,12 @@ function penalty = get_penalty(x, optGoal)
     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 penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
@@ -274,10 +295,29 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
                     (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;
     end
+    
+    %% add modulation penalty
+    if true
+        max_modulation = 5;
+        mod_pen_weight = 1.0e10;
+        % calc the penalty
+        mod_excess = max(0, x-5*mean(x));
+        mod_pen = mod_pen_weight*(sum(mod_excess) + numel(mod_excess)* any(mod_excess));
+        penalty = penalty + mod_pen;
+    end
 end
 
 % ---- MAKE ROI ROBUST ----
@@ -293,62 +333,52 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % Y - Y>0 moves image down
     % Z - in/out.
     
-    shift_mag = 1; % vox of shift
+    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 ####====----
     
 %     sss_scene_list={[0,0,0]};
-    sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
+    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]};
 %     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
 %         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
 
     % ----====#### CHANGE ROBUSTNESS HERE ####====----
-
-    rgs_scene_list={[0,0,0]};
-
 %     [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;
@@ -358,7 +388,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
         end
     end
 end
-% ------ supp: RO case SSS ------
+%% ------ supp: RO case SSS ------
 function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
     % translate the target image 
     

WiscPlanPhotonkV125/matlab_frontend/dicomrt/dicomrt2geometry.m

@@ -123,8 +123,8 @@ if downsample_flag == true
         ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);
         ct_ymesh = linspace(ct_ymesh(1)+dy/2, ct_ymesh(end)-dy/2, numel(ct_ymesh)*downsample_factor);
     else
-        downsample_factor = 0.125;
-        downsample_factor = 1/8;
+        downsample_factor = 0.25;
+        downsample_factor = 1/2;
         dx = ct_xmesh(2)-ct_xmesh(1);
         dy = ct_ymesh(2)-ct_ymesh(1);
         ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);

WiscPlanPhotonkV125/matlab_frontend/get_beamlets.m

@@ -1,8 +1,16 @@
 
 
-function [beamlets, numBeamlet] = get_beamlets(Geometry, patient_dir)
+function [beamlets, numBeamlet] = get_beamlets(Geometry, patient_dir, OptGoals)
 % 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);
+
+%% add idx lists for each goal!
+
 beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
 beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
 
@@ -15,7 +23,7 @@ else
     batchSize = 200;
 end
 
-beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
+beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
 % --- join beamlets into beams
 % load([patient_dir '\all_beams.mat'])
 % beamletOrigin=[0 0 0];
@@ -44,14 +52,20 @@ beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
 end
 
 % ---- GET BEAMLETS ----
-function beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
+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 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);
@@ -59,7 +73,7 @@ function beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
         end
 
 %         beamlet_batch(idx, beam_i_temp) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
-        beamlet_batch(idx, beam_i_temp) = 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

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -80,6 +80,8 @@ handles.uitable1.ColumnFormat{5}{6} = 'max_exp';
 handles.uitable1.ColumnFormat{5}{7} = 'LeastSquare';
 handles.uitable1.ColumnFormat{5}{8} = 'min_perc_Volume';
 handles.uitable1.ColumnFormat{5}{9} = 'max_perc_Volume';
+handles.uitable1.ColumnFormat{5}{8} = 'min_sq_voxwgt';
+handles.uitable1.ColumnFormat{5}{9} = 'max_sq_voxwgt';
 
 % == populate the first entry ==
 handles.uitable1.Data = handles.uitable1.Data(1,:);
@@ -192,13 +194,35 @@ for i = 1 : size(handles.uitable1.Data, 1)
     elseif strcmp (handles.uitable1.Data{i,3}, 'Dose map')
 %         warning('Works only for NRRD (for now)')
         dose_in = nrrdread(handles.uitable1.Data{i,4});
-        colorwash(Geometry.data-1000, dose_in, [-500, 500], [0, 1.5*(max(dose_in(:)))], 1, goal_i.name)
+        colorwash(Geometry.data-1000, dose_in, [-500, 500], [0, 1.0*(max(dose_in(:)))], 1, goal_i.name)
         goal_i.D_final = dose_in(goal_i.ROI_idx);
         % loading of data comes here
         % matrix size checks
     end
     % == goal penalty function
     goal_i.function = handles.uitable1.Data{i,5};
+    
+    % == goal penalty function
+    if strcmp (handles.uitable1.Data{i,5}, 'min_sq_voxwgt')
+        handles.uitable1.Data{i,2};
+        voxwgt_in = nrrdread(handles.uitable1.Data{i, 7});
+%         load(handles.uitable1.Data{i, 7});
+%         voxwgt_in = miTLM_3;
+        goal_i.wgt_map = voxwgt_in(goal_i.ROI_idx);
+        goal_i.wgt_map = goal_i.wgt_map / mean(goal_i.wgt_map(:)); % normalize
+        
+    elseif strcmp (handles.uitable1.Data{i,5}, 'max_sq_voxwgt')
+        voxwgt_in = nrrdread(handles.uitable1.Data{i, 7});
+%         load(handles.uitable1.Data{i, 7});
+%         voxwgt_in = miTLM_3;
+        goal_i.wgt_map = voxwgt_in(goal_i.ROI_idx);
+        goal_i.wgt_map = goal_i.wgt_map / mean(goal_i.wgt_map(:)); % normalize
+    else
+        if isfield (goal_i, 'wgt_map')
+            goal_i = rmfield(goal_i, 'wgt_map');
+        end
+    end
+    
     % == goal weight
     goal_i.opt_weight = handles.uitable1.Data{i,6} / numel(goal_i.ROI_idx); % normalize to volume of target area
     % == optional
@@ -262,6 +286,22 @@ if eventdata.Indices(2) == 3;
     end
 end
 
+% --- update DP or fixed dose goals
+if eventdata.Indices(2) == 5;
+    mode = handles.uitable1.Data{eventdata.Indices(1), eventdata.Indices(2)}
+    switch mode
+        case 'min_sq_voxwgt'
+            [FileName,PathName,FilterIndex] = uigetfile('F:\021_WiscPlan_data\FET_repeat_005_1\matlab_files\*.mat', 'Select voxel weight map');
+            handles.uitable1.Data{eventdata.Indices(1), 7} = [PathName, FileName];
+            
+        case 'max_sq_voxwgt'
+            [FileName,PathName,FilterIndex] = uigetfile('F:\021_WiscPlan_data\FET_repeat_005_1\matlab_files\*.mat', 'Select voxel weight map');
+            handles.uitable1.Data{eventdata.Indices(1), 7} = [PathName, FileName];
+            
+        otherwise
+    end
+end
+
 % --- update with number of voxels
 if eventdata.Indices(2) == 2;
     for i = 1:size(handles.Data.Geometry.ROIS,2)

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -52,7 +52,11 @@ ypmax = 1.25;
 % y-prime points in the z-direction in the CT coordinate system
 
 % ============================================= End of user-supplied inputs
-executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';
+% executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';
+
+executable_path = mfilename('fullpath');
+executable_path = [executable_path(1:end-37), 'WiscPlanEXE\RyanCsphoton.x86.exe']
+
 kernel_file = 'Kernels.mat';
 geometry_file = fullfile(patient_dir, 'matlab_files\Geometry.mat');
 
@@ -278,10 +282,12 @@ end
 all_beams = horzcat(batches{:});
 num_beams_per_batch = ceil(numel(all_beams)/num_batches);
 batches = cell(num_batches,1);
-for k = 1:(num_batches-1)
-    batches{k} = all_beams(1+num_beams_per_batch*(k-1):num_beams_per_batch*k);
+for k = 1:(num_batches)
+    beams_idx = 1+num_beams_per_batch*(k-1):num_beams_per_batch*k;
+    beams_idx (beams_idx>numel(all_beams)) = [];
+    batches{k} = all_beams(beams_idx);
 end
-batches{num_batches} = all_beams(1+num_beams_per_batch*(k):end);
+% batches{num_batches} = all_beams(1+num_beams_per_batch*(k):end);
 
 
 % Everything else in this file is related to saving the batches in a

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -61,8 +61,9 @@ classdef XTPS < handle
             obj.handles.hSaveOptimResultsMenu = uimenu(obj.handles.hFileMenu, 'Label', '3a.2 Save Optimization Results');
             obj.handles.hGetFullDose          = uimenu(obj.handles.hFileMenu, 'Label', '3a.3 Get Full Dose');
             
-            obj.handles.hRO_GoalSpec          = uimenu(obj.handles.hFileMenu, 'Label', '3b. RO Goal Spec');
+            obj.handles.hRO_GoalSpec          = uimenu(obj.handles.hFileMenu, 'Label', '3a. RO Goal Spec');
             obj.handles.hRO_Optimize          = uimenu(obj.handles.hFileMenu, 'Label', '3b. RO Optimization');
+            obj.handles.hRO_fullDose          = uimenu(obj.handles.hFileMenu, 'Label', '3c. Full dose calc');
             
             
             obj.handles.hLoadGeometryMenu     = uimenu(obj.handles.hFileMenu, 'Label', 'Load Geometry', 'Separator','on');
@@ -93,6 +94,7 @@ classdef XTPS < handle
             
             set(obj.handles.hRO_GoalSpec, 'Callback',   @obj.RO_GoalSpec);
             set(obj.handles.hRO_Optimize, 'Callback',    @obj.RO_Optimize);
+            set(obj.handles.hRO_fullDose, 'Callback',    @obj.RO_fullDose);
             
             % Grozomah
             set(obj.handles.hMergeBeamlets, 'Callback',      @obj.MergeBeamlets_Callback);
@@ -481,7 +483,14 @@ classdef XTPS < handle
         %-------------------------------------------------------------------------------
 % Grozomah
         function MergeBeamlets_Callback(obj, src, evt)
-            merge_beamlets(obj.handles.hSVPS.Geometry.num_batches, obj.handles.hSVPS.Geometry.patient_dir);
+            if isfield(obj.handles.hSVPS.Geometry, 'num_batches')
+                merge_beamlets(obj.handles.hSVPS.Geometry.num_batches, obj.handles.hSVPS.Geometry.patient_dir);
+            else
+%                 num_batches = 1;
+                str = input('Enter num of beam batches: ','s');
+                merge_beamlets(str2double(str), obj.handles.hSVPS.Geometry.patient_dir);
+            end
+            
             disp('Beamlets merged!')
         end
         %-------------------------------------------------------------------------------
@@ -563,6 +572,19 @@ classdef XTPS < handle
 %             obj.handles.hSVPS.optResults = optResults;
             
         end
+        
+        function RO_fullDose(obj, src, evt)
+            disp('Full dose calculation called ...')
+            
+            Pat_path = [obj.handles.hSVPS.Geometry.patient_dir];
+            [NLP_file,NLP_path,indx] = uigetfile([obj.handles.hSVPS.Geometry.patient_dir '\matlab_files\*.mat'], 'Select NLP file' );
+
+            NLP_getFullDose(NLP_path, NLP_file);
+            
+            disp('Full dose calculation completed!')
+            
+        end
+        
 
 %-------------------------------------------------------------------------------
         function LoadGeometryMenu_Callback(obj, src, evt)

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -8,6 +8,7 @@ for k = 1:num_batches
     batch_doses = read_ryan_beamlets([pat_dir, '\batch_dose' num2str(k-1) '.bin']);
     B = [B batch_doses];
 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
 f = waitbar(0, 'Correcting beamlet shift');
@@ -21,8 +22,8 @@ for beam_i = 1:numel(B)
     shifted_img = imtranslate(tabula, [-0.5,-0.5,-0.5]);
     nonzero_idx = find(shifted_img>0.005*max(shifted_img(:)));
 
-    B{1, beam_i}.non_zero_indices = nonzero_idx;
-    B{1, beam_i}.non_zero_values = shifted_img(nonzero_idx);
+    B{beam_i}.non_zero_indices = nonzero_idx;
+    B{beam_i}.non_zero_values = shifted_img(nonzero_idx);
 end
 close(f)
 

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -65,7 +65,7 @@ while true
     fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target. New start size: ' num2str(N_svox) '\n'])
     
     if converge_factor < 0.01
-        orthoslice(superMask)
+%         orthoslice(superMask)
         break
     end
 end 

planEvaluation/plot_DVH.m

@@ -26,7 +26,7 @@ function plot_DVH(Geometry, dose)
 %         / num_vox;
     title(['NLP optimized DVH'])
     legend(names)
-    axis([0 70 0 100])
+    axis([0 110 0 100])
 end
 function plot_DVH_robust(dose, optGoal, optGoal_idx)
     % this function plots the DVHs of the given dose