Fast generation of points inside triangulated objects obtained by cross-sectional contours

Michael Lahanas and Kostas Karouzakis

In medical applications three-dimensional (3D) anatomical models are usually produced from contours from cross-sectional images. Conventional algorithms are described which are used to generate points inside these objects which are required for example for volumetric purposes or for the calculation of the dose distributions in radiotherapy or brachytherapy [LB00], see Fig. 1. These algorithms use triangles obtained from a triangulation of the contours. Using a preprocessing method we generate from the contours so called 2D virtual polygons. This enables us to reduce the 3D point in a polyhedron test into a 2D point in a virtual polygon test. Additional we use optimized bounding boxes [LK00] which further improves the efficiency of this method.

Fig. 1 100000 sampling Points in a part of the body which includes a rib tumor with 50000 sampling points used for volumetric and dosimetric purposes.


Test of if a point is within a 2D contour

We consider a 2D contour as a closed polygon made up of N vertices (xi,yi) i, Î 1,..., N. There are several methods which can prove if a point P=(xp,yp) is within such a polygon but the most efficient method is the topology-based method. In order to determine if a point P is inside the polygon we consider a horizontal ray emanating from P towards the right, see Fig. 1.

The Intersection number test.

In this method we count the number of intersections of a ray with all edges of the polygon. If the number of edges of the polygon intersected by the ray is even then the point P is outside the polygon, and if the number of intersections is odd then the point P lies inside the polygon. This is a consequence of the Jordan curve theorem in topology and a proof is given by Courant et al. [CR41] There are some special cases when a ray passes exactly through a vertex. These situations can be considered using tie breaking rules or by a method such as described by Edelsbrunner et al. [ED90]

Fig 2 Schematic diagram illustrating a generalized point P in relation to a 2D cross-section through a polygon (a) with and (b) without a critical structure inside the cross-section. The number of intersections is shown for each point.

The winding number test in 2D.

An alternative test for whether a point is inside or outside a 2D contour can be achieved using the property that if a point P is inside a polygon such as that described by the contour, then the sum of all angles at P which are made between the adjacent lines from each vertex to the point P equals 2p otherwise it is 0. This method is algorithmic more complex and therefore more time consuming than the intersection number method.

Test of if a point is within a 3D Polyhedron

We also assume that the volume is described as a polyhedron. As for the 2D case there exist a similar intersection and winding number test.

The Intersection number test with the Polyhedron surface.

We have used a variant of the intersection algorithm described by O’Rourke:

Algorithm: Point in Polyhedron

Compute bounding radius R

loop forever

r0 = random ray of length R

r = q + r0.

crossings = 0.

for each face F of polyhedron P do


if degenerate intersection

then Go back to loop

else increment crossings appropriately.

if crossings odd

then q is inside P.

else q is outside P.


O’Rourke Point in Polyhedron Algorithm, see [OR98].

This algorithm avoids degenerate cases such as if the ray crosses an edge or a polyhedron vertex. If such a case occurs a new ray is generated and the test is repeated. Usually it is common to triangulate the facets, so that a line triangle intersection routine is used which is simpler than a general polygon line intersection test.

B. The winding number test.

The winding number test can be extended in three dimensions[CA95]:

Project each face of the polyhedron onto a unit sphere of center P and compute the signed area of the spherical polygon thus determined. The sign is positive if the spherical polygon has counterclockwise orientation and negative otherwise.

The sum S of the signed areas of the projections of all facets of the simple polyhedron P onto the unit sphere p is necessarily 0, 4p or 4p.

Fig. 3 Spherical Triangle from the projection from a point P of a triangular polyhedron facet on a unit sphere.

If the facets of P have clockwise or counterclockwise orientation and the point P is inside the polyhedron then S = -4p or 4p respectively. If the point is outside the polyhedron then regardless of the orientation S = 0.

The area of a spherical triangle can be obtained by Girard s formula:

S = A + B + C = 2p

Where A,B,C are the spherical angles at each vertex. The spherical angle a at a given vertex A is the angle determined by the tangents to the two great circles corresponding to the sides that cross at A. a is also the angle formed by the planes containing those two great circles. Thus a can be determined from the normal vectors to each of these planes, which can be computer without projecting the vertices onto the sphere.

Generation of the intermediate polygons.

For the generation of the polygons we have to consider that the triangles are not necessary stored in the appropriate order. We take a random triangle and calculate the intersection point P1 of this triangle with the plane between the two consecutive slices. Then we take the next side of the triangle and the next intersection point P2. Then we search the triangle which has the same edge A2A3. This is the triangle number 2 with the intersection point P3. We continue this procedure until we arrive again to the point P1. The intersection points P1,P2,P3 define the polygon.

Fig. 4 Generation of a polygon P1P2 from the triangulated surface between two slices (dotted line).

Generation of the virtual polygons.

We have developed a method which avoids the generation of polygons from the intersection of the triangulated surface and the plane through the point to be tested which is computation complex. We use a preprocessing method where so called virtual polygons are generated and stored.

For each consecutive pair of actual contours a virtual polygon is generated. The vertices of this polygon are functions of z, i.e. the ith vertex is given as Pi(x,y) = (fi(z), gi(z)). A polygon is therefore a function of z and we call therefore these polygons virtual polygons. We use the same method as shown in Fig. 4. but this time we do not store the actual x,y values for each polygon vertex but the function (fi(z), gi(z)).

Fig. 5 Example of an object defined by three contours and the two corresponding virtual polygons shown in red color generated for some z coordinates.

Fig. 6 Extension of point in 3D Polyhedron considering objects inside the polyhedron.


[LB00] M. Lahanas, D. Baltas, N. Milickovic, S. Giannouli, and N. Zamboglou,Generation of uniformly distributed dose points for anatomy-based-three-dimensional dose optimization in brachytherapy, Med. Phys. 27, 1034-1046, 2000. Abstract

[LK00] M. Lahanas, T. Kemmerer, N. Milickovic, D. Baltas, N. Zamboglou," Optimized bounding boxes for three-dimensional treatment planning in brachytherapy," accepted for publication in Med. Phys. 2000.

[CR41] R. Courant and H. Robbins, "What is Mathematics?," Oxford University Press, (1941).

[ED90] H. Edelsbrunner and E. P. Mücke, "Simulation of Simplicity: A technique to cope with degenerate cases in geometric algorithms," ACM Transactions on Graphics 9, 66-104 (1990).

[CA95] P. C. P. Carvaiho, P. R. Cavalcanti, Point in Polyhedron Testing Using Spherical Polygons, Graphics Gems V edited by A. W. Paeth, Academic Press 1995

[OR98]J. O’Rourke, Computational Geometry in C , Cambridge, Second Edition, Academic Press 1998