Subject: Summary: Points inside/outside of a polyhedron

From: Robert Fraczkiewicz <[email-masked]>
Date: Mon, 22 Apr 1996 11:23:15 +0059 (CDT)

Dear Netters,

Thank you all who replied to my posting. In a few days I have received 
more great ideas than I would have generated in months! Once again, 
it confirms how important for us is to keep the CCL alive.

Cheers,
Robert Fraczkiewicz

=========================================================================
-------------------------Original posting--------------------------------

Suppose you have given N arbitrary points (in 3D space) that are vertices
of a polyhedron: (x1,y1,z1),.....,(xN,yN,zN). What is the fastest 
algorithm for determining if a given point (x,y,z) is inside or outside of
this polyhedron? (Calculating distances to every vertex is not a good 
idea). I will post summary to the net.

Best regards,

Robert Fraczkiewicz
UT Medical Branch
--------------------------------------------------------------------------

================================Summary===================================
--------------------------------------------------------------------------
Dear Dr. Fraczkiewicz,

I guess the first thing I'd do is see whether the point lay within
a sphere just barely inerior to the polyhedron, by calculating the
square radius of the point (rsq=x^2+y^2+z^2) and seeing whether it
is less than R^2, where R is the radius of the sphere. All other
points may be inside or outside, but at the least this will decrease
the number of points to be tested substantially. You could also
construct a sphere which is just barely exterior with radius rho 
to see whether rsq > rho^2 . All remaining points will lie
within a thin layer. These remaining points can be subjected to
some (more numerically intensive) clever test, which you may 
come up with or which someone
else from the CCL may suggest.
Note that I did not calculate r, in order to avoid the sqrt function.

Hope this is helpful,
rqt

************************************************************************
Check out the Molecular Monte Carlo home page!
http://www.cooper/edu/engineering/chemechem/monte.html
************************************************************************
Robert Q. Topper                       email:   [email-masked]
Department of Chemistry                phone:   (212) 353-4378
The Cooper Union                       FAX:     (212) 353-4341 
51 Astor Place                         subway:  take the 6 to Astor Place 
New York, NY 10003 USA                          or the N/R to 8th St/NYU
http://www.cooper/edu/engineering/chemechem/depts_info/topper.html
************************************************************************

--------------------------------------------------------------------------

    Robert,

      A method I have used is to pass a line through the test point that is 
  parallel to one of the coordinate axes (say X), and then count the number of 
  those intersections of the line with polyhedron faces that have greater X-
  coordinate than the test point. If the number of such intersections is odd,
  than the test point is in the interior, if even it is in the exterior. 

      If you save the min. and max. y and z coordinates of every vertex in an 
  array, you can quickly eliminate most faces from consideration for a given
  test point - for those faces that can intersect the line, you apply a test
  (appropriate for the number of sides of the face) to determine if the line
  through the test point intersects the interior of the face.

       Randy

All opinions expressed here are mine, not my employer's

///////////////////////////////////////////////////////////////////////// 
\\ Randy J. Zauhar, PhD             | E-mail: [email-masked]        //
\\ Tripos, Inc.                     |       : [email-masked]           //
\\ 1699 S. Hanley Rd., Suite 303    |  Phone: (314) 647-1099 Ext. 3382 //
\\ St. Louis, MO 63144              |                                  //
/////////////////////////////////////////////////////////////////////////
**                                                                     **
**  "If you have conceptions of things that you can have no conception **
**   of, then the conception and the thing appear to co-incide."       **
**   --- C.G. Jung                                                     **
*************************************************************************

--------------------------------------------------------------------------

Hi Robert,
  There's a really nice convex hull algorithm that runs in ~nlog(n) time from
the Univ. of Minn. Geometry Center.  You can get it over the WWW via:
http://www.geom.umn.edu/software/download/
                                             --Mike Colvin
                                               Sandia National Labs.

---------------------------------------------------------------------------

Dear Robert,

The algorithm I would employ would do the following, starting with the 
set of points Pi(x,y,z) (x = 1 to n) and O(x,y) is the point we wnt to test.
1 - define average A(x,z,y) of set of points Pi(x,z,y).
2 - define vector AO from A through O.
3 - find three closest points to this vector that are on the positive 
end (test with dot product of vectors).  Let's call them Pa, Pb, and Pc.
4 - calculate point of intersection of the plane defined by Pa,Pb&Pc and 
the vector AO.  Lets call it P.
5 - now compare distances:
	if |AO| < |AP| then O is within the polyhedron.

Perhaps there are better ways, but this is the easiest I can think of 
without resorting to pen and paper.

Cheers, Craig

"Imagine if every Thursday your shoes exploded if you tied them the
 usual way.  This happens to us all the time with computers, and nobody
 thinks of complaining."
                -- Jeff Raskin, interviewed in Doctor Dobb's Journal

Craig Taverner
Structural Chemistry, University of the Witwatersrand, South Africa
tel:   +27-11-716-2290                fax: +27-11-716-3826
email: [email-masked]
www:   http://www.gh.wits.ac.za/craig

------------------------------------------------------------------------------

I would be tempted to recommend that you calculate the following:
  say X-x1,..., X-xN but not sqrt ( (X-x1) **2 + (Y-y1) **2 + (Z-z1)**2)

which would allow you to quickly select which points are inside as opposed
to outside.
Should any of these preliminary calculations have a hint that they may indeed
be inside, then you would do a full distance calculation.
In other words, testing differences are much faster than calculating
the full distance
I hope this helps
Carlos

            **********************************************
            *  Dr. Carlos H. Faerman                     *
            *  Biochemistry, Molecular and Cell Biology  *
            *  Cornell University                        *
            *  221 Biotechnology Building                *
            *  Ithaca, NY 14853                          *
            *  USA                                       *
            *  ........................................  *
	    *  [email-masked]           *
	    *  Phone: (607)255-8432                      *
            **********************************************

-------------------------------------------------------------------------------

Dear Robert,

  Use the Gauss theorem:
1) find all faces and normals
2) calculate an integral of  dS/R^2 for each face (R is a radius-vector from
your point to an surface element of the face. This integral can be calculated
for any point and an oriented triangle.
3) sum contributions for all the faces
4) if the sum is 0 it is outside, otherwise should be 4PI which means that it
is inside.

Ruben Abagyan

-------------------------------------------------------------------------------

I am sure that the graphics people have long found a solution
to this problem, but just for fun let's try to make up something:

I'll assume that "fast" refers to repeating the operation for many
different test points and a fixed polyhedron. If not, forget this
idea.

First, make a list of all the sides of the polyhedron. Generate the
normal form of the planes they are in with the normal vectors pointing
inward. Once this is done, all you have to do for testing a specific
point is to insert the coordinates of your test point into all the
equations (three multiplications and three additions per side) and
check if all the results are positive. If so, you are inside the
polyhedron, otherwise outside.

If that's not fast enough, there are two ways to reduce the number
of sides to be checked:

1) Find a sphere that encloses the polyhedron. All test points that
   are outside this sphere don't have to be checked in detail.

2) Divide the polyhedron into two halves by a suitable dividing
   plane. Check first in which half you are, and then check onyl
   the sides in this half.

-------------------------------------------------------------------------------
Konrad Hinsen                     | E-Mail: [email-masked]
Departement de chimie             | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. Centre-Ville     | Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase experimentale)
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

This is a common problem in computer graphics so you may want to check
some of the standard texts, such as Hearn and Baker.

If your test points might be anywhere you can quickly eliminate some
by checking if they are outside a bounding sphere.  Depending on the
shape of the polyhedron some other space partitioning may be useful.

Dave Heisterberg

-------------------------------------------------------------------------------

Just an idea. It's very efficient, but not sure it's right.

First transform the current coordinate system to one with the given point 
(x,y,z) as the origin. I.e. subtrate (xi,yi,zi), where i=1,...,N, by 
(x,y,z). This can be done very fast.

Second look at the sign of the transformed coordinates of all vertices. 
If all xi's, where i=1,...,N, have the same sign, the point is outside 
of the polyhedron. The statement holds for all yi's or all zi's. So if any 
one of the three conditions holds, the point will be outside. Otherwise if 
none of the conditions hold, it's inside.

Good Luck!

Rui Luo
CARB, UMBI
04/15

------------------------------------------------------------------------------

From: [email-masked]
Subject: Computational Geometry in C
To: robert

This is a good book by Joseph O'Rourke - provides algorithms for 2D and
3D geometry calculations - ftp site is provided.  Id suggest contacting
O'Rourke at [email-masked] if you cant find the book 
(Cambridge Univ. Press, 1994, isbn 0-521-445922 paper, 346 pages).

------------------------------------------------------------------------------

> Just an idea. It's very efficient, but not sure it's right.
> 
> First transform the current coordinate system to one with the given point 
> (x,y,z) as the origin. I.e. subtrate (xi,yi,zi), where i=1,...,N, by 
> (x,y,z). This can be done very fast.
> 
> Second look at the sign of the transformed coordinates of all vertices. 
> If all xi's, where i=1,...,N, have the same sign, the point is outside 
> of the polyhedron. The statement holds for all yi's or all zi's. So if any 
> one of the three conditions holds, the point will be outside. Otherwise if 
> none of the conditions hold, it's inside.

What you are effectively checking is whether the point is inside or
outside a box enclosing the polyhedron. But there are points inside
that box which are outside the polyhedron itself.

-------------------------------------------------------------------------------
Konrad Hinsen                     | E-Mail: [email-masked]
Departement de chimie             | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. Centre-Ville     | Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase experimentale)
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------

>>>>> "Rui" == Rui Luo <[email-masked]> writes:

    Rui> Second look at the sign of the transformed coordinates of all
    Rui> vertices.  If all xi's, where i=1,...,N, have the same sign,
    Rui> the point is outside of the polyhedron. The statement holds
    Rui> for all yi's or all zi's. So if any one of the three
    Rui> conditions holds, the point will be outside. Otherwise if
    Rui> none of the conditions hold, it's inside.

Your last statement is incorrect.  Consider the 2D case:

    ----
   /   /
 +/   /
 /   /
/   /
---

If the point '+' is the origin, then all of the Xi's do not have the
same sign, nor do all of the Yi's, so none of the above conditions
hold, yet the point is still outside the polygon.  Your method is
another implementation of a bounding box.  Usually, however, it is
done by simply finding the maximum and minimum value of each
coordinate, x, y, and z, and checking to see whether the point's
coordinates lie within the extrema of the polyhedrons coordinates.  If
most of the points of interest likely lie outside the polyhedron this
is a quick initial screening method.

-Harry
______________________________________________________________________
Avant!			 Harry C. Johnson IV      [email-masked]
1101 Slater Rd		  Software Engineer       Phone: (919)941-6726
RTP, NC 27703		                            FAX: (919)941-6700
______________________________________________________________________

------------------------------------------------------------------------

On Mon, 15 Apr 1996, Rui Luo wrote:
> 
> Just an idea. It's very efficient, but not sure it's right.
> 
> First transform the current coordinate system to one with the given point 
> (x,y,z) as the origin. I.e. subtrate (xi,yi,zi), where i=1,...,N, by 
> (x,y,z). This can be done very fast.
> 
> Second look at the sign of the transformed coordinates of all vertices. 
> If all xi's, where i=1,...,N, have the same sign, the point is outside 
> of the polyhedron. The statement holds for all yi's or all zi's. So if any 
> one of the three conditions holds, the point will be outside. Otherwise if 
> none of the conditions hold, it's inside.
> 
> 
It seems to me that this test will be reliable if the polyhedron is known 
to be convex, but other wise might be easily fooled:

       ___
       \  \
        \  \
         \  \
        *|  |
         /  /
        /  /
       /  /
      ----

Steve Williams
Chemistry, ASU
Boone, NC USA

snip

----------------------------------------------------------------------------

> Just an idea. It's very efficient, but not sure it's right.
> 
> First transform the current coordinate system to one with the given point 
> (x,y,z) as the origin. I.e. subtrate (xi,yi,zi), where i=1,...,N, by 
> (x,y,z). This can be done very fast.
> 
> Second look at the sign of the transformed coordinates of all vertices. 
> If all xi's, where i=1,...,N, have the same sign, the point is outside 
> of the polyhedron. The statement holds for all yi's or all zi's. So if any 
> one of the three conditions holds, the point will be outside. Otherwise if 
> none of the conditions hold, it's inside.
> 
> 
> Good Luck!
> 
> Rui Luo
> 
While this is a necessary condition for the point to be inside
the polyhedron, it isn't sufficient. A polyhedron that has
a concave portion may span both side of the point without
actually enclosing it. 

=============================================================================
Carlos L. Simmerling, Ph.D.                 
Department of Pharmaceutical Chemistry       Phone: (415) 476-7985 
Box 0446, Room S-926                          Fax:  (415) 476-0688
University of California 
San Francisco, CA 94143                    E-mail: [email-masked]
=============================================================================

-----------------------------------------------------------------------------

Dear Dr. Fraczkiewicz :

I will suppose that your polyhedron is convex, and made of points (Ai),i=1,n.

You can first determine each face of the polyhedron. Any given triad (Ai,Aj,Ak)
of points defines a plane for which a normal vector is the vectorial product
Nijk=AiAj^AiAk. If for any other point Al the scalar product Nijk.AiAl has
always the same sign, this means that all other points Al (l<>i,j,k) are on the
same side of the plane defined by (Ai,Aj,Ak), and thus these three points define
a facet of the convex polyhedron. Knowing the sign of Nijk.AiAl, you can easily
define a normal vector N'ijk=Nijk or -Nijk pointing out of the polyhedron.

Then, if you want to determine if a point B is in or out of the polyhedron,
simply compute the scalar product AiB.N'ijk for all faces (Ai,Aj,Ak). If
this product is negative for any given face (Ai,Aj,Ak), B is inside.
Otherwise you can stop the algorithm whenever this product becomes positive :
B is outside.

If your polyhedron is not convex, the solution is likely to be more complex...

Best regards,

Nicolas Froloff.

------------------------------------------------------------------------------

Take any ray beginning at the given point. Count the number of
points of intersection of this ray with polyhedron surface. If it is
even, then the point is outside the polyhedron, if it is odd, then inside.
One precaution: the ray may slide on some edges of the polyhedron.

Best regards,
Alexander Golbraikh

------------------------------------------------------------------------
Dr. Alexander Golbraikh                 e-mail: [email-masked]
Latvian Institute of Organic Synthesis
Aizkraukles st. 21                      Phone:  +(371)-7-552-354
Riga, LV-1006                           FAX:    +(371)-7-821-038
LATVIA
------------------------------------------------------------------------

------------------------------------------------------------------------------

Alexander Golbraikh wrote:
: 
: 
: 
: Take any ray beginning at the given point. Count the number of
: points of intersection of this ray with polyhedron surface. If it is
: even, then the point is outside the polyhedron, if it is odd, then inside.
: One precaution: the ray may slide on some edges of the polyhedron.
: 
: Best regards,
: Alexander Golbraikh
: 
: 
       ^ /
       |/
       <     As you say, the ray may slide. The "o" point at the left is
       |\    outside the "polyhedron" yet passes through the polyhedron
       | \   surface only once if you send the ray only through the vertices
       o  \  defined by... say the Calphas of a protein.

	-Brian
-- 
=============================================================================
|   .---------.| Brian W. Beck      |    E-mail Addresses:                  |
|/\ |         || Biochem/Biophysics |        [email-masked]        |
|| \\     WSU || Washington St. Univ|   [email-masked]                  |
|\  -        *|| 206 Fulmer         |  URL  http://elmo.chem.wsu.edu/~brian |
| |           || Pullman, WA, USA   |    VOICE    (509) 335-4083            |
| \___________||       99164-4660   |      FAX    (509) 335-9688            |
=============================================================================

------------------------------------------------------------------------------

There are a few quick checks:

If the x coords of all vertices are all less than or all greater than the
x coord of the pt, then the point must be outside.  Similarly for y and z.

Find the center of an enclosing sphere for the polyhedron.  If the point
is outside the enclosing sphere, then it is outside the polyhedron.

Find the center of a maximal enclosed sphere for the polyhedron.  If the
point is inside the enclosed sphere, then it must be inside the polyhedron.

There is also the scheme of shooting out a ray from the point.  If it
crosses the polyhedron surface an odd number of times, then it is inside,
if an even number of times, then it is outside.
-Todd Wipke
 Molecular Engineering Laboratory
 University of California
 Santa Cruz, CA
 [email-masked]

------------------------------------------------------------------------------

You can find a long answer in the comp.graphics.algorithm FAQ archived at:

 http://www.cis.ohio-state.edu/hypertext/faq/usenet/graphics/algorithms-faq/faq.html

In any case:

 1) If your polyhedron is convex, it is very efficient to have determined
    the Ax+By+Cz+D=0 plane equation for each face with the sign convention,
    for instance, that all points lying in the 'inside' hemi-space
    are positive (i.e. A x0 + B y0 + C z0 + D > 0). In that case, checking
    if a new point is inside or outside is just a question of substituting
    the cartesian coordinates of the point in the plane equations of
    every face and testing whether all signs are positive or there is
    one or more negative.

 2) No matter if your polyhedron is convex or not, an inside point is such that
    every ray that begins in the point will intersect the polyhedron faces
    an odd number of times. Given a new point, you can always decide
    an arbitrary direction and check for the intersections of the line
    with the polyhedron faces. If there is zero or an even number of
    intersections, your point is outside the polyhedron, otherwise it is
    inside.

Hope this helps you.

    Sincerely,
          Victor Lua~na

--
    HomePage   http://www.uniovi.es/~quimica.fisica/qcg/vlc/luana.html
+--------------------------------------------+  +---^---/    /
!               Victor Lua~na                !  |   ~       / Just in case
! Departamento de Quimica Fisica y Analitica !  |           | you don't
! Universidad de Oviedo, 33006-Oviedo, Spain ! <            / remember
!  e-mail: [email-masked]  !  |          /  where Oviedo
!          phone: (34)-8-5103491             !  |____  ___/   is  ;-)
!           fax: (34)-8-5103125              !       \/
+--------------------------------------------+ 

------------------------------------------------------------------------------

OK, here is another one:

1) Make a Delone tesselation for the polyhedron vortices.

If it is convex, all Delone simplexes (tetrahedra) are inside
the polyhedron.
In this case if your point is inside of any of the simplexes (simplex is
a tetrahedron - so check is simple (center of tetrahedron is always
inside, connect center with your point and see wether the border is crossed)
it is inside polyhedron.

If polyhedron is not convex, there will be Delone simplexes connecting
vortices which are entirely outside the polyhedron, and they can be excluded

I believe it should work
Mike

Michael Kotelyanskii PhD
Department of Chemical Engineering
University of Delaware                          [email-masked]
Newark DE 19716                                    

------------------------------------------------------------------------------

If I remember correctly, this is a standard computer science algorithm
which is described in Robert Sedgewick's book "Algorithms in C."
Normally I'd check first before sending this, but I'm at home, I'm leaving town
tomorrow, and the book is currently at my office.

Wilfred Tang
([email-masked])

------------------------------------------------------------------------------

On Tue, 16 Apr 1996, Brian W. Beck wrote:

> Alexander Golbraikh wrote:
> :
> :
> :
> : Take any ray beginning at the given point. Count the number of
> : points of intersection of this ray with polyhedron surface. If it is
> : even, then the point is outside the polyhedron, if it is odd, then inside.
> : One precaution: the ray may slide on some edges of the polyhedron.
> :
> : Best regards,
> : Alexander Golbraikh
> :
> :
>        ^ /
>        |/
>        <     As you say, the ray may slide. The "o" point at the left is
>        |\    outside the "polyhedron" yet passes through the polyhedron
>        | \   surface only once if you send the ray only through the vertices
>        o  \  defined by... say the Calphas of a protein.
>
> 	-Brian
> --

Thank you for your remark. I mean the point must be taken into
account in case the ray crosses the surface of the polyhedron, i.e. when
any whatever small segment of the ray containing this point
contains also points both in the inside and in the outside of the
polyhedron. Both cases could be easily discriminated. You may also ask
what to do if the given point belongs to the polyhedron surface. In this
case all depends on which polyhedron we consider: the closed or the open
one. (If we consider that the border (the surface) of the polyhedron
belongs to the polyhedron, it is closed). If it is closed, the point
belongs to the polyhedron (but isn't its inner point), if it is open, the
point doesn't belong to the polyhedron.

Alexander

+----------------------------------------------------------------------+
+  Dr. Alexander Golbraikh                   Email: [email-masked]  +
+  Group of Molecular Biophysics                                       +
+  Latvian Institute of Organic Synthesis    FAX:   + (371) 7-821-038  +
+  21 Aizkraukles st.                                                  +
+  Riga LV-1006                              Phone: + (371) 7-552-354  +
+  LATVIA                                                              +
+----------------------------------------------------------------------+

------------------------------------------------------------------------------

  Dear Netters,
  
  Shouldn't be noticed that if you have n points in space, it isn't
 determined what kind of polyhedron you actually have? This is the
 non-convex case. Imagine 4 points in 2D:
 
                           .4
                          3.
                       1.      2.
                       
  You can go 1-2-4-3-1, or 1-3-4-2-1, or 1-3-2-4-1. So the problem is
 in this case quite complicated.
 
  For the "convex" case, staying by 2D (generalization to 3D is 
 straightforward, hopefully :-) I have an idea first to find a gravity
 centre of the polygon. This is the point we are sure it is inside. 
 Then you join g.c. with the point in question- and determine whether
 it crosses the border (point outside) or not (inside).
 
  Andrzej Szymoszek
  [email-masked]  

------------------------------------------------------------------------------

> Just an idea. It's very efficient, but not sure it's right.

Sorry, but it won't work.

> First transform the current coordinate system to one with the given point 
> (x,y,z) as the origin. I.e. subtrate (xi,yi,zi), where i=1,...,N, by 
> (x,y,z). This can be done very fast.
> 
> Second look at the sign of the transformed coordinates of all vertices. 
> If all xi's, where i=1,...,N, have the same sign, the point is outside 
> of the polyhedron. The statement holds for all yi's or all zi's. So if any 
> one of the three conditions holds, the point will be outside. Otherwise if 
> none of the conditions hold, it's inside.

If any one of your test conditions is true (eg, all x positive or all x 
negative) then yes the point is outside the polyhedron, however, if any 
number of tests is false, the point can still be outside the polyhedron.  
The tests only work one way.

Just consider a square on a 2-d surface with all sides at 45 degrees to 
the axes.  Place the origin just outside one side, and notice that your 
have points with both positive and negative values for both x and y axes.

Cheers, Craig

"Imagine if every Thursday your shoes exploded if you tied them the
 usual way.  This happens to us all the time with computers, and nobody
 thinks of complaining."
                -- Jeff Raskin, interviewed in Doctor Dobb's Journal

Craig Taverner
Structural Chemistry, University of the Witwatersrand, South Africa
tel:   +27-11-716-2290                fax: +27-11-716-3826
email: [email-masked]
www:   http://www.gh.wits.ac.za/craig

------------------------------------------------------------------------------

Sorry, but can't a voxel algorhythm be instrumental here?

-- Eugene

------------------------------------------------------------------------------
-----------------------End of summary-----------------------------------------