matRad_constraintJacobian

 matRad IPOPT callback: jacobian function for inverse planning 
 supporting max dose constraint, min dose constraint, min mean dose constraint, 
 max mean dose constraint, min EUD constraint, max EUD constraint, max DVH 
 constraint, min DVH constraint 
 
 call
   jacob = matRad_jacobFunc(optiProb,w,dij,cst)

 input
   optiProb: option struct defining the type of optimization
   w:    bixel weight vector
   dij:  dose influence matrix
   cst:  matRad cst struct

 output
   jacob: jacobian of constraint function

 References

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

0001 function jacob = matRad_constraintJacobian(optiProb,w,dij,cst)
0002 % matRad IPOPT callback: jacobian function for inverse planning
0003 % supporting max dose constraint, min dose constraint, min mean dose constraint,
0004 % max mean dose constraint, min EUD constraint, max EUD constraint, max DVH
0005 % constraint, min DVH constraint
0006 %
0007 % call
0008 %   jacob = matRad_jacobFunc(optiProb,w,dij,cst)
0009 %
0010 % input
0011 %   optiProb: option struct defining the type of optimization
0012 %   w:    bixel weight vector
0013 %   dij:  dose influence matrix
0014 %   cst:  matRad cst struct
0015 %
0016 % output
0017 %   jacob: jacobian of constraint function
0018 %
0019 % References
0020 %
0021 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 
0023 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 %
0025 % Copyright 2016 the matRad development team.
0026 %
0027 % This file is part of the matRad project. It is subject to the license
0028 % terms in the LICENSE file found in the top-level directory of this
0029 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0030 % of the matRad project, including this file, may be copied, modified,
0031 % propagated, or distributed except according to the terms contained in the
0032 % LICENSE file.
0033 %
0034 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0035 
0036 % get current dose / effect / RBExDose vector
0037 %d = matRad_backProjection(w,dij,optiProb);
0038 %d = optiProb.matRad_backProjection(w,dij);
0039 optiProb.BP = optiProb.BP.compute(dij,w);
0040 d = optiProb.BP.GetResult();
0041 
0042 % initialize jacobian (only single scenario supported in optimization)
0043 jacob = sparse([]);
0044 
0045 % initialize projection matrices and id containers
0046 DoseProjection{1}          = sparse([]);
0047 mAlphaDoseProjection{1}    = sparse([]);
0048 mSqrtBetaDoseProjection{1} = sparse([]);
0049 voxelID                     = [];
0050 constraintID                = [];
0051 scenID2                     = [];
0052 
0053 % compute objective function for every VOI.
0054 for i = 1:size(cst,1)
0055     
0056     % Only take OAR or target VOI.
0057     if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') )
0058         
0059         
0060         % loop over the number of constraints for the current VOI
0061         for j = 1:numel(cst{i,6})
0062             
0063             obj = cst{i,6}{j}; %Get the Optimization Object
0064             
0065             % only perform computations for constraints
0066             %if ~isempty(strfind(obj.type,'constraint'))
0067             if isa(obj,'DoseConstraints.matRad_DoseConstraint')
0068                 
0069                 % compute reference
0070                 if (~isequal(obj.name, 'max dose constraint')      && ~isequal(obj.name, 'min dose constraint')      &&...
0071                         ~isequal(obj.name, 'max mean dose constraint') && ~isequal(obj.name, 'min mean dose constraint') && ...
0072                         ~isequal(obj.name, 'min EUD constraint')       && ~isequal(obj.name, 'max EUD constraint'))      && ...
0073                         (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection'))
0074                     
0075                     doses = obj.getDoseParameters();
0076                     
0077                     effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2;
0078                     
0079                     obj = obj.setDoseParameters(effect);
0080                 end
0081                 
0082                 % if conventional opt: just add constraints of nominal dose
0083                 %if strcmp(cst{i,6}{j}.robustness,'none')
0084                 
0085                 d_i = d{1}(cst{i,4}{1});
0086                 
0087                 %jacobVec =  matRad_jacobFunc(d_i,cst{i,6}{j},d_ref);
0088                 jacobSub = obj.computeDoseConstraintJacobian(d_i);
0089                 
0090                 nConst = size(jacobSub,2);
0091                 
0092                 %tic;
0093                 
0094                 %Iterate through columns of the sub-jacobian
0095                 %TODO: Maybe this could all be function of the projection
0096                 %Objects???
0097                 if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobSub) || isa(optiProb.BP,'matRad_ConstantRBEProjection')
0098                     
0099                     startIx = size(DoseProjection{1},2) + 1;
0100                     %First append the Projection matrix with sparse zeros
0101                     DoseProjection{1}          = [DoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)];
0102                     
0103                     %Now directly write the jacobian in there
0104                     DoseProjection{1}(cst{i,4}{1},startIx:end) = jacobSub;
0105                     
0106                 elseif isa(optiProb.BP,'matRad_EffectProjection') && ~isempty(jacobSub)
0107                     
0108                     if isa(optiProb.BP,'matRad_VariableRBEProjection')
0109                         scaledEffect = (dij.gamma(cst{i,4}{1}) + d_i);
0110                         jacobSub = jacobSub./(2*dij.bx(cst{i,4}{1}) .* scaledEffect);
0111                     end
0112                     
0113                     startIx = size(mAlphaDoseProjection{1},2) + 1;
0114                     
0115                     %First append the alphaDose matrix with sparse
0116                     %zeros then insert
0117                     mAlphaDoseProjection{1}    = [mAlphaDoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)];
0118                     mAlphaDoseProjection{1}(cst{i,4}{1},startIx:end) = jacobSub;
0119                     
0120                     %The betadose has a different structure due to the
0121                     %quadratic transformation, but in principle the
0122                     %same as above
0123                     mSqrtBetaDoseProjection{1} =  [mSqrtBetaDoseProjection{1}, sparse(repmat(cst{i,4}{1},nConst,1),repmat(1:numel(cst{i,4}{1}),1,nConst),2*reshape(jacobSub',[],1),dij.doseGrid.numOfVoxels,nConst*numel(cst{i,4}{1}))];
0124                     
0125                     if isempty(constraintID)
0126                         newID = 1;
0127                     else
0128                         newID = constraintID(end)+1;
0129                     end
0130                     
0131                     voxelID = [voxelID;repmat(cst{i,4}{1},nConst,1)];                         %Keep track of voxels for organizing the sqrt(beta)Dose projection later
0132                     constraintID = [constraintID, ...
0133                         reshape(ones(numel(cst{i,4}{1}),1)*[newID:newID+nConst-1],[1 nConst*numel(cst{i,4}{1})])];  %Keep track of constraints for organizing the sqrt(beta)Dose projection later
0134                     
0135                 end
0136                 
0137                       
0138                 %Old implementation with for loop
0139                 %{
0140                 for c = 1:size(jacobSub,2)
0141                     jacobVec = jacobSub(:,c);
0142                     
0143                     if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobVec) || isa(optiProb.BP,'matRad_ConstantRBEProjection')
0144                         
0145                         DoseProjection{1}          = [DoseProjection{1},sparse(cst{i,4}{1},1,jacobVec,dij.doseGrid.numOfVoxels,1)];
0146                         
0147                     elseif isa(optiProb.BP,'matRad_EffectProjection') && ~isempty(jacobVec)
0148                         
0149                         if isa(optiProb.BP,'matRad_VariableRBEProjection')
0150                             scaledEffect = (dij.gamma(cst{i,4}{1}) + d_i);
0151                             jacobVec = jacobVec./(2*dij.bx(cst{i,4}{1}).*scaledEffect);
0152                         end
0153                         
0154                         mAlphaDoseProjection{1}    = [mAlphaDoseProjection{1},sparse(cst{i,4}{1},1,jacobVec,dij.doseGrid.numOfVoxels,1)];
0155                         mSqrtBetaDoseProjection{1} = [mSqrtBetaDoseProjection{1},...
0156                             sparse(cst{i,4}{1},1:numel(cst{i,4}{1}),2*jacobVec,dij.doseGrid.numOfVoxels,numel(cst{i,4}{1}))];
0157                         
0158                         voxelID                 = [voxelID ;cst{i,4}{1}];   %list of voxels relevant for constraints to enable faster computations
0159                         
0160                         if isempty(constraintID)
0161                             lastID = 0;
0162                         else
0163                             lastID = constraintID(end);
0164                         end
0165                         constraintID            = [constraintID, repmat(1 + lastID,1,numel(cst{i,4}{1}))]; %Maps constraints to voxels
0166                         
0167                     end
0168                 end
0169                 %}
0170             end
0171             
0172         end
0173         
0174     end
0175     
0176 end
0177 
0178 % enter if statement also for protons using a constant RBE
0179 if isa(optiProb.BP,'matRad_DoseProjection')
0180     
0181     if ~isempty(DoseProjection{1})
0182         jacob = DoseProjection{1}' * dij.physicalDose{1};
0183     end
0184     
0185 elseif isa(optiProb.BP,'matRad_ConstantRBEProjection')
0186     
0187     if ~isempty(DoseProjection{1})
0188         jacob = DoseProjection{1}' * dij.RBE * dij.physicalDose{1};
0189     end
0190     
0191 elseif isa(optiProb.BP,'matRad_EffectProjection')
0192     
0193     if ~isempty(mSqrtBetaDoseProjection{1}) && ~isempty(mAlphaDoseProjection{1})
0194         mSqrtBetaDoseProjection{1} = mSqrtBetaDoseProjection{1}' * dij.mSqrtBetaDose{1} * w;
0195         mSqrtBetaDoseProjection{1} = sparse(voxelID,constraintID,mSqrtBetaDoseProjection{1},...
0196             size(mAlphaDoseProjection{1},1),size(mAlphaDoseProjection{1},2));
0197         
0198         jacob   = mAlphaDoseProjection{1}' * dij.mAlphaDose{1} +...
0199             mSqrtBetaDoseProjection{1}' * dij.mSqrtBetaDose{1};
0200         
0201     end
0202 end
0203 end