Method IO Graphics PART05.SH

# ! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 5 (of 5)."
# Contents:  AALines/FastMatMul.c FitCurves.c GGVecLib.c NearestPoint.c
# Wrapped by craig@weedeater on Wed Dec 12 20:45:16 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'AALines/FastMatMul.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'AALines/FastMatMul.c'\"
else
echo shar: Extracting \"'AALines/FastMatMul.c'\" \(12335 characters\)
sed "s/^X//" >'AALines/FastMatMul.c' <<'END_OF_FILE'
X/*  FILENAME: FastMatMul.c  [revised 18 AUG 90]
X
X    AUTHOR:  Kelvin Thompson
X
X    DESCRIPTION:  Routines to multiply different kinds of 4x4
X      matrices as fast as possible.  Based on ideas on pages 454,
X      460-461, and 646 of _Graphics_Gems_.
X
X	These routines offer two advantages over the standard
X      V3MatMul in the _Graphics_Gems_ vector library GGVecLib.c.
X      First, the routines are faster.  Second, they can handle
X      input matrices that are the same as the output matrix.
X      The routines have the disadvantage of taking more code
X      space (from unwound loops).
X
X        The routines are about as fast as you can get for
X      general-purpose multiplication.  If you have special
X      knowledge about your system, you may be able to improve
X      them a little more:
X
X        [1]  If you know that your input and output matrices will
X      never overlap: remove the tests near the beginning and end of
X      each routine, and just #define 'mptr' to 'result'.  (The
X      standard library's V3MatMul makes this assumption.)
X
X	[2]  If you know that your compiler supports more than
X      three pointer-to-double registers in a subroutine: make
X      'result' in each routine a register variable.  You might
X      also make the 'usetemp' boolean a register.
X
X	[3]  If you have limited stack space, or your system can
X      access global memory faster than stack:  make each routine's
X      'tempx' a static, or let all routines share a global 'tempx'.
X      (This is useless if assumption [1] holds.)
X*/
X
X/* definitions from "GraphicsGems.h" */
Xtypedef struct Matrix3Struct {	/* 3-by-3 matrix */
X	double element[3][3];
X	} Matrix3;
Xtypedef struct Matrix4Struct {	/* 4-by-4 matrix */
X	double element[4][4];
X	} Matrix4;
X
X/* routines in this file */
XMatrix3 *V2MatMul();     /* general 3x3 matrix multiplier */
XMatrix4 *V3MatMul();     /* general 4x4 matrix multiplier */
XMatrix4 *V3AffMatMul();  /* affine 4x4 matrix multiplier */
XMatrix4 *V3LinMatMul();  /* linear 4x4 matrix multiplier */
X
X/* macro to access matrix element */
X#define MVAL(mptr,row,col)  ((mptr)->element[row][col])
X
X
X
X
X
X/*  ************************************
X    *                                  *
X    *           V2MatMul               *
X    *                                  *
X    ************************************
X
X    DESCRIPTION:  Multiply two general 3x3 matricies.  If one of
X      the input matrices is the same as the output, write the
X      result to a temporary matrix during multiplication, then
X      copy to the output matrix.
X
X    ENTRY:
X      a -- pointer to left matrix
X      b -- pointer to right matrix
X      result -- result matrix
X
X    EXIT:  returns 'result'
X*/
X
XMatrix3 *V2MatMul ( a, b, result )
Xregister Matrix3 *a,*b;
XMatrix3 *result;
X{
Xregister Matrix3 *mptr;
Xint usetemp;  /* boolean */
XMatrix3 tempx;
X
X/* decide where intermediate result goes */
Xusetemp = ( a == result  ||  b == result );
Xif ( usetemp )
X  mptr = & tempx;
Xelse
X  mptr = result;
X
XMVAL(mptr,0,0) =
X     MVAL(a,0,0) * MVAL(b,0,0)
X  +  MVAL(a,0,1) * MVAL(b,1,0)
X  +  MVAL(a,0,2) * MVAL(b,2,0);
X
XMVAL(mptr,0,1) =
X     MVAL(a,0,0) * MVAL(b,0,1)
X  +  MVAL(a,0,1) * MVAL(b,1,1)
X  +  MVAL(a,0,2) * MVAL(b,2,1);
X
XMVAL(mptr,0,2) =
X     MVAL(a,0,0) * MVAL(b,0,2)
X  +  MVAL(a,0,1) * MVAL(b,1,2)
X  +  MVAL(a,0,2) * MVAL(b,2,2);
X
XMVAL(mptr,1,0) =
X     MVAL(a,1,0) * MVAL(b,0,0)
X  +  MVAL(a,1,1) * MVAL(b,1,0)
X  +  MVAL(a,1,2) * MVAL(b,2,0);
X
XMVAL(mptr,1,1) =
X     MVAL(a,1,0) * MVAL(b,0,1)
X  +  MVAL(a,1,1) * MVAL(b,1,1)
X  +  MVAL(a,1,2) * MVAL(b,2,1);
X
XMVAL(mptr,1,2) =
X     MVAL(a,1,0) * MVAL(b,0,2)
X  +  MVAL(a,1,1) * MVAL(b,1,2)
X  +  MVAL(a,1,2) * MVAL(b,2,2);
X
XMVAL(mptr,2,0) =
X     MVAL(a,2,0) * MVAL(b,0,0)
X  +  MVAL(a,2,1) * MVAL(b,1,0)
X  +  MVAL(a,2,2) * MVAL(b,2,0);
X
XMVAL(mptr,2,1) =
X     MVAL(a,2,0) * MVAL(b,0,1)
X  +  MVAL(a,2,1) * MVAL(b,1,1)
X  +  MVAL(a,2,2) * MVAL(b,2,1);
X
XMVAL(mptr,2,2) =
X     MVAL(a,2,0) * MVAL(b,0,2)
X  +  MVAL(a,2,1) * MVAL(b,1,2)
X  +  MVAL(a,2,2) * MVAL(b,2,2);
X
X/* copy temp matrix to result if needed */
Xif ( usetemp )
X  *result = *mptr;
X
Xreturn result;
X}
X
X
X
X
X/*  ************************************
X    *                                  *
X    *           V3MatMul               *
X    *                                  *
X    ************************************
X
X    DESCRIPTION:  Multiply two general 4x4 matricies.  If one of
X      the input matrices is the same as the output, write the
X      result to a temporary matrix during multiplication, then
X      copy to the output matrix.
X
X    ENTRY:
X      a -- pointer to left matrix
X      b -- pointer to right matrix
X      result -- result matrix
X
X    EXIT:  returns 'result'
X*/
X
XMatrix4 *V3MatMul ( a, b, result )
Xregister Matrix4 *a,*b;
XMatrix4 *result;
X{
Xregister Matrix4 *mptr;
Xint usetemp;  /* boolean */
XMatrix4 tempx;
X
X/* decide where intermediate result goes */
Xusetemp = ( a == result  ||  b == result );
Xif ( usetemp )
X  mptr = & tempx;
Xelse
X  mptr = result;
X
XMVAL(mptr,0,0) =
X     MVAL(a,0,0) * MVAL(b,0,0)
X  +  MVAL(a,0,1) * MVAL(b,1,0)
X  +  MVAL(a,0,2) * MVAL(b,2,0)
X  +  MVAL(a,0,3) * MVAL(b,3,0);
X
XMVAL(mptr,0,1) =
X     MVAL(a,0,0) * MVAL(b,0,1)
X  +  MVAL(a,0,1) * MVAL(b,1,1)
X  +  MVAL(a,0,2) * MVAL(b,2,1)
X  +  MVAL(a,0,3) * MVAL(b,3,1);
X
XMVAL(mptr,0,2) =
X     MVAL(a,0,0) * MVAL(b,0,2)
X  +  MVAL(a,0,1) * MVAL(b,1,2)
X  +  MVAL(a,0,2) * MVAL(b,2,2)
X  +  MVAL(a,0,3) * MVAL(b,3,2);
X
XMVAL(mptr,0,3) =
X     MVAL(a,0,0) * MVAL(b,0,3)
X  +  MVAL(a,0,1) * MVAL(b,1,3)
X  +  MVAL(a,0,2) * MVAL(b,2,3)
X  +  MVAL(a,0,3) * MVAL(b,3,3);
X
XMVAL(mptr,1,0) =
X     MVAL(a,1,0) * MVAL(b,0,0)
X  +  MVAL(a,1,1) * MVAL(b,1,0)
X  +  MVAL(a,1,2) * MVAL(b,2,0)
X  +  MVAL(a,1,3) * MVAL(b,3,0);
X
XMVAL(mptr,1,1) =
X     MVAL(a,1,0) * MVAL(b,0,1)
X  +  MVAL(a,1,1) * MVAL(b,1,1)
X  +  MVAL(a,1,2) * MVAL(b,2,1)
X  +  MVAL(a,1,3) * MVAL(b,3,1);
X
XMVAL(mptr,1,2) =
X     MVAL(a,1,0) * MVAL(b,0,2)
X  +  MVAL(a,1,1) * MVAL(b,1,2)
X  +  MVAL(a,1,2) * MVAL(b,2,2)
X  +  MVAL(a,1,3) * MVAL(b,3,2);
X
XMVAL(mptr,1,3) =
X     MVAL(a,1,0) * MVAL(b,0,3)
X  +  MVAL(a,1,1) * MVAL(b,1,3)
X  +  MVAL(a,1,2) * MVAL(b,2,3)
X  +  MVAL(a,1,3) * MVAL(b,3,3);
X
XMVAL(mptr,2,0) =
X     MVAL(a,2,0) * MVAL(b,0,0)
X  +  MVAL(a,2,1) * MVAL(b,1,0)
X  +  MVAL(a,2,2) * MVAL(b,2,0)
X  +  MVAL(a,2,3) * MVAL(b,3,0);
X
XMVAL(mptr,2,1) =
X     MVAL(a,2,0) * MVAL(b,0,1)
X  +  MVAL(a,2,1) * MVAL(b,1,1)
X  +  MVAL(a,2,2) * MVAL(b,2,1)
X  +  MVAL(a,2,3) * MVAL(b,3,1);
X
XMVAL(mptr,2,2) =
X     MVAL(a,2,0) * MVAL(b,0,2)
X  +  MVAL(a,2,1) * MVAL(b,1,2)
X  +  MVAL(a,2,2) * MVAL(b,2,2)
X  +  MVAL(a,2,3) * MVAL(b,3,2);
X
XMVAL(mptr,2,3) =
X     MVAL(a,2,0) * MVAL(b,0,3)
X  +  MVAL(a,2,1) * MVAL(b,1,3)
X  +  MVAL(a,2,2) * MVAL(b,2,3)
X  +  MVAL(a,2,3) * MVAL(b,3,3);
X
XMVAL(mptr,3,0) =
X     MVAL(a,3,0) * MVAL(b,0,0)
X  +  MVAL(a,3,1) * MVAL(b,1,0)
X  +  MVAL(a,3,2) * MVAL(b,2,0)
X  +  MVAL(a,3,3) * MVAL(b,3,0);
X
XMVAL(mptr,3,1) =
X     MVAL(a,3,0) * MVAL(b,0,1)
X  +  MVAL(a,3,1) * MVAL(b,1,1)
X  +  MVAL(a,3,2) * MVAL(b,2,1)
X  +  MVAL(a,3,3) * MVAL(b,3,1);
X
XMVAL(mptr,3,2) =
X     MVAL(a,3,0) * MVAL(b,0,2)
X  +  MVAL(a,3,1) * MVAL(b,1,2)
X  +  MVAL(a,3,2) * MVAL(b,2,2)
X  +  MVAL(a,3,3) * MVAL(b,3,2);
X
XMVAL(mptr,3,3) =
X     MVAL(a,3,0) * MVAL(b,0,3)
X  +  MVAL(a,3,1) * MVAL(b,1,3)
X  +  MVAL(a,3,2) * MVAL(b,2,3)
X  +  MVAL(a,3,3) * MVAL(b,3,3);
X
X/* copy temp matrix to result if needed */
Xif ( usetemp )
X  *result = *mptr;
X
Xreturn result;
X}
X
X
X
X
X
X/*  ************************************
X    *                                  *
X    *        V3AffMatMul               *
X    *                                  *
X    ************************************
X
X    DESCRIPTION:  Multiply two affine 4x4 matricies.  The
X      routine assumes the rightmost column of each input
X      matrix is [0 0 0 1].  The output matrix will have the
X      same property.
X    
X      If one of the input matrices is the same as the output,
X      write the result to a temporary matrix during multiplication,
X      then copy to the output matrix.
X
X    ENTRY:
X      a -- pointer to left matrix
X      b -- pointer to right matrix
X      result -- result matrix
X
X    EXIT:  returns 'result'
X*/
X
XMatrix4 *V3AffMatMul ( a, b, result )
Xregister Matrix4 *a,*b;
XMatrix4 *result;
X{
Xregister Matrix4 *mptr;
Xint usetemp;  /* boolean */
XMatrix4 tempx;
X
X/* decide where intermediate result goes */
Xusetemp = ( a == result  ||  b == result );
Xif ( usetemp )
X  mptr = & tempx;
Xelse
X  mptr = result;
X
XMVAL(mptr,0,0) =
X     MVAL(a,0,0) * MVAL(b,0,0)
X  +  MVAL(a,0,1) * MVAL(b,1,0)
X  +  MVAL(a,0,2) * MVAL(b,2,0);
X
XMVAL(mptr,0,1) =
X     MVAL(a,0,0) * MVAL(b,0,1)
X  +  MVAL(a,0,1) * MVAL(b,1,1)
X  +  MVAL(a,0,2) * MVAL(b,2,1);
X
XMVAL(mptr,0,2) =
X     MVAL(a,0,0) * MVAL(b,0,2)
X  +  MVAL(a,0,1) * MVAL(b,1,2)
X  +  MVAL(a,0,2) * MVAL(b,2,2);
X
XMVAL(mptr,0,3) = 0.0;
X
XMVAL(mptr,1,0) =
X     MVAL(a,1,0) * MVAL(b,0,0)
X  +  MVAL(a,1,1) * MVAL(b,1,0)
X  +  MVAL(a,1,2) * MVAL(b,2,0);
X
XMVAL(mptr,1,1) =
X     MVAL(a,1,0) * MVAL(b,0,1)
X  +  MVAL(a,1,1) * MVAL(b,1,1)
X  +  MVAL(a,1,2) * MVAL(b,2,1);
X
XMVAL(mptr,1,2) =
X     MVAL(a,1,0) * MVAL(b,0,2)
X  +  MVAL(a,1,1) * MVAL(b,1,2)
X  +  MVAL(a,1,2) * MVAL(b,2,2);
X
XMVAL(mptr,1,3) = 0.0;
X
XMVAL(mptr,2,0) =
X     MVAL(a,2,0) * MVAL(b,0,0)
X  +  MVAL(a,2,1) * MVAL(b,1,0)
X  +  MVAL(a,2,2) * MVAL(b,2,0);
X
XMVAL(mptr,2,1) =
X     MVAL(a,2,0) * MVAL(b,0,1)
X  +  MVAL(a,2,1) * MVAL(b,1,1)
X  +  MVAL(a,2,2) * MVAL(b,2,1);
X
XMVAL(mptr,2,2) =
X     MVAL(a,2,0) * MVAL(b,0,2)
X  +  MVAL(a,2,1) * MVAL(b,1,2)
X  +  MVAL(a,2,2) * MVAL(b,2,2);
X
XMVAL(mptr,2,3) = 0.0;
X
XMVAL(mptr,3,0) =
X     MVAL(a,3,0) * MVAL(b,0,0)
X  +  MVAL(a,3,1) * MVAL(b,1,0)
X  +  MVAL(a,3,2) * MVAL(b,2,0)
X  +                MVAL(b,3,0);
X
XMVAL(mptr,3,1) =
X     MVAL(a,3,0) * MVAL(b,0,1)
X  +  MVAL(a,3,1) * MVAL(b,1,1)
X  +  MVAL(a,3,2) * MVAL(b,2,1)
X  +                MVAL(b,3,1);
X
XMVAL(mptr,3,2) =
X     MVAL(a,3,0) * MVAL(b,0,2)
X  +  MVAL(a,3,1) * MVAL(b,1,2)
X  +  MVAL(a,3,2) * MVAL(b,2,2)
X  +                MVAL(b,3,2);
X
XMVAL(mptr,3,3) = 1.0;
X
X/* copy temp matrix to result if needed */
Xif ( usetemp )
X  *result = *mptr;
X
Xreturn result;
X}
X
X
X
X
X
X/*  ************************************
X    *                                  *
X    *        V3LinMatMul               *
X    *                                  *
X    ************************************
X
X    DESCRIPTION:  Multiply two affine 4x4 matricies.  The
X      routine assumes the right column and bottom line
X      of each input matrix is [0 0 0 1].  The output matrix
X      will have the same property.  This is pretty much the
X      same thing as multiplying two 3x3 matrices.
X    
X      If one of the input matrices is the same as the output,
X      write the result to a temporary matrix during multiplication,
X      then copy to the output matrix.
X
X    ENTRY:
X      a -- pointer to left matrix
X      b -- pointer to right matrix
X      result -- result matrix
X
X    EXIT:  returns 'result'
X*/
X
XMatrix4 *V3LinMatMul ( a, b, result )
Xregister Matrix4 *a,*b;
XMatrix4 *result;
X{
Xregister Matrix4 *mptr;
Xint usetemp;  /* boolean */
XMatrix4 tempx;
X
X/* decide where intermediate result goes */
Xusetemp = ( a == result  ||  b == result );
Xif ( usetemp )
X  mptr = & tempx;
Xelse
X  mptr = result;
X
XMVAL(mptr,0,0) =
X     MVAL(a,0,0) * MVAL(b,0,0)
X  +  MVAL(a,0,1) * MVAL(b,1,0)
X  +  MVAL(a,0,2) * MVAL(b,2,0);
X
XMVAL(mptr,0,1) =
X     MVAL(a,0,0) * MVAL(b,0,1)
X  +  MVAL(a,0,1) * MVAL(b,1,1)
X  +  MVAL(a,0,2) * MVAL(b,2,1);
X
XMVAL(mptr,0,2) =
X     MVAL(a,0,0) * MVAL(b,0,2)
X  +  MVAL(a,0,1) * MVAL(b,1,2)
X  +  MVAL(a,0,2) * MVAL(b,2,2);
X
XMVAL(mptr,0,3) = 0.0;
X
XMVAL(mptr,1,0) =
X     MVAL(a,1,0) * MVAL(b,0,0)
X  +  MVAL(a,1,1) * MVAL(b,1,0)
X  +  MVAL(a,1,2) * MVAL(b,2,0);
X
XMVAL(mptr,1,1) =
X     MVAL(a,1,0) * MVAL(b,0,1)
X  +  MVAL(a,1,1) * MVAL(b,1,1)
X  +  MVAL(a,1,2) * MVAL(b,2,1);
X
XMVAL(mptr,1,2) =
X     MVAL(a,1,0) * MVAL(b,0,2)
X  +  MVAL(a,1,1) * MVAL(b,1,2)
X  +  MVAL(a,1,2) * MVAL(b,2,2);
X
XMVAL(mptr,1,3) = 0.0;
X
XMVAL(mptr,2,0) =
X     MVAL(a,2,0) * MVAL(b,0,0)
X  +  MVAL(a,2,1) * MVAL(b,1,0)
X  +  MVAL(a,2,2) * MVAL(b,2,0);
X
XMVAL(mptr,2,1) =
X     MVAL(a,2,0) * MVAL(b,0,1)
X  +  MVAL(a,2,1) * MVAL(b,1,1)
X  +  MVAL(a,2,2) * MVAL(b,2,1);
X
XMVAL(mptr,2,2) =
X     MVAL(a,2,0) * MVAL(b,0,2)
X  +  MVAL(a,2,1) * MVAL(b,1,2)
X  +  MVAL(a,2,2) * MVAL(b,2,2);
X
XMVAL(mptr,2,3) = 0.0;
X
XMVAL(mptr,3,0) = 0.0;
XMVAL(mptr,3,1) = 0.0;
XMVAL(mptr,3,2) = 0.0;
XMVAL(mptr,3,3) = 1.0;
X
X/* copy temp matrix to result if needed */
Xif ( usetemp )
X  *result = *mptr;
X
Xreturn result;
X}
END_OF_FILE
if test 12335 -ne `wc -c <'AALines/FastMatMul.c'`; then
    echo shar: \"'AALines/FastMatMul.c'\" unpacked with wrong size!
fi
# end of 'AALines/FastMatMul.c'
fi
if test -f 'FitCurves.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'FitCurves.c'\"
else
echo shar: Extracting \"'FitCurves.c'\" \(14420 characters\)
sed "s/^X//" >'FitCurves.c' <<'END_OF_FILE'
X/*
XAn Algorithm for Automatically Fitting Digitized Curves
Xby Philip J. Schneider
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#define TESTMODE
X
X/*  fit_cubic.c	*/									
X/*	Piecewise cubic fitting code	*/
X
X#include "GraphicsGems.h"					
X#include <stdio.h>
X#include <malloc.h>
X#include <math.h>
X
Xtypedef Point2 *BezierCurve;
X
X/* Forward declarations */
Xvoid		FitCurve();
Xstatic	void		FitCubic();
Xstatic	double		*Reparameterize();
Xstatic	double		NewtonRaphsonRootFind();
Xstatic	Point2		Bezier();
Xstatic	double 		B0(), B1(), B2(), B3();
Xstatic	Vector2		ComputeLeftTangent();
Xstatic	Vector2		ComputeRightTangent();
Xstatic	Vector2		ComputeCenterTangent();
Xstatic	double		ComputeMaxError();
Xstatic	double		*ChordLengthParameterize();
Xstatic	BezierCurve	GenerateBezier();
Xstatic	Vector2		V2AddII();
Xstatic	Vector2		V2ScaleII();
Xstatic	Vector2		V2SubII();
X
X#define MAXPOINTS	1000		/* The most points you can have */
X
X#ifdef TESTMODE
X
XDrawBezierCurve(n, curve)
Xint n;
XBezierCurve curve;
X{
X	/* You'll have to write this yourself. */
X}
X
X/*
X *  main:
X *	Example of how to use the curve-fitting code.  Given an array
X *   of points and a tolerance (squared error between points and 
X *	fitted curve), the algorithm will generate a piecewise
X *	cubic Bezier representation that approximates the points.
X *	When a cubic is generated, the routine "DrawBezierCurve"
X *	is called, which outputs the Bezier curve just created
X *	(arguments are the degree and the control points, respectively).
X *	Users will have to implement this function themselves 	
X *   ascii output, etc. 
X *
X */
Xmain()
X{
X    static Point2 d[7] = {	/*  Digitized points */
X	{ 0.0, 0.0 },
X	{ 0.0, 0.5 },
X	{ 1.1, 1.4 },
X	{ 2.1, 1.6 },
X	{ 3.2, 1.1 },
X	{ 4.0, 0.2 },
X	{ 4.0, 0.0 },
X    };
X    double	error = 4.0;		/*  Squared error */
X    FitCurve(d, 7, error);		/*  Fit the Bezier curves */
X}
X#endif						 /* TESTMODE */
X
X/*
X *  FitCurve :
X *  	Fit a Bezier curve to a set of digitized points 
X */
Xvoid FitCurve(d, nPts, error)
X    Point2	*d;			/*  Array of digitized points	*/
X    int		nPts;		/*  Number of digitized points	*/
X    double	error;		/*  User-defined error squared	*/
X{
X    Vector2	tHat1, tHat2;	/*  Unit tangent vectors at endpoints */
X
X    tHat1 = ComputeLeftTangent(d, 0);
X    tHat2 = ComputeRightTangent(d, nPts - 1);
X    FitCubic(d, 0, nPts - 1, tHat1, tHat2, error);
X}
X
X
X
X/*
X *  FitCubic :
X *  	Fit a Bezier curve to a (sub)set of digitized points
X */
Xstatic void FitCubic(d, first, last, tHat1, tHat2, error)
X    Point2	*d;			/*  Array of digitized points */
X    int		first, last;	/* Indices of first and last pts in region */
X    Vector2	tHat1, tHat2;	/* Unit tangent vectors at endpoints */
X    double	error;		/*  User-defined error squared	   */
X{
X    BezierCurve	bezCurve; /*Control points of fitted Bezier curve*/
X    double	*u;		/*  Parameter values for point  */
X    double	*uPrime;	/*  Improved parameter values */
X    double	maxError;	/*  Maximum fitting error	 */
X    int		splitPoint;	/*  Point to split point set at	 */
X    int		nPts;		/*  Number of points in subset  */
X    double	iterationError; /*Error below which you try iterating  */
X    int		maxIterations = 4; /*  Max times to try iterating  */
X    Vector2	tHatCenter;   	/* Unit tangent vector at splitPoint */
X    int		i;		
X
X    iterationError = error * error;
X    nPts = last - first + 1;
X
X    /*  Use heuristic if region only has two points in it */
X    if (nPts == 2) {
X	    double dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0;
X
X		bezCurve = (Point2 *)malloc(4 * sizeof(Point2));
X		bezCurve[0] = d[first];
X		bezCurve[3] = d[last];
X		V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]);
X		V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]);
X		DrawBezierCurve(3, bezCurve);
X		return;
X    }
X
X    /*  Parameterize points, and attempt to fit curve */
X    u = ChordLengthParameterize(d, first, last);
X    bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2);
X
X    /*  Find max deviation of points to fitted curve */
X    maxError = ComputeMaxError(d, first, last, bezCurve, u, &splitPoint);
X    if (maxError < error) {
X		DrawBezierCurve(3, bezCurve);
X		return;
X    }
X
X
X    /*  If error not too large, try some reparameterization  */
X    /*  and iteration */
X    if (maxError < iterationError) {
X		for (i = 0; i < maxIterations; i++) {
X	    	uPrime = Reparameterize(d, first, last, u, bezCurve);
X	    	bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2);
X	    	maxError = ComputeMaxError(d, first, last,
X				       bezCurve, uPrime, &splitPoint);
X	    	if (maxError < error) {
X			DrawBezierCurve(3, bezCurve);
X			return;
X	    }
X	    free((char *)u);
X	    u = uPrime;
X	}
X}
X
X    /* Fitting failed -- split at max error point and fit recursively */
X    tHatCenter = ComputeCenterTangent(d, splitPoint);
X    FitCubic(d, first, splitPoint, tHat1, tHatCenter, error);
X    V2Negate(&tHatCenter);
X    FitCubic(d, splitPoint, last, tHatCenter, tHat2, error);
X}
X
X
X/*
X *  GenerateBezier :
X *  Use least-squares method to find Bezier control points for region.
X *
X */
Xstatic BezierCurve  GenerateBezier(d, first, last, uPrime, tHat1, tHat2)
X    Point2	*d;			/*  Array of digitized points	*/
X    int		first, last;		/*  Indices defining region	*/
X    double	*uPrime;		/*  Parameter values for region */
X    Vector2	tHat1, tHat2;	/*  Unit tangents at endpoints	*/
X{
X    int 	i;
X    Vector2 	A[MAXPOINTS][2];	/* Precomputed rhs for eqn	*/
X    int 	nPts;			/* Number of pts in sub-curve */
X    double 	C[2][2];			/* Matrix C		*/
X    double 	X[2];			/* Matrix X			*/
X    double 	det_C0_C1,		/* Determinants of matrices	*/
X    	   	det_C0_X,
X	   		det_X_C1;
X    double 	alpha_l,		/* Alpha values, left and right	*/
X    	   	alpha_r;
X    Vector2 	tmp;			/* Utility variable		*/
X    BezierCurve	bezCurve;	/* RETURN bezier curve ctl pts	*/
X
X    bezCurve = (Point2 *)malloc(4 * sizeof(Point2));
X    nPts = last - first + 1;
X
X 
X    /* Compute the A's	*/
X    for (i = 0; i < nPts; i++) {
X		Vector2		v1, v2;
X		v1 = tHat1;
X		v2 = tHat2;
X		V2Scale(&v1, B1(uPrime[i]));
X		V2Scale(&v2, B2(uPrime[i]));
X		A[i][0] = v1;
X		A[i][1] = v2;
X    }
X
X    /* Create the C and X matrices	*/
X    C[0][0] = 0.0;
X    C[0][1] = 0.0;
X    C[1][0] = 0.0;
X    C[1][1] = 0.0;
X    X[0]    = 0.0;
X    X[1]    = 0.0;
X
X    for (i = 0; i < nPts; i++) {
X        C[0][0] += V2Dot(&A[i][0], &A[i][0]);
X		C[0][1] += V2Dot(&A[i][0], &A[i][1]);
X/*					C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/	
X		C[1][0] = C[0][1];
X		C[1][1] += V2Dot(&A[i][1], &A[i][1]);
X
X		tmp = V2SubII(d[first + i],
X	        V2AddII(
X	          V2ScaleII(d[first], B0(uPrime[i])),
X		    	V2AddII(
X		      		V2ScaleII(d[first], B1(uPrime[i])),
X		        			V2AddII(
X	                  		V2ScaleII(d[last], B2(uPrime[i])),
X	                    		V2ScaleII(d[last], B3(uPrime[i]))))));
X	
X
X	X[0] += V2Dot(&A[i][0], &tmp);
X	X[1] += V2Dot(&A[i][1], &tmp);
X    }
X
X    /* Compute the determinants of C and X	*/
X    det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
X    det_C0_X  = C[0][0] * X[1]    - C[0][1] * X[0];
X    det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];
X
X    /* Finally, derive alpha values	*/
X    if (det_C0_C1 == 0.0) {
X		det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
X    }
X    alpha_l = det_X_C1 / det_C0_C1;
X    alpha_r = det_C0_X / det_C0_C1;
X
X
X    /*  If alpha negative, use the Wu/Barsky heuristic (see text) */
X    if (alpha_l < 0.0 || alpha_r < 0.0) {
X		double	dist = V2DistanceBetween2Points(&d[last], &d[first]) /
X					3.0;
X
X		bezCurve[0] = d[first];
X		bezCurve[3] = d[last];
X		V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]);
X		V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]);
X		return (bezCurve);
X    }
X
X    /*  First and last control points of the Bezier curve are */
X    /*  positioned exactly at the first and last data points */
X    /*  Control points 1 and 2 are positioned an alpha distance out */
X    /*  on the tangent vectors, left and right, respectively */
X    bezCurve[0] = d[first];
X    bezCurve[3] = d[last];
X    V2Add(&bezCurve[0], V2Scale(&tHat1, alpha_l), &bezCurve[1]);
X    V2Add(&bezCurve[3], V2Scale(&tHat2, alpha_r), &bezCurve[2]);
X    return (bezCurve);
X}
X
X
X/*
X *  Reparameterize:
X *	Given set of points and their parameterization, try to find
X *   a better parameterization.
X *
X */
Xstatic double *Reparameterize(d, first, last, u, bezCurve)
X    Point2	*d;			/*  Array of digitized points	*/
X    int		first, last;		/*  Indices defining region	*/
X    double	*u;			/*  Current parameter values	*/
X    BezierCurve	bezCurve;	/*  Current fitted curve	*/
X{
X    int 	nPts = last-first+1;	
X    int 	i;
X    double	*uPrime;		/*  New parameter values	*/
X
X    uPrime = (double *)malloc(nPts * sizeof(double));
X    for (i = first; i <= last; i++) {
X		uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-
X					first]);
X    }
X    return (uPrime);
X}
X
X
X
X/*
X *  NewtonRaphsonRootFind :
X *	Use Newton-Raphson iteration to find better root.
X */
Xstatic double NewtonRaphsonRootFind(Q, P, u)
X    BezierCurve	Q;			/*  Current fitted curve	*/
X    Point2 		P;		/*  Digitized point		*/
X    double 		u;		/*  Parameter value for "P"	*/
X{
X    double 		numerator, denominator;
X    Point2 		Q1[3], Q2[2];	/*  Q' and Q''			*/
X    Point2		Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q''	*/
X    double 		uPrime;		/*  Improved u			*/
X    int 		i;
X    
X    /* Compute Q(u)	*/
X    Q_u = Bezier(3, Q, u);
X    
X    /* Generate control vertices for Q'	*/
X    for (i = 0; i <= 2; i++) {
X		Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0;
X		Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0;
X    }
X    
X    /* Generate control vertices for Q'' */
X    for (i = 0; i <= 1; i++) {
X		Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0;
X		Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0;
X    }
X    
X    /* Compute Q'(u) and Q''(u)	*/
X    Q1_u = Bezier(2, Q1, u);
X    Q2_u = Bezier(1, Q2, u);
X    
X    /* Compute f(u)/f'(u) */
X    numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y);
X    denominator = (Q1_u.x) * (Q1_u.x) + (Q1_u.y) * (Q1_u.y) +
X		      	  (Q_u.x - P.x) * (Q2_u.x) + (Q_u.y - P.y) * (Q2_u.y);
X    
X    /* u = u - f(u)/f'(u) */
X    uPrime = u - (numerator/denominator);
X    return (uPrime);
X}
X
X	
X		       
X/*
X *  Bezier :
X *  	Evaluate a Bezier curve at a particular parameter value
X * 
X */
Xstatic Point2 Bezier(degree, V, t)
X    int		degree;		/* The degree of the bezier curve	*/
X    Point2 	*V;		/* Array of control points		*/
X    double 	t;		/* Parametric value to find point for	*/
X{
X    int 	i, j;		
X    Point2 	Q;	        /* Point on curve at parameter t	*/
X    Point2 	*Vtemp;		/* Local copy of control points		*/
X
X    /* Copy array	*/
X    Vtemp = (Point2 *)malloc((unsigned)((degree+1) 
X				* sizeof (Point2)));
X    for (i = 0; i <= degree; i++) {
X		Vtemp[i] = V[i];
X    }
X
X    /* Triangle computation	*/
X    for (i = 1; i <= degree; i++) {	
X		for (j = 0; j <= degree-i; j++) {
X	    	Vtemp[j].x = (1.0 - t) * Vtemp[j].x + t * Vtemp[j+1].x;
X	    	Vtemp[j].y = (1.0 - t) * Vtemp[j].y + t * Vtemp[j+1].y;
X		}
X    }
X
X    Q = Vtemp[0];
X    free((char *)Vtemp);
X    return Q;
X}
X
X
X/*
X *  B0, B1, B2, B3 :
X *	Bezier multipliers
X */
Xstatic double B0(u)
X    double	u;
X{
X    double tmp = 1.0 - u;
X    return (tmp * tmp * tmp);
X}
X
X
Xstatic double B1(u)
X    double	u;
X{
X    double tmp = 1.0 - u;
X    return (3 * u * (tmp * tmp));
X}
X
Xstatic double B2(u)
X    double	u;
X{
X    double tmp = 1.0 - u;
X    return (3 * u * u * tmp);
X}
X
Xstatic double B3(u)
X    double	u;
X{
X    return (u * u * u);
X}
X
X
X
X/*
X * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
X *Approximate unit tangents at endpoints and "center" of digitized curve
X */
Xstatic Vector2 ComputeLeftTangent(d, end)
X    Point2	*d;			/*  Digitized points*/
X    int		end;		/*  Index to "left" end of region */
X{
X    Vector2	tHat1;
X    tHat1 = V2SubII(d[end+1], d[end]);
X    tHat1 = *V2Normalize(&tHat1);
X    return tHat1;
X}
X
Xstatic Vector2 ComputeRightTangent(d, end)
X    Point2	*d;			/*  Digitized points		*/
X    int		end;		/*  Index to "right" end of region */
X{
X    Vector2	tHat2;
X    tHat2 = V2SubII(d[end-1], d[end]);
X    tHat2 = *V2Normalize(&tHat2);
X    return tHat2;
X}
X
X
Xstatic Vector2 ComputeCenterTangent(d, center)
X    Point2	*d;			/*  Digitized points			*/
X    int		center;		/*  Index to point inside region	*/
X{
X    Vector2	V1, V2, tHatCenter;
X
X    V1 = V2SubII(d[center-1], d[center]);
X    V2 = V2SubII(d[center], d[center+1]);
X    tHatCenter.x = (V1.x + V2.x)/2.0;
X    tHatCenter.y = (V1.y + V2.y)/2.0;
X    tHatCenter = *V2Normalize(&tHatCenter);
X    return tHatCenter;
X}
X
X
X/*
X *  ChordLengthParameterize :
X *	Assign parameter values to digitized points 
X *	using relative distances between points.
X */
Xstatic double *ChordLengthParameterize(d, first, last)
X    Point2	*d;			/* Array of digitized points */
X    int		first, last;		/*  Indices defining region	*/
X{
X    int		i;	
X    double	*u;			/*  Parameterization		*/
X
X    u = (double *)malloc((unsigned)(last-first+1) * sizeof(double));
X
X    u[0] = 0.0;
X    for (i = first+1; i <= last; i++) {
X		u[i-first] = u[i-first-1] +
X	  			V2DistanceBetween2Points(&d[i], &d[i-1]);
X    }
X
X    for (i = first + 1; i <= last; i++) {
X		u[i-first] = u[i-first] / u[last-first];
X    }
X
X    return(u);
X}
X
X
X
X
X/*
X *  ComputeMaxError :
X *	Find the maximum squared distance of digitized points
X *	to fitted curve.
X*/
Xstatic double ComputeMaxError(d, first, last, bezCurve, u, splitPoint)
X    Point2	*d;			/*  Array of digitized points	*/
X    int		first, last;		/*  Indices defining region	*/
X    BezierCurve	bezCurve;		/*  Fitted Bezier curve		*/
X    double	*u;			/*  Parameterization of points	*/
X    int		*splitPoint;		/*  Point of maximum error	*/
X{
X    int		i;
X    double	maxDist;		/*  Maximum error		*/
X    double	dist;		/*  Current error		*/
X    Point2	P;			/*  Point on curve		*/
X    Vector2	v;			/*  Vector from point to curve	*/
X
X    *splitPoint = (last - first + 1)/2;
X    maxDist = 0.0;
X    for (i = first + 1; i < last; i++) {
X		P = Bezier(3, bezCurve, u[i-first]);
X		v = V2SubII(P, d[i]);
X		dist = V2SquaredLength(&v);
X		if (dist >= maxDist) {
X	    	maxDist = dist;
X	    	*splitPoint = i;
X		}
X    }
X    return (maxDist);
X}
Xstatic Vector2 V2AddII(a, b)
X    Vector2 a, b;
X{
X    Vector2	c;
X    c.x = a.x + b.x;  c.y = a.y + b.y;
X    return (c);
X}
Xstatic Vector2 V2ScaleII(v, s)
X    Vector2	v;
X    double	s;
X{
X    Vector2 result;
X    result.x = v.x * s; result.y = v.y * s;
X    return (result);
X}
X
Xstatic Vector2 V2SubII(a, b)
X    Vector2	a, b;
X{
X    Vector2	c;
X    c.x = a.x - b.x; c.y = a.y - b.y;
X    return (c);
X}
END_OF_FILE
if test 14420 -ne `wc -c <'FitCurves.c'`; then
    echo shar: \"'FitCurves.c'\" unpacked with wrong size!
fi
# end of 'FitCurves.c'
fi
if test -f 'GGVecLib.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'GGVecLib.c'\"
else
echo shar: Extracting \"'GGVecLib.c'\" \(10162 characters\)
sed "s/^X//" >'GGVecLib.c' <<'END_OF_FILE'
X/* 
X2d and 3d Vector C Library 
Xby Andrew Glassner
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#include <math.h>
X#include "GraphicsGems.h"
X
X/******************/
X/*   2d Library   */
X/******************/
X
X/* returns squared length of input vector */	
Xdouble V2SquaredLength(a) 
XVector2 *a;
X{	return((a->x * a->x)+(a->y * a->y));
X	};
X	
X/* returns length of input vector */
Xdouble V2Length(a) 
XVector2 *a;
X{
X	return(sqrt(V2SquaredLength(a)));
X	};
X	
X/* negates the input vector and returns it */
XVector2 *V2Negate(v) 
XVector2 *v;
X{
X	v->x = -v->x;  v->y = -v->y;
X	return(v);
X	};
X
X/* normalizes the input vector and returns it */
XVector2 *V2Normalize(v) 
XVector2 *v;
X{
Xdouble len = V2Length(v);
X	if (len != 0.0) { v->x /= len;  v->y /= len; };
X	return(v);
X	};
X
X
X/* scales the input vector to the new length and returns it */
XVector2 *V2Scale(v, newlen) 
XVector2 *v;
Xdouble newlen;
X{
Xdouble len = V2Length(v);
X	if (len != 0.0) { v->x *= newlen/len;   v->y *= newlen/len; };
X	return(v);
X	};
X
X/* return vector sum c = a+b */
XVector2 *V2Add(a, b, c)
XVector2 *a, *b, *c;
X{
X	c->x = a->x+b->x;  c->y = a->y+b->y;
X	return(c);
X	};
X	
X/* return vector difference c = a-b */
XVector2 *V2Sub(a, b, c)
XVector2 *a, *b, *c;
X{
X	c->x = a->x-b->x;  c->y = a->y-b->y;
X	return(c);
X	};
X
X/* return the dot product of vectors a and b */
Xdouble V2Dot(a, b) 
XVector2 *a, *b; 
X{
X	return((a->x*b->x)+(a->y*b->y));
X	};
X
X/* linearly interpolate between vectors by an amount alpha */
X/* and return the resulting vector. */
X/* When alpha=0, result=lo.  When alpha=1, result=hi. */
XVector2 *V2Lerp(lo, hi, alpha, result) 
XVector2 *lo, *hi, *result; 
Xdouble alpha;
X{
X	result->x = LERP(alpha, lo->x, hi->x);
X	result->y = LERP(alpha, lo->y, hi->y);
X	return(result);
X	};
X
X
X/* make a linear combination of two vectors and return the result. */
X/* result = (a * ascl) + (b * bscl) */
XVector2 *V2Combine (a, b, result, ascl, bscl) 
XVector2 *a, *b, *result;
Xdouble ascl, bscl;
X{
X	result->x = (ascl * a->x) + (bscl * b->x);
X	result->y = (ascl * a->y) + (bscl * b->y);
X	return(result);
X	};
X
X/* multiply two vectors together component-wise */
XVector2 *V2Mul (a, b, result) 
XVector2 *a, *b, *result;
X{
X	result->x = a->x * b->x;
X	result->y = a->y * b->y;
X	return(result);
X	};
X
X/* return the distance between two points */
Xdouble V2DistanceBetween2Points(a, b)
XPoint2 *a, *b;
X{
Xdouble dx = a->x - b->x;
Xdouble dy = a->y - b->y;
X	return(sqrt((dx*dx)+(dy*dy)));
X	};
X
X/* return the vector perpendicular to the input vector a */
XVector2 *V2MakePerpendicular(a, ap)
XVector2 *a, *ap;
X{
X	ap->x = -a->y;
X	ap->y = a->x;
X	return(ap);
X	};
X
X/* create, initialize, and return a new vector */
XVector2 *V2New(x, y)
Xdouble x, y;
X{
XVector2 *v = NEWTYPE(Vector2);
X	v->x = x;  v->y = y; 
X	return(v);
X	};
X	
X
X/* create, initialize, and return a duplicate vector */
XVector2 *V2Duplicate(a)
XVector2 *a;
X{
XVector2 *v = NEWTYPE(Vector2);
X	v->x = a->x;  v->y = a->y; 
X	return(v);
X	};
X	
X/* multiply a point by a matrix and return the transformed point */
XPoint2 *V2MulPointByMatrix(p, m)
XPoint2 *p;
XMatrix3 *m;
X{
Xdouble w;
XPoint2 ptmp;
X	ptmp.x = (p->x * m->element[0][0]) + 
X             (p->y * m->element[1][0]) + m->element[2][0];
X	ptmp.y = (p->x * m->element[0][1]) + 
X             (p->y * m->element[1][1]) + m->element[2][1];
X	w    = (p->x * m->element[0][2]) + 
X             (p->y * m->element[1][2]) + m->element[2][2];
X	if (w != 0.0) { ptmp.x /= w;  ptmp.y /= w; }
X	*p = ptmp;
X	return(p);
X	};
X
X/* multiply together matrices c = ab */
X/* note that c must not point to either of the input matrices */
XMatrix3 *V2MatMul(a, b, c)
XMatrix3 *a, *b, *c;
X{
Xint i, j, k;
X	for (i=0; i<3; i++) {
X		for (j=0; j<3; j++) {
X			c->element[i][j] = 0;
X		for (k=0; k<3; k++) c->element[i][j] += 
X				a->element[i][k] * b->element[k][j];
X			};
X		};
X	return(c);
X	};
X
X
X
X
X/******************/
X/*   3d Library   */
X/******************/
X	
X/* returns squared length of input vector */	
Xdouble V3SquaredLength(a) 
XVector3 *a;
X{
X	return((a->x * a->x)+(a->y * a->y)+(a->z * a->z));
X	};
X
X/* returns length of input vector */
Xdouble V3Length(a) 
XVector3 *a;
X{
X	return(sqrt(V3SquaredLength(a)));
X	};
X
X/* negates the input vector and returns it */
XVector3 *V3Negate(v) 
XVector3 *v;
X{
X	v->x = -v->x;  v->y = -v->y;  v->z = -v->z;
X	return(v);
X	};
X
X/* normalizes the input vector and returns it */
XVector3 *V3Normalize(v) 
XVector3 *v;
X{
Xdouble len = V3Length(v);
X	if (len != 0.0) { v->x /= len;  v->y /= len; v->z /= len; };
X	return(v);
X	};
X
X/* scales the input vector to the new length and returns it */
XVector3 *V3Scale(v, newlen) 
XVector3 *v;
Xdouble newlen;
X{
Xdouble len = V3Length(v);
X	if (len != 0.0) {
X	v->x *= newlen/len;   v->y *= newlen/len;  v->z *= newlen/len;
X	};
X	return(v);
X	};
X
X
X/* return vector sum c = a+b */
XVector3 *V3Add(a, b, c)
XVector3 *a, *b, *c;
X{
X	c->x = a->x+b->x;  c->y = a->y+b->y;  c->z = a->z+b->z;
X	return(c);
X	};
X	
X/* return vector difference c = a-b */
XVector3 *V3Sub(a, b, c)
XVector3 *a, *b, *c;
X{
X	c->x = a->x-b->x;  c->y = a->y-b->y;  c->z = a->z-b->z;
X	return(c);
X	};
X
X/* return the dot product of vectors a and b */
Xdouble V3Dot(a, b) 
XVector3 *a, *b; 
X{
X	return((a->x*b->x)+(a->y*b->y)+(a->z*b->z));
X	};
X
X/* linearly interpolate between vectors by an amount alpha */
X/* and return the resulting vector. */
X/* When alpha=0, result=lo.  When alpha=1, result=hi. */
XVector3 *V3Lerp(lo, hi, alpha, result) 
XVector3 *lo, *hi, *result; 
Xdouble alpha;
X{
X	result->x = LERP(alpha, lo->x, hi->x);
X	result->y = LERP(alpha, lo->y, hi->y);
X	result->z = LERP(alpha, lo->z, hi->z);
X	return(result);
X	};
X
X/* make a linear combination of two vectors and return the result. */
X/* result = (a * ascl) + (b * bscl) */
XVector3 *V3Combine (a, b, result, ascl, bscl) 
XVector3 *a, *b, *result;
Xdouble ascl, bscl;
X{
X	result->x = (ascl * a->x) + (bscl * b->x);
X	result->y = (ascl * a->y) + (bscl * b->y);
X	result->y = (ascl * a->z) + (bscl * b->z);
X	return(result);
X	};
X
X
X/* multiply two vectors together component-wise and return the result */
XVector3 *V3Mul (a, b, result) 
XVector3 *a, *b, *result;
X{
X	result->x = a->x * b->x;
X	result->y = a->y * b->y;
X	result->z = a->z * b->z;
X	return(result);
X	};
X
X/* return the distance between two points */
Xdouble V3DistanceBetween2Points(a, b)
XPoint3 *a, *b;
X{
Xdouble dx = a->x - b->x;
Xdouble dy = a->y - b->y;
Xdouble dz = a->z - b->z;
X	return(sqrt((dx*dx)+(dy*dy)+(dz*dz)));
X	};
X
X/* return the cross product c = a cross b */
XVector3 *V3Cross(a, b, c)
XVector3 *a, *b, *c;
X{
X	c->x = (a->y*b->z) - (a->z*b->y);
X	c->y = (a->z*b->x) - (a->x*b->z);
X	c->z = (a->x*b->y) - (a->y*b->x);
X	return(c);
X	};
X
X/* create, initialize, and return a new vector */
XVector3 *V3New(x, y, z)
Xdouble x, y, z;
X{
XVector3 *v = NEWTYPE(Vector3);
X	v->x = x;  v->y = y;  v->z = z;
X	return(v);
X	};
X
X/* create, initialize, and return a duplicate vector */
XVector3 *V3Duplicate(a)
XVector3 *a;
X{
XVector3 *v = NEWTYPE(Vector3);
X	v->x = a->x;  v->y = a->y;  v->z = a->z;
X	return(v);
X	};
X
X	
X/* multiply a point by a matrix and return the transformed point */
XPoint3 *V3MulPointByMatrix(p, m)
XPoint3 *p;
XMatrix4 *m;
X{
Xdouble w;
XPoint3 ptmp;
X	ptmp.x = (p->x * m->element[0][0]) + (p->y * m->element[1][0]) + 
X		 (p->z * m->element[2][0]) + m->element[3][0];
X	ptmp.y = (p->x * m->element[0][1]) + (p->y * m->element[1][1]) + 
X		 (p->z * m->element[2][1]) + m->element[3][1];
X	ptmp.z = (p->x * m->element[0][2]) + (p->y * m->element[1][2]) + 
X		 (p->z * m->element[2][2]) + m->element[3][2];
X	w =    (p->x * m->element[0][3]) + (p->y * m->element[1][3]) + 
X		 (p->z * m->element[2][3]) + m->element[3][3];
X	if (w != 0.0) { ptmp.x /= w;  ptmp.y /= w;  ptmp.z /= w; }
X	*p = ptmp;
X	return(p);
X	};
X
X/* multiply together matrices c = ab */
X/* note that c must not point to either of the input matrices */
XMatrix4 *V3MatMul(a, b, c)
XMatrix4 *a, *b, *c;
X{
Xint i, j, k;
X	for (i=0; i<4; i++) {
X		for (j=0; j<4; j++) {
X			c->element[i][j] = 0;
X			for (k=0; k<4; k++) c->element[i][j] += 
X				a->element[i][k] * b->element[k][j];
X			};
X		};
X	return(c);
X	};
X
X/* binary greatest common divisor by Silver and Terzian.  See Knuth */
X/* both inputs must be >= 0 */
Xgcd(u, v)
Xint u, v;
X{
Xint k, t, f;
X	if ((u<0) || (v<0)) return(1); /* error if u<0 or v<0 */
X	k = 0;  f = 1;
X	while ((0 == (u%2)) && (0 == (v%2))) {
X		k++;  u>>=1;  v>>=1,  f*=2;
X		};
X	if (u&01) { t = -v;  goto B4; } else { t = u; }
X	B3: if (t > 0) { t >>= 1; } else { t = -((-t) >> 1); };
X	B4: if (0 == (t%2)) goto B3;
X
X	if (t > 0) u = t; else v = -t;
X	if (0 != (t = u - v)) goto B3;
X	return(u*f);
X	};	
X
X/***********************/
X/*   Useful Routines   */
X/***********************/
X
X/* return roots of ax^2+bx+c */
X/* stable algebra derived from Numerical Recipes by Press et al.*/
Xint quadraticRoots(a, b, c, roots)
Xdouble a, b, c, *roots;
X{
Xdouble d, q;
Xint count = 0;
X	d = (b*b)-(4*a*c);
X	if (d < 0.0) { *roots = *(roots+1) = 0.0;  return(0); };
X	q =  -0.5 * (b + (SGN(b)*sqrt(d)));
X	if (a != 0.0)  { *roots++ = q/a; count++; }
X	if (q != 0.0) { *roots++ = c/q; count++; }
X	return(count);
X	};
X
X
X/* generic 1d regula-falsi step.  f is function to evaluate */
X/* interval known to contain root is given in left, right */
X/* returns new estimate */
Xdouble RegulaFalsi(f, left, right)
Xdouble (*f)(), left, right;
X{
Xdouble d = (*f)(right) - (*f)(left);
X	if (d != 0.0) return (right - (*f)(right)*(right-left)/d);
X	return((left+right)/2.0);
X	};
X
X/* generic 1d Newton-Raphson step. f is function, df is derivative */
X/* x is current best guess for root location. Returns new estimate */
Xdouble NewtonRaphson(f, df, x)
Xdouble (*f)(), (*df)(), x;
X{
Xdouble d = (*df)(x);
X	if (d != 0.0) return (x-((*f)(x)/d));
X	return(x-1.0);
X	};
X
X
X/* hybrid 1d Newton-Raphson/Regula Falsi root finder. */
X/* input function f and its derivative df, an interval */
X/* left, right known to contain the root, and an error tolerance */
X/* Based on Blinn */
Xdouble findroot(left, right, tolerance, f, df)
Xdouble left, right, tolerance;
Xdouble (*f)(), (*df)();
X{
Xdouble newx = left;
X	while (ABS((*f)(newx)) > tolerance) {
X		newx = NewtonRaphson(f, df, newx);
X		if (newx < left || newx > right) 
X			newx = RegulaFalsi(f, left, right);
X		if ((*f)(newx) * (*f)(left) <= 0.0) right = newx;  
X			else left = newx;
X		};
X	return(newx);
X	}; 
END_OF_FILE
if test 10162 -ne `wc -c <'GGVecLib.c'`; then
    echo shar: \"'GGVecLib.c'\" unpacked with wrong size!
fi
# end of 'GGVecLib.c'
fi
if test -f 'NearestPoint.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'NearestPoint.c'\"
else
echo shar: Extracting \"'NearestPoint.c'\" \(12269 characters\)
sed "s/^X//" >'NearestPoint.c' <<'END_OF_FILE'
X/*
XSolving the Nearest Point-on-Curve Problem 
Xand
XA Bezier Curve-Based Root-Finder
Xby Philip J. Schneider
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X /*	point_on_curve.c	*/		
X									
X#include <stdio.h>
X#include <malloc.h>
X#include <math.h>
X#include "GraphicsGems.h"
X
X#define TESTMODE
X
X/*
X *  Forward declarations
X */
XPoint2  NearestPointOnCurve();
Xstatic	int	FindRoots();
Xstatic	Point2	*ConvertToBezierForm();
Xstatic	double	ComputeXIntercept();
Xstatic	int	ControlPolygonFlatEnough();
Xstatic	int	CrossingCount();
Xstatic	Point2	Bezier();
Xstatic	Vector2	V2ScaleII();
X
Xint		MAXDEPTH = 64;	/*  Maximum depth for recursion */
X
X#define	EPSILON	(ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */
X#define	DEGREE	3			/*  Cubic Bezier curve		*/
X#define	W_DEGREE 5			/*  Degree of eqn to find roots of */
X
X#ifdef TESTMODE
X/*
X *  main :
X *	Given a cubic Bezier curve (i.e., its control points), and some
X *	arbitrary point in the plane, find the point on the curve
X *	closest to that arbitrary point.
X */
Xmain()
X{
X   
X static Point2 bezCurve[4] = {	/*  A cubic Bezier curve	*/
X	{ 0.0, 0.0 },
X	{ 1.0, 2.0 },
X	{ 3.0, 3.0 },
X	{ 4.0, 2.0 },
X    };
X    static Point2 arbPoint = { 3.5, 2.0 }; /*Some arbitrary point*/
X    Point2	pointOnCurve;		 /*  Nearest point on the curve */
X
X    /*  Find the closest point */
X    pointOnCurve = NearestPointOnCurve(arbPoint, bezCurve);
X    printf("pointOnCurve : (%4.4f, %4.4f)\n", pointOnCurve.x,
X		pointOnCurve.y);
X}
X#endif /* TESTMODE */
X
X
X/*
X *  NearestPointOnCurve :
X *  	Compute the parameter value of the point on a Bezier
X *		curve segment closest to some arbtitrary, user-input point.
X *		Return the point on the curve at that parameter value.
X *
X */
XPoint2 NearestPointOnCurve(P, V)
X    Point2 	P;			/* The user-supplied point	  */
X    Point2 	*V;			/* Control points of cubic Bezier */
X{
X    Point2	*w;			/* Ctl pts for 5th-degree eqn	*/
X    double 	t_candidate[W_DEGREE];	/* Possible roots		*/     
X    int 	n_solutions;		/* Number of roots found	*/
X    double	t;			/* Parameter value of closest pt*/
X
X    /*  Convert problem to 5th-degree Bezier form	*/
X    w = ConvertToBezierForm(P, V);
X
X    /* Find all possible roots of 5th-degree equation */
X    n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0);
X    free((char *)w);
X
X    /* Compare distances of P to all candidates, and to t=0, and t=1 */
X    {
X		double 	dist, new_dist;
X		Point2 	p;
X		Vector2  v;
X		int		i;
X
X	
X	/* Check distance to beginning of curve, where t = 0	*/
X		dist = V2SquaredLength(V2Sub(&P, &V[0], &v));
X        	t = 0.0;
X
X	/* Find distances for candidate points	*/
X        for (i = 0; i < n_solutions; i++) {
X	    	p = Bezier(V, DEGREE, t_candidate[i], NULL, NULL);
X	    	new_dist = V2SquaredLength(V2Sub(&P, &p, &v));
X	    	if (new_dist < dist) {
X                	dist = new_dist;
X	        		t = t_candidate[i];
X    	    }
X        }
X
X	/* Finally, look at distance to end point, where t = 1.0 */
X		new_dist = V2SquaredLength(V2Sub(&P, &V[DEGREE], &v));
X        	if (new_dist < dist) {
X            	dist = new_dist;
X	    	t = 1.0;
X        }
X    }
X
X    /*  Return the point on the curve at parameter value t */
X    printf("t : %4.12f\n", t);
X    return (Bezier(V, DEGREE, t, NULL, NULL));
X}
X
X
X/*
X *  ConvertToBezierForm :
X *		Given a point and a Bezier curve, generate a 5th-degree
X *		Bezier-format equation whose solution finds the point on the
X *      curve nearest the user-defined point.
X */
Xstatic Point2 *ConvertToBezierForm(P, V)
X    Point2 	P;			/* The point to find t for	*/
X    Point2 	*V;			/* The control points		*/
X{
X    int 	i, j, k, m, n, ub, lb;	
X    double 	t;			/* Value of t for point P	*/
X    int 	row, column;		/* Table indices		*/
X    Vector2 	c[DEGREE+1];		/* V(i)'s - P			*/
X    Vector2 	d[DEGREE];		/* V(i+1) - V(i)		*/
X    Point2 	*w;			/* Ctl pts of 5th-degree curve  */
X    double 	cdTable[3][4];		/* Dot product of c, d		*/
X    static double z[3][4] = {	/* Precomputed "z" for cubics	*/
X	{1.0, 0.6, 0.3, 0.1},
X	{0.4, 0.6, 0.6, 0.4},
X	{0.1, 0.3, 0.6, 1.0},
X    };
X
X
X    /*Determine the c's -- these are vectors created by subtracting*/
X    /* point P from each of the control points				*/
X    for (i = 0; i <= DEGREE; i++) {
X		V2Sub(&V[i], &P, &c[i]);
X    }
X    /* Determine the d's -- these are vectors created by subtracting*/
X    /* each control point from the next					*/
X    for (i = 0; i <= DEGREE - 1; i++) { 
X		d[i] = V2ScaleII(V2Sub(&V[i+1], &V[i], &d[i]), 3.0);
X    }
X
X    /* Create the c,d table -- this is a table of dot products of the */
X    /* c's and d's							*/
X    for (row = 0; row <= DEGREE - 1; row++) {
X		for (column = 0; column <= DEGREE; column++) {
X	    	cdTable[row][column] = V2Dot(&d[row], &c[column]);
X		}
X    }
X
X    /* Now, apply the z's to the dot products, on the skew diagonal*/
X    /* Also, set up the x-values, making these "points"		*/
X    w = (Point2 *)malloc((unsigned)(W_DEGREE+1) * sizeof(Point2));
X    for (i = 0; i <= W_DEGREE; i++) {
X		w[i].y = 0.0;
X		w[i].x = (double)(i) / W_DEGREE;
X    }
X
X    n = DEGREE;
X    m = DEGREE-1;
X    for (k = 0; k <= n + m; k++) {
X		lb = MAX(0, k - m);
X		ub = MIN(k, n);
X		for (i = lb; i <= ub; i++) {
X	    	j = k - i;
X	    	w[i+j].y += cdTable[j][i] * z[j][i];
X		}
X    }
X
X    return (w);
X}
X
X
X/*
X *  FindRoots :
X *	Given a 5th-degree equation in Bernstein-Bezier form, find
X *	all of the roots in the interval [0, 1].  Return the number
X *	of roots found.
X */
Xstatic int FindRoots(w, degree, t, depth)
X    Point2 	*w;			/* The control points		*/
X    int 	degree;		/* The degree of the polynomial	*/
X    double 	*t;			/* RETURN candidate t-values	*/
X    int 	depth;		/* The depth of the recursion	*/
X{  
X    int 	i;
X    Point2 	Left[W_DEGREE+1],	/* New left and right 		*/
X    	  	Right[W_DEGREE+1];	/* control polygons		*/
X    int 	left_count,		/* Solution count from		*/
X		right_count;		/* children			*/
X    double 	left_t[W_DEGREE+1],	/* Solutions from kids		*/
X	   		right_t[W_DEGREE+1];
X
X    switch (CrossingCount(w, degree)) {
X       	case 0 : {	/* No solutions here	*/
X	     return 0;	
X	     break;
X	}
X	case 1 : {	/* Unique solution	*/
X	    /* Stop recursion when the tree is deep enough	*/
X	    /* if deep enough, return 1 solution at midpoint 	*/
X	    if (depth >= MAXDEPTH) {
X			t[0] = (w[0].x + w[W_DEGREE].x) / 2.0;
X			return 1;
X	    }
X	    if (ControlPolygonFlatEnough(w, degree)) {
X			t[0] = ComputeXIntercept(w, degree);
X			return 1;
X	    }
X	    break;
X	}
X}
X
X    /* Otherwise, solve recursively after	*/
X    /* subdividing control polygon		*/
X    Bezier(w, degree, 0.5, Left, Right);
X    left_count  = FindRoots(Left,  degree, left_t, depth+1);
X    right_count = FindRoots(Right, degree, right_t, depth+1);
X
X
X    /* Gather solutions together	*/
X    for (i = 0; i < left_count; i++) {
X        t[i] = left_t[i];
X    }
X    for (i = 0; i < right_count; i++) {
X 		t[i+left_count] = right_t[i];
X    }
X
X    /* Send back total number of solutions	*/
X    return (left_count+right_count);
X}
X
X
X/*
X * CrossingCount :
X *	Count the number of times a Bezier control polygon 
X *	crosses the 0-axis. This number is >= the number of roots.
X *
X */
Xstatic int CrossingCount(V, degree)
X    Point2	*V;			/*  Control pts of Bezier curve	*/
X    int		degree;			/*  Degreee of Bezier curve 	*/
X{
X    int 	i;	
X    int 	n_crossings = 0;	/*  Number of zero-crossings	*/
X    int		sign, old_sign;		/*  Sign of coefficients	*/
X
X    sign = old_sign = SGN(V[0].y);
X    for (i = 1; i <= degree; i++) {
X		sign = SGN(V[i].y);
X		if (sign != old_sign) n_crossings++;
X		old_sign = sign;
X    }
X    return n_crossings;
X}
X
X
X
X/*
X *  ControlPolygonFlatEnough :
X *	Check if the control polygon of a Bezier curve is flat enough
X *	for recursive subdivision to bottom out.
X *
X */
Xstatic int ControlPolygonFlatEnough(V, degree)
X    Point2	*V;		/* Control points	*/
X    int 	degree;		/* Degree of polynomial	*/
X{
X    int 	i;			/* Index variable		*/
X    double 	*distance;		/* Distances from pts to line	*/
X    double 	max_distance_above;	/* maximum of these		*/
X    double 	max_distance_below;
X    double 	error;			/* Precision of root		*/
X    Vector2 	t;			/* Vector from V[0] to V[degree]*/
X    double 	intercept_1,
X    	   	intercept_2,
X	   		left_intercept,
X		   	right_intercept;
X    double 	a, b, c;		/* Coefficients of implicit	*/
X    					/* eqn for line from V[0]-V[deg]*/
X
X    /* Find the  perpendicular distance		*/
X    /* from each interior control point to 	*/
X    /* line connecting V[0] and V[degree]	*/
X    distance = (double *)malloc((unsigned)(degree + 1) * 					sizeof(double));
X    {
X	double	abSquared;
X
X	/* Derive the implicit equation for line connecting first *'
X    /*  and last control points */
X	a = V[0].y - V[degree].y;
X	b = V[degree].x - V[0].x;
X	c = V[0].x * V[degree].y - V[degree].x * V[0].y;
X
X	abSquared = (a * a) + (b * b);
X
X        for (i = 1; i < degree; i++) {
X	    /* Compute distance from each of the points to that line	*/
X	    	distance[i] = a * V[i].x + b * V[i].y + c;
X	    	if (distance[i] > 0.0) {
X				distance[i] = (distance[i] * distance[i]) / abSquared;
X	    	}
X	    	if (distance[i] < 0.0) {
X				distance[i] = -((distance[i] * distance[i]) / 						abSquared);
X	    	}
X		}
X    }
X
X
X    /* Find the largest distance	*/
X    max_distance_above = 0.0;
X    max_distance_below = 0.0;
X    for (i = 1; i < degree; i++) {
X		if (distance[i] < 0.0) {
X	    	max_distance_below = MIN(max_distance_below, distance[i]);
X		};
X		if (distance[i] > 0.0) {
X	    	max_distance_above = MAX(max_distance_above, distance[i]);
X		}
X    }
X    free((char *)distance);
X
X    {
X	double	det, dInv;
X	double	a1, b1, c1, a2, b2, c2;
X
X	/*  Implicit equation for zero line */
X	a1 = 0.0;
X	b1 = 1.0;
X	c1 = 0.0;
X
X	/*  Implicit equation for "above" line */
X	a2 = a;
X	b2 = b;
X	c2 = c + max_distance_above;
X
X	det = a1 * b2 - a2 * b1;
X	dInv = 1.0/det;
X	
X	intercept_1 = (b1 * c2 - b2 * c1) * dInv;
X
X	/*  Implicit equation for "below" line */
X	a2 = a;
X	b2 = b;
X	c2 = c + max_distance_below;
X	
X	det = a1 * b2 - a2 * b1;
X	dInv = 1.0/det;
X	
X	intercept_2 = (b1 * c2 - b2 * c1) * dInv;
X    }
X
X    /* Compute intercepts of bounding box	*/
X    left_intercept = MIN(intercept_1, intercept_2);
X    right_intercept = MAX(intercept_1, intercept_2);
X
X    error = 0.5 * (right_intercept-left_intercept);    
X    if (error < EPSILON) {
X		return 1;
X    }
X    else {
X		return 0;
X    }
X}
X
X
X
X/*
X *  ComputeXIntercept :
X *	Compute intersection of chord from first control point to last
X *  	with 0-axis.
X * 
X */
Xstatic double ComputeXIntercept(V, degree)
X    Point2 	*V;			/*  Control points	*/
X    int		degree; 		/*  Degree of curve	*/
X{
X    double	XLK, YLK, XNM, YNM, XMK, YMK;
X    double	det, detInv;
X    double	S, T;
X    double	X, Y;
X
X    XLK = 1.0 - 0.0;
X    YLK = 0.0 - 0.0;
X    XNM = V[degree].x - V[0].x;
X    YNM = V[degree].y - V[0].y;
X    XMK = V[0].x - 0.0;
X    YMK = V[0].y - 0.0;
X
X    det = XNM*YLK - YNM*XLK;
X    detInv = 1.0/det;
X
X    S = (XNM*YMK - YNM*XMK) * detInv;
X    T = (XLK*YMK - YLK*XMK) * detInv;
X    
X    X = 0.0 + XLK * S;
X    Y = 0.0 + YLK * S;
X
X    return X;
X}
X
X
X/*
X *  Bezier : 
X *	Evaluate a Bezier curve at a particular parameter value
X *      Fill in control points for resulting sub-curves if "Left" and
X *	"Right" are non-null.
X * 
X */
Xstatic Point2 Bezier(V, degree, t, Left, Right)
X    int 	degree;		/* Degree of bezier curve	*/
X    Point2 	*V;			/* Control pts			*/
X    double 	t;			/* Parameter value		*/
X    Point2 	*Left;		/* RETURN left half ctl pts	*/
X    Point2 	*Right;		/* RETURN right half ctl pts	*/
X{
X    int 	i, j;		/* Index variables	*/
X    Point2 	Vtemp[W_DEGREE+1][W_DEGREE+1];
X
X
X    /* Copy control points	*/
X    for (j =0; j <= degree; j++) {
X		Vtemp[0][j] = V[j];
X    }
X
X    /* Triangle computation	*/
X    for (i = 1; i <= degree; i++) {	
X		for (j =0 ; j <= degree - i; j++) {
X	    	Vtemp[i][j].x =
X	      		(1.0 - t) * Vtemp[i-1][j].x + t * Vtemp[i-1][j+1].x;
X	    	Vtemp[i][j].y =
X	      		(1.0 - t) * Vtemp[i-1][j].y + t * Vtemp[i-1][j+1].y;
X		}
X    }
X    
X    if (Left != NULL) {
X		for (j = 0; j <= degree; j++) {
X	    	Left[j]  = Vtemp[j][0];
X		}
X    }
X    if (Right != NULL) {
X		for (j = 0; j <= degree; j++) {
X	    	Right[j] = Vtemp[degree-j][j];
X		}
X    }
X
X    return (Vtemp[degree][0]);
X}
X
Xstatic Vector2 V2ScaleII(v, s)
X    Vector2	*v;
X    double	s;
X{
X    Vector2 result;
X
X    result.x = v->x * s; result.y = v->y * s;
X    return (result);
X}
END_OF_FILE
if test 12269 -ne `wc -c <'NearestPoint.c'`; then
    echo shar: \"'NearestPoint.c'\" unpacked with wrong size!
fi
# end of 'NearestPoint.c'
fi
echo shar: End of archive 5 \(of 5\).
cp /dev/null ark5isdone
MISSING=""
for I in 1 2 3 4 5 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 5 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0

  .