Added GUI for goal definition and a better interaction with the optimizer.

PF/Simple_DP_eval.m

@@ -1,43 +0,0 @@
-
-
-pat_path='C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\PatientData_dog5_HD';
-
-% geometry
-load([pat_path '\matlab_files\Geometry.mat'])
-load([pat_path '\matlab_files\optResults.mat'])
-
-
-CTin=Geometry.data;
-Din = optResults.dose{end};
-colorwash(CTin, Din, [-500, 500], [0,60])
-
-roi_max = size(Geometry.ROIS,2);
-
-for i= 1:roi_max
-    ROI = logical(zeros(size(CTin)));
-    ROI(Geometry.ROIS{1, 1}.ind) = 1;
-    orthoslice(ROI, [0,1])
-
-    [dvh, dosebins] = dvhist (Din, ROI);
-    figure
-    plot(dosebins, dvh, ...
-    'Color', [0.5,0.0,0.1] , ...
-    'LineStyle', '-', ...
-    'DisplayName', 'test');
-end
-
-
-
-for roi_idx = 1:numel(obj.handles.hSVPS.Geometry.ROIS)
-    if obj.handles.hSVPS.Geometry.ROIS{roi_idx}.visible == true % display == ON
-        fprintf('%s\n', obj.handles.hSVPS.Geometry.ROIS{roi_idx}.name);
-%                     [dvh dosebins] = dvhist(obj.handles.hSVPS.optResults.dose{end}, obj.handles.hSVPS.Geometry.ROIS{roi_idx}.ind);
-        % temp solution for karthik dosimetry data
-        [dvh dosebins] = dvhist(obj.handles.hSVPS.optResults.dose{end}, ...
-            obj.handles.hSVPS.Geometry, ...
-            roi_idx, ...
-            nfrac);                    
-        
-    end
-end
-

PF/make_DP_maps_Dog.m

@@ -1,46 +0,0 @@
-
-
-if true
-    %% load up NRRD
-
-    [petIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_PET.nrrd');
-    [segIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_SEG.nrrd');
-
-    segIdx = find(segIn);
-    
-    matrix = 0.4*60 + 0.6*(70/0.03)*petIn;
-%     matrix = 40 + 0*petIn;
-    
-    [X,Y,Z] = meshgrid(1:size(petIn, 1),1:size(petIn, 2),1:size(petIn, 3));
-    Y_vox = Y(segIdx);
-    matrix = 60*ones(size(Y));
-    matrix(segIdx) = 60 + 5*(mean(Y_vox)-Y(segIdx));
-    
-%     orthoslice(matrix)
-    
-    % save nrrds
-    filename = ['C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_DP_target.nrrd'];
-    pixelspacing = meta.spacedirections;
-    origin = meta.spaceorigin;
-    encoding ='raw'; % 'raw', 'ascii' or 'gzip'
-
-    ok = nrrdWriter(filename, matrix, pixelspacing, origin, encoding)
-
-end
-if false
-    
-    [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
-    
-    orthoslice(matrix)
-    orthoslice(D_full)
-    
-    colorwash(Geometry.data, matrix, [500, 1500], [0,70])
-    colorwash(Geometry.data, D_full, [500, 1500], [0,70])
-    
-    slice_i = 24;
-    
-    D_diff = D_full - matrix;
-    orthoslice(D_diff)
-    
-    
-end

PF/make_DP_maps_Tomo.m

@@ -1,54 +0,0 @@
-
-
-if true
-    %% load up NRRD
-
-    
-    load('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\matlab_files\Geometry.mat')
-    
-    CTin = Geometry.data;
-%     orthoslice(CTin)
-    
-    ROImap=zeros(size(CTin));
-    ROImap(Geometry.ROIS{1, 2}.ind)=1;
-    
-    doseMap = zeros(size(ROImap));
-    for z = 1:size(ROImap, 3)
-        doseMap(:,:,z) = bwdist(1-ROImap(:,:,z));
-    end
-%     orthoslice(doseMap)
-    
-    matrix = 50* ones(size(CTin))+2* doseMap;
-    
-%     matrix = 50* ones(size(CTin));
-%     temp = 50* ones(size(CTin)) + 2*(max(doseMap(:))-doseMap);
-%     matrix(Geometry.ROIS{1, 2}.ind)=temp(Geometry.ROIS{1, 2}.ind);
-%     orthoslice(matrix)
-    
-    
-    % save nrrds
-    filename = ['C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd'];
-    pixelspacing = Geometry.voxel_size;
-    origin = Geometry.start;
-    encoding ='raw'; % 'raw', 'ascii' or 'gzip'
-
-    ok = nrrdWriter(filename, matrix, pixelspacing, origin, encoding)
-
-end
-if false
-    
-    [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
-    
-    orthoslice(matrix)
-    orthoslice(D_full)
-    
-    colorwash(Geometry.data, matrix, [500, 1500], [0,70])
-    colorwash(Geometry.data, D_full, [500, 1500], [0,70])
-    
-    slice_i = 24;
-    
-    D_diff = D_full - matrix;
-    orthoslice(D_diff)
-    
-    
-end

ROI_goals_prep.m

@@ -7,10 +7,19 @@ function ROI_goals_prep
 % you want the plan to be!
 
 %% ---=== INPUT PARAMS ===---
-patient = 'avastin_009';
-outName = 'ROI_goals_DP';
+patient = 'medivation_01';
+outName = 'ROI_goals_simple';
 
 switch patient
+    case 'medivation_01'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2\';
+        paths.CT_in = ['CT_crop_ 1'];
+        paths.target_bin_in  = ['CTV_crop_ 1'];
+        paths.target_fzy_in  = ['CTV_crop_ 1'];
+        paths.body_bin_in    = ['CTV_crop_ 1'];
+        paths.wiscplan = 'F:\21_WiscPlan_data\Medivation_0391811_crop1_2';
+        
+        
     case 'avastin_009'
         paths.in = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
         paths.CT_in = ['AV009b_1ct_resized'];
@@ -84,6 +93,12 @@ fprintf('and beamlets.\n')
 
 %% ---=== GET OPTGOAL ===---
 switch patient
+    case 'medivation_01'
+        ROI_goals.optGoal = make_ROI_goals_medivation_011(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_medivation_011(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+        
     case 'avastin_009'
         ROI_goals.optGoal = make_ROI_goals_avastin_009(Geometry, beamlets);
         ROI_goals.optGoal_beam = make_ROI_goals_avastin_009(Geometry, beamlets_joined);
@@ -135,6 +150,63 @@ fprintf('Writing ROI_goals...')
 save([paths.in '\RODP_files\' outName '.mat'], 'ROI_goals')
 fprintf(' done!\n')
 end
+
+function optGoal = make_ROI_goals_medivation_011(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'PTV_min';
+    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
+    goal_1.ROI_idx = ROI_idx;
+    goal_1.imgDim = size(Geometry.data);
+    goal_1.D_final = 60;
+    goal_1.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target   = ones(numel(ROI_idx), 1) * 60;
+%     goal_1.target   = minDose(ROI_idx);
+    goal_1.opt_weight = 77 / numel(ROI_idx); % normalize to volume of target area
+    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_1;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_2.name = 'PTV_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
+    goal_2.ROI_idx = ROI_idx;
+    goal_2.imgDim = size(Geometry.data);
+    goal_2.D_final = 63;
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target   = ones(numel(ROI_idx), 1) * 63;
+%     goal_2.target   = maxDose(ROI_idx);
+    goal_2.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_2;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_3.name = 'ring_max';
+    goal_3.ROI_name = Geometry.ROIS{1, 3}.name;
+    ROI_idx = Geometry.ROIS{1, 3}.ind;
+    goal_3.ROI_idx = ROI_idx;
+    goal_3.imgDim = size(Geometry.data);
+    goal_3.D_final = 20;
+    goal_3.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target   = ones(numel(ROI_idx), 1) * 20;
+    goal_3.opt_weight = 5 / numel(ROI_idx); % normalize to volume of target area
+    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_3;
+    % -- END DEFINITION OF GOAL --
+   
+end
+
+
 function optGoal = make_ROI_goals_avastin_009_DP(Geometry, beamlets, minDose, maxDose)
     optGoal={};
     

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -1,11 +1,14 @@
 
 
-
+% 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
+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
@@ -15,31 +18,41 @@ function [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;
-patient = 'temp';
-goalsName= 'ROI_goals_DP';
-
-switch patient
-    case 'temp'
-        patient_dir = 'C:\000-temp\DELETEME';
-        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis'
-    case 'avastin_009'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
-        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
-    case 'gbm_005'
-%         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
-        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
-    case 'gbm_015'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
-        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
-    case 'gbm_022'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
-        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
-    case 'avastin_009'
-    otherwise
+% 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
 
@@ -51,18 +64,19 @@ end
 fprintf('starting NLP optimization process... \n')
 
 % % -- LOAD GEOMETRY AND BEAMLETS --
-load([patient_dir '\matlab_files\Geometry.mat'])
+load([paths.patient_dir '\matlab_files\Geometry.mat'])
 
 
 %% -- OPTIMIZATION TARGETS --
-load([DP_dir '\RODP_files\' goalsName '.mat']);
+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, patient_dir);
+[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, paths.patient_dir);
 
 
 % -- make them robust --
@@ -134,7 +148,7 @@ D_full = reshape(beamlets * w_fin, size(Geometry.data));
 %% save outputs
 NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
-save([patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
+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]);
@@ -153,8 +167,8 @@ function penalty = get_penalty(x, optGoal)
     
     for nrs_i = 1:optGoal{1}.NbrRandScenarios 
         for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
-            for rrs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
-                fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rrs_i);
+            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
@@ -163,7 +177,7 @@ function penalty = get_penalty(x, optGoal)
     penalty=max(fobj);
 end
 % ------ supp: penalty for single scenario ------
-function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
+function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
     penalty = 0;
     % for each condition
     for goal_i = 1:numel(optGoal)
@@ -172,41 +186,41 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
             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}.rrs{rrs_i}.target) -...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)));
+                    (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}.rrs{rrs_i}.beamlets_pruned * x)-...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target)));
+                    (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}.rrs{rrs_i}.beamlets_pruned * x)-...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target));
+                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}.rrs{rrs_i}.beamlets_pruned * x)-...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target));
+                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}.rrs{rrs_i}.beamlets_pruned * x) - ...
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target).^2);
+                    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}.rrs{rrs_i}.target) -...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) > 0)) ...
-                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
+                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}.rrs{rrs_i}.target) -...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) < 0)) ...
-                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
+                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
@@ -222,7 +236,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % RO_params - should have the information below
     % nrs - random scenarios
     % sss - system setup scenarios
-    % rrs - random range scenarios
+    % rgs - random range scenarios
     
     % X - X>0 moves image right
     % Y - Y>0 moves image down
@@ -241,16 +255,16 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 
     % ----====#### CHANGE ROBUSTNESS HERE ####====----
 
-    rrs_scene_list={[0,0,0]};
+    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\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(rrs_scene_list);
+        optGoal{i}.NbrRangeScenarios    =numel(rgs_scene_list);
     end
     
 
@@ -271,7 +285,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
                 [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
                 % beamlets stay the same
                 
-                for rrs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rrs
+                for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rgs
                     % modify target and beamlets
                     targetImg4=targetImg3;
                     % beamlets stay the same
@@ -293,9 +307,9 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
                     beamlets_pruned = beamlets(ROI_idx, :);
                     
                     % save to optGoal output
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.ROI_idx            = ROI_idx;
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned    = beamlets_pruned;
-                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target             = target;
+                    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

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v2.m

@@ -0,0 +1,311 @@
+
+
+% 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
+% Inputs: (Pat_path, path2goal) or none
+%   If paths are used they should be passed as strings.
+% Outputs: full dose image dose: [D_full, w_fin, Geometry, optGoal]
+% 
+% [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
+% 
+% 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};
+end
+
+N_fcallback1 = 10000;
+N_fcallback2 = 200000;
+
+%% 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)
+    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, :));
+end
+
+% optGoal_idx = ROI_goals.optGoal_idx;
+% targetMinMax_idx = ROI_goals.targetMinMax_idx;
+
+% -- 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}.D_final(optGoal{1}.ROI_idx) ./ (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([Pat_path, '\matlab_files\NLP_result.mat'], 'NLP_result');
+
+warning('this part needs modification')
+plot_DVH(D_full, optGoal)
+colorwash(Geometry.data, D_full, [-500, 500], [0, 60]);
+% 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);
+                    
+                    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/dicomrt2geometry.m

@@ -16,13 +16,25 @@ function Geometry = dicomrt2geometry(dicomrt_dir)
 % See also readPinnGeometry.m
 %
 % Author: Xiaohu Mo
+% edits: Peter Ferjancic 2019/01
 
 Geometry = [];
 
 if nargin < 1
-    dicomrt_dir = uifile('getdir', 'Select the directory contains the DICOM-RT images');
+%     dicomrt_dir = uifile('getdir', 'Select the directory contains the DICOM-RT images');
 % !!! Grozomah
-%     dicomrt_dir = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\CT with Contours\DICOM\TomoPhantPlan01';
+    load('WiscPlan_preferences.mat')
+    if isstring(WiscPlan_preferences.inDataPath) 
+    else
+        WiscPlan_preferences.inDataPath = 'C://';
+    end
+    dicomrt_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the directory contains the DICOM-RT images');
+    
+    WiscPlan_preferences.inDataPath = dicomrt_dir;
+    thisDir = mfilename('fullpath');
+    idcs   = strfind(thisDir,'\');
+    prefsdir = thisDir(1:idcs(end-1)-1);
+    save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences')
 end
 
 if exist(dicomrt_dir, 'dir') ~= 7
@@ -45,7 +57,8 @@ filenames = cell(0);
 % remove directories
 for i = 1:numel(dirlist)
     [tilde, tilde, ext] = fileparts(dirlist(i).name);
-    if ~dirlist(i).isdir && strcmpi(ext, '.dcm')
+%     if ~dirlist(i).isdir && strcmpi(ext, '.dcm') 
+    if ~dirlist(i).isdir && ~strcmpi(ext, '.txt')   % <--- Grozomah: Fix this back to default!
         filenames{end+1} = dirlist(i).name;
     end
 end
@@ -106,6 +119,7 @@ if downsample_flag == true
         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/4;
         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/dicomrt/dicomrt_loadct.m

@@ -72,8 +72,9 @@ end
 
 % Now get CT images and create 3D Volume
 
-% Set parameters
-CToffset=1000;
+% Set parameters 
+CToffset=1000;  % CT offset!
+% CToffset=0;  % CT offset!
 nct=0;
 
 % Define cell array to store info

WiscPlanPhotonkV125/matlab_frontend/dicomrt/dicomrt_rasterizeVoi.m

@@ -29,6 +29,11 @@ for kk = 1:length(voi2use)
     rasterVOI = zeros(xdim, ydim, zdim, 'int8');
     num_curves = numel(VOI{kk,2});
     for k = 1:num_curves
+        
+%         if k == 6 %% manual fix for hippocampus
+%             continue
+%         end
+        
         % the N x 3 curve matrix
         curve = VOI{kk,2}{k};
         % skip empty curves
@@ -48,8 +53,8 @@ for kk = 1:length(voi2use)
             error('Closest slice found at %f cm', min(abs(ct_zmesh-voi_z)));
         elseif numel(slice) > 1
             warning('Multiple matching slices found');
-        elseif dist_z(slice) >= abs(ct_zmesh(1)-ct_zmesh(2))
-            error('Closest slice found at %f cm', min(abs(ct_zmesh-voi_z)));
+%         elseif dist_z(slice) >= abs(ct_zmesh(1)-ct_zmesh(2)) % grozomah
+%             error('Closest slice found at %f cm', min(abs(ct_zmesh-voi_z)));
         end
         % binary mask of current curve
         if exist('InPolygon', 'file') == 3

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -0,0 +1,257 @@
+function varargout = goal_def_UI(varargin)
+% GOAL_DEF_TEST MATLAB code for goal_def_test.fig
+%      GOAL_DEF_TEST, by itself, creates a new GOAL_DEF_TEST or raises the existing
+%      singleton*.
+%
+%      H = GOAL_DEF_TEST returns the handle to a new GOAL_DEF_TEST or the handle to
+%      the existing singleton*.
+%
+%      GOAL_DEF_TEST('CALLBACK',hObject,eventData,handles,...) calls the local
+%      function named CALLBACK in GOAL_DEF_TEST.M with the given input arguments.
+%
+%      GOAL_DEF_TEST('Property','Value',...) creates a new GOAL_DEF_TEST or raises the
+%      existing singleton*.  Starting from the left, property value pairs are
+%      applied to the GUI before goal_def_test_OpeningFcn gets called.  An
+%      unrecognized property name or invalid value makes property application
+%      stop.  All inputs are passed to goal_def_test_OpeningFcn via varargin.
+%
+%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
+%      instance to run (singleton)".
+%
+% See also: GUIDE, GUIDATA, GUIHANDLES
+
+% Edit the above text to modify the response to help goal_def_test
+
+% Last Modified by GUIDE v2.5 28-Mar-2019 16:16:14
+
+% Begin initialization code - DO NOT EDIT
+gui_Singleton = 1;
+gui_State = struct('gui_Name',       mfilename, ...
+                   'gui_Singleton',  gui_Singleton, ...
+                   'gui_OpeningFcn', @goal_def_test_OpeningFcn, ...
+                   'gui_OutputFcn',  @goal_def_test_OutputFcn, ...
+                   'gui_LayoutFcn',  [] , ...
+                   'gui_Callback',   []);
+if nargin && ischar(varargin{1})
+    gui_State.gui_Callback = str2func(varargin{1});
+end
+
+if nargout
+    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
+else
+    gui_mainfcn(gui_State, varargin{:});
+end
+% End initialization code - DO NOT EDIT
+end
+
+% --- Executes just before goal_def_test is made visible.
+function goal_def_test_OpeningFcn(hObject, eventdata, handles, varargin)
+% This function has no output args, see OutputFcn.
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+% varargin   command line arguments to goal_def_test (see VARARGIN)
+
+% Choose default command line output for goal_def_test
+handles.output = hObject;
+
+%% GET THE INITIAL GEOMETRY DATA
+disp('varargin')
+if numel(varargin) == 0
+    load('WiscPlan_preferences.mat')
+    [handles.Data.Geo_fileName,handles.Data.Geo_path,FilterIndex] = uigetfile([WiscPlan_preferences.patientDataPath '\matlab_files\*.mat'], 'Select Geometry file')
+    load([handles.Data.Geo_path, handles.Data.Geo_fileName])
+end
+
+
+%% POPULATE THE TABLE 
+% == populate ROI options ==
+for i = 1: numel(Geometry.ROIS)
+    handles.uitable1.ColumnFormat{2}{i} = Geometry.ROIS{i}.name;
+end
+% == populate function options ==
+handles.uitable1.ColumnFormat{5}{1} = 'min';
+handles.uitable1.ColumnFormat{5}{2} = 'max';
+handles.uitable1.ColumnFormat{5}{3} = 'min_sq';
+handles.uitable1.ColumnFormat{5}{4} = 'max_sq';
+handles.uitable1.ColumnFormat{5}{5} = 'LeastSquare';
+handles.uitable1.ColumnFormat{5}{6} = 'min_perc_Volume';
+handles.uitable1.ColumnFormat{5}{7} = 'max_perc_Volume';
+
+% == populate the first entry ==
+handles.uitable1.Data = handles.uitable1.Data(1,:);
+handles.uitable1.Data{1} = 'Goal 1';
+handles.uitable1.Data{2} = 'Target';
+handles.uitable1.Data{3} = 'Fixed dose';
+handles.uitable1.Data{4} = '60';
+handles.uitable1.Data{5} = 'min_sq';
+handles.uitable1.Data{6} = 100;
+handles.uitable1.Data{7} = 'null';
+
+% Update handles structure
+guidata(hObject, handles);
+
+% UIWAIT makes goal_def_test wait for user response (see UIRESUME)
+% uiwait(handles.figure1);
+end
+
+% --- Outputs from this function are returned to the command line.
+function varargout = goal_def_test_OutputFcn(hObject, eventdata, handles) 
+% varargout  cell array for returning output args (see VARARGOUT);
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Get default command line output from handles structure
+varargout{1} = handles.output;
+end
+
+% --- Executes on button press in addRow.
+function addRow_Callback(hObject, eventdata, handles)
+% adds new empty row to table
+handles.uitable1.Data{end+1, 1} = ['Goal ' num2str(size(handles.uitable1.Data,1)+1)];
+end
+% --- Executes on button press in removeRow.
+function removeRow_Callback(hObject, eventdata, handles)
+% Removes row number selected by Remove_row_N
+newTable = handles.uitable1.Data;
+rowToDelete = str2double(handles.Remove_row_N.String);
+if rowToDelete>size(newTable, 2)
+    disp(['No row ' num2str(rowToDelete)])
+else
+    newTable(rowToDelete,:) = [];
+    handles.uitable1.Data = newTable;
+end
+end
+function Remove_row_N_Callback(hObject, eventdata, handles)
+% Identifies what is the number for removing a row and checks validty of
+% input.
+num = str2double(handles.Remove_row_N.String);
+if isnan(num)
+    disp('You''re a silly goose. Enter a number.')
+elseif num<1
+    disp('Row number must be greater than 0')
+    num = 1;
+    handles.Remove_row_N.String = 1;
+else
+    r = rem(num, 1);
+    if r ~= 0
+        disp('Enter an integer, duh.')
+        num = round(num);
+        handles.Remove_row_N.String = round(num);
+    end
+    disp(num)
+end
+end
+
+% --- Executes on button press in Save_goal_button.
+function Save_goal_button_Callback(hObject, eventdata, handles)
+% hObject    handle to Save_goal_button (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+Save_goal_function(handles);
+end
+function Save_goal_function(handles)
+% this function creates and outputs the OptGoals file
+
+[file,path,indx] = uiputfile([handles.Data.Geo_path, '*.mat'],'Choose where to save OptGoals');
+OptGoals.geometryPath = [handles.Data.Geo_path handles.Data.Geo_fileName];
+
+disp('Loading Geometry ...')
+load([OptGoals.geometryPath])
+% disp('Loading beamlets ...')
+% [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, Geometry.patient_dir);
+
+OptGoals.data={};
+for i = 1 : size(handles.uitable1.Data, 1)
+    % -- START DEFINITION OF GOAL --
+    % == goal name
+    goal_i.name = handles.uitable1.Data{i,1};
+    % == ROI name
+    goal_i.ROI_name = handles.uitable1.Data{i,2};
+    % == ROI index
+    for j = 1:size(Geometry.ROIS,2)
+        if strcmp(Geometry.ROIS{j}.name, handles.uitable1.Data{i,2})
+            goal_i.ROI_idx = Geometry.ROIS{j}.ind;
+            break
+        end
+        if j == size(Geometry.ROIS,2)
+            error(['Invalid ROI name selected in goal ' num2str(i)])
+        end
+    end
+    % == image dimensions
+    goal_i.imgDim = size(Geometry.data);
+    % == goal dose
+    if strcmp (handles.uitable1.Data{i,3}, 'Fixed dose')
+        goal_i.D_final = str2double(handles.uitable1.Data{i,4}) * ones(size(Geometry.data));
+    elseif strcmp (handles.uitable1.Data{i,3}, 'Dose map')
+        error('havent written this part yet')
+        % loading of data comes here
+        % matrix size checks
+    end
+    % == goal penalty function
+    goal_i.function = handles.uitable1.Data{i,5};
+    % == goal weight
+    goal_i.opt_weight = handles.uitable1.Data{i,6} / numel(goal_i.ROI_idx); % normalize to volume of target area
+    % == beamlets dose selection
+%     goal_i.beamlets_pruned = beamlets(goal_i.ROI_idx, :);
+    goal_i.optionalParam = handles.uitable1.Data{i,7};
+    % -- END DEFINITION OF GOAL -- 
+    goal_i.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    OptGoals.data{end+1}=goal_i;
+
+end
+disp('-- all OptGoals exported!')
+
+OptGoals.goals = [handles.uitable1.Data];
+
+save( [path, file], 'OptGoals')
+end
+
+
+% --- Executes on button press in Load_goal_button.
+function Load_goal_button_Callback(hObject, eventdata, handles)
+% hObject    handle to Load_goal_button (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+[file,path,indx] = uigetfile([handles.Data.Geo_path, '*.mat'],'Choose OptGoals to load');
+load( [path, file])
+disp('load goal')
+
+handles.uitable1.Data = OptGoals.goals;
+
+end
+
+
+% --- Executes when entered data in editable cell(s) in uitable1.
+function uitable1_CellEditCallback(hObject, eventdata, handles)
+% hObject    handle to uitable1 (see GCBO)
+% eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
+%	Indices: row and column indices of the cell(s) edited
+%	PreviousData: previous data for the cell(s) edited
+%	EditData: string(s) entered by the user
+%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
+%	Error: error string when failed to convert EditData to appropriate value for Data
+% handles    structure with handles and user data (see GUIDATA)
+
+if eventdata.Indices(2) == 3;
+    mode = handles.uitable1.Data{eventdata.Indices(1), eventdata.Indices(2)}
+    switch mode
+        case 'Dose map'
+            [FileName,PathName,FilterIndex] = uigetfile('F:\021_WiscPlan_data\FET_repeat_005_1\matlab_files\*.mat', 'Select dose reference file')
+            warning('loading of said profile comes here')
+            handles.uitable1.Data{eventdata.Indices(1), 4} = [PathName, FileName];
+            
+        case 'Fixed dose'
+            answer = inputdlg('Enter desired goal dose');
+            handles.uitable1.Data{eventdata.Indices(1), 4} = answer{1};
+            
+        otherwise
+            error('This doesnt work')
+    end
+end
+
+
+
+end

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -23,12 +23,20 @@
 function num_batches = helicalDosecalcSetup7(patient_dir)
 % -- INPUT:
 % patient_dir: specify where kernel and [geometry] files are located
-num_batches = 3; % number of cores you want to run the beam calculation
 
 iso = [0 0 0]; % Point about which gantry rotations begin
 SAD = 85;  % source-axis distance for the x-ray source ##
 pitch = 0.86; % fraction of beam with couch translates per rotation
 
+%% --- make the figure prompt for number of angles and beamlets
+str = inputdlg({'Enter number of calc cores', 'Enter number of angles (51 default)', ...
+    'Enter number of beamlets (64 default)'}, 'input', [1,35], {'3', '51', '64'});
+num_batches = str2double(str{1}); % number of cores you want to run the beam calculation 
+% -- (3 for a 4-core comp to prevent lockdown)
+N_angles = str2double(str{2}); % 51 for full resolution
+Mxp = str2double(str{3}); % Mxp = 64;  number of MLC leaves;
+Nyp = 1;  % always 1 for Tomo due to binary mlc
+
 % define the overall beam field size for each beam angle
 % xpmin = -0.5;     % -field width / 2
 % xpmax = 0.5;      % +field width / 2
@@ -43,11 +51,6 @@ ypmax = 0.3125;
 
 % y-prime points in the z-direction in the CT coordinate system
 
-% Number of beamlets in the BEV for each direction
-Mxp = 32;  % Mxp = 64;  number of MLC leaves;
-Nyp = 1;  % always 1 for Tomo due to binary mlc
-
-
 % ============================================= End of user-supplied inputs
 executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';
 kernel_file = 'Kernels.mat';
@@ -70,8 +73,8 @@ targetMask(Geometry.ROIS{target_idx}.ind) = 1;
 % the target
 
 targetMaskZ = sum(sum(targetMask,1),2);
-zBow = find(targetMaskZ>0, 1, 'first')*Geometry.voxel_size(3) + Geometry.start(3) + ypmin;
-zStern = find(targetMaskZ>0, 1, 'last')*Geometry.voxel_size(3) + Geometry.start(3) + ypmax;
+zBow = (find(targetMaskZ>0, 1, 'first')-1)*Geometry.voxel_size(3) + Geometry.start(3) + ypmin;
+zStern = (find(targetMaskZ>0, 1, 'last')+1)*Geometry.voxel_size(3) + Geometry.start(3) + ypmax;
 [subi, subj, subk] = ind2sub(size(Geometry.data), Geometry.ROIS{target_idx}.ind);
 iso = [Geometry.start(1)+Geometry.voxel_size(1)*mean(subi) ...
     Geometry.start(2)+Geometry.voxel_size(2)*mean(subj) 0];
@@ -88,7 +91,7 @@ Nrot = ceil(abs(zBow - zStern)/(pitch*fieldWidth));
 
 
 % Nphi = Nrot*51;  % number of angles used in the calculation
-Nphi = Nrot*21;  % Grozomah
+Nphi = Nrot * N_angles;  % Grozomah
 
 % define the limits of the angles that will be used for the calculation
 % ##phimin = 0;  % starting angle in radians

WiscPlanPhotonkV125/matlab_frontend/lab/RDXTPS_geometry_setup_wizard_ASP.m

@@ -41,84 +41,94 @@ end
 
 
 %% -----===== BTV and Ring creation =====-----
-% waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
-% ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
-% 
-% [target_idx okay] = listdlg('ListString', ROI_names, ...
-%     'SelectionMode', 'single', 'Name', 'Target Selection', ...
-%     'PromptString', 'Please select the target ROI. ');
-% if okay ~= 1
-%     msgbox('Plan creation aborted');
-%     delete(hWaitbar);
-%     return;
-% end
-% 
-% [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
-%     'BTV margin specification', 1, {'0.6'});
-% if isempty(BTV_margin_answer)
-%     BTV_margin = 0.6;
-% else
-%     BTV_margin = str2double(BTV_margin_answer{1});
-% end
-% 
-% % target_idx and BTV_margin are set. Expand PTV to BTV
-% PTVmask = false(size(Geometry.rhomw));
-% % for target_idx = target_indices
-%         PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
-% % end
-% if BTV_margin > 0
-%     if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
-%         bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-%         Geometry.BTV = bwD <= BTV_margin;
-%     end
-% end
-% 
-% % Create btv
-% Geometry.ROIS{end+1} = Geometry.ROIS{end};
-% Geometry.ROIS{end}.name = 'BTV';
-% Geometry.ROIS{end}.num_curves = 0;
-% Geometry.ROIS{end}.curves = {};
-% Geometry.ROIS{end}.ind = find(Geometry.BTV);
-% Geometry.ROIS{end}.visible = false;
-% 
-% % Create ring
-% [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
-%     'Ring margin specification', 1, {'1'});
-% if isempty(ring_margin_answer)
-%     ring_margin = 1;
-% else
-%     ring_margin = str2double(ring_margin_answer{1});
-% end
-% bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-% Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
-% Geometry.Ring = xor(Geometry.Ring, PTVmask);
-% Geometry.ROIS{end+1} = Geometry.ROIS{end};
-% Geometry.ROIS{end}.name = 'Ring';
-% Geometry.ROIS{end}.num_curves = 0;
-% Geometry.ROIS{end}.curves = {};
-% Geometry.ROIS{end}.ind = find(Geometry.Ring);
-% Geometry.ROIS{end}.visible = false;
-
-% % Create high comformity ring structure for adaptive step size
-% [HCRing_margin_answer] = inputdlg({sprintf('Please enter the HC ring margin (cm):\n(default 1.2 cm or 2 sigma, enter 0 to skip)')}, ...
-%     'HCRing margin specification', 1, {'1.2'});
-% if isempty(HCRing_margin_answer)
-%     HCRing_margin = 1.2;
-% else
-%     HCRing_margin = str2double(HCRing_margin_answer{1});
-% end
-% 
-% if HCRing_margin > 0
-%     bwD = bwdistsc(~Geometry.BTV, Geometry.voxel_size);
-%     % shrinked BTV from original BTV, such that original_BTV > target_VOI > shrinked_BTV
-%     Geometry.SBTV = bwD >= HCRing_margin;
-%     Geometry.HCRing = xor(Geometry.BTV, Geometry.SBTV);
+waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
+ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
+
+[target_idx okay] = listdlg('ListString', ROI_names, ...
+    'SelectionMode', 'single', 'Name', 'Target Selection', ...
+    'PromptString', 'Please select the target ROI. ');
+if okay ~= 1
+    msgbox('Plan creation aborted');
+    delete(hWaitbar);
+    return;
+end
+
+[BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
+    'BTV margin specification', 1, {'0.6'});
+if isempty(BTV_margin_answer)
+    BTV_margin = 0.6;
+else
+    BTV_margin = str2double(BTV_margin_answer{1});
+end
+
+% target_idx and BTV_margin are set. Expand PTV to BTV
+PTVmask = false(size(Geometry.rhomw));
+% for target_idx = target_indices
+        PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
 % end
+if BTV_margin > 0
+    if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
+        bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+        Geometry.BTV = bwD <= BTV_margin;
+    end
+end
+
+% Create btv
+Geometry.ROIS{end+1} = Geometry.ROIS{end};
+Geometry.ROIS{end}.name = 'BTV';
+Geometry.ROIS{end}.num_curves = 0;
+Geometry.ROIS{end}.curves = {};
+Geometry.ROIS{end}.ind = find(Geometry.BTV);
+Geometry.ROIS{end}.visible = false;
+
+% Create ring
+[ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
+    'Ring margin specification', 1, {'1'});
+if isempty(ring_margin_answer)
+    ring_margin = 1;
+else
+    ring_margin = str2double(ring_margin_answer{1});
+end
+bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
+Geometry.Ring = xor(Geometry.Ring, PTVmask);
+Geometry.ROIS{end+1} = Geometry.ROIS{end};
+Geometry.ROIS{end}.name = 'Ring';
+Geometry.ROIS{end}.num_curves = 0;
+Geometry.ROIS{end}.curves = {};
+Geometry.ROIS{end}.ind = find(Geometry.Ring);
+Geometry.ROIS{end}.visible = false;
+
+% Create high comformity ring structure for adaptive step size
+[HCRing_margin_answer] = inputdlg({sprintf('Please enter the HC ring margin (cm):\n(default 1.2 cm or 2 sigma, enter 0 to skip)')}, ...
+    'HCRing margin specification', 1, {'1.2'});
+if isempty(HCRing_margin_answer)
+    HCRing_margin = 1.2;
+else
+    HCRing_margin = str2double(HCRing_margin_answer{1});
+end
+
+if HCRing_margin > 0
+    bwD = bwdistsc(~Geometry.BTV, Geometry.voxel_size);
+    % shrinked BTV from original BTV, such that original_BTV > target_VOI > shrinked_BTV
+    Geometry.SBTV = bwD >= HCRing_margin;
+    Geometry.HCRing = xor(Geometry.BTV, Geometry.SBTV);
+end
 
 %% -----===== Save geometry files =====-----
 waitbar(2/total_num_steps, hWaitbar, 'Step 3: Save geometry files');
 
-patient_dir = uifile('getdir', 'Save the patient data to directory');
+load('WiscPlan_preferences.mat')
+patient_dir = uigetdir(WiscPlan_preferences.patientDataPath, 'Save the patient data to directory');
+
+WiscPlan_preferences.patientDataPath = patient_dir;
+thisDir = mfilename('fullpath');
+idcs   = strfind(thisDir,'\');
+prefsdir = thisDir(1:idcs(end-1)-1);
+save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences');
+
+
+% patient_dir = uifile('getdir', 'Save the patient data to directory');
 % !!! Grozomah
 % patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
 

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -49,22 +49,21 @@ classdef XTPS < handle
 
             % File menu
             obj.handles.hFileMenu             = uimenu(obj.handles.hMainFigure, 'Label','File');
-            obj.handles.hImportDicomRTMenu    = uimenu(obj.handles.hFileMenu, 'Label', '1. Import DicomRT Geometry');
-
-            % ---- Grozomah: This part edited by Peter
-            obj.handles.hNRRD1                = uimenu(obj.handles.hFileMenu, 'Label', '1a. Geometry wizard (am/nrrd)');
+            obj.handles.hImportDicomRTMenu    = uimenu(obj.handles.hFileMenu, 'Label', '1a. Import DicomRT Geometry');
+            obj.handles.hNRRD1                = uimenu(obj.handles.hFileMenu, 'Label', '1b. Geometry wizard (am/nrrd)');
             obj.handles.hNRRD2                = uimenu(obj.handles.hFileMenu, 'Label', '1b. Add Geometry ROI');
-            obj.handles.hImgPerturb           = uimenu(obj.handles.hFileMenu, 'Label', '1e. Implement image pertrubations');
+            obj.handles.hImgDownsample        = uimenu(obj.handles.hFileMenu, 'Label', '1(c). Downsample data');
+            
+            obj.handles.hDosecalcMenu         = uimenu(obj.handles.hFileMenu, 'Label', '2.1 Beamlet Dose Calculation');
+            obj.handles.hMergeBeamlets        = uimenu(obj.handles.hFileMenu, 'Label', '2.2 Merge Beamlets');
             
-            obj.handles.hDosecalcMenu         = uimenu(obj.handles.hFileMenu, 'Label', '2. Beamlet Dose Calculation');
+            obj.handles.hOptimMenu            = uimenu(obj.handles.hFileMenu, 'Label', '3a.1 Classical Optimization');
+            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');
             
-            % ---- Grozomah: This part edited by Peter
-            obj.handles.hMergeBeamlets                = uimenu(obj.handles.hFileMenu, 'Label', '2a. Merge Beamlets');
-            obj.handles.hTest2                = uimenu(obj.handles.hFileMenu, 'Label', '2b. Test 2');
+            obj.handles.hRO_GoalSpec          = uimenu(obj.handles.hFileMenu, 'Label', '3b. RO Goal Spec');
+            obj.handles.hRO_Optimize          = uimenu(obj.handles.hFileMenu, 'Label', '3b. RO Optimization');
             
-            obj.handles.hOptimMenu            = uimenu(obj.handles.hFileMenu, 'Label', '3. Plan Optimization');
-            obj.handles.hSaveOptimResultsMenu = uimenu(obj.handles.hFileMenu, 'Label', '4.a Save Optimization Results');
-            obj.handles.hGetFullDose          = uimenu(obj.handles.hFileMenu, 'Label', '4.b Get Full Dose');
             
             obj.handles.hLoadGeometryMenu     = uimenu(obj.handles.hFileMenu, 'Label', 'Load Geometry', 'Separator','on');
             obj.handles.hLoadOptresultsMenu   = uimenu(obj.handles.hFileMenu, 'Label', 'Load OptResults');
@@ -90,12 +89,13 @@ classdef XTPS < handle
             % Grozomah
             set(obj.handles.hNRRD1, 'Callback',         @obj.NRRD1_Callback);
             set(obj.handles.hNRRD2, 'Callback',         @obj.NRRD2_Callback);
-            set(obj.handles.hImgPerturb, 'Callback',    @obj.ImgPerturb_Callback);
+            set(obj.handles.hImgDownsample, 'Callback',    @obj.ImgDownsample_Callback);
+            
+            set(obj.handles.hRO_GoalSpec, 'Callback',   @obj.RO_GoalSpec);
+            set(obj.handles.hRO_Optimize, 'Callback',    @obj.RO_Optimize);
             
-             
             % Grozomah
             set(obj.handles.hMergeBeamlets, 'Callback',      @obj.MergeBeamlets_Callback);
-            set(obj.handles.hTest2, 'Callback',              @obj.Test2_Callback);
             set(obj.handles.hGetFullDose, 'Callback',        @obj.GetFullDose_Callback);
             
             set(obj.handles.hSaveOptimResultsMenu, 'Callback',      @obj.SaveOptimResultsMenu_Callback);
@@ -410,6 +410,7 @@ classdef XTPS < handle
         end
 %-------------------------------------------------------------------------------
         function ImportDicomRTMenu_Callback(obj, src, evt)
+            
             [obj.handles.hSVPS.Geometry obj.patient_dir] = RDXTPS_geometry_setup_wizard_ASP();
             obj.load_geometry();
             
@@ -441,9 +442,11 @@ classdef XTPS < handle
         end
         
 %-------------------------------------------------------------------------------
-        % Grozomah
-        function ImgPerturb_Callback(obj, src, evt)
-            disp('Perturb callback implemented')
+% Grozomah
+        function ImgDownsample_Callback(obj, src, evt)
+            disp('Downsampling data ... to be implemented.')
+            
+            
         end
         
         
@@ -473,25 +476,19 @@ classdef XTPS < handle
 %             system(execline);
 
             obj.handles.hSVPS.Geometry.num_batches = helicalDosecalcSetup7(obj.handles.hSVPS.Geometry.patient_dir)
-
-
         end
         %-------------------------------------------------------------------------------
-        % Grozomah
+% Grozomah
         function MergeBeamlets_Callback(obj, src, evt)
             merge_beamlets(obj.handles.hSVPS.Geometry.num_batches, obj.handles.hSVPS.Geometry.patient_dir);
             disp('Beamlets merged!')
         end
         %-------------------------------------------------------------------------------
-        % Grozomah
-        function Test2_Callback(obj, src, evt)
-            disp('Test 2 works!')
-        end
-%-------------------------------------------------------------------------------
         function OptimMenu_Callback(obj, src, evt)
             % TODO fix patient_dir system
 %             if isempty(obj.patient_dir)
-                obj.patient_dir = uifile('getdir', 'Select the patient directory');
+            load('WiscPlan_preferences.mat')
+            obj.patient_dir = uigetdir([WiscPlan_preferences.patientDataPath], 'Select the patient directory');
 %             end
 
             source_filename = fullfile(obj.patient_dir, 'batch_dose.bin');
@@ -531,19 +528,37 @@ classdef XTPS < handle
         function SaveOptimResultsMenu_Callback(obj, src, evt)
             % TODO fix patient_dir system
 %             if isempty(obj.patient_dir)
-                obj.patient_dir = uifile('getdir', 'Select the patient directory');
+            load('WiscPlan_preferences.mat')
+            obj.patient_dir = uigetdir([WiscPlan_preferences.patientDataPath], 'Select the patient directory');
+%                 obj.patient_dir = uifile('getdir', 'Select the patient directory');
 %             end
 
             [tilde, obj.handles.hSVPS.optResults] = loadOptResults(obj.patient_dir, 'optInput.txt');
             obj.load_optResults();
         end
 %-------------------------------------------------------------------------------
-        % Grozomah
+% Grozomah
         function GetFullDose_Callback(obj, src, evt)
             get_full_dose([obj.patient_dir])
             disp('Full dose calculated.')
         end
-        
+%-------------------------------------------------------------------------------
+% Grozomah
+        function RO_GoalSpec(obj, src, evt)
+            disp('RO Goal spec called!')
+            goal_def_UI
+        end
+%-------------------------------------------------------------------------------
+% Grozomah
+        function RO_Optimize(obj, src, evt)
+            disp('RO Optimization called!')
+            
+            Pat_path = [obj.handles.hSVPS.Geometry.patient_dir];
+            [Obj_file,Obj_path,indx] = uigetfile([obj.handles.hSVPS.Geometry.patient_dir '\matlab_files\*.mat'], 'Select OptGoal file' );
+            path2goal = [Obj_path, Obj_file];
+            [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v2(Pat_path, path2goal);
+        end
+
 %-------------------------------------------------------------------------------
         function LoadGeometryMenu_Callback(obj, src, evt)
             [filename pathname] = uifile('get', '*.mat', 'Choose Geometry file to load');
@@ -554,11 +569,13 @@ classdef XTPS < handle
         end
 %-------------------------------------------------------------------------------
         function LoadOptresultsMenu_Callback(obj, src, evt)
-            [filename pathname] = uifile('get', '*.mat', 'Choose OptResults file to load');
+            load('WiscPlan_preferences.mat')
+            [FileName,PathName,FilterIndex] = uigetfile([WiscPlan_preferences.patientDataPath '\*.mat'], 'Choose OptResults file to load');
+%             [filename pathname] = uifile('get', '*.mat', 'Choose OptResults file to load');
             % do nothing if cancelled
-            if isscalar(filename) && filename == 0; return; end
+            if isscalar(FileName) && FileName == 0; return; end
 
-            obj.load_optResults(fullfile(pathname, filename));
+            obj.load_optResults(fullfile(PathName, FileName));
             
 %             % Extra dose scaling. Temporary code for comparison with Tomo plans
 %             answer = inputdlg({'Dose', '% Volume'}, 'Scale final dose');

WiscPlanPhotonkV125/matlab_frontend/load2geometry.m

@@ -3,7 +3,25 @@
 function [Geometry] = load2geometry
     
     %% -----===== Get CT files =====-----
-    Geometry.data_dir = uigetdir('C:\Box Sync\Optimal Stopping\Data\CDP', 'Select folder with patient data');
+    
+    load('WiscPlan_preferences.mat')
+    if ~isfield(WiscPlan_preferences,'inDataPath')
+        Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
+    elseif WiscPlan_preferences.inDataPath ~= 0
+        Geometry.data_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the directory with patient data');
+    else
+        Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
+    end
+    WiscPlan_preferences.inDataPath = Geometry.data_dir;
+    thisDir = mfilename('fullpath');
+    idcs   = strfind(thisDir,'\');
+    prefsdir = thisDir(1:idcs(end)-1);
+    save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences')
+    
+%     Geometry.data_dir = uigetdir('C:\Box Sync\Optimal Stopping\Data\CDP', 'Select folder with patient data');
+%     Geometry.data_dir = uigetdir('\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2', 'Select folder with patient data');
+    
+
     [IMG_in, IMG_path, filterIdx] = uigetfile([{'*.nrrd'; '*.am'}], 'Select CT image', Geometry.data_dir);
 %     'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\'
     
@@ -31,15 +49,34 @@ function [Geometry] = load2geometry
             return 
         case 1
             [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
+            
+            % offset the CT, because that is what WiscPlan wants
+            Geometry.data = Geometry.data + 1000;
+            
             disp('NRRD file loaded!')
             Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
             Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
             [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
-            Geometry.voxel_size = meta.spacedirections;
+            
+            disp(['current voxel size: ' num2str(meta.spacedirections(1)) ' ' ...
+                num2str(meta.spacedirections(2)) ' ' num2str(meta.spacedirections(3)) ' cm.'])
+            
+            x = input('Is this correct? (y/n)', 's')
+            switch x
+                case 'n'
+                    rescale_f = input('specify resampling factor')
+                case 'y'
+                    rescale_f = 1;
+            end
+            
+            Geometry.voxel_size = meta.spacedirections * rescale_f;
             Geometry.start=meta.spaceorigin;
         case 2
             imgIn=am2mat([IMG_path, IMG_in]);
             Geometry.data = permute(imgIn.data, [2,1,3]);
+            % offset the CT, because that is what WiscPlan wants
+            Geometry.data = Geometry.data + 1000;
+            
             disp('NRRD file loaded!')
             Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
             Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
@@ -97,63 +134,66 @@ function [Geometry] = load2geometry
     mkdir(fullfile(Geometry.patient_dir, 'opt_output'));
 
     %% -----===== BTV and Ring creation =====-----
-%     waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
-%     ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
-% 
-%     [target_idx okay] = listdlg('ListString', ROI_names, ...
-%         'SelectionMode', 'single', 'Name', 'Target Selection', ...
-%         'PromptString', 'Please select the target ROI. ');
-%     if okay ~= 1
-%         msgbox('Plan creation aborted');
-%         delete(hWaitbar);
-%         return;
-%     end
-% 
-%     [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
-%         'BTV margin specification', 1, {'0.6'});
-%     if isempty(BTV_margin_answer)
-%         BTV_margin = 0.6;
-%     else
-%         BTV_margin = str2double(BTV_margin_answer{1});
-%     end
-% 
-%     % target_idx and BTV_margin are set. Expand PTV to BTV
-%     PTVmask = false(size(Geometry.rhomw));
-%     % for target_idx = target_indices
-%             PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
-%     % end
-%     if BTV_margin > 0
-%         if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
-%             bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-%             Geometry.BTV = bwD <= BTV_margin;
-%         end
-%     end
-% 
-%     % Create btv
-%     Geometry.ROIS{end+1} = Geometry.ROIS{end};
-%     Geometry.ROIS{end}.name = 'BTV';
-%     Geometry.ROIS{end}.num_curves = 0;
-%     Geometry.ROIS{end}.curves = {};
-%     Geometry.ROIS{end}.ind = find(Geometry.BTV);
-%     Geometry.ROIS{end}.visible = false;
-% 
-%     % Create ring
-%     [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
-%         'Ring margin specification', 1, {'1'});
-%     if isempty(ring_margin_answer)
-%         ring_margin = 1;
-%     else
-%         ring_margin = str2double(ring_margin_answer{1});
-%     end
-%     bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-%     Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
-%     Geometry.Ring = xor(Geometry.Ring, PTVmask);
-%     Geometry.ROIS{end+1} = Geometry.ROIS{end};
-%     Geometry.ROIS{end}.name = 'Ring';
-%     Geometry.ROIS{end}.num_curves = 0;
-%     Geometry.ROIS{end}.curves = {};
-%     Geometry.ROIS{end}.ind = find(Geometry.Ring);
-%     Geometry.ROIS{end}.visible = false;
+    waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
+    ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
+
+    [target_idx okay] = listdlg('ListString', ROI_names, ...
+        'SelectionMode', 'single', 'Name', 'Target Selection', ...
+        'PromptString', 'Please select the target ROI. ');
+    if okay ~= 1
+        msgbox('Plan creation aborted');
+        delete(hWaitbar);
+        return;
+    end
+
+    [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
+        'BTV margin specification', 1, {'0.6'});
+    if isempty(BTV_margin_answer)
+        BTV_margin = 0.6;
+    else
+        BTV_margin = str2double(BTV_margin_answer{1});
+    end
+
+    % target_idx and BTV_margin are set. Expand PTV to BTV
+    PTVmask = false(size(Geometry.rhomw));
+    % for target_idx = target_indices
+            PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
+    % end
+    if BTV_margin > 0
+        if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
+            bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+            Geometry.BTV = bwD <= BTV_margin;
+        elseif BTV_margin < min(Geometry.voxel_size)
+            warning('BTV margin too small!')
+            Geometry.BTV = zeros(size(Geometry.data));
+        end
+    end
+
+    % Create btv
+    Geometry.ROIS{end+1} = Geometry.ROIS{end};
+    Geometry.ROIS{end}.name = 'BTV';
+    Geometry.ROIS{end}.num_curves = 0;
+    Geometry.ROIS{end}.curves = {};
+    Geometry.ROIS{end}.ind = find(Geometry.BTV);
+    Geometry.ROIS{end}.visible = false;
+
+    % Create ring
+    [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
+        'Ring margin specification', 1, {'1'});
+    if isempty(ring_margin_answer)
+        ring_margin = 1;
+    else
+        ring_margin = str2double(ring_margin_answer{1});
+    end
+    bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+    Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
+    Geometry.Ring = xor(Geometry.Ring, PTVmask);
+    Geometry.ROIS{end+1} = Geometry.ROIS{end};
+    Geometry.ROIS{end}.name = 'Ring';
+    Geometry.ROIS{end}.num_curves = 0;
+    Geometry.ROIS{end}.curves = {};
+    Geometry.ROIS{end}.ind = find(Geometry.Ring);
+    Geometry.ROIS{end}.visible = false;
 
     
     %% -----===== Save the matlab files =====-----

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -9,6 +9,23 @@ for k = 1:num_batches
     B = [B batch_doses];
 end
 
+% apply the shift (0.5 vox up and left) to make dose match geometry
+f = waitbar(0, 'Correcting beamlet shift');
+for beam_i = 1:numel(B)
+    new_dose_list = B{beam_i};
+    waitbar(beam_i/numel(B),f)
+    
+    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{1, beam_i}.non_zero_indices = nonzero_idx;
+    B{1, beam_i}.non_zero_values = shifted_img(nonzero_idx);
+end
+close(f)
+
 for k = 1:numel(B)
     B{k}.num = k-1;
 end

WiscPlanPhotonkV125/matlab_frontend/nrrd2ROI.m

@@ -24,9 +24,12 @@ function Geometry=nrrd2ROI(Geometry, patient_dir)
                         error('You should never see this.')
                 end
                 
+                WiscPlan_preferences.patientDataPath = Geometry.patient_dir;
+                save(which('WiscPlan_preferences.mat'), 'WiscPlan_preferences');
+                
                 Geometry.ROIS{end+1}.ind = find(ROI_img>0);
-                str = input('Enter ROI name: \n','s');
-                Geometry.ROIS{end}.name= str;
+                str = inputdlg('Enter ROI name:')
+                Geometry.ROIS{end}.name= str{1};
                 save(fullfile(Geometry.patient_dir, 'matlab_files', 'Geometry.mat'), 'Geometry');
 
             case 'No'

data/TomoPhantPlan01/ctlist.txt

@@ -1,40 +0,0 @@
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.1.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.10.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.11.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.12.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.13.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.14.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.15.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.16.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.17.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.18.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.19.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.2.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.20.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.21.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.22.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.23.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.24.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.25.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.26.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.27.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.28.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.29.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.3.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.30.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.31.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.32.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.33.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.34.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.35.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.36.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.37.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.38.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.39.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.4.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.40.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.5.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.6.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.7.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.8.dcm
-CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.9.dcm

data/TomoPhantPlan01/ctlist.txt.sort.txt

@@ -1,40 +0,0 @@
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.40.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.39.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.38.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.37.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.36.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.35.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.34.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.33.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.32.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.31.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.30.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.29.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.28.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.27.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.26.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.25.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.24.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.23.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.22.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.21.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.20.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.19.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.18.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.17.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.16.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.15.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.14.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.13.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.12.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.11.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.10.dcm
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.9.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.8.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.7.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.6.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.5.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.4.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.3.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.2.dcm 
-C:\010-work\003_localGit\WiscPlan_v2\data\TomoPhantPlan01\CT.1.2.826.0.1.3680043.2.200.1428810724.6.92314.144.1.dcm 

optGoals/make_ROI_goals.m

@@ -1,443 +0,0 @@
-
-
-function make_ROI_goals (Geometry, beamlets, beamlets_joined, patient)
-% this is a support function for NLP_beamlet_optimizer. It creates the
-% goals planned for specific patients.
-    
-switch patient
-    case 'patient'
-        ROI_goals.optGoal = make_ROI_goals_Pat(Geometry, beamlets);
-        ROI_goals.optGoal_beam = make_ROI_goals_Pat(Geometry, beamlets_joined);
-        
-        ROI_goals.optGoal_idx=[1,2]; 
-        ROI_goals.targetMinMax_idx=[1,3];
-        
-    case 'tomoPhantom'   
-        ROI_goals.optGoal = make_ROI_goals_Pat2(Geometry, beamlets);
-        ROI_goals.optGoal_beam = make_ROI_goals_Pat2(Geometry, beamlets_joined);
-        
-        ROI_goals.optGoal_idx=[1,3];
-        ROI_goals.targetMinMax_idx=[1,2];
-        
-    case 'phantom_HD'
-        ROI_goals.optGoal = make_ROI_goals(Geometry, beamlets);
-        ROI_goals.optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
-        
-        ROI_goals.optGoal_idx=[1,3];
-        ROI_goals.targetMinMax_idx=[1,2];
-        
-    case 'doggo'
-%         optGoal = make_ROI_goals_DOG(Geometry, beamlets);
-%         optGoal_beam = make_ROI_goals_DOG(Geometry, beamlets_joined);
-
-        ROI_goals.optGoal = make_ROI_goals_DOG_2(Geometry, beamlets);
-        ROI_goals.optGoal_beam = make_ROI_goals_DOG_2(Geometry, beamlets_joined);
-        
-        ROI_goals.optGoal_idx=[1,2];
-        ROI_goals.targetMinMax_idx=[1,3];
-    case 'gbm_005'
-        ROI_goals.optGoal = make_ROI_goals_gbm_005(Geometry, beamlets);
-        ROI_goals.optGoal_beam = make_ROI_goals_gbm_005(Geometry, beamlets_joined);
-        
-        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
-        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
-    otherwise
-        error('invalid case')
-    end
-    
-    save(['C:\010-work\003_localGit\WiscPlan_v2\optGoals\' patient, '.mat'], 'ROI_goals')
-end
-
-% ---- MAKE ROI GOALS ----
-% min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
-function optGoal = make_ROI_goals_Pat(Geometry, beamlets)
-    optGoal={};
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'Target_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 60;
-    goal_1.function = 'min';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
-    goal_1.opt_weight = 10 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'Target_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 65;
-    goal_2.function = 'max';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 0.8 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'ROI 1_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 3}.name;
-    ROI_idx = Geometry.ROIS{1, 3}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 10;
-    goal_3.function = 'max';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_Pat2(Geometry, beamlets)
-    optGoal={};
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'Target_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 62;
-    goal_1.function = 'min';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
-    goal_1.opt_weight = 10 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'Target_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 63;
-    goal_2.function = 'max';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.05;
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'ROI 1_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 3}.name;
-    ROI_idx = Geometry.ROIS{1, 3}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 10;
-    goal_3.function = 'max';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Ring_max';
-    goal_4.ROI_name = Geometry.ROIS{1, 6}.name;
-    ROI_idx = Geometry.ROIS{1, 6}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 56;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_DOG(Geometry, beamlets)
-    optGoal={};
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'Target_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 62;
-    goal_1.function = 'min_sq';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
-    goal_1.opt_weight = 10 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    goal_2.function = 'max_sq';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_sq';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'Doggo_max2';
-    goal_3.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 50;
-    goal_3.function = 'max_sq';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_DOG_2(Geometry, beamlets)
-    optGoal={};
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'Target_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 62;
-    goal_1.function = 'min_sq';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
-    goal_1.opt_weight = 10 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    goal_2.function = 'max_sq';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_sq';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target_alpha = 1.05;
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Doggo_max2';
-    goal_4.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 50;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_DOG_3(Geometry, beamlets)
-    optGoal={};
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'Target_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 62;
-    goal_1.function = 'min_perc_Volume';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
-    goal_1.opt_weight = 10 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    goal_2.function = 'max_sq';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.targetImg   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_perc_Volume';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target_alpha = 1.05;
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 0.5 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Doggo_max2';
-    goal_4.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 50;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets)
-    optGoal={};
-    
-    pathIn = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\';
-    [minDose, minDose_meta] = nrrdread([pathIn 'RODP_files\FET_FGL005_B1_DP_minDose.nrrd']);
-    [maxDose, maxDose_meta] = nrrdread([pathIn 'RODP_files\FET_FGL005_B1_DP_maxDose.nrrd']);
-    minDose = double(minDose);
-    maxDose = double(maxDose);
-    
-    % -- START DEFINITION OF GOAL --
-    goal_1.name = 'CTV_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_1.ROI_idx = ROI_idx;
-    goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = minDose(ROI_idx);
-    goal_1.function = 'min';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-%     goal_1.target_alpha = 1.00;
-    goal_1.target   = 60 % minDose(ROI_idx);
-    goal_1.opt_weight = 40 / numel(ROI_idx); % normalize to volume of target area
-    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_1;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_2.name = 'CTV_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = maxDose(ROI_idx);
-    goal_2.function = 'max_sq';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-%     goal_1.target_alpha = 1.00;
-    goal_2.target   = 60 % maxDose(ROI_idx);
-    goal_2.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_2;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'head_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 0.1;
-    goal_3.function = 'max_sq';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 7 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-   
-end

planEvaluation/plot_DVH.m

@@ -1,6 +1,9 @@
 % ---- DVH PLOT FUNCTION ----
-function plot_DVH(dose, optGoal, optGoal_idx, targetMinMax_idx)
+function plot_DVH(dose, optGoal)
     % this function plots the DVHs of the given dose
+    
+    
+    
     nfrac=1;
     
     figure