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@inra.fr 0053 % Created: 2007-01-27 0054 0055 % HISTORY 0056 % 2010-10-21 change cylinder argument convention, add bounds check and doc 0057 % 2010-10-21 add check for points on cylinders, update doc 0058 0059 0060 %% Parse input arguments 0061 0062 % default arguments 0063 checkBounds = true; 0064 % type of cylinder, one of {'closed', 'open', 'infinite'} 0065 type = 'closed'; 0066 0067 % parse inputs 0068 while length(varargin)>1 0069 var = varargin{1}; 0070 if strcmpi(var, 'checkbounds') 0071 checkBounds = varargin{2}; 0072 elseif strcmpi(var, 'type') 0073 type = varargin{2}; 0074 else 0075 error(['Unkown argument: ' var]); 0076 end 0077 varargin(1:2) = []; 0078 end 0079 0080 0081 %% Parse cylinder parameters 0082 0083 % Starting point of the line 0084 l0 = line(1:3); 0085 0086 % Direction vector of the line 0087 dl = line(4:6); 0088 0089 % position of cylinder extremities 0090 c1 = cylinder(1:3); 0091 c2 = cylinder(4:6); 0092 0093 % Direction vector of the cylinder 0094 dc = c2 - c1; 0095 0096 % Radius of the cylinder 0097 r = cylinder(7); 0098 0099 0100 %% Resolution of a quadratic equation to find the increment 0101 0102 % normalisation coefficient corresponding to direction of vector 0103 coef = dc / dot(dc, dc); 0104 0105 % Substitution of parameters 0106 e = dl - dot(dl,dc) * coef; 0107 f = (l0-c1) - dot(l0-c1, dc) * coef; 0108 0109 % Coefficients of 2-nd order equation 0110 A = dot(e, e); 0111 B = 2 * dot(e,f); 0112 C = dot(f,f) - r^2; 0113 0114 % compute discriminant 0115 delta = B^2 - 4*A*C; 0116 0117 % check existence of solution(s) 0118 if delta < 0 0119 points = zeros(0, 3); 0120 return; 0121 end 0122 0123 % extract roots 0124 pos1 = (-B + sqrt(delta)) / (2*A); 0125 pos2 = (-B - sqrt(delta)) / (2*A); 0126 posList = [pos1;pos2]; 0127 0128 0129 %% Estimation of point positions 0130 0131 % process the smallest position 0132 pos1 = min(posList); 0133 0134 % Point on the line: l0 + x*dl = p 0135 point1 = l0 + pos1 * dl; 0136 0137 % process the greatest position 0138 pos2 = max(posList); 0139 0140 % Point on the line: l0 + x*dl = p 0141 point2 = l0 + pos2 * dl; 0142 0143 % Format result 0144 points = [point1 ; point2]; 0145 0146 0147 %% Check if points are located between bounds 0148 0149 % if checkBounds option is not set, we can simply skip the rest 0150 if ~checkBounds || strncmpi(type, 'infinite', 1) 0151 return; 0152 end 0153 0154 % compute cylinder axis 0155 axis = [c1 dc]; 0156 0157 % compute position on axis 0158 ts = linePosition3d(points, axis); 0159 0160 % check bounds for open cylinder 0161 % (keep only intersection points whose projection is between the two 0162 % cylinder extremities) 0163 if strncmpi(type, 'open', 1) 0164 ind = ts>=0 & ts<=1; 0165 points = points(ind, :); 0166 return; 0167 end 0168 0169 % which intersection fall before and after bounds 0170 ind1 = find(ts < 0); 0171 ind2 = find(ts > 1); 0172 0173 % case of both intersection on the same side -> no intersection 0174 if length(ind1) == 2 || length(ind2) == 2 0175 points = zeros(0, 3); 0176 return; 0177 end 0178 0179 % Process the remaining case of closed cylinder 0180 % -> compute eventual intersection(s) with end faces 0181 if ~isempty(ind1) 0182 plane = createPlane(c1, dc); 0183 points(ind1, :) = intersectLinePlane(line, plane); 0184 end 0185 if ~isempty(ind2) 0186 plane = createPlane(c2, dc); 0187 points(ind2, :) = intersectLinePlane(line, plane); 0188 end 0189 0190