INTERSECTLINECYLINDER Compute intersection points between a line and a cylinder.
POINTS = intersectLineCylinder(LINE, CYLINDER)
Returns intersection points between a line and a cylinder.
Input parameters:
LINE = [x0 y0 z0 dx dy dz]
CYLINDER = [x1 y1 z1 x2 y2 z2 R]
Output:
POINTS = [x1 y1 z1 ; x2 y2 z2]
POINTS = intersectLineCylinder(LINE, CYLINDER, 'checkBounds', B)
Where B is a boolean (TRUE by default), check if the points are within
the bounds defined by the two extreme points. If B is false, the
cylinder is considered to be infinite.
Example
% Compute intersection between simple vertical cylinder and line
line = [60 60 60 1 2 3];
cylinder = [20 50 50 80 50 50 30];
points = intersectLineCylinder(line, cylinder);
% Display the different shapes
figure;
drawCylinder(cylinder);
hold on; light;
axis([0 100 0 100 0 100]);
drawLine3d(line);
drawPoint3d(points, 'ko');
% Compute intersections when one of the points is outside the
% cylinder
line = [80 60 60 1 2 3];
cylinder = [20 50 50 80 50 50 30];
intersectLineCylinder(line, cylinder)
ans =
67.8690 35.7380 23.6069
See also
lines3d, intersectLinePlane, drawCylinder, cylinderSurfaceArea
References
See the link:
http://www.gamedev.net/community/forums/topic.asp?topic_id=4677890001 function points = intersectLineCylinder(line, cylinder, varargin) 0002 %INTERSECTLINECYLINDER Compute intersection points between a line and a cylinder. 0003 % 0004 % POINTS = intersectLineCylinder(LINE, CYLINDER) 0005 % Returns intersection points between a line and a cylinder. 0006 % 0007 % Input parameters: 0008 % LINE = [x0 y0 z0 dx dy dz] 0009 % CYLINDER = [x1 y1 z1 x2 y2 z2 R] 0010 % 0011 % Output: 0012 % POINTS = [x1 y1 z1 ; x2 y2 z2] 0013 % 0014 % POINTS = intersectLineCylinder(LINE, CYLINDER, 'checkBounds', B) 0015 % Where B is a boolean (TRUE by default), check if the points are within 0016 % the bounds defined by the two extreme points. If B is false, the 0017 % cylinder is considered to be infinite. 0018 % 0019 % Example 0020 % % Compute intersection between simple vertical cylinder and line 0021 % line = [60 60 60 1 2 3]; 0022 % cylinder = [20 50 50 80 50 50 30]; 0023 % points = intersectLineCylinder(line, cylinder); 0024 % % Display the different shapes 0025 % figure; 0026 % drawCylinder(cylinder); 0027 % hold on; light; 0028 % axis([0 100 0 100 0 100]); 0029 % drawLine3d(line); 0030 % drawPoint3d(points, 'ko'); 0031 % 0032 % 0033 % % Compute intersections when one of the points is outside the 0034 % % cylinder 0035 % line = [80 60 60 1 2 3]; 0036 % cylinder = [20 50 50 80 50 50 30]; 0037 % intersectLineCylinder(line, cylinder) 0038 % ans = 0039 % 67.8690 35.7380 23.6069 0040 % 0041 % 0042 % See also 0043 % lines3d, intersectLinePlane, drawCylinder, cylinderSurfaceArea 0044 % 0045 % References 0046 % See the link: 0047 % http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 0048 % 0049 0050 % ------ 0051 % Author: David Legland, from a file written by Daniel Trauth (RWTH) 0052 % E-mail: david.legland@inrae.fr 0053 % Created: 2007-01-27 0054 % Copyright 2007-2024 0055 0056 %% Parse input arguments 0057 0058 % default arguments 0059 checkBounds = true; 0060 % type of cylinder, one of {'closed', 'open', 'infinite'} 0061 type = 'closed'; 0062 0063 % parse inputs 0064 while length(varargin)>1 0065 var = varargin{1}; 0066 if strcmpi(var, 'checkbounds') 0067 checkBounds = varargin{2}; 0068 elseif strcmpi(var, 'type') 0069 type = varargin{2}; 0070 else 0071 error(['Unkown argument: ' var]); 0072 end 0073 varargin(1:2) = []; 0074 end 0075 0076 0077 %% Parse cylinder parameters 0078 0079 % Starting point of the line 0080 l0 = line(1:3); 0081 0082 % Direction vector of the line 0083 dl = line(4:6); 0084 0085 % position of cylinder extremities 0086 c1 = cylinder(1:3); 0087 c2 = cylinder(4:6); 0088 0089 % Direction vector of the cylinder 0090 dc = c2 - c1; 0091 0092 % Radius of the cylinder 0093 r = cylinder(7); 0094 0095 0096 %% Resolution of a quadratic equation to find the increment 0097 0098 % normalisation coefficient corresponding to direction of vector 0099 coef = dc / dot(dc, dc); 0100 0101 % Substitution of parameters 0102 e = dl - dot(dl,dc) * coef; 0103 f = (l0-c1) - dot(l0-c1, dc) * coef; 0104 0105 % Coefficients of 2-nd order equation 0106 A = dot(e, e); 0107 B = 2 * dot(e,f); 0108 C = dot(f,f) - r^2; 0109 0110 % compute discriminant 0111 delta = B^2 - 4*A*C; 0112 0113 % check existence of solution(s) 0114 if delta < 0 0115 points = zeros(0, 3); 0116 return; 0117 end 0118 0119 % extract roots 0120 pos1 = (-B + sqrt(delta)) / (2*A); 0121 pos2 = (-B - sqrt(delta)) / (2*A); 0122 posList = [pos1;pos2]; 0123 0124 0125 %% Estimation of point positions 0126 0127 % process the smallest position 0128 pos1 = min(posList); 0129 0130 % Point on the line: l0 + x*dl = p 0131 point1 = l0 + pos1 * dl; 0132 0133 % process the greatest position 0134 pos2 = max(posList); 0135 0136 % Point on the line: l0 + x*dl = p 0137 point2 = l0 + pos2 * dl; 0138 0139 % Format result 0140 points = [point1 ; point2]; 0141 0142 0143 %% Check if points are located between bounds 0144 0145 % if checkBounds option is not set, we can simply skip the rest 0146 if ~checkBounds || strncmpi(type, 'infinite', 1) 0147 return; 0148 end 0149 0150 % compute cylinder axis 0151 axis = [c1 dc]; 0152 0153 % compute position on axis 0154 ts = line3dPosition(points, axis); 0155 0156 % check bounds for open cylinder 0157 % (keep only intersection points whose projection is between the two 0158 % cylinder extremities) 0159 if strncmpi(type, 'open', 1) 0160 ind = ts>=0 & ts<=1; 0161 points = points(ind, :); 0162 return; 0163 end 0164 0165 % which intersection fall before and after bounds 0166 ind1 = find(ts < 0); 0167 ind2 = find(ts > 1); 0168 0169 % case of both intersection on the same side -> no intersection 0170 if length(ind1) == 2 || length(ind2) == 2 0171 points = zeros(0, 3); 0172 return; 0173 end 0174 0175 % Process the remaining case of closed cylinder 0176 % -> compute eventual intersection(s) with end faces 0177 if ~isempty(ind1) 0178 plane = createPlane(c1, dc); 0179 points(ind1, :) = intersectLinePlane(line, plane); 0180 end 0181 if ~isempty(ind2) 0182 plane = createPlane(c2, dc); 0183 points(ind2, :) = intersectLinePlane(line, plane); 0184 end 0185 0186