Skip to content

Vectorized Raycasts #528

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 3 commits into
base: dev
Choose a base branch
from
Draft

Vectorized Raycasts #528

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
57c02e3
vectorized the siddon raytracer and put it in a class
wahln Sep 17, 2021
37e6a1b
fix for single ray tracing
wahln Sep 17, 2021
ee1b8f0
Merge branch 'dev' of https://github.com/e0404/matRad into dev_vector…
wahln Feb 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
139 changes: 97 additions & 42 deletions matRad_rayTracing.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
function [radDepthV, radDepthCube] = matRad_rayTracing(stf,ct,V,rot_coordsV,lateralCutoff)
function [radDepthV, radDepthCube] = matRad_rayTracing(stfElement,ct,V,rot_coordsV,lateralCutoff)
% matRad visualization of two-dimensional dose distributions on ct including
% segmentation
%
%
% call
% [radDepthV, radDepthCube] = matRad_rayTracing(stf,ct,V,rot_coordsV,lateralCutoff)
%
% input
% stf: matRad steering information struct of one beam
% stfElement: matRad steering information struct of single(!) 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
% lateralCutoff: lateral cut off used for ray tracing (optional)

%
% output
Expand All @@ -22,23 +22,77 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 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
% Copyright 2015 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
matRad_cfg = MatRad_Config.instance();
matRad_cfg.dispDeprecationWarning('Calls to matRad_rayTracing will be replaced by usage of the matRad_RayTracer class in the future!');

% At the moment we only implement the Siddon ray-tracer in [1]
hTracer = matRad_RayTracerSiddon();

if nargin >= 5
hTracer.lateralCutOff = lateralCutoff;
end

if nargin >= 4
[radDepthV,radDepthCube] = hTracer.traceCube(stfElement,ct,V,rot_coordsV);
elseif nargin >= 3
[radDepthV,radDepthCube] = hTracer.traceCube(stfElement,ct,V);
else
[radDepthV,radDepthCube] = hTracer.traceCube(stfElement,ct);
end

%{

%If no subset of voxels is specified, take all of them
if nargin < 3
V = transpose(1:prod(ct.cubeDim));
end

%if we don't provide rotated patient coordinates we can compute
%them here on our own
if nargin < 4
[yCoordsV_vox, xCoordsV_vox, zCoordsV_vox] = ind2sub(ct.cubeDim,V);
xCoordsV = xCoordsV_vox(:)*ct.resolution.x-stfElement.isoCenter(1);
yCoordsV = yCoordsV_vox(:)*ct.resolution.y-stfElement.isoCenter(2);
zCoordsV = zCoordsV_vox(:)*ct.resolution.z-stfElement.isoCenter(3);
coordsV = [xCoordsV yCoordsV zCoordsV];

% Get Rotation Matrix
% Do not transpose matrix since we usage of row vectors &
% transformation of the coordinate system need double transpose

rotMat_system_T = matRad_getRotationMatrix(stfElement.gantryAngle,stfElement.couchAngle);

% Rotate coordinates (1st couch around Y axis, 2nd gantry movement)
rot_coordsV = coordsV*rotMat_system_T;

%
rot_coordsV(:,1) = rot_coordsV(:,1)-stfElement.sourcePoint_bev(1);
rot_coordsV(:,2) = rot_coordsV(:,2)-stfElement.sourcePoint_bev(2);
rot_coordsV(:,3) = rot_coordsV(:,3)-stfElement.sourcePoint_bev(3);
end

%If we don't provide a lateral cutoff
if nargin < 5
matRad_cfg = MatRad_Config.instance();
lateralCutoff = matRad_cfg.propDoseCalc.defaultGeometricCutOff;
end

% set up rad depth cube for results
radDepthCube = repmat({NaN*ones(ct.cubeDim)},ct.numOfCtScen);

% set up ray matrix direct behind last voxel
rayMx_bev_y = max(rot_coordsV(:,2)) + max([ct.resolution.x ct.resolution.y ct.resolution.z]);
rayMx_bev_y = rayMx_bev_y + stf.sourcePoint_bev(2);
rayMx_bev_y = rayMx_bev_y + stfElement.sourcePoint_bev(2);

% set up list with bev coordinates for calculation of radiological depth
coords = zeros(prod(ct.cubeDim),3);
Expand All @@ -55,28 +109,28 @@
[candidateRaysCoords_X,candidateRaysCoords_Z] = meshgrid(rayMxSpacing*[floor(-500/rayMxSpacing):ceil(500/rayMxSpacing)]);

% check which rays should be used
for i = 1:stf.numOfRays

ix = (candidateRaysCoords_X(:) - (1+rayMx_bev_y/stf.SAD)*stf.ray(i).rayPos_bev(1)).^2 + ...
(candidateRaysCoords_Z(:) - (1+rayMx_bev_y/stf.SAD)*stf.ray(i).rayPos_bev(3)).^2 ...
<= lateralCutoff^2;
for i = 1:stfElement.numOfRays
ix = (candidateRaysCoords_X(:) - (1+rayMx_bev_y/stfElement.SAD)*stfElement.ray(i).rayPos_bev(1)).^2 + ...
(candidateRaysCoords_Z(:) - (1+rayMx_bev_y/stfElement.SAD)*stfElement.ray(i).rayPos_bev(3)).^2 ...
<= lateralCutoff^2;

candidateRayMx(ix) = 1;

end

% set up ray matrix
rayMx_bev = [candidateRaysCoords_X(logical(candidateRayMx(:))) ...
rayMx_bev_y*ones(sum(candidateRayMx(:)),1) ...
candidateRaysCoords_Z(logical(candidateRayMx(:)))];
rayMx_bev_y*ones(sum(candidateRayMx(:)),1) ...
candidateRaysCoords_Z(logical(candidateRayMx(:)))];

% figure,
% for jj = 1:length(rayMx_bev)
% plot(rayMx_bev(jj,1),rayMx_bev(jj,3),'rx'),hold on
% plot(rayMx_bev(jj,1),rayMx_bev(jj,3),'rx'),hold on
% end

% Rotation matrix. Transposed because of row vectors
rotMat_vectors_T = transpose(matRad_getRotationMatrix(stf.gantryAngle,stf.couchAngle));
rotMat_vectors_T = transpose(matRad_getRotationMatrix(stfElement.gantryAngle,stfElement.couchAngle));

% rotate ray matrix from bev to world coordinates
rayMx_world = rayMx_bev * rotMat_vectors_T;
Expand All @@ -86,43 +140,43 @@

% perform ray tracing over all rays
for i = 1:size(rayMx_world,1)

% run siddon ray tracing algorithm
[~,l,rho,~,ixHitVoxel] = matRad_siddonRayTracer(stf.isoCenter, ...
ct.resolution, ...
stf.sourcePoint, ...
rayMx_world(i,:), ...
ct.cube);

% find voxels for which we should remember this tracing because this is
[~,l,rho,~,ixHitVoxel] = matRad_siddonRayTracer(stfElement.isoCenter, ...
ct.resolution, ...
stfElement.sourcePoint, ...
rayMx_world(i,:), ...
ct.cube);
% find voxels for which we should remember this tracing because this is
% the closest ray by projecting the voxel coordinates to the
% intersection points with the ray matrix and checking if the distance
% intersection points with the ray matrix and checking if the distance
% in x and z direction is smaller than the resolution of the ray matrix
scale_factor = (rayMx_bev_y - stf.sourcePoint_bev(2)) ./ ...
coords(ixHitVoxel,2);

scale_factor = (rayMx_bev_y - stfElement.sourcePoint_bev(2)) ./ ...
coords(ixHitVoxel,2);
x_dist = coords(ixHitVoxel,1).*scale_factor - rayMx_bev(i,1);
z_dist = coords(ixHitVoxel,3).*scale_factor - rayMx_bev(i,3);

ixRememberFromCurrTracing = x_dist > -raySelection & x_dist <= raySelection ...
& z_dist > -raySelection & z_dist <= raySelection;

if any(ixRememberFromCurrTracing) > 0

& z_dist > -raySelection & z_dist <= raySelection;
if any(ixRememberFromCurrTracing)
for j = 1:ct.numOfCtScen
% calc radiological depths

% eq 14
% It multiply voxel intersections with \rho values.
d = l .* rho{j}; %Note. It is not a number "one"; it is the letter "l"

% Calculate accumulated d sum.
dCum = cumsum(d)-d/2;

% write radiological depth for voxel which we want to remember
radDepthCube{j}(ixHitVoxel(ixRememberFromCurrTracing)) = dCum(ixRememberFromCurrTracing);
end
end
end

end

Expand All @@ -131,4 +185,5 @@
radDepthV{i} = radDepthCube{i}(V);
end

%}

13 changes: 13 additions & 0 deletions matRad_siddonRayTracer.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,17 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


matRad_cfg = MatRad_Config.instance();
matRad_cfg.dispDeprecationWarning('Calls to matRad_siddonRayTracer will be replaced by usage of the matRad_RayTracerSiddon class in the future!');

% At the moment we only implement the Siddon ray-tracer in [1]
hTracer = matRad_RayTracerSiddon();

[alphas,l,rho,d12,ix] = hTracer.traceRay(isocenter,resolution,sourcePoint,targetPoint,cubes);

%{

% Add isocenter to source and target point. Because the algorithm does not
% works with negatives values. This put (resolution.x,resolution.y,resolution.z)
% in the center of first voxel
Expand Down Expand Up @@ -217,3 +228,5 @@
for i = 1:numel(cubes)
rho{i} = cubes{i}(ix);
end

%}
Loading