Updated version with external goal specification.

ROI_goals_prep.m

@@ -0,0 +1,128 @@
+
+
+
+function ROI_goals_prep
+% this function creates and saves dose painting plans based on Geometry
+% and beamlets created in WiscPlan. Manually change hardcode of whatever
+% you want the plan to be!
+
+
+%% ---=== LOAD DATA ===---
+patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
+
+% load Geometry
+load([patient_dir '\matlab_files\Geometry.mat']);
+fprintf('Loaded geometry ')
+
+% load beamlets
+[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, patient_dir);
+fprintf('and beamlets.\n')
+
+%% ---=== GET OPTGOAL ===---
+ROI_goals.optGoal = make_ROI_goals_gbm_005(Geometry, beamlets);
+ROI_goals.optGoal_beam = make_ROI_goals_gbm_005(Geometry, beamlets_joined);
+
+ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+
+%% ---=== SAVE OPTGOAL ===---
+fprintf('Writing ROI_goals...')
+DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+save([DP_dir '\RODP_files\ROI_goals.mat'], 'ROI_goals')
+fprintf(' done!\n')
+end
+
+
+function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_DP_minDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_DP_maxDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'CTV_min';
+    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
+    goal_1.ROI_idx = ROI_idx;
+    goal_1.imgDim = size(Geometry.data);
+    goal_1.D_final = minDose(ROI_idx);
+    goal_1.function = 'min';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target   = minDose(ROI_idx); % minDose(ROI_idx);
+    goal_1.opt_weight = 40 / numel(ROI_idx); % normalize to volume of target area
+    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_1;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_2.name = 'CTV_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
+    goal_2.ROI_idx = ROI_idx;
+    goal_2.imgDim = size(Geometry.data);
+    goal_2.D_final = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+    goal_2.opt_weight = 2 / numel(ROI_idx); % normalize to volume of target area
+    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_2;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_3.name = 'head_max';
+    goal_3.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_3.ROI_idx = ROI_idx;
+    goal_3.imgDim = size(Geometry.data);
+    goal_3.D_final = 20;
+    goal_3.function = 'max_sq';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 7 / numel(ROI_idx); % normalize to volume of target area
+    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_3;
+    % -- END DEFINITION OF GOAL --
+   
+end
+
+function beamlets = get_beamlets(beamlet_cell_array, numVox);
+    wbar1 = waitbar(0, 'Creating beamlet array');
+    numBeam = size(beamlet_cell_array,2);
+    batchSize=100;
+    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;
+        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
+

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -21,71 +21,74 @@ N_fcallback2 = 200000;
 patient = 'gbm_005';
 
 switch patient
-    case 'patient'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
-    case 'tomoPhantom'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';        
-    case 'phantom_HD'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
-    case 'doggo'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\archive\PatientData_dog5_3';
     case 'gbm_005'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
+        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
     otherwise
         error('invalid case')
 end
 
-merge_beamlets(4, patient_dir);
+% merge_beamlets(4, patient_dir);
 
 
 %% PROGRAM STARTS HERE
 % - no tocar lo que hay debajo -
-fprintf('starting NLP optimization process... ')
+fprintf('starting NLP optimization process... \n')
 
-% -- LOAD GEOMETRY AND BEAMLETS --
+% % -- LOAD GEOMETRY AND BEAMLETS --
 load([patient_dir '\matlab_files\Geometry.mat'])
-% beamlet_batch_filename = [patient_dir '\beamlet_batch_files\' 'beamletbatch0.bin'];
-beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
-beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
-fprintf('\n loaded beamlets...')
-
-% -- SET INITIAL PARAMETERS --
-numVox  = numel(Geometry.data);
-numBeamlet = size(beamlet_cell_array,2);
-
-%% -- BEAMLET DOSE DELIVERY --
-beamlets = get_beamlets(beamlet_cell_array, numVox);
-% show_joint_beamlets(beamlets, size(Geometry.data), 7:9)
-fprintf('\n beamlet array constructed...')
-% - merge 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
-beamlets_joined=[];
-target_joined=[];
-wbar2 = waitbar(0, 'merging beamlets into beams');
-numBeam=numel(unique(beam_i_list));
-for beam_i = unique(beam_i_list)
-    beam_full = sum(beamlets(:,beam_i_list == beam_i), 2); 
-    beamlets_joined(:,end+1) = beam_full;
-    waitbar(beam_i/numBeam, wbar2)
-end
-close(wbar2)
+% % beamlet_batch_filename = [patient_dir '\beamlet_batch_files\' 'beamletbatch0.bin'];
+% beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
+% beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
+% fprintf('\n loaded beamlets...')
+% 
+% % -- SET INITIAL PARAMETERS --
+% numVox  = numel(Geometry.data);
+% numBeamlet = size(beamlet_cell_array,2);
+% 
+% %% -- BEAMLET DOSE DELIVERY --
+% beamlets = get_beamlets(beamlet_cell_array, numVox);
+% % show_joint_beamlets(beamlets, size(Geometry.data), 7:9)
+% fprintf('\n beamlet array constructed...')
+% % - merge 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
+% beamlets_joined=[];
+% target_joined=[];
+% wbar2 = waitbar(0, 'merging beamlets into beams');
+% numBeam=numel(unique(beam_i_list));
+% for beam_i = unique(beam_i_list)
+%     beam_full = sum(beamlets(:,beam_i_list == beam_i), 2); 
+%     beamlets_joined(:,end+1) = beam_full;
+%     waitbar(beam_i/numBeam, wbar2)
+% end
+% close(wbar2)
 
 
 %% -- OPTIMIZATION TARGETS --
-make_ROI_goals(Geometry, beamlets, beamlets_joined, patient);
+% make_ROI_goals(Geometry, beamlets, beamlets_joined, patient);
+
+% [optGoal, optGoal_beam, optGoal_idx, targetMinMax_idx] = get_ROI_goals(patient);
+
+load([DP_dir '\RODP_files\ROI_goals.mat']);
+
+optGoal = ROI_goals.optGoal;
+optGoal_beam = ROI_goals.optGoal_beam;
+optGoal_idx = ROI_goals.optGoal_idx;
+targetMinMax_idx = ROI_goals.targetMinMax_idx;
+
+[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, patient_dir);
 
-[optGoal, optGoal_beam, optGoal_idx, targetMinMax_idx] = get_ROI_goals(patient);
 
 % -- make them robust --
 RO_params=0;
@@ -236,52 +239,6 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
     end
 end
 
-% ---- GET BEAMLETS ----
-function beamlets = get_beamlets(beamlet_cell_array, numVox);
-    wbar1 = waitbar(0, 'Creating beamlet array');
-    numBeam = size(beamlet_cell_array,2);
-    batchSize=100;
-    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;
-        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
-
 
 % ---- MAKE ROI ROBUST ----
 function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);

WiscPlanPhotonkV125/matlab_frontend/get_beam_lets.m

@@ -0,0 +1,81 @@
+
+
+function [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(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);
+beamlets = get_beamlets(beamlet_cell_array, numVox);
+% 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
+beamlets_joined=[];
+target_joined=[];
+wbar2 = waitbar(0, 'merging beamlets into beams');
+numBeam=numel(unique(beam_i_list));
+for beam_i = unique(beam_i_list)
+    beam_full = sum(beamlets(:,beam_i_list == beam_i), 2); 
+    beamlets_joined(:,end+1) = beam_full;
+    waitbar(beam_i/numBeam, wbar2)
+end
+close(wbar2)
+end
+
+% ---- GET BEAMLETS ----
+function beamlets = get_beamlets(beamlet_cell_array, numVox);
+    wbar1 = waitbar(0, 'Creating beamlet array');
+    numBeam = size(beamlet_cell_array,2);
+    batchSize=100;
+    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;
+        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

vol_prep.m

@@ -0,0 +1,249 @@
+
+
+function vol_prep
+% this function calls three other functions:
+% - vol_prep_bin
+% - vol_prep_fuzzy
+% - vol_prep_DPmap
+% that are needed for expanding the initial segmentations done in the
+% segmentation GUI
+
+fprintf('Starting binary volume expansion...')
+vol_prep_bin
+fprintf(' done!\n')
+
+fprintf('Starting fuzzy volume expansion...')
+vol_prep_fuzzy
+fprintf(' done!\n')
+
+fprintf('Creating dose plans ...')
+vol_prep_DPmap
+fprintf(' done!\n')
+end
+
+
+function vol_prep_bin
+% this function takes in the segmentations and creates PTV/GTV/CTV binary volumes,
+% as well as appropriately cropped healthy tissue (body) segmentation.
+
+CTV_margin = 0.5; % GTV->CTV margin, in cm
+PTV_margin = 0.5; % GTV->CTV margin, in cm
+body_margin = 7.0; % how far from tumor do we still include body, in cm
+
+pathIn = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\';
+
+[CT_in, CT_meta] = nrrdread([pathIn 'FET_FGL005_B1_CT2FET.nrrd']);
+[GTV_in, GTV_meta] = nrrdread([pathIn 'B1_seg\FET_FGL005_B1_seg_thr2.0.nrrd']);
+[body_in, body_meta] = nrrdread([pathIn 'B1_seg\FET_FGL005_B1_head.nrrd']);
+
+
+%% identify target voxels - these are used when identifiying valid beamlets
+% voxels that need to be considered (above the threshold)
+p_threshold = 0;
+GTV_binary = zeros(size(GTV_in));
+GTV_binary(GTV_in>p_threshold) = 1;
+
+% Account for infiltration (GTV -> CTV)
+CTV = zeros(size(GTV_binary));
+bwD = bwdistsc(GTV_binary, GTV_meta.spacedirections);
+CTV(bwD<CTV_margin) = 1;
+
+% expand to PTV (CTV -> PTV)
+PTV = zeros(size(GTV_binary));
+bwD = bwdistsc(GTV_binary, GTV_meta.spacedirections);
+PTV(bwD<PTV_margin) = 1;
+
+%% identify body volume of interest
+body = body_in;
+bwD = bwdistsc(GTV_binary, GTV_meta.spacedirections);
+body(PTV>0) = 0;
+body(bwD>body_margin) = 0; % how far from the tumor are we still interested in the body
+
+
+%% save outputs
+filename = [pathIn 'B1_seg\FET_FGL005_B1_seg_thr2.0_binPTV.nrrd'];
+matrix = single(PTV);
+pixelspacing = GTV_meta.spacedirections;
+origin = GTV_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
+filename = [pathIn 'B1_seg\FET_FGL005_B1_head_crop_bin.nrrd'];
+matrix = single(body);
+pixelspacing = body_meta.spacedirections;
+origin = body_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
+
+colorwash(CT_in, PTV,  [-500, 500], [0, 1.2]);
+colorwash(CT_in, body, [-500, 500], [0, 1.2]);
+
+end
+
+function vol_prep_fuzzy
+% this function takes in the segmentations and creates fuzzy (probabilistic)
+% PTV/GTV/CTV volumes
+
+infilt_range = 1.2; % range of infiltration, in cm. Should be value of sigma (or 1/3 decay range). Assuming normal infil.
+
+pathIn = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\';
+
+[CT_in, CT_meta] = nrrdread([pathIn 'FET_FGL005_B1_CT2FET.nrrd']);
+[GTV_in, GTV_meta] = nrrdread([pathIn 'B1_seg\FET_FGL005_B1_seg_combined.nrrd']);
+[body_in, body_meta] = nrrdread([pathIn 'B1_seg\FET_FGL005_B1_head.nrrd']);
+
+% colorwash(CT_in, GTV_in,  [-500, 500], [0, 1.2]);
+
+%% Account for infiltration (GTV -> CTV)
+
+% -- get distance from central dot --
+ir_a = 1.5; % infiltration range factor
+
+CTV = zeros(size(GTV_in));
+bwD = bwdistsc(GTV_in>0.05, GTV_meta.spacedirections);
+CTV(bwD< (ir_a * infilt_range) ) = 1;
+
+% for each voxel in the area of interest, figure out the highest
+% probability based on range from tumor * p that tumor is there
+GTV_idx_lst = find(GTV_in   >0.05);
+CTV_idx_lst = find(CTV      >0.0);
+prob_CTV = zeros(size(GTV_in));
+
+f = waitbar(0,'calculating expansion');
+
+% define area of interest
+Ksize = 1+2*ceil((ir_a * infilt_range)/min(GTV_meta.spacedirections)); % get kernel for convolution.
+kernel_base = zeros(Ksize,Ksize,Ksize); 
+kernel_base(ceil(Ksize/2),ceil(Ksize/2),ceil(Ksize/2))=1;
+bwD = bwdistsc(kernel_base, GTV_meta.spacedirections);
+kernel = zeros(size(kernel_base));
+kernel(bwD< (ir_a * infilt_range)) = 1;
+
+CTV_idx_lst_num = numel(CTV_idx_lst);
+prob_CTV_solutions = zeros(1,CTV_idx_lst_num);
+
+
+tic 
+for i = 1:CTV_idx_lst_num
+    i_idx = CTV_idx_lst(i);
+    [y1, x1, z1] = ind2sub(size(GTV_in), i_idx);
+    
+    waitbar(i/numel(CTV_idx_lst),f);
+    
+    % get the crop of nearby voxels in tumor
+%     tabula = zeros(size(GTV_in));
+%     tabula(i_idx) = 1;
+%     tabula2 = convn(tabula, kernel, 'same');
+%     idx_nearby = find(tabula2>0);
+    
+    idx_nearby = f_neighbors(i_idx, size(GTV_in), ir_a, infilt_range, GTV_meta.spacedirections);
+    
+    % overlap
+    idx_GTV_nearby = intersect (idx_nearby, GTV_idx_lst);
+    
+    if numel(idx_GTV_nearby) <1
+        prob_CTV_solutions(i) = 0;
+        error('this is fucky. no voxels within range')
+    end
+    
+    p_list=zeros(1, numel(idx_GTV_nearby));
+    p_list_max = 0;
+    
+    for j = 1:numel(idx_GTV_nearby)
+        j_idx = idx_GTV_nearby(j);
+        if GTV_in(j_idx) < p_list_max
+            continue
+        end
+        
+        [y2, x2, z2] = ind2sub(size(GTV_in), j_idx);
+        d = sqrt((((y2-y1)*GTV_meta.spacedirections(1))^2 + ...
+            ((x2-x1)*GTV_meta.spacedirections(2))^2 + ((z2-z1)*GTV_meta.spacedirections(3))^2));
+        % get probability of invisible infiltration at given voxel, based
+        % on single voxel from tumor
+        p_list(j) = GTV_in(j_idx) * exp(-(d*d)/(infilt_range*infilt_range));
+        
+        p_list_max = max(p_list_max, p_list(j));
+        
+        if p_list(j) == 1
+            prob_CTV_solutions(i) = 1;
+            break
+        end
+    end
+    prob_CTV_solutions(i) = max(p_list);
+end
+toc
+prob_CTV(CTV_idx_lst) = prob_CTV_solutions;
+close(f)
+
+colorwash(CT_in, GTV_in,  [-500, 500], [0, 1.2]);
+colorwash(CT_in, prob_CTV,  [-500, 500], [0, 1.2]);
+
+
+%% save outputs
+filename = [pathIn 'RODP_files\FET_FGL005_B1_seg_fuzzyCTV.nrrd'];
+matrix = single(prob_CTV);
+pixelspacing = GTV_meta.spacedirections;
+origin = GTV_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
+end
+function neigh_mat_idx = f_neighbors(idx, mat_size, ir_a, infilt_range, voxsize);
+% this function returns indeces of the box nearby around 
+    [y, x, z] = ind2sub(mat_size, idx);
+    
+    y1 = max(1, floor(y - ir_a*infilt_range/(voxsize(1))));
+    y2 = min(mat_size(1), ceil(y + ir_a*infilt_range/(voxsize(1))));
+    
+    x1 = max(1, floor(x - ir_a*infilt_range/(voxsize(2))));
+    x2 = min(mat_size(1), ceil(x + ir_a*infilt_range/(voxsize(2))));
+    
+    z1 = max(1, floor(z - ir_a*infilt_range/(voxsize(3))));
+    z2 = min(mat_size(1), ceil(z + ir_a*infilt_range/(voxsize(3))));
+    
+    tabula = zeros(mat_size);
+    tabula(y1:y2, x1:x2, z1:z2)=1;
+    
+    neigh_mat_idx = find( tabula>0);
+
+end
+
+function vol_prep_DPmap
+% this function takes in the segmentations and creates fuzzy (probabilistic)
+% PTV/GTV/CTV volumes
+
+D_min = 80;
+D_max = 83;
+
+pathIn = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\';
+
+[CT_in, CT_meta] = nrrdread([pathIn 'FET_FGL005_B1_CT2FET.nrrd']);
+[CTV_fuzzy, CTV_meta] = nrrdread([pathIn 'RODP_files\FET_FGL005_B1_seg_fuzzyCTV.nrrd']);
+
+dose_min = f_dose_min(CTV_fuzzy, D_min);
+dose_max = f_dose_max(CTV_fuzzy, D_max);
+
+orthoslice(dose_min, [0, D_max+10]) % plot min dose boundary
+orthoslice(dose_max, [0, D_max+10]) % plot max dose boundary
+
+%% save outputs
+filename = [pathIn 'RODP_files\FET_FGL005_B1_DP_minDose.nrrd'];
+matrix = single(dose_min);
+pixelspacing = CTV_meta.spacedirections;
+origin = CTV_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
+filename = [pathIn 'RODP_files\FET_FGL005_B1_DP_maxDose.nrrd'];
+matrix = single(dose_max);
+pixelspacing = CTV_meta.spacedirections;
+origin = CTV_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
+end
+function dose_min = f_dose_min(CTV_fuzzy, D_min)
+% get low boundary for ideal dose plan
+dose_min = D_min* min(1, 2* max(0, (CTV_fuzzy - 0.25)));
+% orthoslice(dose_min)
+end
+function dose_max = f_dose_max(CTV_fuzzy, D_max)
+% get hihg boundary for ideal dose plan
+dose_max = D_max* min(1, 2* max(0, (CTV_fuzzy - 0.0)));
+end