WiscPlan_olderMV/opt/linlsqOptSetup.m

%% linlsqOptSetup.m

% Stephen R Bowen April 2009
%
% Runs in MATLAB or writes setup files necessary to carry out an
% optimization of photon beamlet or proton spot weight distributions by
% iteratitively minimizing a linear least squares objective function.  The
% objective function has the following components:
% 
% -Tissue relative importance alpha/Ntissue (alpha binary)
% -Overdose and underdose relative penalties in each tissue Bplus & Bminus
% -Maximum and minimum doses at each voxel within a tissue dplus & dminus
% -Overdose-volume and underdose-volume relative penalties BVplus & BVminus
% -Maximum and minimum doses dVplus & dVminus evaluated at each relative
%  tissue maximum and minimum volume percentatges Vplus & Vminus 

%%

%%%%% Start of user supplied inputs %%%%%

% specify relative location of Geometry file including ROI masks
% geometryFile = '../plan/HN002b/HN002b_geometry.mat';

% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV003f\AV003f_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004c\AV004c_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004c\AV004c_cont_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004c\AV004c_PVC_cont_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004c\AV004c_thresh_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004c\AV004c_PVC_thresh_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\HN002b\HN002b_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV001c\AV001c_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV001f\AV001f_geometry.mat';
geometryFile = 'D:\FTP\Dose Painting\AAPM09\HN003c\HN003c_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV004f\AV004f_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV002f\AV002f_geometry.mat';
% geometryFile = 'D:\FTP\Dose Painting\AAPM09\AV005c\AV005c_geometry.mat';

% define optimization settings
optSettings = [];
optSettings.optInfo.Niterations = 1000;                                     % total number of iterations of beamlet weight optimization
optSettings.optInfo.Nperbatch = 250;                                        % iteration interval between writing planned dose and weights to file

optimizing = 1;
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_3MM_continuous_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_5MM_continuous_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_7MM_continuous_';

% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_3MM_SUV70_threshold_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_5MM_SUV70_threshold_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_7MM_SUV70_threshold_';

% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_3MM_SUV40_CTV54_GTV72_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_5MM_SUV40_CTV54_GTV72_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_7MM_SUV40_CTV54_GTV72_';

% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_5MM_SUV40_CTV54_GTV94_';
% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_7MM_SUV40_CTV54_GTV99_';

% optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_5MM_SUV40_CTV54_GTV65_';
optSettings.optInfo.doseBatchBaseName = 'HN003c_U01_7MM_SUV40_CTV54_GTV63_';



switch optimizing
    case 1
        optSettings.optInfo.optRun = 'yes';                                          % flag to run optimization in MATLAB
        optSettings.optInfo.saveOptInfo = 'no';                                    % save opt settings in optinput.txt file located in opt folder
        optSettings.beamletInfo.saveBeamletHeader = 'no';                           % save beamlet info in beamlet_header.txt file located in beamlet folder
        optSettings.prescInfo.savePresc = 'no';                                    % save presc info in prescription.txt file located in input folder
    case 0
        optSettings.optInfo.optRun = 'no';                                          % flag to run optimization in MATLAB
        optSettings.optInfo.saveOptInfo = 'yes';                                    % save opt settings in optinput.txt file located in opt folder
        optSettings.beamletInfo.saveBeamletHeader = 'yes';                           % save beamlet info in beamlet_header.txt file located in beamlet folder
        optSettings.prescInfo.savePresc = 'yes';                                    % save presc info in prescription.txt file located in input folder
end

% optSettings.optInfo.optRun = 'yes';                                          % flag to run optimization in MATLAB
optSettings.optInfo.waitForResults = 'no';                                  % flag to load final dose distribution and DVH into workspace
% optSettings.optInfo.saveOptInfo = 'no';                                    % save opt settings in optinput.txt file located in opt folder
optSettings.beamletInfo.saveBeamlets = 'no';                                % save beamlets when dose calc performed in MATLAB and not Condor
% optSettings.beamletInfo.saveBeamletHeader = 'no';                           % save beamlet info in beamlet_header.txt file located in beamlet folder
% % specify location of beamlet folder relative to opt folder

% optSettings.beamletInfo.beamletFolder = 'AV003fbeamlets'; 
% optSettings.beamletInfo.beamletFolder = 'AV004cbeamlets'; 
% optSettings.beamletInfo.beamletFolder = 'HN002bbeamlets';
% optSettings.beamletInfo.beamletFolder = 'AV001cbeamlets'; 
% optSettings.beamletInfo.beamletFolder = 'AV001fbeamlets'; 
optSettings.beamletInfo.beamletFolder = '../../plan/HN003cbeamlets';                         % specify location of beamlet folder relative to opt folder
% optSettings.beamletInfo.beamletFolder = '../../plan/AV004fbeamlets';                         % specify location of beamlet folder relative to opt folder
% optSettings.beamletInfo.beamletFolder = '../../plan/AV002fbeamlets';
% optSettings.beamletInfo.beamletFolder = '../../plan/AV005cbeamlets';

optSettings.optInfo.optExeName = [optSettings.optInfo.doseBatchBaseName 'linlsqOpt.exe'];

optSettings.beamletInfo.Nbeamletbatches = 222;                              % total number of beamlet batches
optSettings.beamletInfo.Nbeamlets = 12592;                                  % total number of non-zero weighted beamets
% optSettings.prescInfo.savePresc = 'no';                                    % save presc info in prescription.txt file located in input folder
optSettings.optInfo.optFolder = 'C:\WiscPlan\opt\linlsqOpt';                % directory where opt executable is located must be absolute path
optSettings.optInfo.modFactor = 6;                                        % mod factor truncating maximum beamlet weight relative to average weight

%   Prescription table
%   (code name)     (Pinnacle tissue name) (alpha)  (Bplus) (dplus)  (Bminus)  (dminus) (BVplus) (dVplus) (Vplus) (BVminus)  (dVminus)  (Vminus)        
prescTable = {      
    'cord'               'cord'              1         1        0        0         0        0        0        0       0         0         0         
    'GTV'                'GTV'               1      1000       70     1000        70        0        0        0       0         0         0        
%     'cont_base'          'cont_base'         1       1000       60      1000        60        0        0        0       0         0         0        
%     'cont_boost'         'cont_boost'        1       2000       80      2000        80        0        0        0       0         0         0        
%     'PVC_cont_base'      'PVC_cont_base'     1      1000       60       1000        60        0        0        0       0         0         0        
%     'PVC_cont_boost'     'PVC_cont_boost'    1      2000       80       2000        80        0        0        0       0         0         0        
%     'thresh_base'        'thresh_base'       1      1000       50       1000        50        0        0        0       0         0         0        
%     'thresh_boost'       'thresh_boost'      1      2000       70       2000        70        0        0        0       0         0         0        
%     'PVC_thresh_base'    'PVC_thresh_base'   1      1000      50       1000        50        0        0        0       0         0         0        
%     'PVC_thresh_boost'   'PVC_thresh_boost'  1      2000      70       2000        70        0        0        0       0         0         0        
%     'CTV60'              'CTV 60'            0         0        0        0         0        0        0        0       0         0         0                
%     'CTV54'              'CTV 54'            0         0        0        0         0        0        0        0       0         0         0          
%     'CTV-rt-scv-5'       'CTV rt scv 54'     0         0        0        0         0        0        0        0       0         0         0                  
%     'CTV-lt-scv-50'      'CTV lt scv 50'     0         0        0        0         0        0        0        0       0         0         0                  
%     'larynx'             'larynx'            0         1        0        0         0        0        0        0       0         0         0                
%     'opticApparatus'     'optic apparatus'   0         1        0        0         0        0        0        0       0         0         0             
    'left_parotid'       'left_parotid'      1         1        0        0         0        0        0        0       0         0         0         
    'normal_tissue'      'normal_tissue'     1        10        0        0         0        0        0        0       0         0         0 
    'oral_cavity'        'oral_cavity'       1         1        0        0         0        0        0        0       0         0         0        
    'right_Parotid'      'right_parotid'     1         1        0        0         0        0        0        0       0         0         0         
    

%     'Brainstem'          'Brainstem'         0         1        0        0         0        0        0        0       0         0         0         
%     'eyes'               'eyes'              0         1        0        0         0        0        0        0       0         0         0        
%     'couch'              'couch'             0         0        0        0         0        0        0        0       0         0         0        
%     'external'           'external'          0         0        0        0         0        0        0        0       0         0         0        
%     'ext+1cm'            'ext + 1 cm'        0         0        0        0         0        0        0        0       0         0         0        
%     'ring'               'ring'              0         1        0        0         0        0        0        0       0         0         0 
%     'PTV70'              'PTV  70'           0       200       70      200        70        0        0        0       0         0         0 
%     'PTV60'              'PTV  60'           0         0        0        0         0        0        0        0       0         0         0 
%     'PTV54'              'PTV 54'            0        10        0        0         0        0        0        0       0         0         0 
%     'PTV50'              'PTV 50'            0        10        0        0         0        0        0        0       0         0         0 
%     'LResidParotid'      'L resid parotid'   0         1        0        0         0        0        0        0       0         0         0    

%     'PTV60min70'         'PTV 60 - PTV 70'   0        10       60      200        60        0        0        0       0         0         0 
    }; 

%%%%% End of user supplied inputs %%%%%

%%

% % calc total number of initial non-zero weighted beamlets from dose calc if
% % saving new beamlet header file
if strcmp(optSettings.beamletInfo.saveBeamletHeader,'yes')
    optSettings.beamletInfo.Nbeamlets = beamletNum([optSettings.optInfo.optFolder '\' optSettings.beamletInfo.beamletFolder],optSettings.beamletInfo.Nbeamletbatches);
end

load(geometryFile);  % load geometry file

% match isocenter of coordinate system used to calculate photon beamlets or
% proton spots, defined by the patient geometry
x_iso = Geometry.start(1) + (Geometry.voxel_size(1).*double(size(Geometry.data,1)))./2;
y_iso = Geometry.start(2) + (Geometry.voxel_size(2).*double(size(Geometry.data,2)))./2;
z_iso = Geometry.start(3);

iso = [x_iso y_iso z_iso];

% size of CT grid
[M,N,Q] = size(Geometry.data);

% Create the prescription structure
presc = [];   % erase the old prescription structure, if it exists
presc.siz = [M N Q];     % dimensions of the prescription and tissue masks

ROIFlags = zeros(1,size(prescTable,1));  % flags to ensure all ROIs are found

% assign the parameters to the prescription structure
for j=1:size(prescTable,1)
    % pull out the prescription indices for each tissue
    tissName = prescTable{j,2};  % current tissue name

    % find the current tissue in the Geometry.ROIS cells
    for k=1:length(Geometry.ROIS)
        ind = 0;
        if strcmp(tissName,Geometry.ROIS{k}.name)
            ind = Geometry.ROIS{k}.ind;  % extract ROI indices
            break;
        end
    end
    if ind == 0 
        error([tissName ' can not be found in Geometry or is empty!']);
    end
    
    presc.tissue(j).name = prescTable{j,1};
    presc.tissue(j).ind = ind;
    presc.tissue(j).alpha = prescTable{j,3};
    presc.tissue(j).betaPlus = prescTable{j,4};
    presc.tissue(j).dPlus = single(zeros(size(ind)) + prescTable{j,5});
    presc.tissue(j).betaMinus = prescTable{j,6};
    presc.tissue(j).dMinus = single(zeros(size(ind)) + prescTable{j,7});
    presc.tissue(j).betaVPlus = prescTable{j,8};
    presc.tissue(j).dVPlus = prescTable{j,9};
    presc.tissue(j).vPlus = prescTable{j,10};
    presc.tissue(j).betaVMinus = prescTable{j,11};
    presc.tissue(j).dVMinus = prescTable{j,12};
    presc.tissue(j).vMinus = prescTable{j,13};
    ROIFlags(j) = 1;
end

% ensure that all of the ROIS were found
if sum(ROIFlags) ~= length(ROIFlags)
    error('Missed one or more ROIS.\n');
else
    fprintf('All ROIs located successfully for downsampled prescription.\n');
end

optSettings.prescInfo.presc = presc;  % set prescription

% Save prescription and beamlet information, which is necessary for initial
% guess calculation.
linlsqOptimization(optSettings);  

% calculate the initial guess if flagged
[w0,dose0] = calcNormalizedInitialGuess([optSettings.optInfo.optFolder '/optInput.txt']);
optSettings.initialGuessInfo.w0 = w0;
optSettings.initialGuessInfo.saveInitialGuess = 'yes';                     
linlsqOptimization(optSettings);

if (optSettings.optInfo.optRun == 'no')
    fprintf('Optimization setup complete.\n');
else
    [weights, dose] = linlsqOptimization(optSettings);
    fprintf('Optimization run complete.\n');
    
    if (optSettings.optInfo.waitForResults == 'yes')
    [optSettings,optResults] = loadOptResults([optSettings.optInfo.optFolder '/optInput.txt'],'last');
    frpintf('Final dose distribution and DVHs loaded.\n');
    end
    
end
clear all;