Algorithms







What follows is a chaotic collection of some maths and explanations of
mainly (say 98 %) computer graphics algorithms....

Most stuff are my personal
saved items from the comp.graphics.algorithms newsgroup.

See what you can do with it.













.... 3d rotations ...

!!Gimbal lock?!!

Don't let it get you down!

In the following text, we will take the example
of transforming the points data of a spaceship.


We've got our spaceship aligned down our z.
First of all, we want to roll the spaceship...
If we do it now, it saves about 5000% calculation.
Thats serious cpu.
So... we do our standard rotate about the Z...
I'm SURE you've got the calculation for that...
But hell, we're all in the same boat, here it is:

This is a matrix. The lines are to help me.
+-------+-------+-------+
| cos t |-sin t |   0   |
+-------+-------+-------+
| sin t | cos t |   0   |
+-------+-------+-------+
|   0   |   0   |   1   |
+-------+-------+-------+
Now, being a sane rational person, you'll look at that,
see all the 0 and 1, and notice that you have four multiplies
to do.
Remember that: '4' four (4)

This comes in later

This method implies that our roll is an angle...
not a bad assumption at all.

Now, we have some vector somewhere which is the orientation of the plane..
the trickiest bit we'll encounter here, is actually calculating the axis...
We're REALLY gonna need a LOT of bits in fix to store this,
'cos this is the most (in fact the only thing) prone to deterioration...

We take this to be some axis that the user can control...
The real hard work starts because we have to use rotate this
by whatever, and it'll be quite complex to do interactively...

However, suppose we do have this vector, as we almost certainly will
have...
We can take this vector, and call it 'd', because its the [d]irection the
plane
is pointing...
Now, we want to map the z onto d...
So, we rotate about our arbitrary axis...
Uh oh!

First of all, we must calculate the cross product of the two vectors...
Uh oh! six multiplies! But... its only a preprocess step..

HOLD it right there! ALL of those six multiplies will be by either 0 or 1
Remember, Z is always {0,0,1}
So...
Lets quote the cross product:

aXb= ay*bz - by*az
     bx*az - ax*bz
     ax*by - bx*ay

Now, lets say that b is the z...

So, we can insert the values to yield:

aXz = ay*1 - 0*az
      0*az - ax*1
      ax*0 - 0*ay

Well, this quickly simplifies to:

aXz = ay
      -ax
      0

This makes things a LOT easier..
Two variable assignments, and the whole cross product is done...
Now, we've preprocessed this, we've still got that arbitrary axis
rotation to do...
But how many degrees to rotate it around?!??!?!?

Well, we know that:    A.B = |A||B|cos angle_between
A and B are both unit vectors, so:

A.B=cos angle between...

But B=z
So... az=cos angle between

You can quickly verify this for yourself...

So, to recap:
We want to rotate about {ay,-ax,0} by angle arc-cosine(az)
Now... Here's the matrix for a rotation about an arbitrary axis...
I calculated this with a little help from a maths tutor, who found
some info on quarternions in "the mathematical gazette" (don't ask ;-)
and then simplified it into english....
I found out how to express an axis in quarternions, and then how
to rotate a quarternion...
The inns and outs were rather complex, and involved seven A4 sides of
working...  The things I will do for fun.

Anyway, it ended up as a huge expansion, which concurs with Foley
van Dam, Feiner, Hughes, et al... I shall quote theirs, as it uses
some trigonometric identities, which make it more eloquent:

We'll say our axis={x,y,z} for typing's sake
I just hope my squared symbol shows up...
Just a check:

x² = x*x
that work?

Also, t=theta

Matrix:
x²+cos t(1-x²)       xy(1-cos t)-zsin t  xz(1-cos t)+ysin t
xy(1-cos t)+zsin t   y²+cos t(1-y²)      yz(1-cos t)-xsin t
xz(1-cos t)-ysin t   yz(1-cos t)+xsin t  z²+cos t(1-z²)

Hopefully, you can make sense of that...
However, we can simplify it.
Remember that the z component of our arbitrary axis was 0...
That means that z in the above = 0...
Well then, this yields:

x²+cos t(1-x²)       xy(1-cos t)                     ysin t
xy(1-cos t)          y²+cos t(1-y²)                  xsin t
           -ysin t               xsin t     cos t


(I have merely overwritten the *z parts with ' ')
Straightening this up yields:


x²+cos t(1-x²)  xy(1-cos t)     ysin t
xy(1-cos t)     y²+cos t(1-y²)  xsin t
-ysin t         xsin t          cos t


Now, this is very neat!
This gives us nine multiplies per point...
add that nine to the four from above...
eleven multiplies to produce an aircraft at any orientation!

Some people (me especially) hate stuff presented in matrix form,
so i'll multiply out:

for point {a,b,c}
and axis {x,y,z}
and angle t

a'=a(x²+cost(1-x²))  +  b(xy(1-cos t))   - c(ysin t)
b'=a(xy(1-cost))     +  b(y²+cost(1-y²)) + c(xsin t)
c'=a(ysin t)         -  b(xsin t)        + c(cos t)


Et voila!


However, you will note that we had two matrix form transformations,
one for the roll, and another for the arbitrary.
It +is+ possible to multiply these, but I haven't, just so that
you can derive some mathematical enjoyment from this "experience"
(also, i can't be bothered to multiply them :-)

Source:
I've provided some sample source.
This is ONLY so that you can get a further idea of how to implement,
it's relatively simple, and optimised towards simplifying the
instructions used...
It just offers a couple of functions... Don't expect this to work...
It might, it might not... (it's not actually part of my implementation,
as that was in ARM assembler - little use to most ppl here!)

So, here it is:
[
// +-------------------------------+
// | SpaceShip transformation code |
// |  Code to rotate a 3d object   |
// |  with 'behaviour' similar to  |
// |  a spaceship or aircraft      |
// +-------------------------------+

//  This code specifies that the spaceship has been designed
// as relative to the z axis...
//  i.e. that the point it rotates around is at the origin,
// and its nose lies on the z axis.
/*    Y   -Z
       | /
     /\|/____X
     \  \
    / /\/
   |_/
  Z
Sorry, but its hard (for me) to do a spaceship in ascii!
*/


// Usage of functions:
//	void	Preprocess_CalculateMatrix(float x,float y,float z)
// Use this to calculate the matrix needed to transform all the points
// pass it the NORMALISED components of the vector describing
// the direction of the nose from the origin.
// (i.e. to point it straight towards you, send {0,0,1})

//	void	SetRoll(float angle)
// Use this to set up the angle of roll (in radians, because thats what
math.h wants)

//	Point	Rotate_This_Point(Point)
// use this to rotate a point.

// I've just 'invented' a type here
typedef float Point[3];
// this is just to help with parameter passing...

// you don't need to use ANY of these variables, they're just for internal
working

float Rot_Matrix[3][3];
float ax,ay;
float st,ct;
float zs,zc;
float x,y,z;
// sorry these aren't very descriptive, but they get used a lot,
// and it looks horrible otherwise

void
Preprocess_CalculateMatrix(float x, float y, float z)
{
	// right then...
	// first off, sorry for using c++ comments
	// secondarily, the input values:

		// x,y,z are the *****NORMALIZED****** components of the vector
		// pointing in the direction which you want your object to now be aligned
		// IMPORTANT: Initially, the object is aligned down the Z axis

	ax=y;
	ay=-x;

	ct=z;
	st=sin(acs(ct));
	// WARNING this line seems to imply a limitation that the this cannot face
	// backwards... there ought to be some form of test which would fix this,
	// but my trig seems to have gone to sleep.
	// anyway... not too hard to invert sign of all z components ;-)


	Rot_Matrix[0][0]=((ax*ax)+(ct*(1-(ax*ax))));
	Rot_Matrix[1][0]=((ax*ay*(1-ct)));
	Rot_Matrix[2][0]=(ay*st);

	Rot_Matrix[0][1]=((ax*ay*(1-ct)));
	Rot_Matrix[1][1]=((ay*ay)+(ct*(1-(ay*ay))));
	Rot_Matrix[2][1]=-(ax*st);

	Rot_Matrix[0][2]=-(ay*st);
	Rot_Matrix[1][2]=(ax*st);
	Rot_Matrix[2][2]=(ct);
}

void
SetRoll(float angle)
{
	zs=sin(angle);
	zc=cos(angle);
}

Point
RotatePoint(Point p)
{

	x=p[0];
	y=p[1];
	z=p[2];

	// roll
	p[0]=(zc*x)+(zs*y);
	p[1]=(zc*y)-(zs*x);
	x=p[0];
	y=p[1];

	p[0]=(x*Rot_Matrix[0][0])+(y*Rot_Matrix[0][1])+(z*Rot_Matrix[0][2]);
	p[1]=(x*Rot_Matrix[1][0])+(y*Rot_Matrix[1][1])+(z*Rot_Matrix[1][2]);
	p[2]=(x*Rot_Matrix[2][0])+(y*Rot_Matrix[2][1])+(z*Rot_Matrix[2][2]);
	return p;
}

// and thats it...
]

Again, this is just a (small) extension to my code, to give you a little
idea of how this
could be implemented... I offer no guarantee that this will work!












I have made a small proggy just to check out how I could do a blur.
I used xor to get a nice looking image, and then calculated the
average-color of the 4 surrounding colors.
Is this the right way to do it ?
It seems to go slow as hell... Can anyone give me any tip on how to
increase the speed of the blur ?
(in pascal or somewhat readable asm...)

Anders Pedersen

Code posted under:

program _xorblur;
uses crt,dos;
TYPE Virtual = Array [1..64000] of byte;  { The size of our Virtual
Screen }
     VirtPtr = ^Virtual;                  { Pointer to the virtual
screen }

var x,y,i:integer;
col:byte;
Virscr,virscr2 : VirtPtr;                     { Our first Virtual screen

}
Vaddr,vaddr2  : word;                        { The segment of our
virtual screen}

procedure flip32(source,dest:Word); assembler;
  { This copies the entire screen at "source" to destination }
asm
  push    ds
  mov     ax, [Dest]
  mov     es, ax
  mov     ax, [Source]
  mov     ds, ax
  xor     si, si
  xor     di, di
  mov     cx, 16000
  db      $F3, $66, $A5
  pop     ds
end;

Procedure Cls32 (Where:word;Col : Byte); assembler;

   { This clears the screen to the specified color }
asm
   push    es
   mov     cx, 16000;
   mov     es,[where]
   xor     di,di
   mov     al,[col]
   mov     ah,al
   mov     dx, ax
   db      $66, $C1, $E0, $10         {shl eax, 16}
   mov     ax, dx
   db      $F3, $66, $AB              {rep stosd}
   pop     es
End;
Procedure SetUpVirtual;
   { This sets up the memory needed for the virtual screen }
BEGIN
  GetMem (VirScr,64000);
  vaddr := seg (virscr^);   {The segment for page 1}
  GetMem (VirScr2,64000);
  vaddr2 := seg (virscr2^); {The segment for page 2}
END; {Setupvirtual}

Procedure ShutDown;
   { This frees the memory used by the virtual screen }
BEGIN
  FreeMem (VirScr,64000);
  FreeMem (VirScr2,64000);
END;  {Shutdown}

procedure init;assembler;
asm
 mov ax,13h
 int 10h
end;

procedure setpal(nr,r,g,b:byte);
begin;
 port[$3c8]:=nr;
 port[$3c9]:=r;
 port[$3c9]:=g;
 port[$3c9]:=b;
end;

procedure putpixel(x,y:integer;col:byte);
begin;
 mem[vaddr:x+(y*320)]:=col;
end;


procedure blur;
var c,c1,c2,c3,c4:integer;
color:single;
begin;
for y:=1 to 191 do
    for x:=1 to 320 do
         begin;
        c:=mem[$A000:x+(320*y)];
        c1:=mem[$A000:(x-1)+(320*y)];
        c2:=mem[$A000:(x+1)+(320*y)];
        c3:=mem[$A000:x+(320*(y-1))];
        c4:=mem[$A000:x+(320*(y+1))];
        color:=round((c1+c2+c3+c4)/4);
        mem[vaddr:x+(320*y)]:=round(color);
end;
end;


begin;
init;
setupvirtual;
for i:=1 to 256 do
setpal(i,0,i,i div 2);
cls32(vaddr,0);
for y:=1 to 191 do
 for x:=1 to 319 do
  putpixel(x,y,(x xor y));
flip32(vaddr,sega000);
repeat;
cls32(vaddr,0);
blur;
flip32(vaddr,sega000);
until keypressed;
shutdown;
asm mov ax,3h;int 10h;end;
end.
















I my ray tracer I transform the direction of the starting ray like
this:

Vector vpn( Normalize( world->camera->target - world->camera->position
) );
Vector vup( 0,-1,0 );
vup = Normalize( vup- vpn*DotProd( vup, vpn ) );
Vector vrg( CrossProd( vup, vpn ) );

Ray ray( world->camera->position, vpn + vrg*(x-160)/160 +
vup*(y-120)/120 );

But with this method, I get a few problems.
1. I can't conrtol the distance to the viewingplane.
2. How do I find the screen x- and y- coordinates for an arbitary
3d-coordinate. I need this for post-rendering of lens flares.

If I always had my camera placed in (0,0,0), looking at (x-160, y-120,
250), there would be no problem. One solution is to, instead of
transform the ray with the above algorithm, rotate the world
coordinates with vrg,vup,vpn, but this doesn't work when an object is
defined by a center and a rotation angle.















--------------------------------------------------------------------------------------------------------
COMPUTER GRAPHICS FAQ                                Table of Contents
---------------------------------------------------------------------------------------------------------

0. General Information
   0.01: Charter of comp.graphics.algorithms
   0.02: Are the postings to comp.graphics.algorithms archived?
   0.03: How can I get this FAQ?
   0.04: What are some must-have books on graphics algorithms?
   0.05: Are there any online references?
   0.06: Are there other graphics related FAQs?
   0.07: Where is all the source?


 
   I have left out Section 0 since it contained no useful algorithms



1. 2D Computations: Points, Segments, Circles, Etc.
   1.01: How do I rotate a 2D point?
   1.02: How do I find the distance from a point to a line?
   1.03: How do I find intersections of 2 2D line segments?
   1.04: How do I generate a circle through three points?
   1.05: How can the smallest circle enclosing a set of points be found?
   1.06: Where can I find graph layout algorithms?

2. 2D Polygon Computations
   2.01: How do I find the area of a polygon?
   2.02: How can the centroid of a polygon be computed?
   2.03: How do I find if a point lies within a polygon?
   2.04: How do I find the intersection of two convex polygons?
   2.05: How do I do a hidden surface test (backface culling) with 2d
points?
   2.06: How do I find a single point inside a simple polygon?
   2.07: How do I find the orientation of a simple polygon?

3. 2D Image/Pixel Computations
   3.01: How do I rotate a bitmap?
   3.02: How do I display a 24 bit image in 8 bits?
   3.03: How do I fill the area of an arbitrary shape?
   3.04: How do I find the 'edges' in a bitmap?
   3.05: How do I enlarge/sharpen/fuzz a bitmap?
   3.06: How do I map a texture on to a shape?
   3.07: How do I detect a 'corner' in a collection of points?
   3.08: Where do I get source to display (raster font format)?
   3.09: What is morphing/how is it done?
   3.10: How do I quickly draw a filled triangle?
   3.11: D Noise functions and turbulence in Solid texturing.

   3.12: How do I generate realistic sythetic textures?
   3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)?

4. Curve Computations
   4.01: How do I generate a bezier curve that is parallel to another
bezier?
   4.02: How do I split a bezier at a specific value for t?
   4.03: How do I find a t value at a specific point on a bezier?
   4.04: How do I fit a bezier curve to a circle?

5. 3D computations
   5.01: How do I rotate a 3D point?
   5.02: What is ARCBALL and where is the source?
   5.03: How do I clip a polygon against a rectangle?
   5.04: How do I clip a polygon against another polygon?
   5.05: How do I find the intersection of a line and a plane?
   5.06: How do I determine the intersection between a ray and a polygon?
   5.07: How do I determine the intersection between a ray and a sphere?
   5.08: How do I find the intersection of a ray and a bezier surface?
   5.09: How do I ray trace caustics?
   5.10: What is the marching cubes algorithm?
   5.11: What is the status of the patent on the "marching cubes"
algorithm?
   5.12: How do I do a hidden surface test (backface culling) with 3d
points?
   5.13: Where can I find algorithms for 3D collision detection?
   5.14: How do I perform basic viewing in 3d?
   5.15: How do I optimize a 3D polygon mesh?
   5.16: How can I perform volume rendering?
   5.17: Where can I get the spline description of the famous teapot etc.?
   5.18: How can the distance between two lines in space be computed?
   5.19: How can I compute the volume of a polyhedron?

6. Geometric Structures and Mathematics
   6.01: Where can I get source for Voronoi/Delaunay triangulation?
   6.02: Where do I get source for convex hull?
   6.03: Where do I get source for halfspace intersection?
   6.04: What are barycentric coordinates?
   6.05: How do I generate a random point inside a triangle?
   6.06: How do I evenly distribute N points on (tesselate) a sphere?
   6.07: What are coordinates for the vertices of an icosohedron?
   6.08: How do I generate random points on the surface of a sphere?

7. Contributors
   7.01: How can you contribute to this FAQ?
   7.02: Contributors.  Who made this all possible.



   I have left out Section 7 also since it contained no useful algorithms




Search e.g. for "Section 6" to find that section.
Search e.g. for "Subject 6.04" to find that item.

----------------------------------------------------------------------
Section 1. 2D Computations: Points, Segments, Circles, Etc.
----------------------------------------------------------------------
Subject 1.01: How do I rotate a 2D point?

    In 2-D, the 2x2 matrix is very simple.  If you want to rotate a
    column vector v by t degrees using matrix M, use

        M = {{cos t, -sin t}, {sin t, cos t}} in M*v.

    If you have a row vector, use the transpose of M (turn rows into
    columns and vice versa).  If you want to combine rotations, in 2-D
    you can just add their angles, but in higher dimensions you must
    multiply their matrices.


----------------------------------------------------------------------
Subject 1.02: How do I find the distance from a point to a line?


    Let the point be C (XC,YC) and the line be AB (XA,YA) to (XB,YB).
    The length of the line segment AB is L:

        L=((XB-XA)**2+(YB-YA)**2)**0.5

    and
            (YA-YC)(YA-YB)-(XA-XC)(XB-XA)
        r = -----------------------------
                        L**2

            (YA-YC)(XB-XA)-(XA-XC)(YB-YA)
        s = -----------------------------
                        L**2

    Let I be the point of perpendicular projection of C onto AB, the

        XI=XA+r(XB-XA)
        YI=YA+r(YB-YA)

    Distance from A to I = r*L
    Distance from C to I = s*L

    If r<0      I is on backward extension of AB
    If r>1      I is on ahead extension of AB
    If 0<=r<=1  I is on AB

    If s<0      C is left of AB (you can just check the numerator)
    If s>0      C is right of AB
    If s=0      C is on AB


----------------------------------------------------------------------
Subject 1.03: How do I find intersections of 2 2D line segments?

    This problem can be extremely easy or extremely difficult depends
    on your applications.  If all you want is the intersection point,
    the following should work:

    Let A,B,C,D be 2-space position vectors.  Then the directed line
    segments AB & CD are given by:

        AB=A+r(B-A), r in [0,1]
        CD=C+s(D-C), s in [0,1]

    If AB & CD intersect, then

        A+r(B-A)=C+s(D-C), or

        XA+r(XB-XA)=XC+s(XD-XC)
        YA+r(YB-YA)=YC+s(YD-YC)  for some r,s in [0,1]

    Solving the above for r and s yields

            (YA-YC)(XD-XC)-(XA-XC)(YD-YC)
        r = -----------------------------  (eqn 1)
            (XB-XA)(YD-YC)-(YB-YA)(XD-XC)

            (YA-YC)(XB-XA)-(XA-XC)(YB-YA)

        s = -----------------------------  (eqn 2)
            (XB-XA)(YD-YC)-(YB-YA)(XD-XC)

    Let I be the position vector of the intersection point, then

        I=A+r(B-A) or

        XI=XA+r(XB-XA)
        YI=YA+r(YB-YA)

    By examining the values of r & s, you can also determine some
    other limiting conditions:

        If 0<=r<=1 & 0<=s<=1, intersection exists
            r<0 or r>1 or s<0 or s>1 line segments do not intersect

        If the denominator in eqn 1 is zero, AB & CD are parallel
        If the numerator in eqn 1 is also zero, AB & CD are coincident

    If the intersection point of the 2 lines are needed (lines in this
    context mean infinite lines) regardless whether the two line
    segments intersect, then

        If r>1, I is located on extension of AB
        If r<0, I is located on extension of BA
        If s>1, I is located on extension of CD
        If s<0, I is located on extension of DC

    Also note that the denominators of eqn 1 & 2 are identical.

    References:

    [O'Rourke] pp. 249-51
    [Gems III] pp. 199-202 "Faster Line Segment Intersection,"


----------------------------------------------------------------------
Subject 1.04: How do I generate a circle through three points?

    Let the three given points be a, b, c.  Use _0 and _1 to represent
    x and y coordinates. The coordinates of the center p=(p_0,p_1) of
    the circle determined by a, b, and c are:

        A = b_0 - a_0;
        B = b_1 - a_1;
        C = c_0 - a_0;
        D = c_1 - a_1;
    
        E = A*(a_0 + b_0) + B*(a_1 + b_1);
        F = C*(a_0 + c_0) + D*(a_1 + c_1);
    
        G = 2.0*(A*(c_1 - b_1)-B*(c_0 - b_0));
    
        p_0 = (D*E - B*F) / G;
        p_1 = (A*F - C*E) / G;
 
    If G is zero then the three points are collinear and no finite-radius
    circle through them exists.  Otherwise, the radius of the circle is:

            r^2 = (a_0 - p_0)^2 + (a_1 - p_1)^2

    Reference:

    [O' Rourke] p. 201. Simplified by Jim Ward.

----------------------------------------------------------------------
Subject 1.05: How can the smallest circle enclosing a set of points be
found?

    This circle is often called the minimum spanning circle.  It can be
    computed in O(n log n) time for n points.  The center lies on
    the furthest-point Voronoi diagram.  Computing the diagram constrains
    the search for the center.  Constructing the diagram can be
accomplished
    by a 3D convex hull algorithm; for that connection, see, e.g.,
    [O'Rourke, p.195ff].  For a direct algorithm, see:
      S. Skyum, "A simple algorithm for computing the smallest enclosing
circle"
      Inform. Process. Lett. 37 (1991) 121--125.

----------------------------------------------------------------------
Subject 1.06: Where can I find graph layout algorithms?

    See the paper by Eades and Tamassia, "Algorithms for Drawing
    Graphs: An Annotated Bibliography," Tech Rep CS-89-09, Dept.  of
    CS, Brown Univ, Feb. 1989.

    A revised version of the annotated bibliography on graph drawing
    algorithms by Giuseppe Di Battista, Peter Eades, and Roberto
    Tamassia is now available from

    ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.tex.gz and
    ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.ps.gz


----------------------------------------------------------------------
Section 2. 2D Polygon Computations
----------------------------------------------------------------------
Subject 2.01: How do I find the area of a polygon?

    The signed area can be computed in linear time by a simple sum.
    The key formula is this:

        If the coordinates of vertex v_i are x_i and y_i,
        twice the signed area of a polygon is given by

        2 A( P ) = sum_{i=0}^{n-1} (x_i y_{i+1} - y_i x_{i+1}).

    Here n is the number of vertices of the polygon.
    References: [O' Rourke] pp. 18-27; [Gems II] pp. 5-6:
    "The Area of a Simple Polygon", Jon Rokne.

    To find the area of a planar polygon not in the x-y plane, use:
    
       2 A(P) = abs(N . (sum_{i=0}^{n-1} (v_i x v_{i+1})))
    
       where N is a unit vector normal to the plane. The `.' represents the
       dot product operator, the `x' represents the cross product operator,
       and abs() is the absolute value function.
    
    [Gems II] pp. 170-171:
    "Area of Planar Polygons and Volume of Polyhedra", Ronald N. Goldman.

----------------------------------------------------------------------
Subject 2.02: How can the centroid of a polygon be computed?

    The centroid (a.k.a. the center of mass, or center of gravity)
    of a polygon can be computed as the weighted sum of the centroids
    of a partition of the polygon into triangles.  The centroid of a
    triangle is simply the average of its three vertices, i.e., it
    has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3.  This 

    suggests first triangulating the polygon, then forming a sum
    of the centroids of each triangle, weighted by the area of
    each triangle, the whole sum normalized by the total polygon area.
    This indeed works, but there is a simpler method:  the triangulation
    need not be a partition, but rather can use positively and
    negatively oriented triangles (with positive and negative areas),
    as is used when computing the area of a polygon.  This leads to
    a very simple algorithm for computing the centroid, based on a
    sum of triangle centroids weighted with their signed area.
    The triangles can be taken to be those formed by one fixed vertex
    v0 of the polygon, and the two endpoints of consecutive edges of
    the polygon: (v1,v2), (v2,v3), etc.  The area of a triangle with
    vertices a, b, c is half of this expression:
                (b[X] - a[X]) * (c[Y] - a[Y]) -
                (c[X] - a[X]) * (b[Y] - a[Y]);
 
    Code available at ftp://grendel.csc.smith.edu/pub/code/centroid.c (3K).
    Reference: [Gems IV] pp.3-6; also includes code.

----------------------------------------------------------------------
Subject 2.03: How do I find if a point lies within a polygon?

    The definitive reference is "Point in Polyon Strategies" by
    Eric Haines [Gems IV]  pp. 24-46.
    The code in the Sedgewick book Algorithms (2nd Edition, p.354) is 
    incorrect.

    The essence of the ray-crossing method is as follows.
    Think of standing inside a field with a fence representing the polygon.

    Then walk north. If you have to jump the fence you know you are now 
    outside the poly. If you have to cross again you know you are now 
    inside again; i.e., if you were inside the field to start with, the
total 
    number of fence jumps you would make will be odd, whereas if you were 
    ouside the jumps will be even.

    The code below is from Wm. Randolph Franklin 
    with some minor modifications for speed.  It returns 1 for strictly
    interior points, 0 for strictly exterior, and 0 or 1 for points on
    the boundary.  The boundary behavior is complex but determined;
    in particular, for a partition of a region into polygons, each point
    is "in" exactly one polygon.  See the references below for more detail.

    int pnpoly(int npol, float *xp, float *yp, float x, float y)
    {
      int i, j, c = 0;
      for (i = 0, j = npol-1; i < npol; j = i++) {
        if ((((yp[i]<=y) && (y