[comp.graphics] Summary point in volume.

thender@uxh.cso.uiuc.edu (10/24/90)

     Hello all.
Some time ago I posed the question of a point in a volume.  I would like
to thank all who responded to my question and apologize to those who 
requested a summary of the responses.  Things have just slowed down enough
to attempt a summary.  First, I would like to restate the problem.

     We needed to determine if a point was inside a volume or not.  The 
volume had several constraints which should simplify the problem.

       1.  The volume was always defined by 8 points, or it could be
           stated as a 6 sided volume.
       2.  We left the definition of the sides open to allow more freedom,
           i.e. we had four points that defined a side, but no equations
           or approximations at this time.
       3.  The sides could be concave or convex.  This was the tough part.

Now, that the problem is again stated, I would like to let all who read
this know that I am knew at this net stuff.  Therefore, if I do not
preform the curtesies appropriate to the net, achnowlegdements, etc.,
please let me know and I will try to correct the error in my ways.  So
at this point I would like to say that the following is other peoples work
and suggestions.  I'm not sure how most follow ups are done, but here goes.
Also, I hope that people who sent me mail do not mind being posted. I felt
I could either post without names or give credit where credit is due.  I feel
that the credit in this case is good therefore no harm, I hope.  If it is
not stated in the quotes I include, everyone is speaking for themselves not
their companies.


    1. The most common response was to check surface normals, this is great
       for convex volumes but not concave.  I didn't make myself clear in
       the first posting that the volume could have concave sides.

    2. A similar approach to the point in a cell routine could be used.

    3. The following is directly from mail I recieved from Dan Gordon at 
       the Computer Science Department at Texas A&M University,
       GIG'EM AGGIES.(I'm a former Aggie)

Quote start:
One simple solution depends on how exactly you get your distorted cube.
Under the following conditions, there is a simple solution:

1. There are 3 "well-behaved" functions f, g, h of 3 variables each,
such that the distorted cube is obtained from the unit cube by applying
the 3 functions. In other words, every point (x,y,z) in the space of the
unit cube is mapped to the point (f(x,y,z), g(x,y,z), h(x,y,z)) in the
space of the distorted cube.

2. The inverses of f, g, h are well-defined (and hopefully, easy to
compute).

Under these conditions, suppose you are given the point (X,Y,Z) in the
space of the distorted cube. Compute the point (f^-1(X,Y,Z), g^-1(X,Y,Z),
h^-1(X,Y,Z)) in the space of the unit cube, and check if that point is
inside the unit cube. The "well-behavedness" of the functions should be 
such that, mathematically, this method should be correct. It all depends 
on what functions you have.

Dan Gordon.
Quote end:


    4. The following was recieved from Sam Uselton.  Some parts
       are left out that were questions to me to clarify the problem.

Quote start:
There are several hacks ("approximations") you can use, but be aware that
they don't always give precisely correct answers.  
(1) One of the measures of CFD grid quality is that the faces of cells
should be close to planar when possible.  Therefore, you can find a best
fit plane and use it and be right most of the time.
(2) You can take the convex hull, which will result in splitting each 
quadrilateral face into two triangles.  Adjacent cells will have a small
overlap in volumes.  
(3) You can make essentially arbitrary (or other heuristic) splits of 
nonplanar faces (into 2 triangles) and propagate the choice to the
adjacent cell that shares that face.
(4) You can use a bilinear interpolation (linked, not two independent
linear interpolations) to "define" the surface.

Sam Uselton		uselton@nas.nasa.gov
employed by CSC		(415) 604-3985
working for NASA
speaking for myself
Quote end:


    4. The following came from Joseph from the University of 
       New Hampshire.  This seemed like one of the best and surest(sp)
       ways of doing this.

Quote start;
	I have come up with a fairly efficient solution to the point
in cube problem you posted in comp.graphics. Following is a copy of
the usenet article I posted that describes the solution. Hope it
works, and hope you find it useful.

I shall describe a solution to the "point in distorted cube"
problem that uses only "point in front of plane" tests
(at most 18 of them). I have not attempted to prove if this
is faster than the general "point in polygonal volume algorithm",
but it seems to be.

I consider only the interpretation that each distorted face is represented
by 2 triangles. Further, I assume the faces do not intersect each other.

In the following, by "point in front of
the plane (ax+by+cz>0) of a triangle" I mean in the direction of
the outward facing normal.

The algorithm considers each of the 6 faces in turn: if 
the point is behind each of these distorted faces then the point is inside
the cube else it it not.

The following pseudocode determines if a point is behind a
distorted face (potentially inside the volume) or in front
(definitely not in the volume):

If (a face is convex)
then
    if(point is BEHIND both the planes of the two triangles that 
       make up the face)
    then
	the point is behind the distorted face
    else
	the point is in front of the distorted face
else /* the face is concave */
    if(point is IN FRONT OF both the planes corresponding to the
       two triangles that make up the face)
    then
	the point is  in front of the distorted face
    else
	the point is behind the distorted face

We have to determine if a face is convex. Consider the following
face:

      A +-------+ B      It is divided into triangles ABC and ACD.
	| .     |        If point B is behind the plane of triangle
	|   .   |        ACD then the face is convex, else it is
	|     . |        concave.
      D +-------+ C

Notes: (1) The key insight is changing the test from "behind the face"
	   to "in front of the face" if the face is concave.
       (2) The algorithm works even if the faces are planar.
       (3) There are at most three "point in front of plane" tests
	   per face, so at most 18 such tests overall.
       (4) Use a bounding box to quickly determine if some points
	   are outside the volume.
       (5) If the problem is to determine if MANY points are inside
	   the SAME distorted cube, the code can be optimized in
	   several fairly obvious ways (for example, the tests for
	   convexity need be done only once per face).
       (6) I learned of this problem just this morning, and haven't
	   tested anything yet, so I could be making some big
	   mistake. Let me know!

Joseph

PS. If the cube is not very distorted, it may be worth it to
compute a "anti-bounding" box, i.e., to compute the smallest
box that lies entirely *within* the  distorted cube, and test using this
box first.
--------------------------------------------------------------------
Joseph M. Joy                        Internet: jmj@unh.edu
University of New Hampshire          Tel:      (603)862-4335 (office)
Dept. of Computer Science                      (603)742-9449 (home)
Durham, NH 03824                                
Quote end:


    5. There was a fifth response that was unique, however I have miss
       placed it and would like to apologize to the sender whoever it is
       out there.

In case anyone is interested our first try at this will probably go something
like this.  This is in two-d but can be extended to 3-d with trilinear
interpolations instead of bilinear.
Define 2 equations:

     x = ax + bx*[xi] + cx*[eta] +dx*[xi]*[eta]
     y = ay + by*[yi] + cy*[eta] +dy*[yi]*[eta]
  where [xi] and [eta] go from 0 to 1.  For 2-d we have four points
and four unknowns, ax,bx,cx,....  Therefore, we know what ax,bx,..
from the corners of the cell.  Now, given a new (x,y) solve for ([xi],[eta]).
If [xi] or [eta] > 1 or [xi] or [eta] < 0 then the point is outside the cell.
This is not exact and could give that some points are in that are not.
However, once we found the point we were going to use the same interpolation
to go to another domain, so we really didn't lose much.  This, I believe
is in line with what was discussed in  #3. above.

Thank you all again for the responses.  And again I hope I haven't made to
many unforgivable net mistakes.

Thanks
Todd Henderson
U. of Illinois
CFD Lab