please dont rip this site

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


file: /Techref/method/io/graphics/Part05.sh, 52KB, , updated: 1999/6/25 11:02, local time: 2022/8/18 00:47,
TOP NEW HELP FIND: 
18.204.56.185:LOG IN

 ©2022 These pages are served without commercial sponsorship. (No popup ads, etc...).Bandwidth abuse increases hosting cost forcing sponsorship or shutdown. This server aggressively defends against automated copying for any reason including offline viewing, duplication, etc... Please respect this requirement and DO NOT RIP THIS SITE. Questions?
Please DO link to this page! Digg it! / MAKE!

<A HREF="http://massmind.org/techref/method/io/graphics/Part05.sh"> method io graphics Part05</A>

Did you find what you needed?

 

Welcome to massmind.org!

 

Welcome to massmind.org!

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  .