matRad_siddonRayTracer

 siddon ray tracing through 3D cube to calculate the radiological depth 
 according to Siddon 1985 Medical Physics
 
 call
   [alphas,l,rho,d12,vis] = matRad_siddonRayTracer(isocenter, ...
                               resolution, ...
                               sourcePoint, ...
                               targetPoint, ...
                               cubes)

 input
   isocenter:      isocenter within cube [voxels]
   resolution:     resolution of the cubes [mm/voxel]
   sourcePoint:    source point of ray tracing
   targetPoint:    target point of ray tracing
   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

 References
   [1] http://www.ncbi.nlm.nih.gov/pubmed/4000088

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 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.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [alphas,l,rho,d12,ix] = matRad_siddonRayTracer(isocenter, ...
0002                                     resolution, ...
0003                                     sourcePoint, ...
0004                                     targetPoint, ...
0005                                     cubes)
0006 % siddon ray tracing through 3D cube to calculate the radiological depth
0007 % according to Siddon 1985 Medical Physics
0008 %
0009 % call
0010 %   [alphas,l,rho,d12,vis] = matRad_siddonRayTracer(isocenter, ...
0011 %                               resolution, ...
0012 %                               sourcePoint, ...
0013 %                               targetPoint, ...
0014 %                               cubes)
0015 %
0016 % input
0017 %   isocenter:      isocenter within cube [voxels]
0018 %   resolution:     resolution of the cubes [mm/voxel]
0019 %   sourcePoint:    source point of ray tracing
0020 %   targetPoint:    target point of ray tracing
0021 %   cubes:          cell array of cubes for ray tracing (it is possible to pass
0022 %                   multiple cubes for ray tracing to save computation time)
0023 %
0024 % output (see Siddon 1985 Medical Physics for a detailed description of the
0025 % variales)
0026 %   alphas          relative distance between start and endpoint for the
0027 %                    intersections with the cube
0028 %   l               lengths of intersestions with cubes
0029 %   rho             densities extracted from cubes
0030 %   d12             distance between start and endpoint of ray tracing
0031 %   ix              indices of hit voxels
0032 %
0033 % References
0034 %   [1] http://www.ncbi.nlm.nih.gov/pubmed/4000088
0035 %
0036 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0037 %
0038 % Copyright 2015 the matRad development team.
0039 %
0040 % This file is part of the matRad project. It is subject to the license
0041 % terms in the LICENSE file found in the top-level directory of this
0042 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0043 % of the matRad project, including this file, may be copied, modified,
0044 % propagated, or distributed except according to the terms contained in the
0045 % LICENSE file.
0046 %
0047 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0048 
0049 % Add isocenter to source and target point. Because the algorithm does not
0050 % works with negatives values. This put (resolution.x,resolution.y,resolution.z)
0051 % in the center of first voxel
0052 
0053 sourcePoint = sourcePoint + isocenter;
0054 targetPoint = targetPoint + isocenter;
0055 
0056 % Save the numbers of planes.
0057 [yNumPlanes, xNumPlanes, zNumPlanes] = size(cubes{1});
0058 xNumPlanes = xNumPlanes + 1;
0059 yNumPlanes = yNumPlanes + 1;
0060 zNumPlanes = zNumPlanes + 1;
0061 
0062 % eq 11
0063 % Calculate the distance from source to target point.
0064 d12 = norm(sourcePoint-targetPoint);
0065 
0066 % eq 3
0067 % Position of first planes in millimeter. 0.5 because the central position
0068 % of the first voxel is at [resolution.x resolution.y resolution.z]
0069 xPlane_1 = .5*resolution.x;
0070 yPlane_1 = .5*resolution.y;
0071 zPlane_1 = .5*resolution.z;
0072 
0073 % Position of last planes in milimiter
0074 xPlane_end = (xNumPlanes - .5)*resolution.x;
0075 yPlane_end = (yNumPlanes - .5)*resolution.y;
0076 zPlane_end = (zNumPlanes - .5)*resolution.z;
0077 
0078 % eq 4
0079 % Calculate parametrics values of \alpha_{min} and \alpha_{max} for every
0080 % axis, intersecting the ray with the sides of the CT.
0081 if targetPoint(1) ~= sourcePoint(1)
0082     aX_1   = (xPlane_1 - sourcePoint(1)) / (targetPoint(1) - sourcePoint(1));
0083     aX_end = (xPlane_end - sourcePoint(1)) / (targetPoint(1) - sourcePoint(1));
0084 else
0085     aX_1   = [];
0086     aX_end = [];
0087 end
0088 if targetPoint(2) ~= sourcePoint(2)
0089     aY_1   = (yPlane_1 - sourcePoint(2)) / (targetPoint(2) - sourcePoint(2));
0090     aY_end = (yPlane_end - sourcePoint(2)) / (targetPoint(2) - sourcePoint(2));
0091 else
0092     aY_1   = [];
0093     aY_end = [];
0094 end
0095 if targetPoint(3) ~= sourcePoint(3)
0096     aZ_1   = (zPlane_1 - sourcePoint(3)) / (targetPoint(3) - sourcePoint(3));
0097     aZ_end = (zPlane_end - sourcePoint(3)) / (targetPoint(3) - sourcePoint(3));
0098 else
0099     aZ_1   = [];
0100     aZ_end = [];
0101 end
0102 
0103 % eq 5
0104 % Compute the \alpha_{min} and \alpha_{max} in terms of parametric values
0105 % given by equation 4.
0106 alpha_min = max([0 min(aX_1,aX_end) min(aY_1,aY_end) min(aZ_1,aZ_end)]);
0107 alpha_max = min([1 max(aX_1,aX_end) max(aY_1,aY_end) max(aZ_1,aZ_end)]);
0108 
0109 % eq 6
0110 % Calculate the range of indeces who gives parametric values for
0111 % intersected planes.
0112 if targetPoint(1) == sourcePoint(1)
0113     i_min = []; i_max = [];
0114 elseif targetPoint(1) > sourcePoint(1)
0115     i_min = xNumPlanes - (xPlane_end - alpha_min * (targetPoint(1) - sourcePoint(1)) - sourcePoint(1))/resolution.x;
0116     i_max = 1          + (sourcePoint(1) + alpha_max * (targetPoint(1) - sourcePoint(1)) - xPlane_1)/resolution.x;
0117     % rounding
0118     i_min = ceil(1/1000*(round(1000*i_min)));
0119     i_max = floor(1/1000*(round(1000*i_max)));
0120 else
0121     i_min = xNumPlanes - (xPlane_end - alpha_max * (targetPoint(1) - sourcePoint(1)) - sourcePoint(1))/resolution.x;
0122     i_max = 1          + (sourcePoint(1) + alpha_min * (targetPoint(1) - sourcePoint(1)) - xPlane_1)/resolution.x;
0123     i_min = ceil(1/1000*(round(1000*i_min)));
0124     i_max = floor(1/1000*(round(1000*i_max)));
0125 end
0126 if targetPoint(2) == sourcePoint(2)
0127     j_min = []; j_max = [];
0128 elseif targetPoint(2) > sourcePoint(2)
0129     j_min = yNumPlanes - (yPlane_end - alpha_min * (targetPoint(2) - sourcePoint(2)) - sourcePoint(2))/resolution.y;
0130     j_max = 1          + (sourcePoint(2) + alpha_max * (targetPoint(2) - sourcePoint(2)) - yPlane_1)/resolution.y;
0131     j_min = ceil(1/1000*(round(1000*j_min)));
0132     j_max = floor(1/1000*(round(1000*j_max)));
0133 else
0134     j_min = yNumPlanes - (yPlane_end - alpha_max * (targetPoint(2) - sourcePoint(2)) - sourcePoint(2))/resolution.y;
0135     j_max = 1          + (sourcePoint(2) + alpha_min * (targetPoint(2) - sourcePoint(2)) - yPlane_1)/resolution.y;
0136     j_min = ceil(1/1000*(round(1000*j_min)));
0137     j_max = floor(1/1000*(round(1000*j_max)));
0138 end
0139 if targetPoint(3) == sourcePoint(3)
0140     k_min = []; k_max = [];
0141 elseif targetPoint(3) >= sourcePoint(3)
0142     k_min = zNumPlanes - (zPlane_end - alpha_min * (targetPoint(3) - sourcePoint(3)) - sourcePoint(3))/resolution.z;
0143     k_max = 1          + (sourcePoint(3) + alpha_max * (targetPoint(3) - sourcePoint(3)) - zPlane_1)/resolution.z;
0144     k_min = ceil(1/1000*(round(1000*k_min)));
0145     k_max = floor(1/1000*(round(1000*k_max)));
0146 else
0147     k_min = zNumPlanes - (zPlane_end - alpha_max * (targetPoint(3) - sourcePoint(3)) - sourcePoint(3))/resolution.z;
0148     k_max = 1          + (sourcePoint(3) + alpha_min * (targetPoint(3) - sourcePoint(3)) - zPlane_1)/resolution.z;
0149     k_min = ceil(1/1000*(round(1000*k_min)));
0150     k_max = floor(1/1000*(round(1000*k_max)));
0151 end
0152 
0153 % eq 7
0154 % For the given range of indices, calculate the paremetrics values who
0155 % represents intersections of the ray with the plane.
0156 if i_min ~= i_max
0157     if targetPoint(1) > sourcePoint(1)
0158         alpha_x = (resolution.x*(i_min:1:i_max)-sourcePoint(1)-.5*resolution.x)/(targetPoint(1)-sourcePoint(1));
0159     else
0160         alpha_x = (resolution.x*(i_max:-1:i_min)-sourcePoint(1)-.5*resolution.x)/(targetPoint(1)-sourcePoint(1));
0161     end
0162 else
0163     alpha_x = [];
0164 end
0165 if j_min ~= j_max
0166     if targetPoint(2) > sourcePoint(2)
0167         alpha_y = (resolution.y*(j_min:1:j_max)-sourcePoint(2)-.5*resolution.y)/(targetPoint(2)-sourcePoint(2));
0168     else
0169         alpha_y = (resolution.y*(j_max:-1:j_min)-sourcePoint(2)-.5*resolution.y)/(targetPoint(2)-sourcePoint(2));
0170     end
0171 else
0172     alpha_y = [];
0173 end
0174 if k_min ~= k_max
0175     if targetPoint(3) > sourcePoint(3)
0176         alpha_z = (resolution.z*(k_min:1:k_max)-sourcePoint(3)-.5*resolution.z)/(targetPoint(3)-sourcePoint(3));
0177     else
0178         alpha_z = (resolution.z*(k_max:-1:k_min)-sourcePoint(3)-.5*resolution.z)/(targetPoint(3)-sourcePoint(3));
0179     end
0180 else
0181     alpha_z = [];
0182 end
0183 
0184 % eq 8
0185 % Merge parametrics sets.
0186 alphas = unique([alpha_min alpha_x alpha_y alpha_z alpha_max]);
0187 
0188 % eq 10
0189 % Calculate the voxel intersection length.
0190 l = d12*diff(alphas);
0191 
0192 % eq 13
0193 % Calculate \alpha_{middle}
0194 alphas_mid = .5*(alphas(1:end-1)+alphas(2:end));
0195 
0196 % eq 12
0197 % Calculate the voxel indices: first convert to physical coords
0198 i_mm = sourcePoint(1) + alphas_mid*(targetPoint(1) - sourcePoint(1));
0199 j_mm = sourcePoint(2) + alphas_mid*(targetPoint(2) - sourcePoint(2));
0200 k_mm = sourcePoint(3) + alphas_mid*(targetPoint(3) - sourcePoint(3));
0201 % then convert to voxel index
0202 i = round(i_mm/resolution.x);
0203 j = round(j_mm/resolution.y);
0204 k = round(k_mm/resolution.z);
0205 
0206 % Handle numerical instabilities at the borders.
0207 i(i<=0) = 1; j(j<=0) = 1; k(k<=0) = 1;
0208 i(i>xNumPlanes-1) = xNumPlanes-1;
0209 j(j>yNumPlanes-1) = yNumPlanes-1;
0210 k(k>zNumPlanes-1) = zNumPlanes-1;
0211 
0212 % Convert to linear indices
0213 ix = j + (i-1)*size(cubes{1},1) + (k-1)*size(cubes{1},1)*size(cubes{1},2); 
0214 
0215 % obtains the values from cubes
0216 rho = cell(numel(cubes),1);
0217 for i = 1:numel(cubes)
0218     rho{i} = cubes{i}(ix);
0219 end