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 endon 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 precomputing 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) = [(x0xc) (y0yc) (z0zc)]
Q(0) = P(0)  (<P(0),U>/U^2)*U
= [(x0xc) (y0yc) (z0zc)]
((x0xc)xu + (y0yc)yu + (z0zc)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, mathtalkga 
Request for Answer Clarification by
tflyga
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
mathtalkga
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 GramSchmidt 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
tflyga
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
mathtalkga
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, mathtalkga

Clarification of Answer by
mathtalkga
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, mathtalkga

Clarification of Answer by
mathtalkga
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, mathtalkga

Request for Answer Clarification by
tflyga
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
tflyga
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
tflyga
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
mathtalkga
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, mathtalkga
