Changed DP map files,

smaller changes

NLP_beamlet_optimizer.m

@@ -19,6 +19,9 @@ switch patient
     case 'patient'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
         blet_in_beam=5; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
+    case 'tomoPhantom'
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
+        blet_in_beam=7; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     case 'phantom_HD'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
         blet_in_beam=20; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
@@ -79,9 +82,13 @@ close(wbar2)
 %% -- OPTIMIZATION TARGETS --
 switch patient
     case 'patient'
-        optGoal = make_ROI_goals(Geometry, beamlets);
-        optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
-        N_beamlets_in_beam = 10;
+%         optGoal = make_ROI_goals(Geometry, beamlets);
+%         optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
+%         N_beamlets_in_beam = 10;
+    case 'tomoPhantom'   
+        optGoal = make_ROI_goals_2(Geometry, beamlets);
+        optGoal_beam = make_ROI_goals_2(Geometry, beamlets_joined);
+        N_beamlets_in_beam = 7;
     case 'phantom_HD'
         optGoal = make_ROI_goals(Geometry, beamlets);
         optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
@@ -159,8 +166,23 @@ NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
 % save('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat', 'NLP_result');
 
-optGoal_idx=[1,2];
-plot_DVH(D_full, optGoal, optGoal_idx)
+switch patient
+    case 'patient'
+        optGoal_idx=[1,2];
+        targetMinMax_idx=[1,3];
+    case 'tomoPhantom'   
+        optGoal_idx=[1,3];
+        targetMinMax_idx=[1,2];
+    case 'phantom_HD'
+
+    case 'doggo'
+        optGoal_idx=[1,2];
+        targetMinMax_idx=[1,3];
+    otherwise
+        error('invalid case')
+end
+
+plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
 % plot_DVH_robust(D_full, optGoal, optGoal_idx)
 end
 
@@ -220,12 +242,19 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
                 d_penalty = 1.0e-1* sum(((...
                     optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) - ...
                     optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target).^2);
-            case 'min_Volume'
+            case 'min_perc_Volume'
                 % penalize by amount of volume under threshold
-                n_vox = find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
-                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) > 0);
+                perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) > 0)) ...
+                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
+                d_penalty = 3.0e5 * min(perc_vox-0.05, 0)
                 
-%                 d_penalty = 1.0e0 * 
+            case 'max_perc_Volume'
+                % penalize by amount of volume under threshold
+                perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) < 0)) ...
+                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
+                d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
                     
         end
         penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
@@ -331,6 +360,76 @@ function optGoal = make_ROI_goals(Geometry, beamlets)
     % -- END DEFINITION OF GOAL --
     
 end
+function optGoal = make_ROI_goals_2(Geometry, beamlets)
+    optGoal={};
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'Target_min';
+    goal_1.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_1.ROI_idx = ROI_idx;
+    goal_1.imgDim = size(Geometry.data);
+    goal_1.D_final = 62;
+    goal_1.function = 'min';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1.00;
+    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
+    goal_1.opt_weight = 10 / 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 = 'Target_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_2.ROI_idx = ROI_idx;
+    goal_2.imgDim = size(Geometry.data);
+    goal_2.D_final = 63;
+    goal_2.function = 'max';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1.05;
+    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
+    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 = 'ROI 1_max';
+    goal_3.ROI_name = Geometry.ROIS{1, 3}.name;
+    ROI_idx = Geometry.ROIS{1, 3}.ind;
+    goal_3.ROI_idx = ROI_idx;
+    goal_3.imgDim = size(Geometry.data);
+    goal_3.D_final = 10;
+    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 = 2 / 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 --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_4.name = 'Ring_max';
+    goal_4.ROI_name = Geometry.ROIS{1, 6}.name;
+    ROI_idx = Geometry.ROIS{1, 6}.ind;
+    goal_4.ROI_idx = ROI_idx;
+    goal_4.imgDim = size(Geometry.data);
+    goal_4.D_final = 56;
+    goal_4.function = 'max_sq';
+    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
+    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_4.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_4;
+    % -- END DEFINITION OF GOAL --
+    
+end
 function optGoal = make_ROI_goals_DOG(Geometry, beamlets)
     optGoal={};
     
@@ -469,6 +568,77 @@ function optGoal = make_ROI_goals_DOG_2(Geometry, beamlets)
     % -- END DEFINITION OF GOAL --
     
 end
+function optGoal = make_ROI_goals_DOG_3(Geometry, beamlets)
+    optGoal={};
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'Target_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 = 62;
+    goal_1.function = 'min_perc_Volume';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1.00;
+    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
+    goal_1.opt_weight = 10 / 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 = 'Doggo_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
+    ROI_idx = Geometry.ROIS{1, 4}.ind;
+    goal_2.ROI_idx = ROI_idx;
+    goal_2.imgDim = size(Geometry.data);
+    goal_2.D_final = 20;
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.targetImg   = ones(numel(ROI_idx), 1) * goal_2.D_final;
+    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
+    goal_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
+    goal_2.dvh_col = [0.2, 0.9, 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 = 'Target_max';
+    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
+    goal_3.ROI_idx = ROI_idx;
+    goal_3.imgDim = size(Geometry.data);
+    goal_3.D_final = 63;
+    goal_3.function = 'max_perc_Volume';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1.05;
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 0.5 / numel(ROI_idx); % normalize to volume of target area
+    goal_3.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_3;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_4.name = 'Doggo_max2';
+    goal_4.ROI_name = Geometry.ROIS{1, 4}.name;
+    ROI_idx = Geometry.ROIS{1, 4}.ind;
+    goal_4.ROI_idx = ROI_idx;
+    goal_4.imgDim = size(Geometry.data);
+    goal_4.D_final = 50;
+    goal_4.function = 'max_sq';
+    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
+    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_4.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_4;
+    % -- END DEFINITION OF GOAL --
+    
+end
 
 % ---- MAKE ROI ROBUST ----
 function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
@@ -489,7 +659,9 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 %     nrs_scene_list={[0,0,0]};
 %     sss_scene_list={[0,0,0], [-2,0,0], [2,0,0], [0,-2,0], [0,2,0]};
 %     rrs_scene_list={[0,0,0]};
+
     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
     
     for i = 1:numel(optGoal)
         optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
@@ -551,7 +723,7 @@ function targetImg3=get_RO_sss(targetImg2, sss_scene_shift);
 end
 
 % ---- DVH PLOT FUNCTION ----
-function plot_DVH(dose, optGoal, optGoal_idx)
+function plot_DVH(dose, optGoal, optGoal_idx, targetMinMax_idx)
     % this function plots the DVHs of the given dose
     nfrac=1;
     
@@ -564,6 +736,14 @@ function plot_DVH(dose, optGoal, optGoal_idx)
         plot(dosebins, dvh,'Color', optGoal{roi_idx}.dvh_col,'LineStyle', '-','DisplayName', optGoal{roi_idx}.ROI_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(1)}.target*1.07 >0)) ...
+        / num_vox;
+    title([num2str(perc_vox_min) ' % below 0.95D_0, ' num2str(perc_vox_max) ' % above D 1.07D_0'])
+    
     
 end
 function plot_DVH_robust(dose, optGoal, optGoal_idx)

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -49,7 +49,7 @@ ypmax = 1.0;    % +jaw width / 2
 % y-prime points in the z-direction in the CT coordinate system
 
 % Number of beamlets in the BEV for each direction
-Mxp = 20;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
+Mxp = 7;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
 Nyp = 1;  % always 1 for Tomo due to binary mlc
 
 
@@ -93,7 +93,7 @@ Nrot = ceil(abs(zBow - zStern)/(pitch*fieldWidth));
 
 
 % Nphi = Nrot*51;  % number of angles used in the calculation
-Nphi = Nrot*51;  % Grozomah
+Nphi = Nrot*21;  % Grozomah
 
 % load(geometry_file);
 

make_DP_maps_Tomo.m

@@ -0,0 +1,54 @@
+
+
+if true
+    %% load up NRRD
+
+    
+    load('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\matlab_files\Geometry.mat')
+    
+    CTin = Geometry.data;
+%     orthoslice(CTin)
+    
+    ROImap=zeros(size(CTin));
+    ROImap(Geometry.ROIS{1, 2}.ind)=1;
+    
+    doseMap = zeros(size(ROImap));
+    for z = 1:size(ROImap, 3)
+        doseMap(:,:,z) = bwdist(1-ROImap(:,:,z));
+    end
+%     orthoslice(doseMap)
+    
+    matrix = 50* ones(size(CTin))+2* doseMap;
+    
+%     matrix = 50* ones(size(CTin));
+%     temp = 50* ones(size(CTin)) + 2*(max(doseMap(:))-doseMap);
+%     matrix(Geometry.ROIS{1, 2}.ind)=temp(Geometry.ROIS{1, 2}.ind);
+%     orthoslice(matrix)
+    
+    
+    % save nrrds
+    filename = ['C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd'];
+    pixelspacing = Geometry.voxel_size;
+    origin = Geometry.start;
+    encoding ='raw'; % 'raw', 'ascii' or 'gzip'
+
+    ok = nrrdWriter(filename, matrix, pixelspacing, origin, encoding)
+
+end
+if false
+    
+    [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
+    
+    orthoslice(matrix)
+    orthoslice(D_full)
+    
+    colorwash(Geometry.data, matrix, [500, 1500], [0,70])
+    colorwash(Geometry.data, D_full, [500, 1500], [0,70])
+    
+    slice_i = 24;
+    
+    D_diff = D_full - matrix;
+    orthoslice(D_diff)
+    
+    
+end

test_ReadBeamlets.m

@@ -14,7 +14,7 @@ beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
 
 % close all
 
-N=2002;
+N=95;
 data=beamlet_cell_array{N};
 tabula=zeros(data.x_count, data.y_count, data.z_count);
 tabula(data.non_zero_indices)=data.non_zero_values;