Updates in beamlet smoothing, preparation to mobe it to Guy

NLP_run_loop.m

@@ -5,36 +5,33 @@
 % [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
 
 Pat_path = 'F:\021_WiscPlan_data\Uulke_prostate_01';
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift0.mat';
-% [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift1.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss0.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift2.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss1.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift3.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss2.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift4.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss3.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift5.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss4.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
-path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v15_m3_1e4_shift6.mat';
+path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_boost_v22_sss5.mat';
 [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
 
 
+
 NLP_path = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\';
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift0.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss0.mat';
 NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift1.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss1.mat';
 NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift2.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss2.mat';
 NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift3.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss3.mat';
 NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift4.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss4.mat';
 NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift5.mat';
-NLP_getFullDose(NLP_path, NLP_file)
-NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift6.mat';
+NLP_file = 'NLP_result_OG2_4_boost_v22_sss5.mat';
 NLP_getFullDose(NLP_path, NLP_file)
 
 
@@ -42,3 +39,4 @@ NLP_getFullDose(NLP_path, NLP_file)
 
 
 
+

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v3.m

@@ -142,17 +142,27 @@ for i_goal = 1:size(OptGoals.goals,1)
     end
 end
 
-
-
 % -- make them robust --
 RO_params=0;
 optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
-
-save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
+% save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
+
+% -- get beamlet indeces --
+load([Pat_path, '\all_beams.mat'])
+Nbeamlets = all_beams{1}.Mxp;   % number of beamlets in a beam - usually 64 or 32
+weightTable = zeros(100,Nbeamlets);
+for ind_bmlt = 1:numel(all_beams)
+    bLet_idx.y(ind_bmlt) = floor(all_beams{1, ind_bmlt}.num/Nbeamlets)+1;
+    bLet_idx.x(ind_bmlt) = rem(all_beams{1, ind_bmlt}.num, Nbeamlets)+1;
+    weightTable(bLet_idx.y(ind_bmlt),bLet_idx.x(ind_bmlt)) = 1;
+end
+bLet_idx.idx = find(weightTable>0);
+bLet_idx.Nbeamlets = Nbeamlets;
+disp('.')
 
 % -- CALLBACK OPTIMIZATION FUNCTION --
-fun1 = @(x) get_penalty(x, optGoal_beam);
-fun2 = @(x) get_penalty(x, optGoal);
+fun1 = @(x) get_penalty(x, optGoal_beam, bLet_idx);
+fun2 = @(x) get_penalty(x, optGoal,     bLet_idx);
 
 % -- OPTIMIZATION PARAMETERS --
 % define optimization parameters
@@ -222,6 +232,8 @@ D_full = reshape(beamlets * w_fin, size(Geometry.data));
 NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
 NLP_result.sss_scene_list = optGoal{1}.sss_scene_list;
+NLP_result.maxModulation = OptGoals.maxModulation;
+NLP_result.BeamSmoothMax = OptGoals.BeamSmoothMax;
 save([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
 
 
@@ -233,7 +245,7 @@ end
 
 %% support functions
 % ---- PENALTY FUNCTION ----
-function penalty = get_penalty(x, optGoal)
+function penalty = get_penalty(x, optGoal, bLet_idx)
     % this function gets called by the optimizer. It checks the penalty for
     % all the robust implementation and returns the worst result.
     
@@ -246,7 +258,7 @@ 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 rgs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
-                fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i);
+                fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx);
                 sc_i = sc_i + 1;
             end
         end
@@ -261,7 +273,7 @@ function penalty = get_penalty(x, optGoal)
     
 end
 % ------ supp: penalty for single scenario ------
-function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
+function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx)
     penalty = 0;
     % for each condition
     for goal_i = 1:numel(optGoal)
@@ -342,10 +354,12 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
         mod_pen_weight = 1.0e10;
         % calc the penalty
         mod_excess = max(0, x-max_modulation*mean(x));
+%         mod_excess = max(0, x-max_modulation*mean(x(x>1)));
         mod_pen1 = mod_pen_weight*(sum(mod_excess) + numel(mod_excess)* any(mod_excess));
         penalty = penalty + mod_pen1;
     end
-    %% add penlty for single off beamlets - version 1
+    
+    %% add penlty for single off beamlets - version 2
     if false
         if isfield(optGoal{goal_i}, 'BeamSmoothMax')
             BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
@@ -356,12 +370,16 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
         mod_pen_weight = BeamSmoothMax; %1.0e6
         x_down = zeros(size(x));
         x_down(2:end) = x(1:end-1);
+        x_down2 = zeros(size(x));
+        x_down2(3:end) = x(1:end-2);
         x_up = zeros(size(x));
         x_up(1:end-1) = x(2:end);
-        mod_pen2 = mod_pen_weight*sum(((x-x_down).*(x-x_up)).^2)/(mean(x)^(2*2)*numel(x));
+        x_up2 = zeros(size(x));
+        x_up2(1:end-2) = x(3:end);
+        mod_pen2 = mod_pen_weight*sum((x_down+x_up-2*x).^2 + 0.5*(x_down2+x_up2-2*x).^2)/(mean(x)^(2)*numel(x));
         penalty = penalty + mod_pen2;
     end
-    %% add penlty for single off beamlets - version 2
+    %% add penlty for single off beamlets - version 3
     if true
         if isfield(optGoal{goal_i}, 'BeamSmoothMax')
             BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
@@ -370,15 +388,30 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
             BeamSmoothMax = 1e6;
         end
         mod_pen_weight = BeamSmoothMax; %1.0e6
-        x_down = zeros(size(x));
-        x_down(2:end) = x(1:end-1);
-        x_down2 = zeros(size(x));
-        x_down2(3:end) = x(1:end-2);
-        x_up = zeros(size(x));
-        x_up(1:end-1) = x(2:end);
-        x_up2 = zeros(size(x));
-        x_up2(1:end-2) = x(3:end);
-        mod_pen2 = mod_pen_weight*sum((x_down+x_up-2*x).^2 + 0.5*(x_down2+x_up2-2*x).^2)/(mean(x)^(2)*numel(x));
+        
+        % make 2D beamletWeight map
+        maxY = max(bLet_idx.y);
+        tabula = zeros(maxY, bLet_idx.Nbeamlets);
+        tabula(bLet_idx.idx) = x;
+        % for each index calculate laplace
+        myWeights = 1/12*[1,3,1; 1,-12,1; 1,3,1];
+        i=2:maxY-1;
+        j=2:bLet_idx.Nbeamlets-1;
+        tabula2(i,j) = ...
+             myWeights(1)*tabula(i-1,j-1) + myWeights(4)*tabula(i-1,j) + myWeights(7)*tabula(i-1,j+1)...
+            +myWeights(2)*tabula(i,j-1)  + myWeights(8)*tabula(i,j+1) ...
+            +myWeights(3)*tabula(i+1,j-1) + myWeights(6)*tabula(i+1,j) + myWeights(9)*tabula(i+1,j+1)...
+            +myWeights(5)*tabula(i,j);
+        tabula2(1,j) = ...
+             myWeights(2)*tabula(1,j-1) + myWeights(8)*tabula(1,j+1)...
+            +myWeights(3)*tabula(2,j-1) + myWeights(6)*tabula(2,j) + myWeights(9)*tabula(2,j+1)...
+            +myWeights(5)*tabula(1,j);
+        tabula2(maxY,j) =...
+             myWeights(1)*tabula(maxY-1,j-1) + myWeights(4)*tabula(maxY-1,j) + myWeights(7)*tabula(maxY-1,j+1)...
+            +myWeights(2)*tabula(maxY,j-1)  + myWeights(8)*tabula(maxY,j+1) ...
+            +myWeights(5)*tabula(maxY,j);
+        % make sum of squares
+        mod_pen2 = mod_pen_weight*sum((tabula2(:)).^2)/(mean(x.^2)*numel(x));
         penalty = penalty + mod_pen2;
     end
 end

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -190,7 +190,7 @@ for i = 1 : size(handles.uitable1.Data, 1)
     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));
+        goal_i.D_final = str2double(handles.uitable1.Data{i,4}) * ones(size(goal_i.ROI_idx));
     elseif strcmp (handles.uitable1.Data{i,3}, 'Dose map')
 %         warning('Works only for NRRD (for now)')
         dose_in = nrrdread(handles.uitable1.Data{i,4});

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -489,6 +489,10 @@ for n=1:numel(batches)
     end
 end
 
+
+all_beams{1}.Mxp = Mxp;
+all_beams{1}.N_angles = N_angles;
+all_beams{1}.num_batches = num_batches;
 save([patient_dir '\all_beams.mat'], 'all_beams');
 % for k = 1:numel(batches)
 %         system([fullfile(patient_dir,sprintf('run%d.cmd',k-1)) ' &']);

WiscPlanPhotonkV125/matlab_frontend/lab/RDXTPS_geometry_setup_wizard_ASP.m

@@ -121,6 +121,8 @@ waitbar(2/total_num_steps, hWaitbar, 'Step 3: Save geometry files');
 load('WiscPlan_preferences.mat')
 patient_dir = uigetdir(WiscPlan_preferences.patientDataPath, 'Save the patient data to directory');
 
+Geometry.data_dir = patient_dir;
+
 WiscPlan_preferences.patientDataPath = patient_dir;
 thisDir = mfilename('fullpath');
 idcs   = strfind(thisDir,'\');