matRad dij sampling function
This function samples.
call
[ixNew,bixelDoseNew] =
matRad_DijSampling(ix,bixelDose,radDepthV,rad_distancesSq,sType,Param)
input
ix: indices of voxels where we want to compute dose influence data
bixelDose: dose at specified locations as linear vector
radDepthV: radiological depth vector
rad_distancesSq: squared radial distance to the central ray
sType: can either be set to 'radius' or 'dose'. These are two different ways
to determine dose values that are keept as they are and dose values used for sampling
Param: In the case of radius based sampling, dose values having a radial
distance below r0 [mm] are keept anyway and sampling is only done beyond r0.
In the case of dose based sampling, dose values having a relative dose greater
the threshold [0...1] are keept and sampling is done for dose values below the relative threshold
output
ixNew: reduced indices of voxels where we want to compute dose influence data
bixelDoseNew reduced dose at specified locations as linear vector
References
[1] http://dx.doi.org/10.1118/1.1469633
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Copyright 2016 the matRad development team.
This file is part of the matRad project. It is subject to the license
terms in the LICENSE file found in the top-level directory of this
distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
of the matRad project, including this file, may be copied, modified,
propagated, or distributed except according to the terms contained in the
LICENSE file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0001 function [ixNew,bixelDoseNew] = matRad_DijSampling(ix,bixelDose,radDepthV,rad_distancesSq,sType,Param) 0002 % matRad dij sampling function 0003 % This function samples. 0004 % 0005 % call 0006 % [ixNew,bixelDoseNew] = 0007 % matRad_DijSampling(ix,bixelDose,radDepthV,rad_distancesSq,sType,Param) 0008 % 0009 % input 0010 % ix: indices of voxels where we want to compute dose influence data 0011 % bixelDose: dose at specified locations as linear vector 0012 % radDepthV: radiological depth vector 0013 % rad_distancesSq: squared radial distance to the central ray 0014 % sType: can either be set to 'radius' or 'dose'. These are two different ways 0015 % to determine dose values that are keept as they are and dose values used for sampling 0016 % Param: In the case of radius based sampling, dose values having a radial 0017 % distance below r0 [mm] are keept anyway and sampling is only done beyond r0. 0018 % In the case of dose based sampling, dose values having a relative dose greater 0019 % the threshold [0...1] are keept and sampling is done for dose values below the relative threshold 0020 % 0021 % output 0022 % ixNew: reduced indices of voxels where we want to compute dose influence data 0023 % bixelDoseNew reduced dose at specified locations as linear vector 0024 % 0025 % References 0026 % [1] http://dx.doi.org/10.1118/1.1469633 0027 % 0028 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0029 % 0030 % Copyright 2016 the matRad development team. 0031 % 0032 % This file is part of the matRad project. It is subject to the license 0033 % terms in the LICENSE file found in the top-level directory of this 0034 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 0035 % of the matRad project, including this file, may be copied, modified, 0036 % propagated, or distributed except according to the terms contained in the 0037 % LICENSE file. 0038 % 0039 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0040 0041 %% define default parameters as a fallback 0042 defaultType = 'radius'; 0043 deltaRadDepth = 5; % step size of radiological depth 0044 defaultLatCutOff = 25; % default lateral cut off 0045 defaultrelDoseThreshold = 0.01; % default relative dose threshold 0046 0047 relDoseThreshold = defaultrelDoseThreshold; 0048 LatCutOff = defaultLatCutOff; 0049 Type = sType; 0050 0051 % if the input index vector is of type logical convert it to linear indices 0052 if islogical(ix) 0053 ix = find(ix); 0054 end 0055 0056 %% parse inputs 0057 if sum(strcmp(sType,{'radius','dose'})) == 0 0058 Type = defaultType; 0059 end 0060 0061 % if an parameter is provided then use it 0062 if nargin>5 0063 if exist('Param','var') 0064 if strcmp(sType,'radius') 0065 LatCutOff = Param; 0066 elseif strcmp(sType,'dose') 0067 relDoseThreshold = Param; 0068 end 0069 end 0070 end 0071 0072 %% remember dose values inside the inner core 0073 switch Type 0074 case {'radius'} 0075 ixCore = rad_distancesSq < LatCutOff^2; % get voxels indices having a smaller radial distance than r0 0076 case {'dose'} 0077 ixCore = bixelDose > relDoseThreshold * max(bixelDose); % get voxels indices having a greater dose than the thresholdDose 0078 end 0079 0080 bixelDoseCore = bixelDose(ixCore); % save dose values that are not affected by sampling 0081 0082 if all(ixCore) 0083 %% all bixels are in the core 0084 %exit function with core dose only 0085 ixNew = ix; 0086 bixelDoseNew = bixelDoseCore; 0087 else 0088 logIxTail = ~ixCore; % get voxels indices beyond r0 0089 linIxTail = find(logIxTail); % convert logical index to linear index 0090 numTail = numel(linIxTail); 0091 bixelDoseTail = bixelDose(linIxTail); % dose values that are going to be reduced by sampling 0092 ixTail = ix(linIxTail); % indices that are going to be reduced by sampling 0093 0094 %% sample for each radiological depth the lateral halo dose 0095 radDepthTail = (radDepthV(linIxTail)); % get radiological depth in the tail 0096 0097 % cluster radiological dephts to reduce computations 0098 B_r = int32(ceil(radDepthTail)); % cluster radiological depths; 0099 maxRadDepth = double(max(B_r)); 0100 C = int32(linspace(0,maxRadDepth,round(maxRadDepth)/deltaRadDepth)); % coarse clustering of rad depths 0101 0102 ixNew = zeros(numTail,1); % inizialize new index vector 0103 bixelDoseNew = zeros(numTail,1); % inizialize new dose vector 0104 linIx = int32(1:1:numTail)'; 0105 IxCnt = 1; 0106 0107 %% loop over clustered radiological depths 0108 for i = 1:numel(C)-1 0109 ixTmp = linIx(B_r >= C(i) & B_r < C(i+1)); % extracting sub indices 0110 if isempty(ixTmp) 0111 continue 0112 end 0113 subDose = bixelDoseTail(ixTmp); % get tail dose in current cluster 0114 subIx = ixTail(ixTmp); % get indices in current cluster 0115 thresholdDose = max(subDose); 0116 r = rand(numel(subDose),1); % get random samples 0117 ixSamp = r<=(subDose/thresholdDose); 0118 NumSamples = sum(ixSamp); 0119 0120 ixNew(IxCnt:IxCnt+NumSamples-1,1) = subIx(ixSamp); % save new indices 0121 bixelDoseNew(IxCnt:IxCnt+NumSamples-1,1) = thresholdDose; % set the dose 0122 IxCnt = IxCnt + NumSamples; 0123 end 0124 0125 0126 % cut new vectors and add inner core values 0127 ixNew = [ix(ixCore); ixNew(1:IxCnt-1)]; 0128 bixelDoseNew = [bixelDoseCore; bixelDoseNew(1:IxCnt-1)]; 0129 end 0130 0131 0132 0133 end 0134 0135