matRad_EUD

0001 classdef matRad_EUD < DoseObjectives.matRad_DoseObjective
0002 % matRad_EUD Implements a penalized equivalent uniform dose objective
0003 %   See matRad_DoseObjective for interface description
0004 %
0005 % References
0006 %   -
0007 %
0008 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009 %
0010 % Copyright 2020 the matRad development team.
0011 %
0012 % This file is part of the matRad project. It is subject to the license
0013 % terms in the LICENSE file found in the top-level directory of this
0014 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0015 % of the matRad project, including this file, may be copied, modified,
0016 % propagated, or distributed except according to the terms contained in the
0017 % LICENSE file.
0018 %
0019 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020     
0021     properties (Constant)
0022         name = 'EUD';
0023         parameterNames = {'EUD^{ref}', 'k'};
0024         parameterTypes = {'dose','numeric'};
0025     end
0026     
0027     properties
0028         parameters = {0, 3.5};
0029         penalty = 1;
0030     end
0031     
0032     methods
0033         function obj = matRad_EUD(penalty,eudRef, eudExponent)
0034             %If we have a struct in first argument
0035             if nargin == 1 && isstruct(penalty)
0036                 inputStruct = penalty;
0037                 initFromStruct = true;
0038             else
0039                 initFromStruct = false;
0040                 inputStruct = [];
0041             end
0042             
0043             %Call Superclass Constructor (for struct initialization)
0044             obj@DoseObjectives.matRad_DoseObjective(inputStruct);
0045             
0046             %now handle initialization from other parameters
0047             if ~initFromStruct
0048                 if nargin >= 3 && isscalar(eudExponent)
0049                     obj.parameters{2} = eudExponent;
0050                 end
0051                 
0052                 if nargin >= 2 && isscalar(eudRef)
0053                     obj.parameters{1} = eudRef;
0054                 end
0055                 
0056                 if nargin >= 1 && isscalar(penalty)
0057                     obj.penalty = penalty;
0058                 end
0059             end
0060         end
0061         
0062         %% Calculates the Objective Function value
0063         function fDose = computeDoseObjectiveFunction(obj,dose)
0064             % get exponent for EUD
0065             k = obj.parameters{2};
0066             
0067             % calculate power sum
0068             powersum = sum(dose.^k);
0069             
0070             
0071             
0072             %Calculate objective
0073             
0074             %This check is not needed since dose is always positive
0075             %if powersum > 0
0076             fDose = obj.penalty * (nthroot(powersum/numel(dose),k) - obj.parameters{1})^2;
0077             %end
0078         end
0079         
0080         %% Calculates the Objective Function gradient
0081         function fDoseGrad  = computeDoseObjectiveGradient(obj,dose)
0082             % get exponent for EUD
0083             k = obj.parameters{2};
0084             
0085             %numerical stability
0086             dose(dose == 0) = 0.001;
0087             
0088             % calculate power sum
0089             powersum = sum(dose.^k);
0090                         
0091             
0092             %This check is not needed since dose is always positive
0093             %if powersum > 0
0094             
0095             %derivatives = nthroot(1/numel(dose),k) * powersum^((1-k)/k) * (dose.^(k-1));
0096             fDoseGrad = 2 * obj.penalty * nthroot(1/numel(dose),k) * powersum^((1-k)/k) * (dose.^(k-1)) .* (nthroot(powersum/numel(dose),k) - obj.parameters{1});
0097             %end
0098             if any(~isfinite(fDoseGrad)) % check for inf and nan for numerical stability
0099                 error(['EUD computation failed. Reduce exponent to resolve numerical problems.']);
0100             end
0101         end
0102     end
0103     
0104 end
0105