matRad visualization of two-dimensional dose distributions on ct including segmentation call radDepthV = matRad_rayTracing(stf,ct,V,rot_coordsV,lateralCutoff) input stf: matRad steering information struct of one beam ct: ct cube V: linear voxel indices e.g. of voxels inside patient. rot_coordsV coordinates in beams eye view inside the patient lateralCutoff: lateral cut off used for ray tracing
0001 function radDepthV = matRad_rayTracing(stf,ct,V,rot_coordsV,lateralCutoff) 0002 % matRad visualization of two-dimensional dose distributions on ct 0003 % including segmentation 0004 % 0005 % call 0006 % radDepthV = matRad_rayTracing(stf,ct,V,rot_coordsV,lateralCutoff) 0007 % 0008 % input 0009 % stf: matRad steering information struct of one beam 0010 % ct: ct cube 0011 % V: linear voxel indices e.g. of voxels inside patient. 0012 % rot_coordsV coordinates in beams eye view inside the patient 0013 % lateralCutoff: lateral cut off used for ray tracing 0014 0015 % 0016 % output 0017 % radDepthV: radiological depth inside the patient 0018 % 0019 % References 0020 % [1] http://www.sciencedirect.com/science/article/pii/S1120179711001359 0021 % 0022 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 % 0024 % Copyright 2015 the matRad development team. 0025 % 0026 % This file is part of the matRad project. It is subject to the license 0027 % terms in the LICENSE file found in the top-level directory of this 0028 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 0029 % of the matRad project, including this file, may be copied, modified, 0030 % propagated, or distributed except according to the terms contained in the 0031 % LICENSE file. 0032 % 0033 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0034 0035 % set up rad depth cube for results 0036 radDepthCube = repmat({NaN*ones(ct.cubeDim)},ct.numOfCtScen); 0037 0038 % set up ray matrix direct behind last voxel 0039 rayMx_bev_y = max(rot_coordsV(:,2)) + max([ct.resolution.x ct.resolution.y ct.resolution.z]); 0040 rayMx_bev_y = rayMx_bev_y + stf.sourcePoint_bev(2); 0041 0042 % set up list with bev coordinates for calculation of radiological depth 0043 coords = zeros(prod(ct.cubeDim),3); 0044 coords(V,:) = rot_coordsV; 0045 0046 % calculate spacing of rays on ray matrix 0047 rayMxSpacing = 1/sqrt(2) * min([ct.resolution.x ct.resolution.y ct.resolution.z]); 0048 0049 % define candidate ray matrix covering 1000x1000mm^2 0050 numOfCandidateRays = 2 * ceil(500/rayMxSpacing) + 1; 0051 candidateRayMx = zeros(numOfCandidateRays); 0052 0053 % define coordinates 0054 [candidateRaysCoords_X,candidateRaysCoords_Z] = meshgrid(rayMxSpacing*[floor(-500/rayMxSpacing):ceil(500/rayMxSpacing)]); 0055 0056 % check which rays should be used 0057 for i = 1:stf.numOfRays 0058 0059 ix = (candidateRaysCoords_X(:) - (1+rayMx_bev_y/stf.SAD)*stf.ray(i).rayPos_bev(1)).^2 + ... 0060 (candidateRaysCoords_Z(:) - (1+rayMx_bev_y/stf.SAD)*stf.ray(i).rayPos_bev(3)).^2 ... 0061 <= lateralCutoff^2; 0062 0063 candidateRayMx(ix) = 1; 0064 0065 end 0066 0067 % set up ray matrix 0068 rayMx_bev = [candidateRaysCoords_X(logical(candidateRayMx(:))) ... 0069 rayMx_bev_y*ones(sum(candidateRayMx(:)),1) ... 0070 candidateRaysCoords_Z(logical(candidateRayMx(:)))]; 0071 0072 % figure, 0073 % for jj = 1:length(rayMx_bev) 0074 % plot(rayMx_bev(jj,1),rayMx_bev(jj,3),'rx'),hold on 0075 % end 0076 0077 % Rotation matrix. Transposed because of row vectors 0078 rotMat_vectors_T = transpose(matRad_getRotationMatrix(stf.gantryAngle,stf.couchAngle)); 0079 0080 % rotate ray matrix from bev to world coordinates 0081 rayMx_world = rayMx_bev * rotMat_vectors_T; 0082 0083 % criterium for ray selection 0084 raySelection = rayMxSpacing/2; 0085 0086 % perform ray tracing over all rays 0087 for i = 1:size(rayMx_world,1) 0088 0089 % run siddon ray tracing algorithm 0090 [~,l,rho,~,ixHitVoxel] = matRad_siddonRayTracer(stf.isoCenter, ... 0091 ct.resolution, ... 0092 stf.sourcePoint, ... 0093 rayMx_world(i,:), ... 0094 ct.cube); 0095 0096 % find voxels for which we should remember this tracing because this is 0097 % the closest ray by projecting the voxel coordinates to the 0098 % intersection points with the ray matrix and checking if the distance 0099 % in x and z direction is smaller than the resolution of the ray matrix 0100 scale_factor = (rayMx_bev_y - stf.sourcePoint_bev(2)) ./ ... 0101 coords(ixHitVoxel,2); 0102 0103 x_dist = coords(ixHitVoxel,1).*scale_factor - rayMx_bev(i,1); 0104 z_dist = coords(ixHitVoxel,3).*scale_factor - rayMx_bev(i,3); 0105 0106 ixRememberFromCurrTracing = x_dist > -raySelection & x_dist <= raySelection ... 0107 & z_dist > -raySelection & z_dist <= raySelection; 0108 0109 if any(ixRememberFromCurrTracing) > 0 0110 0111 for j = 1:ct.numOfCtScen 0112 % calc radiological depths 0113 0114 % eq 14 0115 % It multiply voxel intersections with \rho values. 0116 d = l .* rho{j}; %Note. It is not a number "one"; it is the letter "l" 0117 0118 % Calculate accumulated d sum. 0119 dCum = cumsum(d)-d/2; 0120 0121 % write radiological depth for voxel which we want to remember 0122 radDepthCube{j}(ixHitVoxel(ixRememberFromCurrTracing)) = dCum(ixRememberFromCurrTracing); 0123 end 0124 end 0125 0126 end 0127 0128 % only take voxel inside the patient 0129 for i = 1:ct.numOfCtScen 0130 radDepthV{i} = radDepthCube{i}(V); 0131 end 0132