cylinder equation

((P - C) x D)^2 = r^2

P is a point on the cylinder, C is a point on the cylinder core (the centre of the cylinder, or the base), D is the direction of the cylinder, r is the radius. Pretty much like a sphere equation.

so develop with a ray equation (P = O + S.t), and you get

((O - C) x D + (S x D).t)^2 = r2

((O - C) x D)^2 - r^2 + 2 * ( ((O - C)xD) * (S x D) ).t + (S x D)^2 . t^2 = 0

second order equation in the form

a.t^2 + b.t + c = 0

where

a = (S x D)^2
b = 2 * ( ((O - C) x D) . (S x D) )
c = ((O - C)xD)^2 - r^2

determinant is

d = (b*b - 4.a.c)

if (d < 0), the ray does not intersect

else the solutions are
t0 = (-b-sqrt(d)) / 2.a
t1 = (-b+sqrt(d)) / 2.a

and

tcylmin = min(t0, t1)
tcylmax = max(t0, t1)


now, the intersection points with the bottom and top caps. I assume the cylinder has the point C as the centre of the cylinder, and the cylinder is of length l.

ttop = ((C - O) * D + l/2) / (S * D)
tbot = ((C - O) * D - l/2) / (S * D)

and

tcapmin = min(ttop, tbot)
tcapmax = max(ttop, tbot)


if (tcapmin > tcylmax) || (tcapmax < tcylmin) the ray does not intersect

else the ray intersects at the overlapping interval between [tcylmin, tcylmax] and [tcapmin, tcapmax].

there are optimisations to be made, like you only need to test one cap if the origin of the ray is above or below the cylinder, and no cap if the origin of the ray is ''leveled'' with the cylinder.

another way is to transform the ray into the cylinder space, and find the intersection there, then transform the points back into world space. First, you need to find the quaternion (or 3x3 rotation matrix) to rotate the ray from it''s direction back into an axis (X, Y or Z axis). If you have the cylinder orientation stored as a matrix instead of a directional vector, it''s simply the transpose of the matrix to get the inverse transform.

Cylinder (C, Ori, l)
Ray (P, S)

Q = (P - C) * Ori.Transpose();
T = S * Ori.Transpose();


then with Q and T, you have the cylinder aligned with an axis (as if the cylinder never rotated), and the ray relative to that cylinder. And you can do the cylinder/ray test virtually in 2D.

yeah, I though that was quite standard as a notation. It's the dot product of a vector by itself, or it's squared magnitude.

x is cross product
. is dot product
* is multiplication
- is substraction
+ is addiciton

and I also use '*' as the dot product sometimes. so, in c++

struct Vector{    float x, y, z;    Vector(float _x, float _y, float _z)    : x(_x)    , y(_y)    , z(_z)    {}    Vector& operator += (const Vector& V) { x += V.x; y += V.y; z += V.z; return  *this; }    Vector& operator -= (const Vector& V) { x -= V.x; y -= V.y; z -= V.z; return  *this; }    Vector& operator *= (float k) { x *= k; y *= k; z *= k; return  *this; }    Vector& operator /= (float k) { return (*this) *= 1.0f / k; }    float operator * (const Vector& V) { return x*V.x + y*V.y + z*V.z; } // dot product    Vector operator -(void) const { return Vector(-x, -y, -z); }    Vector& operator ^=(const Vector &V) // Cross product    {        float Tempx	= (y * V.z) - (z * V.y);        float Tempy	= (z * V.x) - (x * V.z);        z		= (x * V.y) - (y * V.x);        x               = Tempx;        y               = Tempy;        return *this;    }    Vector operator + (const Vector& V) const { Vector Temp(*this); return Temp += V; }    Vector operator - (const Vector& V) const { Vector Temp(*this); return Temp -= V; }    Vector operator * (float k) const { Vector Temp(*this); return Temp *= k; }    Vector operator / (float k) const { Vector Temp(*this); return Temp /= k; }    Vector operator ^ (const Vector& V) const {	Vector Temp(*this); return Temp ^= V; }    friend Vector operator * (float k, const Vector& V) { return V * k; }    float MagnitudeSquared(void) const { return (*this) * (*this); }    float Magnitude(void) const { return (float) sqrt(MagnitudeSquared()); }};

[edited by - oliii on December 6, 2003 9:07:40 AM]