Implemented supervoxel grouping and further cleaned up workflow (GUI expansion)

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v2.m

@@ -51,9 +51,9 @@ switch pre_beamWeights
         w_beamlets = NLP_result.weights;
         
         load([Pat_path, '\all_beams.mat'])
-        if numel(all_beams) ~= numel(w_beamlets)
-            error('Provided weight number does not match beamlet number!')
-        end
+%         if numel(all_beams) ~= numel(w_beamlets)
+%             error('Provided weight number does not match beamlet number!')
+%         end
     case 'n'
         disp('Initial beam weights will be calculated.')
 end
@@ -73,10 +73,17 @@ load(path2goal)
 % -- make the optimization optGoal structure --
 
 for i_goal = 1:size(OptGoals.goals,1)
-    answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
     
-    switch answer{1}
-        case '0'
+    if isfield(OptGoals.data{i_goal}, 'SupVox_num')
+        SupVox_num = OptGoals.data{i_goal}.SupVox_num;
+    else
+        answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
+        SupVox_num = str2double(answer{1})
+    end
+    
+    
+    switch SupVox_num
+        case 0
         % if not supervoxel, just select provided ROI_idx
         optGoal{i_goal} = OptGoals.data{i_goal};
         optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
@@ -86,11 +93,10 @@ for i_goal = 1:size(OptGoals.goals,1)
         otherwise
         % -- if supervoxel, merge given columns
         % - make supervoxel map
+        
         mask = zeros(OptGoals.data{i_goal}.imgDim);
         mask(OptGoals.data{i_goal}.ROI_idx) = 1;
-        
-        superMask = superpix_group(mask, str2double(answer{1})); 
-        
+        superMask = superpix_group(mask, SupVox_num); 
         superVoxList = unique(superMask);
         superVoxList = superVoxList(superVoxList>0);
         
@@ -205,7 +211,7 @@ save([Pat_path, '\matlab_files\NLP_result.mat'], 'NLP_result');
 
 
 plot_DVH(Geometry, D_full)
-colorwash(Geometry.data, D_full, [-500, 500], [0, 35]);
+colorwash(Geometry.data, D_full, [500, 1500], [0, 66]);
 
 
 end
@@ -283,7 +289,6 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
     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
@@ -297,7 +302,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % Y - Y>0 moves image down
     % Z - in/out.
     
-    shift_mag = 1; % vox of shift
+    shift_mag = 3; % vox of shift
     nrs_scene_list={[0,0,0]};
 
     

WiscPlanPhotonkV125/matlab_frontend/dicomrt/dicomrt2geometry.m

@@ -25,7 +25,11 @@ if nargin < 1
 % !!! Grozomah
     load('WiscPlan_preferences.mat')
     
-    if isstring(WiscPlan_preferences.inDataPath) 
+    if isfield(WiscPlan_preferences, 'inDataPath')
+        if isstring(WiscPlan_preferences.inDataPath) 
+        else
+            WiscPlan_preferences.inDataPath = 'C://';
+        end
     else
         WiscPlan_preferences.inDataPath = 'C://';
     end
@@ -146,7 +150,7 @@ if downsample_flag == true
     Geometry.data = resamplegrid(Geometry.data, [downsample_factor downsample_factor 1]);
     % Geometry.data = resamplegrid(Geometry.data, 0.5);
     % note that ct_xmesh, ct_ymesh, ct_zmesh is already downsampled
-    Geometry.start(1:2) = Geometry.start(1:2) + Geometry.voxel_size(1:2) / 4;
+    Geometry.start(1:2) = Geometry.start(1:2) + Geometry.voxel_size(1:2) / 2;
 end
 
 Geometry.rhomw = CT2dens(Geometry.data, 'pinn');

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -59,8 +59,9 @@ handles.output = hObject;
 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])
+    [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]);
+    eval('handles.Data.Geometry = Geometry;')
 end
 
 
@@ -87,6 +88,8 @@ handles.uitable1.Data{4} = '60';
 handles.uitable1.Data{5} = 'min_sq';
 handles.uitable1.Data{6} = 100;
 handles.uitable1.Data{7} = 'null';
+handles.uitable1.Data{8} = 0;
+handles.uitable1.Data{9} = 0;
 
 % Update handles structure
 guidata(hObject, handles);
@@ -196,11 +199,12 @@ for i = 1 : size(handles.uitable1.Data, 1)
     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, :);
+    % == optional
     goal_i.optionalParam = handles.uitable1.Data{i,7};
+    % == Number of supervoxels
+    goal_i.SupVox_num = handles.uitable1.Data{i,9};
+    
     % -- 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;
 
@@ -238,6 +242,7 @@ function uitable1_CellEditCallback(hObject, eventdata, handles)
 %	Error: error string when failed to convert EditData to appropriate value for Data
 % handles    structure with handles and user data (see GUIDATA)
 
+% --- update DP or fixed dose goals
 if eventdata.Indices(2) == 3;
     mode = handles.uitable1.Data{eventdata.Indices(1), eventdata.Indices(2)}
     switch mode
@@ -255,6 +260,21 @@ if eventdata.Indices(2) == 3;
     end
 end
 
+% --- update with number of voxels
+if eventdata.Indices(2) == 2;
+    for i = 1:size(handles.Data.Geometry.ROIS,2)
+        if strcmp(handles.Data.Geometry.ROIS{i}.name, handles.uitable1.Data{eventdata.Indices(1),2})
+            ROI_idx_num = numel(handles.Data.Geometry.ROIS{i}.ind);
+            break
+        end
+        if i == size(handles.Data.Geometry.ROIS,2)
+            error(['Invalid ROI name selected in goal ' num2str(i)])
+        end
+    end
+    
+    handles.uitable1.Data{eventdata.Indices(1), 8} = ROI_idx_num;
+    handles.uitable1.Data{eventdata.Indices(1), 9} = 0;
+end
 
 
 end

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -38,23 +38,25 @@ 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 = -20.0;     % -field width / 2
-xpmax = 20.0;      % +field width / 2
-ypmin = -0.3125;  % total jaw width is 0.625 cm
-ypmax = 0.3125;
-% ypmin = -0.5;    % total jaw width is 1 cm
-% ypmax = 0.5;
-
+% beam is 40 cm wide in transverse direction and 1-5 cm (usually 2) in y
+% direction.
+% isocenter is 85 cm from source, ends of jaws are 23 cm from source
+xpmin = -20.0;      % -field width / 2
+xpmax = 20.0;       % +field width / 2
+% ypmin = -0.3125;  % total jaw width is 0.625 cm
+% ypmax = 0.3125;
+ypmin = -0.5;       % total jaw width is 1 cm
+ypmax = 0.5;
+% ypmin = -1.25;    % total jaw width is 2.5 cm
+% ypmax = 1.25;
 % y-prime points in the z-direction in the CT coordinate system
 
 % ============================================= End of user-supplied inputs
 executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';
 kernel_file = 'Kernels.mat';
 geometry_file = fullfile(patient_dir, 'matlab_files\Geometry.mat');
-load(geometry_file);
-Geometry.num_batches = num_batches;
-save(geometry_file, 'Geometry');
 
+load(geometry_file);
 ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
 [target_idx, okay] = listdlg('ListString', ROI_names, ...
     'SelectionMode', 'single', 'Name', 'Target Selection', ...

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -13,7 +13,7 @@ end
 f = waitbar(0, 'Correcting beamlet shift');
 for beam_i = 1:numel(B)
     new_dose_list = B{beam_i};
-    waitbar(beam_i/numel(B),f)
+    waitbar(beam_i/numel(B),f, ['Correcting beamlet shift: ' num2str(beam_i)])
     
     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;

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -1,8 +1,8 @@
 
 
-function superMask = superpix_group(mask, N_svox)
+function superMask = superpix_group(mask, N_svox_in)
 % this function contains supervoxel grouping
-
+fprintf(['\n' 'Creating ' num2str(N_svox_in) ' SupVox:\n' ])
 canvas2 = sqrt(bwdist(1-mask));
 % orthoslice (canvas2)
 
@@ -43,12 +43,23 @@ end
 
 canvas3 = canvas2(ymin:ymax, xmin:xmax, zmin:zmax);
 
-[L,NumLabels] = superpixels3(canvas3,N_svox);
+N_svox = N_svox_in; % number of supervoxels to give as param to parse
+while true
+    [L,NumLabels] = superpixels3(canvas3,N_svox);
 
-superMask = zeros(size(mask));
-superMask(ymin:ymax, xmin:xmax, zmin:zmax) = L;
-superMask(logical(1-mask)) = 0;
+    superMask = zeros(size(mask));
+    superMask(ymin:ymax, xmin:xmax, zmin:zmax) = L;
+    superMask(logical(1-mask)) = 0;
 
-disp([num2str(numel(unique(superMask))-1) ' areas created.' ])
+    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
+        break
+    end
+    N_svox = round(N_svox * N_svox_in/numSupVox);
+    fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target.\n'])
+end 
+fprintf([' -  ok.\n'])
 
 end

planEvaluation/plot_DVH.m

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