I am trying to write a collision detection class for a 3d game engine.
 I need the equation to find out when a moving particle collides with
a cylinder that extends to infinity with both ends.  If you take a
look at my website, http://www.geocities.com/t_F_l_y/, you can see
what I am looking for in the final equations, for example in "Point of
collision Sphere" I got T = (-B - SQRT(B2 - (4)(A)(C) ) / (2A).  The
page that I need help with is at,
http://www.geocities.com/t_F_l_y/PointOfCollisionCylinder.htm, P0 is
the particles initial position and PV is the particles translation
vector, aka velocity, C0 and CV define the cylinder's axis.  C0 is a
point on the cylinder axis and CV is its direction.  You will not have
to worry about the end of the cylinder just imagine it goes off to
infinitely.  The equation I am looking for should be similiar to the
Point of Collision Sphere equation except for it will include some
free variable S.  This may be impossible to do and I might have to
resort to some form of guess and check, if you can prove this then I
will be happy with the answer.

Also the reason for only considering the one root of the quadratic
equation in the "Point of collision Sphere" is that this is the
particle's first contact with the sphere after this point the particle
is in the sphere and the other root occurs when the particle leaves
the sphere, which is not needed.  If you want to write everything out
and show how you figured it out that would be great you can use x^2
for x squared is fine Or even one final equation would be a great
help, Thanks.

---

Request for Question Clarification by mathtalk-ga on 16 Jun 2004 20:28 PDT

Hi, tfly-ga:

Are you working with a collision between a moving point and a
cylinder, or between a moving sphere and a cylinder?

The discussion above would lead me to think that you have "point" P0
moving along a straight trajectory, and you want to find when/where
the point first touches the cylinder (if it ever does).  This is one
problem (pretty straightforward, just a quadratic equation to solve).

However your "PointOfCollisionCylinder.htm" Web page seems to be
describing the contact between a sphere and a cylinder.  This problem
can actually be reduced to the first one, but I'd like to nail down
the scope of your Question before tackling it, so I can write it up
well the first time.

regards, mathtalk-ga

---

Clarification of Question by tfly-ga on 17 Jun 2004 16:36 PDT

I updated that page, and you are right I only want to consider a point
and a cylinder.  However if I wanted to consider a sphere and a
cylinder I would just add the sphere's radius to the cylinder's radius
and solve with a point and a larger cylinder, so disregard any
references to spheres.  I added a pic to my webpage which could be
right or wrong, basically to show how I am trying to solve the
problem.  http://www.geocities.com/t_F_l_y/PointOfCollisionCylinder.htm
basically all I came up with is that some equation should exist and a
minimization function may be necessary.  I am really stumped on this
one.  I read your explanation to focusedny-ga and I understand it, but
I think this problem is much harder.  Thanks Trevor

Hi, Trevor:

I think we can quickly see that the problem "ought" to be an easy one,
and then show how to solve it generally in a few concrete steps.

Imagine that we view the scene from a direction parallel to the
cylinder's axis.  The trajectory of the point then becomes a straight
line, which may or may not intersect the circular outline of the
cylinder.

Even the point's velocity, although reduced in magnitude, remains
"linear" under this end-on projection.  The task is then to develop a
mathematical formulation that accomplishes this.

Let's begin with a summary of the data initially given.  The point has
an initial position:

P(0) = [x0 y0 z0]

and a directional vector:

V = [x' y' z']

so that at time t:

P(t) = P(0) + t*V

The cylinder (assumed to have infinite extent) is defined by a center
axis and a radius r > 0.  The axis itself may be described somewhat
like the point's trajectory, by giving one point and a direction
through that point.  Such is the approach outlined on your linked Web
page:

C = [xc yc zc]

U = [xu yu zu]

Points on the cylinder's axis can then be given in parameterized form:

{ C + s*U | all real s }

Now our first step will be to reduce to the case that C is the origin.
 For this it suffices to subtract the given point C from P(0), thereby
translating all the points of our discussion in like manner.  Hence we
assume:

C = [0 0 0]

Now motion toward or away from the cylinder's axis is attributable to
a component of the point's trajectory perpendicular to the axis (which
now conveniently passes through the origin).  Therefore we
"orthogonalize" the point's trajectory with respect to the axis'
direction U:

Q(t) = P(t) - (<P(t),U>/|U|^2)*U

An explanation of this formula is in order.  First the notation:

< , >

means the inner product (or "dot product") of two vectors, and it
returns a scalar.  This orthogonalization formula would be written a
little more concisely if we assumed U has unit length.  In that case
we'd omit dividing by the length of U squared.

The effect of the above is to project the point's trajectory onto the
plane perpendicular to the cylinder's axis.  Such is the result of
"removing" any component parallel the cylinder's axis by the
calculation shown.

Now the calculation is "linear" in the sense that we can apply the
projection separately to the point's initial position P(0) and its
direction V:

Q(t) = P(t) - (<P(t),U>/|U|^2)*U

     = P(0) - (<P(0),U>/|U|^2)*U + t*(V - (<V,U>/|U|^2)*U)

It remains only to equate the distance of Q(t) from the cylinder's
axis to the desired radius r and determine whether there are two, one,
or no solutions to the resulting quadratic equation.

Note that because of our projection, the nearest point to Q(t) on the
cylinder's axis will be the origin C = [0 0 0].  So the distance is
conveniently just the norm of Q(t).

r^2 = |Q(t)|^2

    = |P(0) - (<P(0),U>/|U|^2)*U + t*(V - (<V,U>/|U|^2)*U)|^2

Since our goal is to find roots t for this quadratic equation, it
makes sense to take advantage of pre-computing the intermediate
vectors:

Q(0) = P(0) - (<P(0),U>/|U|^2)*U

 W   =   V  -  (<V,U>/|U|^2)*U

which are projections of the point's initial position and direction, resp.

If Q(0) = [xq yq zq] and W = [xw yw zw], then:


r^2 = |Q(t)|^2

    = |Q(0) + t*W|^2

    = (xq + t*xw)^2 + (yq + t*yw)^2 + (zq + t*zw)^2

Collecting "like" terms gives:

|W|^2 t^2 + 2 <Q,W> t + |Q(0)|^2 - r^2 = 0

For your purposes the sign of the constant term:

|Q(0)|^2 - r^2

is important, because if this is negative it means the point is
intially inside the cylinder; if it is zero, then on the cylinder's
surface.  The intended application is for the case when the point is
outside the cylinder, and hence when the constant term above is
positive.

Assuming this is so, then there are either two, one, or no positive
real roots.  You will of course identify the least positive real root
t as the time of impact (intersection).  Since you are obvious
familiar with the discriminant of a quadratic equation, I'll not
belabor these cases further.

Let me point out a "degenerate" case, however.  If the point's
trajectory were exactly parallel to the cylinder's axis, then W = [0 0
0] and the point moves neither toward or away from the cylinder,
remaining at the same distance from the cylinder as initially given.

For your convenience I can work out the formulas for Q(0) and W in
terms of the original input data as follows (and at the same time
restoring a general point C to our formulation):

P(0) = [(x0-xc) (y0-yc) (z0-zc)]

Q(0) = P(0) - (<P(0),U>/|U|^2)*U

     =   [(x0-xc) (y0-yc) (z0-zc)]

        ((x0-xc)xu + (y0-yc)yu + (z0-zc)zu)
     -  ----------------------------------- [xu yu zu]
               (xu^2 + yu^2 + zu^2)

 W   =   V  -  (<V,U>/|U|^2)*U

     =  [x' y' z']

         (x'xu + y'yu + z'zu)
     -  ---------------------- [xu yu zu]
         (xu^2 + yu^2 + zu^2)

Summary
=======

Using the original input data for the point's trajectory:

P(0) = [x0 y0 z0]

 V   = [x' y' z']

and the cylinder's axis:

 C   = [xc yc zc]  (any point on axis)

 U   = [xu yu zu]   (axis' direction)

and the cylinder's radius r, then compute vectors:

Q(0)   (projected intial point)

 W     (projected point's direction)

as above, and locate the positive real roots (if any) of:

|W|^2 t^2 + 2 <Q,W> t + |Q(0)|^2 - r^2 = 0

The least positive root t corresponds to the initial time of
intersection.  From t the point of "impact" P(t) = P(0) + t*V can be
found.

One final comment:  I agree that the case of a sphere impacting a
cylinder can be reduced to the case of a point impacting a cylinder by
the device of adding the sphere's radius to that of the cylinder.  The
time of impact is then the same for both geometries, but a little care
is needed to compute the point of impact for the sphere/cylinder as
opposed to the point/cylinder.  In this regard note that the point of
contact between the sphere and cylinder lies along a line from the
moving point/sphere's center perpendicular to the cylinder's axis.

regards, mathtalk-ga

---

Request for Answer Clarification by tfly-ga on 24 Jun 2004 16:31 PDT

Hello,

I have been looking at your answer and have tried to test it.  Can you
please tell me what to plug in for Q in the equation, |W|^2 t^2 + 2
<Q,W> t + |Q(0)|^2 - r^2 = 0.  I Calculated what B aka (Q dot W)
should be, based on A and C for some test points and got 837 and 825
for the two collision points.  So it looks like the equation will work
since A and C are correct due to the 2 values for B being close enough
to be the same value.  But I do not know how to calculate B without t.
 I tried Q(0) for Q ant it was completely wrong.
Trevor

---

Clarification of Answer by mathtalk-ga on 24 Jun 2004 16:54 PDT

Hi, Trevor:

My apology for the oversight/omission.  Actually it should be Q(0) in
the expression for "B":

|W|^2 t^2 + 2 <Q(0),W> t + |Q(0)|^2 - r^2 = 0

Q(0) is the "projection" of P(0) on the orthogonal complement of U,
after offsetting by the point C on the cylinder's axis.  Assuming C
has already been subtracted (or perhaps it was the origin to start
with), then:

Q(0) = P(0) - (<P(0)/U>/|U|^2)*U

Essentially this is one step of a Gram-Schmidt process.  Q(0) is the
"residual component" of P(0) after its component parallel to U is
removed.

Do you want to suggest some test points/lines and a radius to walk
through an example, or would you like me to set some values?

regards, mathtalk

---

Request for Answer Clarification by tfly-ga on 25 Jun 2004 17:59 PDT

Here are some test points, I have worked out T and I am positive the
solution is correct however, I arrived there with a guess and check
method that is unsatisfactory for what I am trying to do in my 3d
program.

POINT1:     (  -5.764,    25.056,    35.162)
POINTV:     ( 0.62525, -0.529917, -0.572931)Normalized

CYLINDER1:  (      0.0,      0.0,       0.0)
CYLINDERV:  ( 0.657522, -0.46255, -0.572931) Normalized

The Radius is 5.0.



Answer:

T1 = 29.991876  	S1 = 23.47
COLLIDE1:   ( 13.028,   9.154,  17.950)

T2 = 36.044133		S2 = 22.41
COLLIDE2:   ( 16.769,   5.991,  14.513)


To obtain T and S I used the following 2 equations however they seem
impossible to solve for T and S.

S = -1.0f*(CVx*(C0x - (P0x+T*PVx)) + CVy*(C0y - (P0y+T*PVy)) +
CVz*(C0z - (P0z+T*PVz)))
    --------------------------------------------------------------------------------
                           (CVx*CVx + CVy*CVy + CVz*CVz);

Radius^2 =((P0x+T*PVx) - (C0x+S*CVx))^2 + ((P0y+T*PVy) -
(C0y+S*CVy))^2 + ((P0z+T*PVz) - (C0z+S*CVz))^2

I could not get your equation to work if you could please try to use
these numbers or create your own numbers if you think these are
inaccurate or too confusing.

Heres the code of your equation:

Ws = (PVx*CVx + PVy*CVy + PVz*CVz) / (CVx*CVx + CVy*CVy + CVz*CVz);
Qs = (P0x*CVx + P0y*CVy + P0z*CVz) / (CVx*CVx + CVy*CVy + CVz*CVz);

Wx = Wx*Ws; Wy = Wy*Ws; Wz = Wz*Ws;
Qx = Qx*Qs; Qy = Qy*Qs; Qz = Qz*Qs;


double A = Wx*Wx + Wy*Wy + Wz*Wz;
double B = Wx*Qx + Wy*Qy + Wz*Qz;
double C = Qx*Qx + Qy*Qy + Qz*Qz - 25.0f;


T*T*A+T*B+C = 676 and 658 These should be zero I think.

---

Clarification of Answer by mathtalk-ga on 25 Jun 2004 19:50 PDT

Hi, Trevor:

I think you may have copied the last coordinate of POINTV into
CYLINDERV by mistake.  As it stands, POINTV is (nearly) normalized,
but not CYLINDERV.  The last two coordinates being equal led me to
suggest that the final coordinate of CYLINDERV should be -0.594737.

It's not necessary to normalize these vectors, so if there was a
simpler form before normalization, I'd be glad to use them.  However
if the above meets with your approval, I'll use that.

regards, mathtalk-ga

---

Clarification of Answer by mathtalk-ga on 25 Jun 2004 23:54 PDT

I worked through the value -0.594737 (for the final coordinate of
CYLINDERV) and found that there is no collision.  It occurs to me that
perhaps the value 0.594737 was instead intended.

Anyway, enough guess... I await your clarification.

regards, mathtalk-ga

---

Clarification of Answer by mathtalk-ga on 26 Jun 2004 09:22 PDT

Hi, Trevor:

It will probably be useful to see the worked example anyway, despite
lack of a collision point.  I think you will be able to understand a
couple of overlooked aspects more quickly this way.

The point C = 0 on the cylinder, so we don't need to offset
the base point of the "point" trajectory.

P(0) = [-5.764,    25.056,    35.162]

Now we find the component Q(0) of P(0) which is orthogonal
to the axial vector:

 V   = [0.657522, -0.46255, -0.594737]
 
Q(0) = P(0) - s*V

where s = <P(0),V>/<V,V>

        = -36.291752002/0.999999782153
        
        = -36.29176

Q(0) = P(0) + 36.29176 * V

     = [18.09863, 8.269246, 13.57795]

Next we make a similar projection of the axis' direction:

 U   = [0.62525, -0.529917, -0.572931]

onto the orthogonal complement of V, obtaining W:

 W   =  U  -  p*V
 
where p = <U,V>/<V,V>

        = 0.996972002997/0.999999782153
        
        = 0.99697222

 W   =  U  - 0.99697222 * V
 
     = [-0.0302812, -0.0687675, 0.0200053]

As a test of our work to this point, we should see that
V is orthogonal (within the limits of our arithmetic) to
both Q(0) and W:

<Q(0),V> = -0.00000159159

 <W,V>  =  -0.0000000401675

While the second inner product seems noticeably closer
to zero than the first, we should note that the magnitude
of W is somewhat less than that of Q(0).  All things
considered these show the vectors are as close to being
perpendicular to V as we should hope.

The reason that W is small is because U and V are close,
i.e. nearly parallel directions.  So the component of U
which is orthogonal to V is correspondingly small.

Now we turn to setting up our quadratic equation:

|W|^2 t^2 + 2 <Q(0),W> t + |Q(0)|^2 - r^2 = 0

where radius r = 5.

|W|^2 = 0.0060461

2 <Q(0),W> = 2 * (-0.84507) = -1.69014

|Q(0)|^2 = 580.30

Consequently:

A = |W|^2 = 0.0060461

B = 2 <Q(0),W> = -1.69014

C = |Q(0)|^2 - r^2 = 580.30 - 25 = 555.30

This gives a "discriminant" of about -10.573, so there
are no real roots.

I can't be sure from the presentation given of your work, but you may
have simply taken the projections of W and Q(0) onto V, rather than
their projections onto the "orthogonal complement" (the residue
obtained by subtracting the projections onto V from W and Q(0) resp.).
 Also it seems you overlooked a factor of 2 in the "B" coefficient.

regards, mathtalk-ga

---

Request for Answer Clarification by tfly-ga on 28 Jun 2004 18:01 PDT

Hello again,

I have tried to get your equation to work but I think it is incorrect.
 I am sorry for messing up my last post I did in fact mix up some of
the numbers.  After looking at what you did and following your numbers
I managed to get the same answer that there was no collision based on
the -10.573 discriminate however, after changing the nubmers to the
ones I originally wanted to try your equation did not seem to work,
meaning it did not satisfy the answer I calculated with guess and
check.  If you equation is correct and I am just missing something I
will give you a $30.00 tip but I am pretty sure it will not solve what
I am trying to solve.

Here are my final numbers please don't use these until I post another
clarifiacation and tell you that they are infact the correct points. 
I am going to first test them again in a 3d plotting program.

Cylinder0 ( 0, 0, 0)
CylinderV (32.884, 23.133, 29.744)

Point0 ( -5.764, 25.056, 35.162)
PointV ( 36.82, -31.206, -33.739)
Point0 + PointV = ( 31.056, -6.15, 1.423)

One Collision at (12.988, 9.1628, 17.979)
Two Collision at (16.7726, 5.9556, 14.5112)

Again please don't respond to this until I have posted another
clarification stating that these number are correct.

Thanks,
Trevor

---

Request for Answer Clarification by tfly-ga on 28 Jun 2004 22:15 PDT

Hello,

I double checked the points again.  They seem to work in the equation
I arrived at and they also work in the 3d plotting program.  When ever
you get a chance can you please take a look at them again.  Maybe your
equation does work and I am just missing something or maybe you did
not understand what I am trying to do, I am going to post a picture
after this and it should hopefully clear everything up, the points in
the picture should match the ones given in the previous post.  Please
just try to run these points though one last time.
Thanks,
Trevor

---

Request for Answer Clarification by tfly-ga on 28 Jun 2004 22:20 PDT

Hello,

Here is the screen shot from the 3d plotting program.

http://www.geocities.com/t_f_l_y/test.jpg

Hope this helps Thanks,
Trevor

---

Clarification of Answer by mathtalk-ga on 29 Jun 2004 21:36 PDT

Hi, Trevor:

I was unable to view the graphic you posted at the link above, but
that's not a problem.  I was able to duplicate your results with my
formulas.

You didn't explicitly say so, but the radius r = 5 is the same as
before.  I could tell by finding the orthogonal projections of the two
collision points on the plane perpendicular to U = CylinderV, the axis
of the cylinder.  Both of these projections are (within the accuracy
shown) 5 units distance from the origin C = Cylinder0.

Recall that in my formulas we have:

P(0) = Point0
 V   = PointV

and we find their "residues" orthogonal to U with these formulas:

Q(0) = P(0) - (<P(0),U>/<U,U>)*U

 W   =  V   -   (<V,U>/<U,U>)*U

Now |U|^2 = <U,U>, the dot of U with itself, is 2501.198681 according
to the Excel spreadsheet I whipped up.

We also have the dot products:

<P(0),U> = 1435.9356

<V,U>  = -514.632334

When the respective scalar multiples of U are subtracted, we get:

Q(0) = [-24.64267071, 11.7753684, 18.08600008]

 W   = [43.58602375, -26.44628624, -27.61904469]

Now the quadratic formula for "time" t to solve is this:

<W,W>t^2 + 2<Q(0),W>t + <Q(0),Q(0)> - r^2 = 0

Since:   <W,W>    =  3361.959152,
     2<Q(0),W>    = -3770.017677,
<Q(0),Q(0)> - r^2 =  1048.023919,

we can use the quadratic formula to determine the values of t that
will give us your two collision points:

t = (3770.017677 +/- SQRT(119378.856))/(2*3361.959152)

  = 0.509302027 _or_ 0.612073191

Adding these multiples of PointV to Point0 give, respectively:

P(0) + 0.509302027*V = [12.98850064, 9.16272094, 17.9786589]

P(0) + 0.612073191*V = [16.77253489, 5.955644006, 14.51126261]

Their agreement with the collision points you found is good, within a
least significant unit in figures cited.

If I can help to Clarify the computations in any way, please let me know!

regards, mathtalk-ga

Great Answer stated very clearly, Thanks for your patience.

For a similar Question and my Comments on it, see:

http://answers.google.com/answers/threadview?id=325446

In that case the issues concerned a point trajectory and a sphere or a
right circular cone.

regards, mathtalk-ga

If the problem is unclear this should help.

Given a Cylinder defined by some line and a radius and another line
for the particles linear path.  Find where the particle hits the
cylinder, the 2 points of intersection, the one point of intersection
or the case of no point of contact. Both the particle's path and the
sphere go to infinity in both directions.

Hope this helps,
TFLY