%% 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;