IMRT optimization interfacing with an external solver - cerr/CERR GitHub Wiki

IMRTP: IMRT optimization interfacing with an external solver

A user may desire to use the IMRTP components of CERR to do IMRT optimization (beamlet/fluence map optimization). Note that this is currently not included with the IMRTP CERR package and requires writing or obtaining optimization code. However, if a user has optimization code (or would like to write it), it can be very valuable to use CERR to calculate inputs to the optimization as well as recalculate dose and evaluate the result.

In this section, we describe the pertinent steps to interface with CERR to do IMRT optimization using an external solver (using the MOSEK optimization toolbox for MATLAB as an example). These steps are:

• Get optimization inputs
• Interface with external solver
• Get optimization outputs

Before we begin, let’s clarify the inputs and outputs of the optimization.

The matrix of the unit-weight beamlet dose distributions is referred to as the ‘influence matrix,’ and is denoted A. Dose values (D_i) are computed by addition of all beamlets weighted by the solution vector w:

D_i(w) = A_i*w  (1)

where the subscript refers to the i’th dose element and the i’th row in A.

Inputs: Influence matrix (A)

Outputs: Beamlet weight vector (w), used to calculate dose (D)

The optimization inputs can be calculated by using the IMRTP GUI. Once the beamlets have been calculated, the inputs can be found inside of the planC data structure as follows. Please note that you must first load planC (and/or have the plan loaded in CERR) and declare

global planC
indexS = planC{end};

before using the following code. For further information on planC, see structure of planC

Influence matrix (The A matrix in eq (1). The size of this matrix is the number of voxels by the number of pencil beams.) For a given IMNumber, potentially defined as follows,

IMNumber = size(planC{indexS.IM},2); %to use the last IM structure created

and for structureNumber which corresponds to the desired structure number in planC, the influence matrix can be found in this way:

influenceM = getSingleGlobalInfluenceM(planC{indexS.IM}(IMNumber).IMDosimetry, structureNumber);

There are several options for solvers to be used for IMRT optimization. For example, the following are viable options:

• MATLAB optimization toolbox
• External solver integrated with MATLAB commands (i.e. MOSEK optimization toolbox for MATLAB)
• External solver which requires input files in a format other than MATLAB code (i.e. MPS, AMPL format as input to a wide variety of solvers such as CPLEX, MINOS, MOSEK)

The first two options involve calling the optimization code directly from MATLAB m-files, while the last option allows for much greater flexibility of solvers but requires the conversion of data from MATLAB to .txt format. In this section, we focus on using MOSEK integrated with MATLAB and allow the reader to extrapolate for his/her particular use.

For further assistance or explanation of anything in this section, please see the MOSEK optimization toolbox for MATLAB manual. Downloads and licenses for MOSEK are available at http://www.mosek.com.

Set up inputs to optimizer

For the MOSEK optimization toolbox for MATLAB, the input to the solver is a struct with different variables included inside it. For example, following is an empty problem input. Please note that this code "as is" will not be useful as input to the solver; it is simply offered as a starting point for one wishes to create their own input based on an optimization formulation.

%create empty problem setup for MOSEK
%constraints
prob.blc = [];
prob.buc = [];
prob.blx = zeros(numPBs,1); %where numPBs is the number of pencil beams
prob.bux = upperBound*ones(numPBs,1);
%objectives
linOptV  = zeros(numPBs,1);
prob.c = linOptV; %linear objective
hessianM = sparse(numPBs,numPBs);
[prob.qosubi,prob.qosubj,prob.qoval] = find(tril(hessianM)); %quadratic objective

Run optimizer

The MOSEK optimizer is run in MATLAB by using the 'mosekopt' command, which is the last line of the following text box. The lines before it are optional parameters that do not need to be included but may be useful to be aware of.

%Optional parameters for the optimization.  There are many more, as outlined here.

param.MSK_IPAR_INTPNT_MAX_ITERATIONS = 1000;
param.MSK_DPAR_INTPNT_TOL_REL_GAP=1e-4; % termination criteria
param.MSK_IPAR_LICENSE_WAIT = 'MSK_ON';
param.MSK_IPAR_LICENSE_PAUSE_TIME = 5000;

%run the actual optimization
[r,res] = mosekopt('minimize info',prob,param,callback); % r is the return code

Get beamlet weight vector

Get output from optimizer. For MOSEK:

wV = res.sol.itr.xx; 

or, depending on the optimization type (see [MOSEK documentation]):

wV = res.sol.bas.xx;

Calculate dose

Calculating the dose according to eq (1),

dose3DM = influenceM * wV;

Display dose

If the CERR viewer is open, dose can be displayed as follows:

showIMDose(dose3DM,'Test dose name');

If the CERR viewer is closed (or open), dose can be added to planC with this code:

IM = planC{indexS.IM}(IMNumber).IMDosimetry;
dose2CERR(dose3DM,[],'Test dose name','CERR test','Test PB distribution.','UniformCT',[],'no', IM.assocScanUID);

Save results into planC

The resulting output can be saved anywhere. Internally in our group, we save the solutions here (for a given step parameter):

% save dose in planC IM structure 
planC{indexS.IM}(IMNumber).IMDosimetry.solutions(step).doseArray = sparse(dose3DM(:));
% save wV into planC
planC{indexS.IM}(IMNumber).IMDosimetry.solutions(step).beamletWeights = wV;

Save planC

Make sure to save planC at the end if the results are useful to keep.

save nameOfFileToSaveTo.mat planC

This concludes the basics of interfacing IMRTP and CERR with an external solver for IMRT optimization.