Wrap up changes for Marcel

WiscPlanPhotonkV125/matlab_frontend/NLP_optimizer_v3.m

@@ -0,0 +1,383 @@
+
+
+% paths.patient_dir
+% paths.Goal_dir (previously called DP_dir)
+% paths.patient
+% paths.goalsName
+
+% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
+% orthoslice(D_full, [0,70])
+
+function [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v2(varargin)
+% This function performs the beamlet optimization
+% [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
+
+if nargin<2
+    load('WiscPlan_preferences.mat')
+    [Pat_path] = uigetdir([WiscPlan_preferences.patientDataPath ], 'Select Patient folder');
+    [Goal_file,Goal_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select OptGoal file');
+    
+    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
+    path2goal = [Goal_path, Goal_file];
+else
+    Pat_path = varargin{1};
+    path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
+    path2goal = varargin{2};
+    [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};
+
+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 -
+fprintf('starting NLP optimization process... \n')
+
+% % -- LOAD GEOMETRY, GOALS, BEAMLETS --
+load(path2geometry)
+load(path2goal)
+[beamlets, numBeamlet] = get_beamlets(Geometry, Pat_path);
+% [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)
+    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, :));
+    
+        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, SupVox_num); 
+        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);
+        
+        h_w1 = waitbar(0, 'merging superboxels');
+        for i_supVox = 1:numel(superVoxList)
+            waitbar(i_supVox/numel(superVoxList), h_w1)
+            supVox_idx = superVoxList(i_supVox);
+            idxList = find(superMask == supVox_idx);
+            optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
+            % -- make new indeces
+            optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
+        end
+        close(h_w1)
+    end
+end
+
+
+
+% -- make them robust --
+RO_params=0;
+optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
+
+% -- CALLBACK OPTIMIZATION FUNCTION --
+fun1 = @(x) get_penalty(x, optGoal_beam);
+fun2 = @(x) get_penalty(x, optGoal);
+
+% -- OPTIMIZATION PARAMETERS --
+% define optimization parameters
+A = [];
+b = [];
+Aeq = [];
+beq = [];
+lb = zeros(1, numBeamlet);
+% lb_beam = zeros(1, numBeam);
+ub = [];
+nonlcon = [];
+
+% define opt limits, and make it fmincon progress
+options = optimoptions('fmincon');
+options.MaxFunctionEvaluations = N_fcallback1;
+options.Display = 'iter';
+options.PlotFcn = 'optimplotfval';
+% options.UseParallel = true;
+options.UseParallel = false;
+% options.OptimalityTolerance = 1e-9;
+
+%% -- 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 = ones(numBeamlet,1);
+        w0 = mean(optGoal{1}.D_final(optGoal{1}.ROI_idx) ./ (optGoal{1}.beamlets_pruned*w0+0.1)) * w0;
+        w_beamlets = double(w0);
+
+        % -- GET BEAM WEIGHTS --
+%         tic
+%         w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb,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
+
+%% FULL OPTIMIZATION
+
+% -- GET FULL BEAMLET WEIGHTS --
+options.MaxFunctionEvaluations = N_fcallback2;
+tic
+fprintf('\n running full optimizer:')
+w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
+fprintf('  done!:')
+t=toc;
+disp(['Optimization time for beamlets = ',num2str(t)]);
+
+
+%% evaluate the results
+D_full = reshape(beamlets * w_fin, size(Geometry.data));
+
+%% save outputs
+NLP_result.dose = D_full;
+NLP_result.weights = w_fin;
+save([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
+
+
+plot_DVH(Geometry, D_full)
+colorwash(Geometry.data, D_full, [500, 1500], [0, 36]);
+
+
+end
+
+%% support functions
+% ---- PENALTY FUNCTION ----
+function penalty = get_penalty(x, optGoal)
+    % this function gets called by the optimizer. It checks the penalty for
+    % all the robust implementation and returns the worst result.
+    
+    NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
+    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
+                fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i);
+                sc_i = sc_i + 1;
+            end
+        end
+    end
+    % take the worst case penalty of evaluated scenarios
+    penalty=max(fobj);
+end
+% ------ supp: penalty for single scenario ------
+function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
+    penalty = 0;
+    % for each condition
+    for goal_i = 1:numel(optGoal)
+        switch optGoal{goal_i}.function
+            % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
+            case 'min'
+                % penalize if achieved dose is lower than target dose
+                d_penalty = 1.0e0 * sum(max(0, ...
+                    (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)));
+            case 'max'
+                % penalize if achieved dose is higher than target dose
+                d_penalty = 1.0e0 * sum(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)));
+            case 'min_sq'
+                % 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.*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.*temp1);
+            case 'min_exp'
+                % 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(exp(temp1));
+            case 'max_exp'
+                % 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(exp(temp1));
+            case 'LeastSquare'
+                % penalize with sum of squares any deviation from target
+                % dose
+                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.0e4 * min(perc_vox-0.05, 0)
+                
+            case 'max_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.0e4 * min(perc_vox-0.05, 0)
+                    
+        end
+        penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
+    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
+    
+    % RO_params - should have the information below
+    % nrs - random scenarios
+    % sss - system setup scenarios
+    % rgs - random range scenarios
+    
+    % X - X>0 moves image right
+    % Y - Y>0 moves image down
+    % Z - in/out.
+    
+    shift_mag = 1; % vox of shift
+    nrs_scene_list={[0,0,0]};
+
+    
+    % ----====#### CHANGE ROBUSTNESS HERE ####====----
+    
+%     sss_scene_list={[0,0,0]};
+    sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
+%     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 ####====----
+
+    rgs_scene_list={[0,0,0]};
+
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
+    
+    for i = 1:numel(optGoal)
+        optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
+        optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
+        optGoal{i}.NbrRangeScenarios    =numel(rgs_scene_list);
+    end
+    
+
+    for goal_i = 1:numel(optGoal)
+        % get target
+        idx=optGoal{goal_i}.ROI_idx;
+        targetImg1=zeros(optGoal{goal_i}.imgDim);
+        targetImg1(idx)=1;
+        % get beamlets
+        
+        for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios          % num. of random scenarios
+            % modify target and beamlets
+            targetImg2=targetImg1;
+            % beamlets stay the same
+            
+            for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
+                % modify target and beamlets
+                [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
+                % beamlets stay the same
+                
+                for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rgs
+                    % modify target and beamlets
+                    targetImg4=targetImg3;
+                    % beamlets stay the same
+                    
+                    %% make new target and beamlets
+                    ROI_idx=[];
+                    ROI_idx=find(targetImg4>0);
+                    
+                    target = optGoal{goal_i}.D_final(idxValid);
+                    
+                    beamlets_pruned = beamlets(ROI_idx, :);
+                    
+                    % save to optGoal output
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx            = ROI_idx;
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned    = beamlets_pruned;
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target             = target;
+                end
+            end
+        end
+    end
+end
+% ------ supp: RO case SSS ------
+function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
+    % translate the target image 
+    
+    targetImg3 = imtranslate(targetImg2,sss_scene_shift);
+    
+    % now we need to figure out if any target voxels fell out during the
+    % shift
+    imgValid = imtranslate(targetImg3,-sss_scene_shift);
+    imgInvalid = (targetImg2-imgValid);
+    
+    idx_1 = find(targetImg2);
+    idx_2 = find(imgInvalid);
+    
+    [idxValid,ia] = setdiff(idx_1,idx_2);
+    
+    [C,ia, ib] = intersect(idx_1,idxValid);
+    
+end
+
+
+
+

WiscPlanPhotonkV125/matlab_frontend/get_beamlets.m

@@ -0,0 +1,90 @@
+
+
+function [beamlets, numBeamlet] = get_beamlets(Geometry, patient_dir)
+% this function loads and returns beamlets and joined beams
+
+beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
+beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
+
+numVox  = numel(Geometry.data);
+numBeamlet = size(beamlet_cell_array,2);
+
+if size(Geometry.data, 1)<129
+    batchSize = 400;
+else
+    batchSize = 200;
+end
+
+beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
+% --- join beamlets into beams
+% load([patient_dir '\all_beams.mat'])
+% beamletOrigin=[0 0 0];
+% beam_i=0;
+% beam_i_list=[];
+% for beamlet_i = 1:numel(all_beams)
+%     newBeamletOrigin = all_beams{beamlet_i}.ip;
+%     if any(newBeamletOrigin ~= beamletOrigin)
+%         beam_i = beam_i+1;
+%         beamletOrigin = newBeamletOrigin;
+%     end
+%     beam_i_list=[beam_i_list, beam_i];
+% end
+% 
+% 
+% wbar2 = waitbar(0, 'merging beamlets into beams');
+% numBeam=numel(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(:,beam_i) = beam_full;
+%     waitbar(beam_i/numBeam, wbar2)
+% end
+% 
+% close(wbar2)
+end
+
+% ---- GET BEAMLETS ----
+function beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize);
+    wbar1 = waitbar(0, 'Creating beamlet array');
+    numBeam = size(beamlet_cell_array,2);
+    beamlets = sparse(0, 0);
+    for beam_i=1:numBeam
+        % for each beam define how much dose it delivers on each voxel
+        idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
+
+        % break the beamlets into multiple batches
+        if rem(beam_i, batchSize)==1;
+            beamlet_batch = sparse(numVox, batchSize);
+            beam_i_temp=1;
+        end
+
+%         beamlet_batch(idx, beam_i_temp) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
+        beamlet_batch(idx, beam_i_temp) = beamlet_cell_array{1, beam_i}.non_zero_values;
+        waitbar(beam_i/numBeam, wbar1, ['Adding beamlet array: #', num2str(beam_i)])
+
+        % add the batch to full set when filled
+        if rem(beam_i, batchSize)==0;
+            beamlets =[beamlets, beamlet_batch];
+        end
+        % crop and add the batch to full set when completed
+        if beam_i==numBeam;
+            beamlet_batch=beamlet_batch(:, 1:beam_i_temp);
+            beamlets =[beamlets, beamlet_batch];
+        end
+        beam_i_temp=beam_i_temp+1;
+
+    end
+    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.
+    
+    beam=sum(beamlets(:, Beam_list), 2);
+    
+    beamImg=reshape(full(beam), IMGsize);
+        
+    orthoslice(beamImg)
+    
+end

WiscPlanPhotonkV125/matlab_frontend/get_full_dose.m

@@ -9,3 +9,7 @@ optResults.dose{end+1} = D;
 optResults.weights{end+1} = optResults.weights{end};
 save(optResultsFile, 'optResults', 'optSettings');
 
+%% save also in alter format
+NLP_result.dose = double(D);
+NLP_result.weights = double(optResults.weights{end});
+save([path_in, '\matlab_files\Classical_opt.mat'], 'NLP_result');

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -45,10 +45,10 @@ 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;
+% 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

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -444,8 +444,9 @@ classdef XTPS < handle
 %-------------------------------------------------------------------------------
 % Grozomah
         function ImgDownsample_Callback(obj, src, evt)
-            disp('Downsampling data ... to be implemented.')
+            disp('Making downsampled geometry')
             
+            make_downsampled_Geometry(obj.handles.hSVPS.Geometry)
             
         end
         
@@ -556,7 +557,7 @@ classdef XTPS < handle
             Pat_path = [obj.handles.hSVPS.Geometry.patient_dir];
             [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);
+            [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(Pat_path, path2goal);
             
             % This part needs work still
 %             obj.handles.hSVPS.optResults = optResults;

WiscPlanPhotonkV125/matlab_frontend/superpix_group.m

@@ -44,6 +44,7 @@ end
 canvas3 = canvas2(ymin:ymax, xmin:xmax, zmin:zmax);
 
 N_svox = N_svox_in; % number of supervoxels to give as param to parse
+converge_factor = 1;
 while true
     [L,NumLabels] = superpixels3(canvas3,N_svox);
 
@@ -54,11 +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
+    if abs(numSupVox-N_svox_in)/N_svox_in < 0.2/converge_factor
+        orthoslice(superMask)
+        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;
+    fprintf([' - not ok. ' num2str(numSupVox/N_svox_in) ' of target. New start size: ' num2str(N_svox) '\n'])
+    
+    if converge_factor < 0.01
+        orthoslice(superMask)
         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'])