Updated files, now work for australian glioma patient #15. Next step: generalizing on the case of patient #22.

ROI_goals_prep.m

@@ -6,39 +6,56 @@ function ROI_goals_prep
 % and beamlets created in WiscPlan. Manually change hardcode of whatever
 % you want the plan to be!
 
+%% ---=== INPUT PARAMS ===---
+patient = 'gbm_015';
+
+switch patient
+    case 'gbm_015'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
+        paths.CT_in = ['FET_FGL015_B1_CT2FET'];
+        paths.target_bin_in  = ['\RODP_files\FET_FGL015_B1_thr2_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\FET_FGL015_B1_thr2_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\FET_FGL015_B1_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
+    otherwise
+        error('invalid case')
+end
 
 %% ---=== LOAD DATA ===---
-patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
 
 % load Geometry
-load([patient_dir '\matlab_files\Geometry.mat']);
+load([paths.wiscplan '\matlab_files\Geometry.mat']);
 fprintf('Loaded geometry ')
 
 % load beamlets
-[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, patient_dir);
+[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, paths.wiscplan);
 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);
+switch patient
+    case 'gbm_015'
+        ROI_goals.optGoal = make_ROI_goals_gbm_015(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_gbm_015(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
+    otherwise
+        error('invalid case')
+end
 
-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')
+save([paths.in '\RODP_files\ROI_goals.mat'], 'ROI_goals')
 fprintf(' done!\n')
 end
 
 
-function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets, minDose, maxDose)
+function optGoal = make_ROI_goals_gbm_015(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']);
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL015_B1_thr2_DP_minDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL015_B1_thr2_DP_maxDose.nrrd']);
     minDose = double(minDose);
     maxDose = double(maxDose);
     
@@ -81,7 +98,7 @@ function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets, minDose, maxDose)
     goal_3.ROI_idx = ROI_idx;
     goal_3.imgDim = size(Geometry.data);
     goal_3.D_final = 20;
-    goal_3.function = 'max_sq';
+    goal_3.function = 'max';
     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

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -18,12 +18,15 @@ function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 
 N_fcallback1 = 5000;
 N_fcallback2 = 200000;
-patient = 'gbm_005';
+patient = 'gbm_015';
 
 switch patient
     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';
+    case 'gbm_015'
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
+        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
     otherwise
         error('invalid case')
 end
@@ -37,49 +40,9 @@ fprintf('starting NLP optimization process... \n')
 
 % % -- 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)
 
 
 %% -- OPTIMIZATION TARGETS --
-% 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;
@@ -162,7 +125,7 @@ NLP_result.weights = w_fin;
 save([patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
 
 plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
-colorwash(Geometry.data, D_full);
+colorwash(Geometry.data, D_full, [-500, 500], [0, 80]);
 % plot_DVH_robust(D_full, optGoal, optGoal_idx)
 end
 

WiscPlanPhotonkV125/matlab_frontend/nrrd2geometry.m

@@ -31,7 +31,7 @@ function [Geometry patient_dir] = nrrd2geometry
     end
     
     % Geometry.rhomw
-y    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
+    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
     Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
     [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
     

vol_prep.m

@@ -1,6 +1,8 @@
 
 
 function vol_prep
+% This series of functions take in segmented PET/CT image paths and outputs
+% Dose Painting (DP) maps, as well as some intermediate results. 
 % this function calls three other functions:
 % - vol_prep_bin
 % - vol_prep_fuzzy
@@ -8,34 +10,47 @@ function vol_prep
 % that are needed for expanding the initial segmentations done in the
 % segmentation GUI
 
+close all
+
+%% ---=== INPUT PARAMS ===---
+paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
+paths.CT_in = ['FET_FGL015_B1_CT2FET'];
+paths.target_bin_in  = ['FET_FGL015_B1_thr2'];
+paths.target_fzy_in  = ['FET_FGL015_B1_thr2'];
+paths.body_bin_in    = ['FET_FGL015_B1_head'];
+
+margins.CTV = 1.0; % GTV->CTV margin, in cm
+margins.PTV = 1.0; % CTV->PTV margin, in cm
+margins.body = 7.0; % how far from tumor do we still include body, in cm
+
+dose.min = 80;
+dose.max = 83;
+
+
+%% ---=== FUNCTION CALLS ===---
 fprintf('Starting binary volume expansion...')
-vol_prep_bin
+vol_prep_bin(paths, margins)
 fprintf(' done!\n')
 
 fprintf('Starting fuzzy volume expansion...')
-vol_prep_fuzzy
+vol_prep_fuzzy(paths, margins)
 fprintf(' done!\n')
 
 fprintf('Creating dose plans ...')
-vol_prep_DPmap
+vol_prep_DPmap(paths, dose)
 fprintf(' done!\n')
 end
 
 
-function vol_prep_bin
+function vol_prep_bin(paths, margins)
 % 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']);
+[CT_in, CT_meta] = nrrdread([paths.in '\' paths.CT_in '.nrrd']);
+[GTV_in, GTV_meta] = nrrdread([paths.in '\' paths.target_bin_in '.nrrd']);
+[body_in, body_meta] = nrrdread([paths.in '\' paths.body_bin_in '.nrrd']);
 
+colorwash(CT_in, GTV_in,  [-500, 500], [0, 1.2]);
 
 %% identify target voxels - these are used when identifiying valid beamlets
 % voxels that need to be considered (above the threshold)
@@ -46,52 +61,45 @@ 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;
+CTV(bwD<margins.CTV) = 1;
 
 % expand to PTV (CTV -> PTV)
 PTV = zeros(size(GTV_binary));
 bwD = bwdistsc(GTV_binary, GTV_meta.spacedirections);
-PTV(bwD<PTV_margin) = 1;
+PTV(bwD<margins.PTV) = 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
-
+body(bwD>margins.body) = 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'];
+filename = [paths.in '\RODP_files\' paths.target_bin_in '_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'];
+filename = [paths.in '\RODP_files\' paths.body_bin_in '_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
+function vol_prep_fuzzy(paths, margins)
 % 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']);
+infilt_range = margins.CTV; % range of infiltration, in cm. Should be value of sigma (or 1/3 decay range). Assuming normal infil.
 
-% colorwash(CT_in, GTV_in,  [-500, 500], [0, 1.2]);
+[CT_in, CT_meta] = nrrdread([paths.in '\' paths.CT_in '.nrrd']);
+[GTV_in, GTV_meta] = nrrdread([paths.in '\' paths.target_bin_in '.nrrd']);
+% [body_in, body_meta] = nrrdread([paths.in '\' paths.body_bin_in '.nrrd']);
 
 %% Account for infiltration (GTV -> CTV)
 
@@ -121,7 +129,6 @@ 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);
@@ -174,12 +181,11 @@ 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'];
+filename = [paths.in '\RODP_files\' paths.target_bin_in '_fuzzyCTV.nrrd'];
 matrix = single(prob_CTV);
 pixelspacing = GTV_meta.spacedirections;
 origin = GTV_meta.spaceorigin;
@@ -206,32 +212,29 @@ function neigh_mat_idx = f_neighbors(idx, mat_size, ir_a, infilt_range, voxsize)
 
 end
 
-function vol_prep_DPmap
+function vol_prep_DPmap(paths, dose)
 % 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']);
+[CT_in, CT_meta] = nrrdread([paths.in '\' paths.CT_in '.nrrd']);
+[CTV_fuzzy, CTV_meta] = nrrdread([paths.in '\RODP_files\' paths.target_bin_in '_fuzzyCTV.nrrd']);
 
-dose_min = f_dose_min(CTV_fuzzy, D_min);
-dose_max = f_dose_max(CTV_fuzzy, D_max);
+dose_min = f_dose_min(CTV_fuzzy, dose.min);
+dose_max = f_dose_max(CTV_fuzzy, dose.max);
 
-orthoslice(dose_min, [0, D_max+10]) % plot min dose boundary
-orthoslice(dose_max, [0, D_max+10]) % plot max dose boundary
+orthoslice(dose_min, [0, dose.max+10]) % plot min dose boundary
+orthoslice(dose_max, [0, dose.max+10]) % plot max dose boundary
 
 %% save outputs
-filename = [pathIn 'RODP_files\FET_FGL005_B1_DP_minDose.nrrd'];
+filename = [paths.in '\RODP_files\' paths.target_bin_in '_DP_minDose.nrrd'];
+% 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'];
+filename = [paths.in '\RODP_files\' paths.target_bin_in '_DP_maxDose.nrrd'];
+% filename = [pathIn 'RODP_files\FET_FGL005_B1_DP_maxDose.nrrd'];
 matrix = single(dose_max);
 pixelspacing = CTV_meta.spacedirections;
 origin = CTV_meta.spaceorigin;