Ray and Object Intersections: Cylinder

by Otto Wilke otto_wilke@hotmail.com

I am vaguely aware that graphics is normally done with vector operations, generic

solids positioned at the origin, and transformation matrices to move rays to and fro.

I thought it would be interesting to use rectangular coordinates and objects located

anywhere in space and oriented in any direction.

This is one of four files covering the plane, the sphere, the cylinder, and the cone.

Let C be a right circular cylinder having radius r and positioned some place in space and oriented in

some direction. Let P1(x1,y1,z1) and P2(x2,y2,z2) be the centers of the circular ends. Let P(x,y,z)

be some point on the cylinder. Let d1 be the distance from P1 to P2, d2 be the distance from P to P1,

and d3 be the distance from P to P2. restart; d1:=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2); d2:=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2); d3:=sqrt((x2-x)^2+(y2-y)^2+(z2-z)^2);

The area of the triangle P P1 P2 is

Area=1/2*r*d1;AreaSquared1:=expand((1/2*r*d1)^2);

Let s be the semi-perimeter of triangle P P1 P2. By Heron's formula the area squared is also

s:=(d1+d2+d3)/2;AreaSquared2:=expand((d1+d2+d3)/2*(s-d1)* (s-d2)*(s-d3));

sort(AreaSquared1-AreaSquared2,[x,y,z])=0;sort(4*%,[x,y,z]);

This last equation gives a point on an infinitely long cylinder, however, the distance from the center of the cylinder to a point on the cylinder must be limited. When P3 is on the cylinder at one end, let the distance from the midpoint of P1 and P2 to P3 be d4.

d4:=sqrt(r^2+(x1-(x1+x2)/2)^2+(y1-(y1+y2)/2)^2+(z1-(z1+z2)/2)^2);

The distance from the midpoint of the cylinder to any P must be less that or equal to d4.

sqrt(r^2+(x-(x1+x2)/2)^2+(y-(y1+y2)/2)^2+(z-(z1+z2)/2)^2)<=d4;

Let L be a line on point P4(x4,y4,z4) and let L have direction numbers a, b, and c. The parametric equations for L are

x=x4+a*t;y=y4+b*t;z=z4+c*t;

Now that we have equations for points on a line and a cylinder and one inequality, we'll ask Maple to solve the simlutaneous equations and inequality.

solve({x=x4+a*t,y=y4+b*t,z=z4+c*t,

sqrt(r^2+(x-(x1+x2)/2)^2+(y-(y1+y2)/2)^2+(z-(z1+z2)/2)^2)<=d4, sort(AreaSquared1-AreaSquared2,[x,y,z])=0},{x,y,z,t});

Maple fails to solve the set. Let's pick P4 as (0,0,0) and a=1, b=2, c=1 and r=1 and P1 as (1,1,0) and P2 as (2,2,3).

x4:=0;y4:=0;z4:=0;a:=1;b:=2;c:=1;x1:=1;y1:=1;z1:=0; x2:=2;y2:=2;z2:=3;r:=1;

Try again.

solve({x=x4+a*t,y=y4+b*t,z=z4+c*t,

sqrt(r^2+(x-(x1+x2)/2)^2+(y-(y1+y2)/2)^2+(z-(z1+z2)/2)^2)<=d4, sort(AreaSquared1-AreaSquared2,[x,y,z])=0},{x,y,z,t});

evalf(allvalues(%));

with(plots):

plot1:=arrow(<0,0,0>, <1.206622805,2.413245610,1.206622805>, width=[0.1, relative], head_length=[0.2, relative], color=red):

plot2:=implicitplot3d(sort(AreaSquared1-AreaSquared2,[x,y,z])=0, x=-1..4, y=-1..4, z=-1..4, grid=[13,13,13]):

plot3:=pointplot3d({[1.206622805,2.413245610,1.206622805] },axes=normal,symbol=box,symbolsize=50,color=red);

display({plot1,plot2,plot3});

Now let's consider the circular ends. The direction numbers of the axis of the cylinder are

A=x2-x1, B=y2-y1, and C=z2-z1. The plane A x + B y + C z +D = 0 is normal to the

axis of the cylinder. P2 is on that plane.

(x2-x1)*x2+(y2-y1)*y2+(z2-z1)*z2 +D2 = 0;

D2:=solve((x2-x1)*x2+(y2-y1)*y2+(z2-z1)*z2 +D2 = 0,D2);

The equation of the plane of the circular end on which P2 lies is

(x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D2 = 0;

The distance d5 from some point P on the circular end to the center of the circle must be less

than or equal to the radius of the cylinder.

d5:=sqrt((x-x2)^2+(y-y2)^2+(z-z2)^2);

Now that we have the equation of the ray, the equation of the circular end at P2, and the

necessary inequality, let's check for an intersection.

solve({x=x4+a*t,y=y4+b*t,z=z4+c*t,

d5<=r, (x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D2 = 0},{x,y,z,t});

No answer. Let's try the other end.

D1:=solve((x2-x1)*x1+(y2-y1)*y1+(z2-z1)*z1 +D1 = 0,D1);

(x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D1 = 0;

d6:=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);

solve({x=x4+a*t,y=y4+b*t,z=z4+c*t,

d6<=r, (x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D1 = 0},{x,y,z,t});

plotPlane1:=implicitplot3d((x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D1 = 0, x=-1..4, y=-1..4, z=-1..4, grid=[13,13,13]): plotPlane2:=implicitplot3d((x2-x1)*x+(y2-y1)*y+(z2-z1)*z +D2 = 0, x=-1..4, y=-1..4, z=-1..4, grid=[13,13,13]): plot4:=pointplot3d({[1/3.,2/3.,1/3.]},axes=normal,symbol=box, symbolsize=50,color=red):

display({plot1,plot2,plot3,plot4,plotPlane1,plotPlane2});