Added loop run, added beamlet smoothing function, added modulation parametrization, improved GUI

NLP_run_loop.m

@@ -0,0 +1,44 @@
+
+
+% Pat_path = 'F:\021_WiscPlan_data\Uulke_prostate_01';
+% path2goal = 'F:\021_WiscPlan_data\Uulke_prostate_01\matlab_files\OG2_4_v13_m3_0e4.mat';
+% [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';
+[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';
+[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';
+[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';
+[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';
+[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';
+[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_getFullDose(NLP_path, NLP_file)
+NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift1.mat';
+NLP_getFullDose(NLP_path, NLP_file)
+NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift2.mat';
+NLP_getFullDose(NLP_path, NLP_file)
+NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift3.mat';
+NLP_getFullDose(NLP_path, NLP_file)
+NLP_file = 'NLP_result_OG2_4_v15_m3_1e4_shift4.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_getFullDose(NLP_path, NLP_file)
+
+
+
+
+
+

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v3.m

@@ -39,11 +39,21 @@ else
     [Goal_path,Goal_file,ext] = fileparts(path2goal);
 end
 
-str = inputdlg({'N of iterations for initial calc', 'N of iterations for full calc', ...
-    'Use pre-existing NLP_result to initiate? (y/n)'}, 'input', [1,35], {'100000', '500000', 'n'});
-N_fcallback1 = str2double(str{1}); % 100000  is a good guesstimate
-N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
-pre_beamWeights = str{3};
+dialogue_box = 'no'
+switch dialogue_box
+    case 'yes'
+        str = inputdlg({'N of iterations for initial calc', 'N of iterations for full calc', ...
+            'Use pre-existing NLP_result to initiate? (y/n)'}, 'input', [1,35], {'100000', '500000', 'n'});
+        N_fcallback1 = str2double(str{1}); % 100000  is a good guesstimate
+        N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
+        pre_beamWeights = str{3};
+    case 'no'
+        disp('dialogue box skipped')
+        N_fcallback1 = 100000;
+        N_fcallback2 = 500000;
+        pre_beamWeights = 'n';
+end
+
 
 switch pre_beamWeights
     case 'y'
@@ -73,16 +83,23 @@ load(path2goal)
 %% -- OPTIMIZATION TARGETS --
 % -- make the optimization optGoal structure --
 for i_goal = 1:size(OptGoals.goals,1)
+        
     if isfield(OptGoals.data{i_goal}, 'SupVox_num')
         SupVox_num = OptGoals.data{i_goal}.SupVox_num;
     else
         answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
         SupVox_num = str2double(answer{1})
     end
+    optGoal{i_goal} = OptGoals.data{i_goal};
+    optGoal{i_goal}.sss_scene_list = OptGoals.sss_scene_list;
+    optGoal{i_goal}.maxModulation = OptGoals.maxModulation;
+    optGoal{i_goal}.BeamSmoothMax = OptGoals.BeamSmoothMax;
+    
+    % modulation
     switch SupVox_num
         case 0
         % if not supervoxel, just select provided ROI_idx
-        optGoal{i_goal} = OptGoals.data{i_goal};
+        
         optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
     
         otherwise
@@ -90,11 +107,12 @@ for i_goal = 1:size(OptGoals.goals,1)
         % - make supervoxel map
         mask = zeros(OptGoals.data{i_goal}.imgDim);
         mask(OptGoals.data{i_goal}.ROI_idx) = 1;
-        superMask = superpix_group(mask, SupVox_num); 
+        % group superpixels
+        superMask = superpix_group(mask, SupVox_num, 'no'); 
+        
         superVoxList = unique(superMask);
         superVoxList = superVoxList(superVoxList>0);
         
-        optGoal{i_goal} = OptGoals.data{i_goal};
         optGoal{i_goal}.ROI_idx_old = optGoal{i_goal}.ROI_idx; % copy old index data
         optGoal{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
         optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
@@ -203,6 +221,7 @@ D_full = reshape(beamlets * w_fin, size(Geometry.data));
 %% save outputs
 NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
+NLP_result.sss_scene_list = optGoal{1}.sss_scene_list;
 save([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
 
 
@@ -222,6 +241,8 @@ function penalty = get_penalty(x, optGoal)
     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
@@ -311,12 +332,54 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
     
     %% add modulation penalty
     if true
-        max_modulation = 5;
+        if isfield(optGoal{goal_i}, 'maxModulation')
+            max_modulation = optGoal{goal_i}.maxModulation;
+        else
+            error('no Max modulation parameter enterd!')
+            max_modulation = 5;
+        end
+        
         mod_pen_weight = 1.0e10;
         % calc the penalty
-        mod_excess = max(0, x-5*mean(x));
-        mod_pen = mod_pen_weight*(sum(mod_excess) + numel(mod_excess)* any(mod_excess));
-        penalty = penalty + mod_pen;
+        mod_excess = max(0, x-max_modulation*mean(x));
+        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
+    if false
+        if isfield(optGoal{goal_i}, 'BeamSmoothMax')
+            BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
+        else
+            error('no Max beam smooth parameter enterd!')
+            BeamSmoothMax = 1e6;
+        end
+        mod_pen_weight = BeamSmoothMax; %1.0e6
+        x_down = zeros(size(x));
+        x_down(2:end) = x(1:end-1);
+        x_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));
+        penalty = penalty + mod_pen2;
+    end
+    %% add penlty for single off beamlets - version 2
+    if true
+        if isfield(optGoal{goal_i}, 'BeamSmoothMax')
+            BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
+        else
+            error('no Max beam smooth parameter enterd!')
+            BeamSmoothMax = 1e6;
+        end
+        mod_pen_weight = BeamSmoothMax; %1.0e6
+        x_down = zeros(size(x));
+        x_down(2:end) = x(1:end-1);
+        x_down2 = zeros(size(x));
+        x_down2(3:end) = x(1:end-2);
+        x_up = zeros(size(x));
+        x_up(1:end-1) = x(2:end);
+        x_up2 = zeros(size(x));
+        x_up2(1:end-2) = x(3:end);
+        mod_pen2 = mod_pen_weight*sum((x_down+x_up-2*x).^2 + 0.5*(x_down2+x_up2-2*x).^2)/(mean(x)^(2)*numel(x));
+        penalty = penalty + mod_pen2;
     end
 end
 
@@ -340,11 +403,13 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 
     
     % ----====#### CHANGE ROBUSTNESS HERE ####====----
-    
+    if isfield(optGoal{1}, 'sss_scene_list')
+        sss_scene_list = optGoal{1}.sss_scene_list;
+    else
+        sss_scene_list={[0,0,0], [-shift_Y,0,0], [shift_Y,0,0], [0,-shift_X,0], [0,shift_X,0], [0,0,-shift_Z], [0,0,shift_Z]};
+        optGoal{1}.sss_scene_list = sss_scene_list;
+    end
 %     sss_scene_list={[0,0,0]};
-    sss_scene_list={[0,0,0], [-shift_Y,0,0], [shift_Y,0,0], [0,-shift_X,0], [0,shift_X,0], [0,0,-shift_Z], [0,0,shift_Z]};
-%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
-%         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
 
     % ----====#### CHANGE ROBUSTNESS HERE ####====----
 %     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -22,7 +22,7 @@ function varargout = goal_def_UI(varargin)
 
 % 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
+% Last Modified by GUIDE v2.5 04-Jun-2020 04:51:32
 
 % Begin initialization code - DO NOT EDIT
 gui_Singleton = 1;
@@ -238,6 +238,25 @@ end
 disp('-- all OptGoals exported!')
 
 OptGoals.goals = [handles.uitable1.Data];
+OptGoals.maxModulation = str2double(handles.maxModulation.String);
+OptGoals.BeamSmoothMax = str2double(handles.BeamSmoothMax.String);
+
+sss_scene_list={[0,0,0]};
+Num = sscanf(handles.sss_Y.String, '%g,').';
+for i=1:numel(Num)
+    sss_scene_list{end+1}=[Num(i),0,0];
+end
+Num = sscanf(handles.sss_X.String, '%g,').';
+for i=1:numel(Num)
+    sss_scene_list{end+1}=[0,Num(i),0];
+end
+Num = sscanf(handles.sss_Z.String, '%g,').';
+for i=1:numel(Num)
+    sss_scene_list{end+1}=[0,0,Num(i)];
+end
+
+OptGoals.sss_scene_list = sss_scene_list;
+
 
 save( [path, file], 'OptGoals')
 end
@@ -253,7 +272,10 @@ load( [path, file])
 disp('load goal')
 
 handles.uitable1.Data = OptGoals.goals;
-
+handles.maxModulation.String = num2str(OptGoals.maxModulation);
+if isfield(OptGoals, 'BeamSmoothMax')
+    handles.MaxBeamSmooth.String = num2str(OptGoals.BeamSmoothMax);
+end
 end
 
 
@@ -320,3 +342,107 @@ end
 
 
 end
+
+
+
+function maxModulation_Callback(hObject, eventdata, handles)
+% hObject    handle to maxModulation (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of maxModulation as text
+%        str2double(get(hObject,'String')) returns contents of maxModulation as a double
+end
+
+% --- Executes during object creation, after setting all properties.
+function maxModulation_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to maxModulation (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+end
+
+
+function sss_Y_Callback(hObject, eventdata, handles)
+% hObject    handle to sss_Y (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of sss_Y as text
+%        str2double(get(hObject,'String')) returns contents of sss_Y as a double
+end
+
+% --- Executes during object creation, after setting all properties.
+function sss_Y_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to sss_Y (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+end
+
+
+function sss_X_Callback(hObject, eventdata, handles)
+% hObject    handle to sss_X (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of sss_X as text
+%        str2double(get(hObject,'String')) returns contents of sss_X as a double
+end
+
+% --- Executes during object creation, after setting all properties.
+function sss_X_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to sss_X (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+end
+
+
+function sss_Z_Callback(hObject, eventdata, handles)
+% hObject    handle to sss_Z (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of sss_Z as text
+%        str2double(get(hObject,'String')) returns contents of sss_Z as a double
+end
+
+% --- Executes during object creation, after setting all properties.
+function sss_Z_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to sss_Z (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+end
+
+
+
+function BeamSmoothMax_Callback(hObject, eventdata, handles)
+% hObject    handle to BeamSmoothMax (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of BeamSmoothMax as text
+%        str2double(get(hObject,'String')) returns contents of BeamSmoothMax as a double
+end

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -1,6 +1,6 @@
 
 
-function superMask = superpix_group(mask, N_svox_in)
+function superMask = superpix_group(mask, N_svox_in, orthoPlot)
 % this function contains supervoxel grouping
 fprintf(['\n' 'Creating ' num2str(N_svox_in) ' SupVox:\n' ])
 canvas2 = sqrt(bwdist(1-mask));
@@ -55,13 +55,19 @@ while true
     numSupVox = numel(unique(superMask))-1; % number of created supervoxels
     fprintf([num2str(numSupVox) ' areas created' ])
     
-    if abs(numSupVox-N_svox_in)/N_svox_in < 0.2/converge_factor
-        orthoslice(superMask)
+    if abs(numSupVox-N_svox_in)/N_svox_in < 0.2/sqrt(converge_factor)
+        
+        switch orthoPlot
+            case 'yes'
+                orthoslice(superMask)
+            case 'no'
+                fprintf('\n supervoxel volume created! (orthoslice skipped)')
+        end
         break
     end
     N_svox2 = N_svox * N_svox_in/numSupVox;
     N_svox = round(converge_factor * N_svox2 + (1-converge_factor)*N_svox);
-    converge_factor = converge_factor * 0.97;
+    converge_factor = converge_factor * 0.95;
     fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target. New start size: ' num2str(N_svox) '\n'])
     
     if converge_factor < 0.01

view_beamlet.m

@@ -0,0 +1,15 @@
+
+
+% data = read_ryan_beamlets('F:\021_WiscPlan_data\Hippocampus_HD_2\batch_dose.bin')
+
+
+idx = 3000;
+
+tabula = zeros([data{idx}.x_count, data{idx}.y_count, data{idx}.z_count]);
+
+tabula(data{idx}.non_zero_indices)   = tabula(data{idx}.non_zero_indices)+data{idx}.non_zero_values;
+tabula(data{idx+1}.non_zero_indices) = tabula(data{idx+1}.non_zero_indices)+data{idx+1}.non_zero_values;
+
+tabula(data{idx+2}.non_zero_indices) = tabula(data{idx+2}.non_zero_indices)+data{idx+2}.non_zero_values;
+
+orthoslice(tabula);