Toolbox for robust and nonrobust IMPT optimization in Matlab Bram - - PowerPoint PPT Presentation
Toolbox for robust and nonrobust IMPT optimization in Matlab Bram L. Gorissen Physics research, department of radiation oncology Massachusetts General Hospital and Harvard Medical School September 17, 2015 Bram L. Gorissen Matlab toolbox for
Bram L. Gorissen
Physics research, department of radiation oncology Massachusetts General Hospital and Harvard Medical School
September 17, 2015
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 1 / 48
Past two years: research on robust optimization of IMPT Created a set of harmonious Matlab scripts Useful as a whole for researchers joining the field Individual scripts may be useful for veterans too
Most scripts can be used independently
Downloadable from 3142.nl
Easy to remember: Dutch π
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 2 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 3 / 48
Three variables encode all patient data
DT Dfilter VVStruct
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 4 / 48
The dose influence matrix D is stored as a sparse matrix
smaller by an order of magnitude compared to a dense matrix for a typical IMPT case
Matlab can efficiently access sparse matrices by column That means filtering on beamlets is fast and filtering by voxels is slow Since we need to access by voxel, we store DT Let f be an index vector that indicates which voxels are in an organ d = D(f,:)*x; % slow DT = D'; clear D; d = (x'*DT(:,f))'; % fast
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 5 / 48
For Robust Optimization, DT is a cell array of filenames DT{1} refers to the D matrix of the nominal scenario DT{2}, DT{3}, etc, refer to D matrices of other scenarios
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 6 / 48
To know which voxels are in organ i, we use “Dfilter{i}” Dfilter{i} is a vector whose length is the number of voxels Dfilter{i}(j) = true iff voxel j is in organ i % to compute the dose for voxels in volume i: d = (x'*DT(:,Dfilter{i}))';
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 7 / 48
VVStruct provides information about organs VVStruct.nVois is the number of organs VVStruct.name{i} is the name of organ i
used in getObjectives and getDoseStatistics
VVStruct.size{i} is the number of voxels in organ i
used to compute objective values such as mean dose
VVStruct.dist vector(i) is the distance from voxel i to the target structure
used for the conformity objective
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 8 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 9 / 48
getPtMetadata.m defines meta data for data sets To add new patient, define:
matlabfolder: folder with DijInfo.mat (which contains VVStruct&Dfilter), and files DT* (written by writesparse) dosegrid: voxel grid, used to export to CERR fields, pencilbeam: used to split up dose into dose per field target structure: used for the conformity objective dijfolder, structures: used for importing data from Astroid sample trampfile, num fractions: used to im-/export tramp files (treatment plans)
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 10 / 48
importDfromAstroid.m imports Astroid data from meta.dijfolder:
read vv.m loads structures.vv read dij.m loads *.dij (opt4d format) sparse3.mex creates DT results (DijInfo.mat and DT*.sparse) are written to meta.matlabfolder
Write your own import script for other data sources
There is example import CORT.m for the CORT data set
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 11 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 12 / 48
cerr 3d dose to vector / cerr dose vector to 3d: convert between a dose vector and a 3d dose cube cerr import plan: export x to CERR (full plan or per field) cerr import dose: export d to CERR cerr show plan: shows a plan from the CERR menu cerr create screenshots: creates screenshots of a series of plans, useful for comparisons cerr create plot: creates vector images of plans, useful for papers compute distance to target: calculates the distance between each voxel and a structure, used to create VVStruct.dist vector cerr get Dfilter for structure: gets the voxel indices of a CERR structure
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 13 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 14 / 48
getDoseStatistics / plotDvh / plotDvhCloud
helps to assess a plan
To remove small elements in a matrix or reclaim space: % remove small elements: DT = removeSmallElements( DT, threshold ); % removing part of a matrix: DT(:,not(Dfilter{1})) = 0; % does not free memory DT = removeSmallElements(DT,0); % reclaims space
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 15 / 48
To load DT/VVStruct/Dfilter for pt=1 (the meaning of ’1’ is defined in getPtMetadata): pt=1; initialize; To filter small elements or reclaim space: % remove small elements: DT = removeSmallElements( DT, threshold ); % removing part of a matrix: DT(:,not(Dfilter{1})) = 0; % does not free memory DT = removeSmallElements(DT,0); % reclaims space
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 16 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 17 / 48
Core of the toolbox: optipopt functions Computes function values and derivatives of objective and constraint functions Supports all types of worst case (composite, objectivewise, voxelwise)
distributionally worst case voxelwise expected value voxelwise distributionally worst case
Robust Optimization currently only in objective
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 18 / 48
case 1 i = i+1; objectives{i} = struct( 'voi', ... 'Skin', 'type', 'maxdose', 'bound', 33 ); i = i+1; objectives{i} = struct( 'voi', ... 'CTV', 'type', 'lbmean', 'bound', 30.6 ); “case 1” defines the objectives id
independent of patient id each set of objectives can be used for any patient (as long as structure names are the same)
each line defines an objective or constraint
constraints are also denoted as “objectives{i}”
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 19 / 48
i = i+1; objectives{i} = struct( 'voi', ... 'Skin', 'type', 'maxdose', 'bound', 33 ); three mandatory properties:
voi: name of the structure type: name of the objective/constraint function bound: numerical value for a constraint, ’objective’ for an objective
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 20 / 48
mean underdose / mean overdose
parameter pi is encoded in ’dose’ field:
mean squared underdose / mean squared overdose
parameter pi is encoded in ’dose’ field
uniform
parameter pi is encoded in ’dose’ field
mindose/maxdose
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 21 / 48
minmean/maxmean (for objectives) lbmean/ubmean (for constraints)
why not specify “minmean” with a bound? that would be a bound on a minmean objective confusing: minmean bounds the dose from above
sumfluence maxfluence lower CVAR / upper CVAR
mean dose in hottest/coldest α% α is encoded in ’alpha’ field (range 0-100)
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 22 / 48
conformity
pi is based on distance to target pi = max{0.5prescribed target dose, prescribed target dose − disti}
pi is the prescribed dose for voxels at the boundary of the target pi does not get below 50% of the prescribed dose
uses the struct field ’prescribed target dose’ and VVStruct.dist vector based on all voxels (set the ’voi’ field to an arbitrary structure) useful for creating conformal plans
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 23 / 48
To define a weighted sum objective:
'type', 'maxdose', 'bound', 'objective' );
'type', 'lbmean', 'bound', 'objective' ); The optional field ’weight’ changes the relative weights (default: 1)
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 24 / 48
Weight vector: vector of length four
First element: nominal weight Second element: expected value weight Third element: worst case weight Fourth element: distributionally worst case weight (for nine objectives, puts 5/9 of the total weight on the worst case and weight 1/9 on the second-fifth worst objectives)
Optional field: ’voxelwise’
Calculates expected / worst dose on a per voxel basis
Optional field: ’group’:
(Distributionally) worst case is determined per group
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 25 / 48
'maxdose', 'bound', 'objective', 'weight', [0 1] );
'lbmean', 'bound', 'objective' ... 'weight', [0 0.5]);
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 26 / 48
'type', 'maxdose', 'bound', 'objective', ... 'weight', [0 0 1], 'group', 1 );
'type', 'lbmean', 'bound', 'objective' ... 'weight', [0 0 0.5], 'group', 1 ); The objectives are in the same ’group’, so they are summed before the worst case is computed.
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 27 / 48
'maxdose', 'bound', 'objective', 'weight', [0 0 1] );
'lbmean', 'bound', 'objective', 'weight', [0 0 0.5] ); Omit group property or set it to a different value for each objective.
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 28 / 48
'type', 'maxdose', 'bound', 'objective', ... 'weight', [0 0 1], 'voxelwise', true );
'type', 'lbmean', 'bound', 'objective', ... 'weight', [0 0 0.5], 'voxelwise', true );
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 29 / 48
'type', 'maxdose', 'bound', 'objective', ... 'weight', [0 0 0 1] );
'type', 'lbmean', 'bound', 'objective', ... 'weight', [0 0 0 0.5] );
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 30 / 48
'type', 'maxdose', 'bound', 'objective', ... 'weight', [0 0 1], 'group', 1 );
'type', 'lbmean', 'bound', 'objective', ... 'weight', [0 0 0.5], 'group', 1 );
'type', 'mindose', 'bound', 'objective', ... 'weight', [0 0 0.5], 'group', 2 );
'type', 'upper CVAR', 'bound', 'objective', ... 'weight', [0 0 0 1], 'pervoxel', true );
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 31 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 32 / 48
Extract opt pkg.zip Open Matlab and go to the opt pkg folder Run compile opt pkg.m Add opt pkg to the Matlab path
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 33 / 48
IPOPT
http://www.coin-or.org/download/binary/Ipopt/ Latest version (as of Sept 14, 2015): Ipopt-3.11.8-linux64mac64win32win64-matlabmexfiles.zip The mex file is all you need
CERR
http://www.cerr.info/gettingstarted.php
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 34 / 48
Modify getPtMetadata.m to define patient data Modify getObjectives.m to define objectives To run 1500 iterations with an infinite time limit: pt=1; initialize; [x,info] = optipopt(DT, getObjectives(VVStruct, ... 1),Dfilter,VVStruct,1500,inf);
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 35 / 48
pt=1; initialize; % press ctrl-c to abort loading DT [x,info] = optipopt(getDfiles(pt), getObjectives( ... VVStruct,1),Dfilter,VVStruct,1500,inf);
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 36 / 48
1500 iterations with BFGS take a while, but results in superior solutions
For faster optimization, use a different stopping criterion, or switch to L-BFGS by commenting out options.ipopt.limited memory max history in optipopt.m
Stopping with ctrl-c makes IPOPT crash
Instead, you can create a text file stop.txt that contains a smaller iteration limit (its existence and contents are checked in each iteration)
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 37 / 48
unfixed bug: readsparse/writesparse do not work on some systems; solution: remove readsparse.mex/writesparse.mex files (.m files take
Matlab crashes when optimizing: first check readsparse/writesparse. Otherwise there may be an error in an IPOPT callback function. Uncomment the line that starts with ”%funcs” in optipopt.m to
file line by line.
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 38 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 39 / 48
Requires function values and derivatives Example: mean squared deviation from prescribed dose
f (x) =
i(Dx − p)2 i
∇f (x) = 2D⊤(Dx − p)
Bottleneck is in computing d = Dx and ∇f (x) = 2D⊤y Same bottleneck for many other objectives f (x) Matlab: single threaded sparse matrix multiplication Toolbox: multithreaded sparse matrix multiplication
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 40 / 48
Multithreaded sparse matrix multiplication Calculations are often done for one organ f = Dfilter{i}; % to compute D(f,:)*x: d = sparseTmatrixvectorproduct(DT,x,f); % to compute D(f,:)'*y: derivative = sparsematrixvectorproduct(DT,y,f);
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 41 / 48
Multithreaded sparse matrix multiplication Calculations are often done for one organ, filtering on voxels is more efficient too f = Dfilter{i}; % to compute D(f,:)*x: d = sparseTmatrixvectorproduct(DT,x,f); % to compute D(f,:)'*y: derivative = sparsematrixvectorproduct(DT,y,f);
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 42 / 48
Matlab’s sparse(i,j,s) has excessive memory overhead Our replacement sparse2 and sparse3 have less overhead sparse3 does not use significantly more memory than to hold DT % method 1: D = sparse(i,j,s); DT = D'; % method 2: D = sparse2(i,j,s); DT = D'; % method 3: DT = sparse3('I.txt', 'J.txt', 'S.txt');
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 43 / 48
Matlab cannot store variables >2GB without compression Matlab uses gzip compression Custom save function uses LZ4 compression
saving takes longer, but files are 15-20% smaller, and decompressing is four times faster
% slow save filename DT; load filename DT; % fast writesparse('filename.sparse', DT); DT = readsparse('filename.sparse');
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 44 / 48
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 45 / 48
Optipopt uses optipopt functions The latter script holds all D matrices in memory Parts corresponding to unused voxels are removed to reduce memory consumption
For nominal and robust scenarios separately
One process has DT in memory and performs multiplications Other processes compute objective values and run IPOPT
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 46 / 48
Virtual beamlets (Unkelbach Med Phys 36, 149, 2009) can be used to reduce memory consumption Currently not supported Relatively simple to hack into the optimization script
Open optipopt functions.m Search for ”DT{” and change:
Lines 18-35: set numBeamlets/numScenarios manually, load single DT Lines 126-143: remove these lines Lines 818-871: implement virtual beamlet logic here
This enables optimization
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 47 / 48
To use virtual beamlets throughout the toolbox:
Create sparse(T)matrixvectorproduct virtual beamlets(DT,x,f,curS) Replace all calls to sparse(T)matrixvectorproduct Fix remaining bugs :-)
Bram L. Gorissen Matlab toolbox for IMPT optimization September 17, 2015 48 / 48