matRad_constraintJacobian

 matRad IPOPT callback: jacobian function for direct aperture optimization
 
 call
   jacob = matRad_daoJacobFunc(optiProb,apertureInfoVec,dij,cst)

 input
   optiProb:        option struct defining the type of optimization
   apertureInfoVec: aperture info vector
   dij:             dose influence matrix
   cst:             matRad cst struct

 output
   jacob:           jacobian of constraint function

 References

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function jacob = matRad_constraintJacobian(optiProb,apertureInfoVec,dij,cst)
0002 % matRad IPOPT callback: jacobian function for direct aperture optimization
0003 %
0004 % call
0005 %   jacob = matRad_daoJacobFunc(optiProb,apertureInfoVec,dij,cst)
0006 %
0007 % input
0008 %   optiProb:        option struct defining the type of optimization
0009 %   apertureInfoVec: aperture info vector
0010 %   dij:             dose influence matrix
0011 %   cst:             matRad cst struct
0012 %
0013 % output
0014 %   jacob:           jacobian of constraint function
0015 %
0016 % References
0017 %
0018 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019 
0020 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 %
0022 % Copyright 2015 the matRad development team.
0023 %
0024 % This file is part of the matRad project. It is subject to the license
0025 % terms in the LICENSE file found in the top-level directory of this
0026 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0027 % of the matRad project, including this file, may be copied, modified,
0028 % propagated, or distributed except according to the terms contained in the
0029 % LICENSE file.
0030 %
0031 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 apertureInfo = optiProb.apertureInfo;
0034 
0035 % update apertureInfo if necessary
0036 if ~isequal(apertureInfoVec,apertureInfo.apertureVector)
0037     apertureInfo = optiProb.matRad_daoVec2ApertureInfo(apertureInfo,apertureInfoVec);
0038 end
0039 
0040 % jacobian of the dao constraints
0041 
0042 % row indices
0043 i = repmat(1:apertureInfo.totalNumOfLeafPairs,1,2);
0044 % column indices
0045 j = [apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs ...
0046      apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs+1:apertureInfo.totalNumOfShapes+2*apertureInfo.totalNumOfLeafPairs];
0047 
0048 % -1 for left leaves, 1 for right leaves
0049 s = [-1*ones(1,apertureInfo.totalNumOfLeafPairs) ones(1,apertureInfo.totalNumOfLeafPairs)];
0050 
0051 jacob_dao = sparse(i,j,s, ...
0052     apertureInfo.totalNumOfLeafPairs, ...
0053     apertureInfo.totalNumOfShapes+2*apertureInfo.totalNumOfLeafPairs, ...
0054     2*apertureInfo.totalNumOfLeafPairs);
0055 
0056 % compute jacobian of dosimetric constrainst
0057 
0058 % dosimetric jacobian in bixel space
0059 jacob_dos_bixel = matRad_constraintJacobian@matRad_OptimizationProblem(optiProb,apertureInfo.bixelWeights,dij,cst);
0060 
0061 % allocate sparse matrix for dosimetric jacobian
0062 jacob_dos = sparse(size(jacob_dos_bixel,1),numel(apertureInfoVec));
0063 
0064 if ~isempty(jacob_dos)
0065     
0066     % 1. calculate jacobian for aperture weights
0067     % loop over all beams
0068     offset = 0;
0069     for i = 1:numel(apertureInfo.beam)
0070 
0071         % get used bixels in beam
0072         ix = ~isnan(apertureInfo.beam(i).bixelIndMap);
0073 
0074         % loop over all shapes and add up the gradients x openingFrac for this shape
0075         for j = 1:apertureInfo.beam(i).numOfShapes            
0076             jacob_dos(:,offset+j) = jacob_dos_bixel(:,apertureInfo.beam(i).bixelIndMap(ix)) ...
0077                                       * apertureInfo.beam(i).shape(j).shapeMap(ix);
0078         end
0079 
0080         % increment offset
0081         offset = offset + apertureInfo.beam(i).numOfShapes;
0082 
0083     end
0084 
0085     % 2. find corresponding bixel to the leaf Positions and aperture
0086     % weights to calculate the jacobian
0087     jacob_dos(:,apertureInfo.totalNumOfShapes+1:end) = ...
0088             ( ones(size(jacob_dos,1),1) * apertureInfoVec(apertureInfo.mappingMx(apertureInfo.totalNumOfShapes+1:end,2))' ) ...
0089          .* jacob_dos_bixel(:,apertureInfo.bixelIndices(apertureInfo.totalNumOfShapes+1:end)) / apertureInfo.bixelWidth;
0090 
0091     % correct the sign for the left leaf positions
0092     jacob_dos(:,apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs) = ...
0093    -jacob_dos(:,apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs);
0094 
0095 end
0096 
0097 % concatenate
0098 jacob = [jacob_dao;jacob_dos];