added dose painting functionality

NLP_beamlet_optimizer.m

@@ -1,5 +1,5 @@
 
-function [D_full, w_fin] = NLP_beamlet_optimizer
+function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 % 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
@@ -11,24 +11,26 @@ function [D_full, w_fin] = NLP_beamlet_optimizer
 % To-do:
 % - Add robusness aspect (+take worst case scenario, see REGGUI)
 
-N_fcallback1 = 10000;
-N_fcallback2 = 50000;
-patient = 'phantom';
+N_fcallback1 = 1000;
+N_fcallback2 = 150000;
+patient = 'doggo';
 
 switch patient
-    case 'phantom'
+    case 'patient'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
         blet_in_beam=5; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     case 'phantom_HD'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
-        blet_in_beam=20; % ???
+        blet_in_beam=20; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     case 'doggo'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_2';
-        blet_in_beam=5; % ???
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_3';
+        blet_in_beam=5; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     otherwise
         error('invalid case')
 end
 
+merge_beamlets(4, patient_dir);
+
 
 %% PROGRAM STARTS HERE
 % - no tocar lo que hay debajo -
@@ -76,7 +78,7 @@ close(wbar2)
 
 %% -- OPTIMIZATION TARGETS --
 switch patient
-    case 'phantom'
+    case 'patient'
         optGoal = make_ROI_goals(Geometry, beamlets);
         optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
         N_beamlets_in_beam = 10;
@@ -84,8 +86,11 @@ switch patient
         optGoal = make_ROI_goals(Geometry, beamlets);
         optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
     case 'doggo'
-        optGoal = make_ROI_goals_DOG(Geometry, beamlets);
-        optGoal_beam = make_ROI_goals_DOG(Geometry, beamlets_joined);
+%         optGoal = make_ROI_goals_DOG(Geometry, beamlets);
+%         optGoal_beam = make_ROI_goals_DOG(Geometry, beamlets_joined);
+
+        optGoal = make_ROI_goals_DOG_2(Geometry, beamlets);
+        optGoal_beam = make_ROI_goals_DOG_2(Geometry, beamlets_joined);
     otherwise
         error('invalid case')
 end
@@ -120,14 +125,16 @@ options = optimoptions('fmincon');
 options.MaxFunctionEvaluations = N_fcallback1;
 options.Display = 'iter';
 options.PlotFcn = 'optimplotfval';
+% options.UseParallel = true;
 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);
-t=toc;
-disp(['Optimization time for beams = ',num2str(t)]);
+% 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
@@ -138,7 +145,7 @@ w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % smal
 
 % -- GET FULL BEAMLET WEIGHTS --
 options.MaxFunctionEvaluations = N_fcallback2;
-tic
+% tic
 w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
 t=toc;
 disp(['Optimization time for beamlets = ',num2str(t)]);
@@ -152,12 +159,13 @@ NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
 % save('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat', 'NLP_result');
 
-optGoal_idx=[1,3];
+optGoal_idx=[1,2];
 plot_DVH(D_full, optGoal, optGoal_idx)
 % plot_DVH_robust(D_full, optGoal, optGoal_idx)
 end
 
 % orthoslice(D_full, [0,70])
+% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
 
 %% support functions
 % ---- PENALTY FUNCTION ----
@@ -170,7 +178,7 @@ function penalty = get_penalty(x, optGoal)
     sc_i = 1;
     
     for nrs_i = 1:optGoal{1}.NbrRandScenarios 
-        for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = ss
+        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);
                 sc_i = sc_i + 1;
@@ -196,12 +204,29 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
                 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)));
+            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));
+                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));
+                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);
+            case 'min_Volume'
+                % penalize by amount of volume under threshold
+                n_vox = 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);
+                
+%                 d_penalty = 1.0e0 * 
+                    
         end
         penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
     end
@@ -241,7 +266,6 @@ function beamlets = get_beamlets(beamlet_cell_array, numVox);
     close(wbar1)
 
 end
-
 function show_joint_beamlets(beamlets, IMGsize, Beam_list)
     % this function overlays and plots multiple beamlets. The goal is to
     % check whether some beamlets are part of the same beam manually.
@@ -268,7 +292,7 @@ function optGoal = make_ROI_goals(Geometry, beamlets)
     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 = 1 / numel(ROI_idx); % normalize to volume of target area
+    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;
@@ -316,11 +340,11 @@ function optGoal = make_ROI_goals_DOG(Geometry, beamlets)
     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';
+    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 = 1 / numel(ROI_idx); % normalize to volume of target area
+    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;
@@ -329,20 +353,123 @@ function optGoal = make_ROI_goals_DOG(Geometry, beamlets)
     % -- 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';
+    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.9, 0.2, 0.2]; % color of the final DVH plot
+    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
+
 % ---- MAKE ROI ROBUST ----
 function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % take regular optimal goal and translate it into several robust cases
@@ -360,9 +487,9 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     rrs_scene_list={[0,0,0]};
     
 %     nrs_scene_list={[0,0,0]};
-    sss_scene_list={[0,0,0],  [-15,0,0], [-10,0,0], [-5,0,0], [5,0,0], [10,0,0], [15,0,0]};
+%     sss_scene_list={[0,0,0], [-2,0,0], [2,0,0], [0,-2,0], [0,2,0]};
 %     rrs_scene_list={[0,0,0]};
-    
+    [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_DP_target.nrrd');
     
     for i = 1:numel(optGoal)
         optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
@@ -383,7 +510,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
             targetImg2=targetImg1;
             % beamlets stay the same
             
-            for sss_i = 1:optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
+            for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
                 % modify target and beamlets
                 targetImg3=get_RO_sss(targetImg2, sss_scene_list{sss_i});
                 % beamlets stay the same
@@ -396,7 +523,16 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
                     %% make new target and beamlets
                     ROI_idx=[];
                     ROI_idx=find(targetImg4>0);
-                    target = ones(numel(ROI_idx), 1) * optGoal{goal_i}.D_final;
+                    
+                    if isfield(optGoal{goal_i}, 'target_alpha')
+                        target = double(optGoal{goal_i}.target_alpha * targetIn(ROI_idx));
+                        disp('exists')
+                    else
+                        target = ones(numel(ROI_idx), 1) * optGoal{goal_i}.D_final;
+                        disp('not exists')
+                    end
+                    
+                    
                     beamlets_pruned = beamlets(ROI_idx, :);
                     
                     % save to optGoal output

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -49,7 +49,7 @@ ypmax = 1.0;    % +jaw width / 2
 % y-prime points in the z-direction in the CT coordinate system
 
 % Number of beamlets in the BEV for each direction
-Mxp = 5;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
+Mxp = 20;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
 Nyp = 1;  % always 1 for Tomo due to binary mlc
 
 
@@ -93,7 +93,7 @@ Nrot = ceil(abs(zBow - zStern)/(pitch*fieldWidth));
 
 
 % Nphi = Nrot*51;  % number of angles used in the calculation
-Nphi = Nrot*21;  % Grozomah
+Nphi = Nrot*51;  % Grozomah
 
 % load(geometry_file);
 

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -1,6 +1,6 @@
 function B = merge_beamlets(num_batches, pat_dir)
 
-num_batches = 4;
+% num_batches = 4;
 
 
 B = [];
@@ -13,3 +13,5 @@ for k = 1:numel(B)
     B{k}.num = k-1;
 end
 write_ryan_beamlets([pat_dir '\batch_dose.bin'],B);
+
+write_ryan_beamlets([pat_dir '\batch_dose_backup.bin'],B);