diff --git a/matRad_rayTracing.m b/matRad_rayTracing.m index cd841017b..0c5cf578d 100644 --- a/matRad_rayTracing.m +++ b/matRad_rayTracing.m @@ -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 @@ -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); @@ -55,11 +109,11 @@ [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; @@ -67,16 +121,16 @@ % 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; @@ -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 @@ -131,4 +185,5 @@ radDepthV{i} = radDepthCube{i}(V); end +%} diff --git a/matRad_siddonRayTracer.m b/matRad_siddonRayTracer.m index 8fde44baf..8ba4007f0 100644 --- a/matRad_siddonRayTracer.m +++ b/matRad_siddonRayTracer.m @@ -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 @@ -217,3 +228,5 @@ for i = 1:numel(cubes) rho{i} = cubes{i}(ix); end + +%} diff --git a/raytracing/matRad_RayTracer.m b/raytracing/matRad_RayTracer.m new file mode 100644 index 000000000..e13a64ebc --- /dev/null +++ b/raytracing/matRad_RayTracer.m @@ -0,0 +1,193 @@ +classdef matRad_RayTracer < handle + %UNTITLED Summary of this class goes here + % Detailed explanation goes here + + properties + lateralCutOff; + end + + + methods + function obj = matRad_RayTracer() + %UNTITLED Construct an instance of this class + % Detailed explanation goes here + matRad_cfg = MatRad_Config.instance(); + obj.lateralCutOff = matRad_cfg.propDoseCalc.defaultGeometricCutOff; + end + + function [alphas,l,rho,d12,ix] = traceRays(this,... + isocenter, ... + resolution, ... + sourcePoints, ... + targetPoints, ... + cubes) + + %Default trivial implementation based on traceRay + nRays = size(targetPoint,1); + nSources = size(sourcePoint,1); + + if nSources ~= nRays && nSources ~= 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Number of source points (%d) needs to be one or equal to number of target points (%d)!',nSources,nRays); + elseif nSources == 1 + sourcePoint = repmat(sourcePoint,nRays,1); + nSources = nRays; + end + + for r = 1:nRays + [alphas{r},l{r},rho{d12},ix{r}] = this.traceRay(isocenter,resolution,sourcePoints(r,:),targetPoints(r,:),cubes); + end + + %pad with NaN values + numval = cellfun(@numel,ix); + maxnumval = max(numval); + + nanpad = @(x) [x(1:end), NaN(maxnumval - length(x),1)]; + + alphas = cellfun(nanpad,alphas,'UniformOutput',false); + l = cellfun(nanpad,l,'UniformOutput',false); + ix = cellfun(nanpad,ix,'UniformOutput',false); + + for c = 1:numel(cubes) + rho{c} = cellfun(nanpad,rho{c},'UniformOutput',false); + end + % Now the output should be consistent + + end + + function [alphas,l,rho,d12,ix] = traceRay(this,... + isocenter, ... + resolution, ... + sourcePoint, ... + targetPoint, ... + cubes) + error('Needs to be implemented!'); + end + + function [radDepthsV,radDepthCube] = traceCube(this,stfElement,ct,V,rot_coordsV) + matRad_cfg = MatRad_Config.instance(); + + if ~isstruct(stfElement) || numel(stfElement) ~= 1 + matRad_cfg.dispError('The RayTracer does not accept stf struct arrays stf and only operates on a single field!'); + end + + %If no subset of voxels is specified, take all of them + if nargin < 4 + 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 < 5 + [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 + + % 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 + stfElement.sourcePoint_bev(2); + + + % set up list with bev coordinates for calculation of radiological depth + coords = zeros(prod(ct.cubeDim),3); + coords(V,:) = rot_coordsV; + + % calculate spacing of rays on ray matrix + rayMxSpacing = 1/sqrt(2) * min([ct.resolution.x ct.resolution.y ct.resolution.z]); + + % define candidate ray matrix covering 1000x1000mm^2 + numOfCandidateRays = 2 * ceil(500/rayMxSpacing) + 1; + candidateRayMx = zeros(numOfCandidateRays); + + % define coordinates + [candidateRaysCoords_X,candidateRaysCoords_Z] = meshgrid(rayMxSpacing*[floor(-500/rayMxSpacing):ceil(500/rayMxSpacing)]); + + % check which rays should be used + for i = 1:stfElement.numOfRays + + ixCandidates = (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 ... + <= this.lateralCutOff^2; + + candidateRayMx(ixCandidates) = 1; + + end + + % set up ray matrix + rayMx_bev = [candidateRaysCoords_X(logical(candidateRayMx(:))) ... + rayMx_bev_y*ones(sum(candidateRayMx(:)),1) ... + candidateRaysCoords_Z(logical(candidateRayMx(:)))]; + + % Rotation matrix. Transposed because of row vectors + 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; + + % criterium for ray selection + raySelection = rayMxSpacing/2; + + %Trace all selected rays + [~,l,rho,~,ixHitVoxel] = this.traceRays(stfElement.isoCenter, ... + ct.resolution, ... + stfElement.sourcePoint, ... + rayMx_world, ... + 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 + % in x and z direction is smaller than the resolution of the ray matrix + scale_factor = NaN(size(ixHitVoxel)); + valid_ix = ~isnan(ixHitVoxel); + scale_factor(valid_ix) = (rayMx_bev_y - stfElement.sourcePoint_bev(2)) ./ ... + coords(ixHitVoxel(valid_ix),2); + + x_dist = NaN(size(ixHitVoxel)); + z_dist = NaN(size(ixHitVoxel)); + + x_dist(valid_ix) = coords(ixHitVoxel(valid_ix),1).*scale_factor(valid_ix); + x_dist = x_dist - rayMx_bev(:,1); + + z_dist(valid_ix) = coords(ixHitVoxel(valid_ix),3).*scale_factor(valid_ix); + z_dist = z_dist - rayMx_bev(:,3); + + % Find indices + ixRememberFromCurrTracing = x_dist > -raySelection & x_dist <= raySelection ... + & z_dist > -raySelection & z_dist <= raySelection; + + % set up rad depth cube for results + radDepthCube = repmat({NaN(ct.cubeDim)},ct.numOfCtScen); + radDepthV = cell(size(radDepthCube)); + + for j = 1:ct.numOfCtScen + rayDistances = l .* rho{j}; + rayWepl = cumsum(rayDistances,2) - rayDistances/2; + + radDepthCube{j}(ixHitVoxel(ixRememberFromCurrTracing)) = rayWepl(ixRememberFromCurrTracing); + + radDepthsV{j} = radDepthCube{j}(V); + end + end + + + end +end + diff --git a/raytracing/matRad_RayTracerSiddon.m b/raytracing/matRad_RayTracerSiddon.m new file mode 100644 index 000000000..6cc775c2c --- /dev/null +++ b/raytracing/matRad_RayTracerSiddon.m @@ -0,0 +1,305 @@ +classdef matRad_RayTracerSiddon < matRad_RayTracer + %UNTITLED2 Summary of this class goes here + % Detailed explanation goes here + + + properties (SetAccess = private) + sourcePoint; + targetPoint; + rayVec; + + xPlanes; + yPlanes; + zPlanes; + numPlanes; + resolution; + cubeDim; + + end + + methods + function this = matRad_RayTracerSiddon() + + end + + function [alphas,l,rho,d12,ix] = traceRay(this,... + isocenter, ... + resolution, ... + sourcePoint, ... + targetPoint, ... + cubes) + + %traceRay Traces an individual ray + % Detailed explanation goes here + + nRays = size(targetPoint,1); + nSources = size(sourcePoint,1); + + if nRays ~= 1 || nSources ~= 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Number of target Points and source points needs to be equal to one! If you want to trace multiple rays at once, use traceRays instead!',nSources,nRays); + end + + [alphas,l,rho,d12,ix] = this.traceRays(isocenter,resolution,sourcePoint,targetPoint,cubes); + end + + function [alphas,l,rho,d12,ix] = traceRays(this,... + isocenter, ... + resolution, ... + sourcePoint, ... + targetPoint, ... + cubes) + + nRays = size(targetPoint,1); + nSources = size(sourcePoint,1); + + if nSources ~= nRays && nSources ~= 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Number of source points (%d) needs to be one or equal to number of target points (%d)!',nSources,nRays); + elseif nSources == 1 + sourcePoint = repmat(sourcePoint,nRays,1); + nSources = nRays; + end + + this.sourcePoint = sourcePoint + isocenter; + this.targetPoint = targetPoint + isocenter; + + this.initializeGeometry(resolution,size(cubes{1})); + + + + % eq 7 & 8 + % Calculate relative distances (alphas) at which intersections + % occur + alphas = this.computeAllAlphas(); + + % eq 11 + % Calculate the distance from source to target point. + d12 = vecnorm(this.rayVec,2,2); + + + % eq 10 + % Calculate the voxel intersection length. + tmpDiff = diff(alphas,1,2); + + l = d12.*tmpDiff; + alphas_mid = alphas(:,1:end-1) + 0.5*tmpDiff; + + + + % eq 12 + % Calculate the voxel indices: first convert to physical coords + % and convert to voxel indices + i = round((this.sourcePoint(:,1) + alphas_mid.*this.rayVec(:,1))./resolution.x); + j = round((this.sourcePoint(:,2) + alphas_mid.*this.rayVec(:,2))./resolution.y); + k = round((this.sourcePoint(:,3) + alphas_mid.*this.rayVec(:,3))./resolution.z); + + + % Handle numerical instabilities at the borders. + i(i<1) = 1; i(i>length(this.xPlanes)-1) = length(this.xPlanes) - 1; + j(j<1) = 1; j(j>length(this.yPlanes)-1) = length(this.yPlanes) - 1; + k(k<1) = 1; k(k>length(this.zPlanes)-1) = length(this.zPlanes) - 1; + + valIx = ~isnan(alphas_mid); + + %In Matlab direct assignment with sub2ind would work, Octave + %however does not like NaN values in the subscripts + ix = NaN(size(valIx)); + ix(valIx) = sub2ind([this.numPlanes(2), this.numPlanes(1), this.numPlanes(3)]-1,j(valIx),i(valIx),k(valIx)); + + for i = 1:numel(cubes) + rho{i} = NaN(size(valIx)); + rho{i}(valIx) = cubes{i}(ix(valIx)); + end + + + end + + + end + + methods (Access = protected) + function alphas = computeAllAlphas(this) + + % Here we setup grids to enable logical indexing when computing + % the alphas along each dimension. All alphas between the + % minimum and maximum index will be computed, with additional + % exclusion of singular plane occurences (max == min) + % All values out of scope will be set to NaN. + + nRays = size(this.rayVec,1); + + % eq 4 + % Calculate parametrics values of \alpha_{min} and \alpha_{max} for every + % axis, intersecting the ray with the sides of the CT. + aX_1 = (this.xPlanes(1) - this.sourcePoint(:,1)) ./ this.rayVec(:,1); + aX_end = (this.xPlanes(end) - this.sourcePoint(:,1)) ./ this.rayVec(:,1); + + tmpIx = this.targetPoint(:,1) == this.sourcePoint(:,1); + aX_1(tmpIx) = NaN; + aX_end(tmpIx) = NaN; + + + aY_1 = (this.yPlanes(1) - this.sourcePoint(:,2)) ./ this.rayVec(:,2); + aY_end = (this.yPlanes(end) - this.sourcePoint(:,2)) ./ this.rayVec(:,2); + + tmpIx = this.targetPoint(:,2) == this.sourcePoint(:,2); + aY_1(tmpIx) = NaN; + aY_end(tmpIx) = NaN; + + aZ_1 = (this.zPlanes(1) - this.sourcePoint(:,3)) ./ this.rayVec(:,3); + aZ_end = (this.zPlanes(end) - this.sourcePoint(:,3)) ./ this.rayVec(:,3); + + tmpIx = this.targetPoint(:,3) == this.sourcePoint(:,3); + aZ_1(tmpIx) = NaN; + aZ_end(tmpIx) = NaN; + + + % eq 5 + % Compute the \alpha_{min} and \alpha_{max} in terms of parametric values + % given by equation 4. + alpha_limits(:,1) = max([zeros(nRays,1) min(aX_1,aX_end) min(aY_1,aY_end) min(aZ_1,aZ_end)],[],2); + alpha_limits(:,2) = min([ones(nRays,1) max(aX_1,aX_end) max(aY_1,aY_end) max(aZ_1,aZ_end)],[],2); + + % eq 6 + % Calculate the range of indeces who gives parametric values for + % intersected planes. + + % get indices + ixTpBiggerSp = this.targetPoint > this.sourcePoint; %this.rayVec > 0 + ixTpSmallerSp = this.targetPoint < this.sourcePoint; %this.rayVec < 0 + + + + alphaTmp = NaN(nRays,2); + alphaTmp(ixTpBiggerSp(:,1),1) = alpha_limits(ixTpBiggerSp(:,1),1); + alphaTmp(ixTpSmallerSp(:,1),1) = alpha_limits(ixTpSmallerSp(:,1),2); + alphaTmp(ixTpBiggerSp(:,1),2) = alpha_limits(ixTpBiggerSp(:,1),2); + alphaTmp(ixTpSmallerSp(:,1),2) = alpha_limits(ixTpSmallerSp(:,1),1); + + i_min = this.numPlanes(1) - (this.xPlanes(end) - alphaTmp(:,1) .* this.rayVec(:,1) - this.sourcePoint(:,1))./this.resolution.x; + i_max = 1 + (this.sourcePoint(:,1) + alphaTmp(:,2) .* this.rayVec(:,1) - this.xPlanes(1))./this.resolution.x; + %rounding + i_min = ceil(1/1000 * (round(1000*i_min))); + i_max = floor(1/1000 * (round(1000*i_max))); + + + %j + alphaTmp(:) = NaN; + alphaTmp(ixTpBiggerSp(:,2),1) = alpha_limits(ixTpBiggerSp(:,2),1); + alphaTmp(ixTpSmallerSp(:,2),1) = alpha_limits(ixTpSmallerSp(:,2),2); + alphaTmp(ixTpBiggerSp(:,2),2) = alpha_limits(ixTpBiggerSp(:,2),2); + alphaTmp(ixTpSmallerSp(:,2),2) = alpha_limits(ixTpSmallerSp(:,2),1); + + j_min = this.numPlanes(2) - (this.yPlanes(end) - alphaTmp(:,1) .* this.rayVec(:,2) - this.sourcePoint(:,2))./this.resolution.y; + j_max = 1 + (this.sourcePoint(:,2) + alphaTmp(:,2) .* this.rayVec(:,2) - this.yPlanes(1))./this.resolution.y; + %rounding + j_min = ceil(1/1000 * (round(1000*j_min))); + j_max = floor(1/1000 * (round(1000*j_max))); + + %k + alphaTmp(:) = NaN; + alphaTmp(ixTpBiggerSp(:,3),1) = alpha_limits(ixTpBiggerSp(:,3),1); + alphaTmp(ixTpSmallerSp(:,3),1) = alpha_limits(ixTpSmallerSp(:,3),2); + alphaTmp(ixTpBiggerSp(:,3),2) = alpha_limits(ixTpBiggerSp(:,3),2); + alphaTmp(ixTpSmallerSp(:,3),2) = alpha_limits(ixTpSmallerSp(:,3),1); + + k_min = this.numPlanes(3) - (this.zPlanes(end) - alphaTmp(:,1) .* this.rayVec(:,3) - this.sourcePoint(:,3))./this.resolution.z; + k_max = 1 + (this.sourcePoint(:,3) + alphaTmp(:,2) .* this.rayVec(:,3) - this.zPlanes(1))./this.resolution.z; + %rounding + k_min = ceil(1/1000 * (round(1000*k_min))); + k_max = floor(1/1000 * (round(1000*k_max))); + + % eq 7 + % For the given range of indices, calculate the paremetrics values who + % represents intersections of the ray with the plane. + + planeIx = 1:length(this.xPlanes); + planeGrid = repmat(this.xPlanes,nRays,1); + planeGrid(planeIx < i_min | planeIx > i_max | (planeIx == i_min & planeIx == i_max) | isnan(i_min) | isnan(i_max)) = NaN; + alpha_x = (planeGrid-this.sourcePoint(:,1))./(this.rayVec(:,1)); + + planeIx = 1:length(this.yPlanes); + planeGrid = repmat(this.yPlanes,nRays,1); + planeGrid(planeIx < j_min | planeIx > j_max | (planeIx == j_min & planeIx == j_max) | isnan(j_min) | isnan(j_max)) = NaN; + alpha_y = (planeGrid-this.sourcePoint(:,2))./(this.rayVec(:,2)); + + planeIx = 1:length(this.zPlanes); + planeGrid = repmat(this.zPlanes,nRays,1); + planeGrid(planeIx < k_min | planeIx > k_max | (planeIx == k_min & planeIx == k_max) | isnan(k_min) | isnan(k_max)) = NaN; + alpha_z = (planeGrid-this.sourcePoint(:,3))./(this.rayVec(:,3)); + + % eq 8 + % Merge parametrics sets. + % The following might look slow but is quite close to Matlab's + % "unique" implementation + alphas = sort([alpha_limits alpha_x alpha_y alpha_z],2); %NaN's are placed at the end when sorting in ascending order + alphas(diff(alphas,1,2) == 0) = NaN; %Remove duplicates + alphas = sort(alphas,2); %Again place NaN's at the end + + %Size Reduction (reduce NaN padding) for further computations + maxNumColumns = max(sum(~isnan(alphas),2)); + alphas = alphas(:,1:maxNumColumns); + end + + + + function this = initializeGeometry(this, ... + resolution, ... + cubeDim) + %initializeGeometry Initializes CT geometry for Tracing + % Allows for vectorized input in source and target points + % input + % resolution: resolution of the cubes [mm/voxel] + % this.sourcePoint: source point(s) of ray tracing. Either an 1x3 + % or Nx3 matrix (with N rays to trace). If + % 1x3 with multiple target points, the one + % this.sourcePoint will be used for all rays. + % this.targetPoint: target point(s) of ray tracing. Is an Nx3 + % matrix where N is the number of rays to + % trace + % cubes: cell array of cubes for ray tracing (it is possible to pass + % multiple cubes for ray tracing to save computation time) + % + % output (see Siddon 1985 Medical Physics for a detailed description of the + % variales) + % alphas relative distance between start and endpoint for the + % intersections with the cube + % l lengths of intersestions with cubes + % rho densities extracted from cubes + % d12 distance between start and endpoint of ray tracing + % ix indices of hit voxels + + % Set private member variables for later use + this.cubeDim = cubeDim; + this.resolution = resolution; + + % Save the numbers of planes. + xNumPlanes = cubeDim(2) + 1; + yNumPlanes = cubeDim(1) + 1; + zNumPlanes = cubeDim(3) + 1; + this.numPlanes = [xNumPlanes,yNumPlanes,zNumPlanes]; + + % eq 3 + % Position of planes in millimeter. 0.5 because the central position + % of the first voxel is at [resolution.x resolution.y resolution.z] + + this.xPlanes = resolution.x*(1:xNumPlanes) - .5*resolution.x; + this.yPlanes = resolution.y*(1:yNumPlanes) - .5*resolution.y; + this.zPlanes = resolution.z*(1:zNumPlanes) - .5*resolution.z; + + + % Construct ray Vectors + this.rayVec = this.targetPoint - this.sourcePoint; + nRays = size(this.rayVec,1); + + + + + + end + + + + end +end \ No newline at end of file