Path: utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!wuarchive!udel!princeton!phoenix!markv
From: markv@phoenix.Princeton.EDU (Mark T Vandewettering)
Newsgroups: comp.graphics
Subject: Re: Point within cylinder problem
Message-ID: <10423@phoenix.Princeton.EDU>
Date: 14 Sep 89 18:49:14 GMT
References: <7967@medusa.cs.purdue.edu>
Reply-To: markv@phoenix.Princeton.EDU (Mark T Vandewettering)
Organization: Princeton University, NJ
Lines: 96

>	Even more helpful would be something that would determine if a
>point is within a cylinder with a hemisphere at each end.  This can be
>done by checking the point against a spere at each endpoint, but if
>there is a way to get this "free" I'd prefer it.

ASCII, and you shall receive.

Formula time it is.

The expanded problem was to determine if a point is inside a cylinder capped
by hemispheres.  We are given the following definitions.

	Point P, the point that we are testing.

	Base and Apex, the points that define the base and apex of the cylinder.

	rad, the radius of the cylinder

I adopt the convention that vectors be capitalized, and scalars be represented
by lowercase names.  We know that point is inside if it is inside either 
sphere at the end, 

	in other words, Length(P - Base) < radius
			or Length(P - Apex) < radius.

	
	Length requires taking a square root, we can square both sides 
	and rewrite the expressions as
			
			DotProduct(P-Base, P-Base) < radius * radius
			and
			DotProduct(P-Apex, P-Apex) < radius * radius

Total maximum cost, assuming that we precompute radius * radius, is
3 subtracts, 3 multiplies and 2 adds for each endpoint, for a total 
of 6 subtracts, 6 multiplies and 4 adds worst case.

The other case is that the point is inside the cylinder.  We can determine
this by projecting the point down onto the base plane, and determining
whether it is within radius of the base point.

We can do some early trimming by precomputing some stuff.  We first of all
need the plane containing Base, and that has a normal along Apex - Base.

Piece of cake, let N = Normalize(Apex - Base), then the equation for the 
plane is 
	a x + b y + c z + d = 0,
		where	N = <a, b, c>
			d = - VecDot(N, Base) ;

The projection P' of P onto the base plane is in the direction of the normal 
to the plane.  We can express this as the set of equations.

	VecDot(P', N) + d = 0
	and P' = P + t N ;

and then attempting to solve for t.  If t is < 0, or t > Length(Apex - Base),
then we can abort the test.  Merging the two steps above, we have

	VecDot(P + t N, N) + d = 0
or	
	VecDot(P, N) + VecDot(t N, N) + d = 0

	VecDot(t N, N) = - (d + VecDot(P, N))

	t VecDot(N, N) = - (d + VecDot(P, N))

	t = - (d + VecDot(P, N))
	      ------------------
	      VecDot(N, N)


Computing t, assuming 1.0 / VecDot(N, N) is precomputed, can be done with
3 adds and 4 multiplies.  Two comparisons are needed to ensure t is between
0 and Length(Apex - Base).   Assuming that it passes, we can then compute 
P' and see if its within radius of Base.

	P' = P + t N ;

	if VecDot(P' - Base) < radius * radius
		then we are inside
	else fail

The above takes 10 adds, 6 multiplies and one comparison.


Total for the whole mess.
	23 adds
	12 multiplies
	5 comparisons


Note:
	This is all off the top of my head, so be careful and check 
	everything.  I do make mistakes on this kind of stuff, although
	I believe that the methodology is sound.