Added superpixel grouping

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v2.m

@@ -10,12 +10,18 @@
 
 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;
 % 
+% Inputs:
+%       () OR
+%       (Pat_path, path2goal) OR
+%       (Pat_path, path2goal, beamlet_weights)
+%   Pat_path, path2goal = strings to patient folder and optimal goals
+%   beamlet_weights = initial beamlet weights
+%   
+% Outputs:
+%       full dose image dose: [D_full, w_fin, Geometry, optGoal]
+% 
 % Made by Peter Ferjancic 1. May 2018
 % Last updated: 1. April 2019
 
@@ -32,8 +38,25 @@ else
     path2goal = varargin{2};
 end
 
-N_fcallback1 = 10000;
-N_fcallback2 = 200000;
+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], {'10000', '500000', 'n'});
+N_fcallback1 = str2double(str{1}); % 10000  is a good guesstimate
+N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
+pre_beamWeights = str{3};
+
+switch pre_beamWeights
+    case 'y'
+        [NLP_file,NLP_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select NLP_result file');
+        load([NLP_path, NLP_file])
+        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
+    case 'n'
+        disp('Initial beam weights will be calculated.')
+end
 
 %% PROGRAM STARTS HERE
 % - no tocar lo que hay debajo -
@@ -44,30 +67,70 @@ 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, :));
+    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 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, :));
+        optGoal_beam{i_goal} = OptGoals.data{i_goal};
+        optGoal_beam{i_goal}.beamlets_pruned = sparse(beamlets_joined(optGoal{i_goal}.ROI_idx, :));
+    
+        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})); 
+        
+        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);
+        
+        optGoal_beam{i_goal} = OptGoals.data{i_goal};
+        optGoal_beam{i_goal}.ROI_idx_old = optGoal_beam{i_goal}.ROI_idx; 
+        optGoal_beam{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
+        optGoal_beam{i_goal}.opt_weight = optGoal_beam{i_goal}.opt_weight * numel(optGoal_beam{i_goal}.ROI_idx_old)/numel(optGoal_beam{i_goal}.ROI_idx);
+        
+        h_w1 = waitbar(0, 'merging superboxels');
+        for i_supVox = 1:numel(superVoxList)
+            supVox_idx = superVoxList(i_supVox);
+            
+            waitbar(i_supVox/numel(superVoxList), h_w1)
+            
+            idxList = find(superMask == supVox_idx);
+            
+            optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
+            optGoal_beam{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets_joined(idxList, :),1));
+            % -- make new indeces
+            optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
+            optGoal_beam{i_goal}.ROI_idx(i_supVox) = idxList(1);
+        end
+        close(h_w1)
+    end
+    
+    
 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);
@@ -88,31 +151,43 @@ options = optimoptions('fmincon');
 options.MaxFunctionEvaluations = N_fcallback1;
 options.Display = 'iter';
 options.PlotFcn = 'optimplotfval';
-options.UseParallel = true;
+% options.UseParallel = true;
+options.UseParallel = false;
 % 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);
+%% -- INITIALIZE BEAMLET WEIGHTS --
+switch pre_beamWeights
+    case 'y'
+        % should have been assigned previously.
+        disp('Provided beamlet weights used for initial comparison')
+    case 'n'
+        % if initial beamlet weights are not provided, get quick estimate
+        fprintf('\n running initial optimizer:')
+        % initialize beamlet weights, OR
+        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 = double(w0_beams);
+
+        % -- GET 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
 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
+
+%% FULL OPTIMIZATION
 
 % -- GET FULL BEAMLET WEIGHTS --
 options.MaxFunctionEvaluations = N_fcallback2;
-% tic
+tic
 fprintf('\n running full optimizer:')
 w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
 fprintf('  done!:')
@@ -128,10 +203,11 @@ 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)
+
+plot_DVH(Geometry, D_full)
+colorwash(Geometry.data, D_full, [-500, 500], [0, 35]);
+
+
 end
 
 %% support functions
@@ -173,27 +249,27 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
                     (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
+                % penalize if achieved dose is lower 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);
+                d_penalty = 1.0e0 * sum(temp1.*temp1);
             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);
+                d_penalty = 1.0e0 * sum(temp1.*temp1);
             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);
+                temp1 = (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 '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)
+                d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
                 
             case 'max_perc_Volume'
                 % penalize by amount of volume under threshold
@@ -221,7 +297,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % Y - Y>0 moves image down
     % Z - in/out.
     
-    shift_mag = 2; % vox of shift
+    shift_mag = 1; % vox of shift
     nrs_scene_list={[0,0,0]};
 
     

WiscPlanPhotonkV125/matlab_frontend/dicomrt/dicomrt2geometry.m

@@ -24,6 +24,7 @@ if nargin < 1
 %     dicomrt_dir = uifile('getdir', 'Select the directory contains the DICOM-RT images');
 % !!! Grozomah
     load('WiscPlan_preferences.mat')
+    
     if isstring(WiscPlan_preferences.inDataPath) 
     else
         WiscPlan_preferences.inDataPath = 'C://';
@@ -119,7 +120,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;
+        downsample_factor = 1/8;
         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/get_beam_lets.m

@@ -22,15 +22,17 @@ for beamlet_i = 1:numel(all_beams)
     end
     beam_i_list=[beam_i_list, beam_i];
 end
-beamlets_joined=[];
-target_joined=[];
+
+
 wbar2 = waitbar(0, 'merging beamlets into beams');
 numBeam=numel(unique(beam_i_list));
-for beam_i = unique(beam_i_list)
+beamlets_joined=sparse(size(beamlets,1), numBeam);
+for beam_i = 1:numel(unique(beam_i_list))
     beam_full = sum(beamlets(:,beam_i_list == beam_i), 2); 
-    beamlets_joined(:,end+1) = beam_full;
+    beamlets_joined(:,beam_i) = beam_full;
     waitbar(beam_i/numBeam, wbar2)
 end
+
 close(wbar2)
 end
 
@@ -38,7 +40,7 @@ end
 function beamlets = get_beamlets(beamlet_cell_array, numVox);
     wbar1 = waitbar(0, 'Creating beamlet array');
     numBeam = size(beamlet_cell_array,2);
-    batchSize=100;
+    batchSize=400;
     beamlets = sparse(0, 0);
     for beam_i=1:numBeam
         % for each beam define how much dose it delivers on each voxel

WiscPlanPhotonkV125/matlab_frontend/goal_def_UI.m

@@ -185,7 +185,10 @@ for i = 1 : size(handles.uitable1.Data, 1)
     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')
+        warning('Works only for NRRD (for now)')
+        
+        dose_in = nrrdread(handles.uitable1.Data{i,4});
+        goal_i.D_final = dose_in(goal_i.ROI_idx);
         % loading of data comes here
         % matrix size checks
     end

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -38,16 +38,12 @@ 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
-% ypmin = -0.5;   % -jaw width / 2
-% ypmax = 0.5;    % +jaw width / 2
-
-% ##
 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;
 
 % y-prime points in the z-direction in the CT coordinate system
 
@@ -55,8 +51,10 @@ ypmax = 0.3125;
 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');
+
 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/lab/XTPS.m

@@ -557,6 +557,10 @@ classdef XTPS < handle
             [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);
+            
+            % This part needs work still
+%             obj.handles.hSVPS.optResults = optResults;
+            
         end
 
 %-------------------------------------------------------------------------------

WiscPlanPhotonkV125/matlab_frontend/merge_beamlets.m

@@ -20,7 +20,7 @@ for beam_i = 1:numel(B)
     
     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

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -0,0 +1,54 @@
+
+
+function superMask = superpix_group(mask, N_svox)
+% this function contains supervoxel grouping
+
+canvas2 = sqrt(bwdist(1-mask));
+% orthoslice (canvas2)
+
+% find box crop of ROI
+ymin = 0;
+ymax = 0;
+for yi = 1:size(mask,1)
+    data = mask(yi, :,:);
+    if max(data(:)) >0
+        if ymin == 0
+            ymin = yi;
+        end
+        ymax = yi;
+    end
+end
+xmin = 0;
+xmax = 0;
+for xi = 1:size(mask,2)
+    data = mask(:,xi,:);
+    if max(data(:)) >0
+        if xmin == 0
+            xmin = xi;
+        end
+        xmax = xi;
+    end
+end
+zmin = 0;
+zmax = 0;
+for zi = 1:size(mask,3)
+    data = mask(:,:,zi);
+    if max(data(:)) >0
+        if zmin == 0
+            zmin = zi;
+        end
+        zmax = zi;
+    end
+end
+
+canvas3 = canvas2(ymin:ymax, xmin:xmax, zmin:zmax);
+
+[L,NumLabels] = superpixels3(canvas3,N_svox);
+
+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.' ])
+
+end

planEvaluation/plot_DVH.m

@@ -1,29 +1,32 @@
 % ---- DVH PLOT FUNCTION ----
-function plot_DVH(dose, optGoal)
+function plot_DVH(Geometry, dose)
     % this function plots the DVHs of the given dose
     
+    nfrac=1;
     
     
-    nfrac=1;
+    colorList = jet(numel(Geometry.ROIS));
+    names = {};
     
     figure
     hold on
-    for roi_idx=optGoal_idx
+    for roi_idx= 1:numel(Geometry.ROIS);
         % plot the histogram
         
-        [dvh dosebins] = get_dvh_data(dose,optGoal{roi_idx}.ROI_idx,nfrac);
-        plot(dosebins, dvh,'Color', optGoal{roi_idx}.dvh_col,'LineStyle', '-','DisplayName', optGoal{roi_idx}.ROI_name);
+        [dvh dosebins] = get_dvh_data(dose, Geometry.ROIS{roi_idx}.ind, nfrac);
+        plot(dosebins, dvh,'Color', colorList(roi_idx, :),'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+        names{roi_idx} = Geometry.ROIS{roi_idx}.name;
     end
     hold off
     % get % of volume above/below threshold
-    num_vox = numel(optGoal{targetMinMax_idx(1)}.target);
-    perc_vox_min = 100* numel(find((dose(optGoal{targetMinMax_idx(1)}.ROI_idx) - optGoal{targetMinMax_idx(1)}.target*0.95) <0)) ...
-        / num_vox;
-    perc_vox_max = 100* numel(find((dose(optGoal{targetMinMax_idx(2)}.ROI_idx) - optGoal{targetMinMax_idx(2)}.target*1.07) >0)) ...
-        / num_vox;
-    title([num2str(perc_vox_min) ' % below 0.95D_0, ' num2str(perc_vox_max) ' % above D 1.07D_0'])
-    
-    axis([0 100 0 100])
+%     num_vox = numel(optGoal{targetMinMax_idx(1)}.target);
+%     perc_vox_min = 100* numel(find((dose(optGoal{targetMinMax_idx(1)}.ROI_idx) - optGoal{targetMinMax_idx(1)}.target*0.95) <0)) ...
+%         / num_vox;
+%     perc_vox_max = 100* numel(find((dose(optGoal{targetMinMax_idx(2)}.ROI_idx) - optGoal{targetMinMax_idx(2)}.target*1.07) >0)) ...
+%         / num_vox;
+    title(['NLP optimized DVH'])
+    legend(names)
+    axis([0 50 0 100])
 end
 function plot_DVH_robust(dose, optGoal, optGoal_idx)
     % this function plots the DVHs of the given dose