HW2 Solutions Page

To create a solution set, I selected the code of a student who I thought did well on the assignment. For allowing me to use their homework, the student receives 1 bonus point added after curving. If you see your homework used, but do not want me to use it, just ask, and I will choose someone else's homework. They, rather than you, will receive the bonus point. Below are the results of this process for this homework:

sphere.txt

R : 0 : 4
0.0f,  -3.0f,    3.0f 
2.0f,  -1.0f,    1.0f
2.0f,   1.0f,    1.0f
0.0f,   3.0f,    3.0f

Point4f.h

#include <assert.h>
#include <float.h>

/*
 CLASS DEFINITIONS
*/

/*-----------------------------------------------*/
struct Point4fS
/*
 PURPOSE: Provides a static data initializer form for Point4f.
*/
{
   float x, y, z, w;
};

/*-----------------------------------------------*/
class Point4f
/*
 PURPOSE: Encapsulates the notions of a 3D point (x,y,z,1) or
          homogenous 3D point (x,y,z,w).
*/
{
public:
   Point4f();
   Point4f( float ix, float iy );
   Point4f( float ix, float iy, float iz );
   Point4f( float ix, float iy, float iz, float iw );
   Point4f( const Point4fS & p );

   inline float getX() const;
   inline float getY() const;
   inline float getZ() const;
   inline float getW() const;
   inline const float * getRawFloats() const;

   void setW( float iw );
   void set( float ix, float iy, float iz );
   void set( float ix, float iy, float iz, float iw );
   void set( const Point4f & ipoint );

   void startRunningMin();
   void startRunningMax();
   void runningMin( const Point4f & ipoint );
   void runningMax( const Point4f & ipoint );

   Point4f operator * ( float alpha ) const;
   Point4f operator + ( const Point4f & b ) const;
   Point4f operator - ( const Point4f & b ) const;

   void cross( const Point4f & ivector );
   void combine( const Point4f & a, float alpha, const Point4f & b, float beta );
   void lerp( const Point4f & a, const Point4f & b, float alpha );
   float distToLine( const Point4f & a, const Point4f & b ) const;
   float length() const;
   float square_length() const;

	 inline void project();

private:
   float x, y, z, w;
};

#include "Point4f.inl"

#endif

HW2.cpp

/******************************************************
 * Copyright (c):   2010, All Rights Reserved.
 * Project:         CS 216 Homework #2
 * File:            bezier.cpp
 * Purpose:         Code to draw surface of revolution based on a spline curve
 * Start date:      Feb. 1, 2010
 * Programmer:      SJSU Students
 *
 * Remarks: Adapted from Hearn and Baker (3rd edition), and sample code by Chris Pollett
 *
 *******************************************************/

#ifdef WIN32
#include 
#endif

#ifdef __APPLE__ //for MAC
/*
 I created my project under xcode. I chose new C++ tool as the kind of project.
 Then I added the two frameworks:
 OpenGL.framework and GLUT.framework. (Frameworks are in /Library/Frameworks)
 
 */
#include  
#include 
#include 
#endif

#ifdef linux // for linux
/*My compile line was:
 g++ -I /usr/X11R6/include -L /usr/X11R6/lib -lglut -lGL \
 -lGLU -lX11 -lXmu -lXi -lm name.cpp -o name
 */
#include 
#include 
#endif

#include 
#include 
#include 
#include 
#include 

using namespace std;

#include "Point4f.h"

/*
 CONSTANT DEFINITIONS
 */

static const char * FILE_NAME = "sphere.txt";
//static const char * FILE_NAME = "vase_1.txt";
//static const char * FILE_NAME = "fig_7_23_a.txt";
//static const char * FILE_NAME = "fig_7_23_b.txt";
//static const char * FILE_NAME = "fig_7_24_a.txt";
//static const char * FILE_NAME = "fig_7_24_b.txt";

// Defines Pi
static const GLfloat _M_PI = 3.1415926535897932384626433832795f;

// The tessellation resolution for curves and surfaces.
static const int SURFACE_TESSELLATION = 16;

// The number of coordinates in a homogenous 3D vector.
static const int HOMOG_3D_COORDS = 4;

// The degree of a cubic polynomial.
static const int CUBIC_DEGREE = 3;
// The order of a cubic polynomial.
static const int CUBIC_ORDER = CUBIC_DEGREE + 1;
// The number of vertices in a cubic path, in the u-axis direction.
static const int CUBIC_RATIONAL_PATCH_VERTS_U = CUBIC_ORDER;
// The number of vertices in a cubic path, in the v-axis direction.
static const int CUBIC_RATIONAL_PATCH_VERTS_V = CUBIC_ORDER;
// The number of control points to increment for each step in a piecewise cubic
// Bezier curve, where the curve shares start and end points of neighboring 
// segments.
static const int PROFILE_CUBIC_CONTROL_POINT_STEP = CUBIC_DEGREE;

// Short alias for unsigned int.
typedef unsigned int uint;

// Alias for cubic rational curve control point array.
typedef Point4f CubicRational3DBezierCurveControlPoints[ CUBIC_RATIONAL_PATCH_VERTS_U ];
// Alias for cubic rational patch control point array.
typedef Point4f CubicRational3DBezierPatchControlPoints[ CUBIC_RATIONAL_PATCH_VERTS_V ][ CUBIC_RATIONAL_PATCH_VERTS_U ];

// Alias for vector of Point4f.
typedef std::vector< Point4f > Point4fArray;
// Alias for vector of float.
typedef std::vector< float > FloatArray;


// Rational cubic Bezier hemicircle x^2+y^2 = 1, y > 0
const CubicRational3DBezierCurveControlPoints cubic_xy_py_circ =
{
	Point4f(    1.0f,    0.0f,         0,    1.0f ),
	Point4f(  1.0f/3,  2.0f/3,         0,  1.0f/3 ), 
	Point4f( -1.0f/3,  2.0f/3,         0,  1.0f/3 ), 
	Point4f(   -1.0f,    0.0f,         0,    1.0f ) 
};

// Rational cubic Bezier hemicircle x^2+y^2 = 1, y < 0
const CubicRational3DBezierCurveControlPoints cubic_xy_ny_circ =
{
	Point4f(   -1.0f,    0.0f,         0,    1.0f ), 
	Point4f( -1.0f/3, -2.0f/3,         0,  1.0f/3 ), 
	Point4f(  1.0f/3, -2.0f/3,         0,  1.0f/3 ), 
	Point4f(    1.0f,    0.0f,         0,    1.0f )
};

const GLfloat RED_COLOR[ 3 ] = { 1.0f, 0.0f, 0.0f };
const GLfloat GREEN_COLOR[ 3 ] = { 0.0f, 1.0f, 0.0f };
const GLfloat YELLOW_COLOR[ 3 ] = { 1.0f, 1.0f, 0.0f };

/*
 PROTOTYPES
 */

void drawCubicRationalBezierPatchSOR( const Point4f profile_rcbez[], const GLfloat knot_vector[],
									 uint num_points, float max_angle, GLenum fill_mode, bool draw_control_cage );
void drawCubicRationalCatmullRomPatchSOR( const Point4f profile_points[],
										 uint num_points, float max_angle, GLenum fill_mode, bool draw_control_cage );
void drawCubicRationalBesselOverhauserPatchSOR( const Point4f profile_points[],
											   const GLfloat knot_vector[], uint num_points, float max_angle,
											   GLenum fill_mode, bool draw_control_cage );

#define POINT4F_INL

/*-----------------------------------------------*/
Point4f::Point4f()
/*
 PURPOSE: Default constructor, performing no initialization.
 REMARKS: This constructor does not initialize the data, so use with caution.
 */
{
}

/*-----------------------------------------------*/
Point4f::Point4f( GLfloat ix, GLfloat iy ) : x( ix ), y( iy ), z( 0.0f ), w( 1.0f )
/*
 PURPOSE: Initialize a 2D homogenous point with given coordinates.
 RECEIVES: 2D Coordinates.
 */
{
}

/*-----------------------------------------------*/
Point4f::Point4f( GLfloat ix, GLfloat iy, GLfloat iz ) : x( ix ), y( iy ), z( iz ), w( 1.0f )
/*
 PURPOSE: Initialize a 3D homogenous point with given coordinates.
 RECEIVES: 3D Coordinates.
 */
{
}

/*-----------------------------------------------*/
Point4f::Point4f( GLfloat ix, GLfloat iy, GLfloat iz, GLfloat iw ) : x( ix ), y( iy ), z( iz ), w( iw )
/*
 PURPOSE: Initialize a 3D homogenous point with given coordinates.
 RECEIVES: 4D Coordinates.
 */
{
}

/*-----------------------------------------------*/
Point4f::Point4f( const Point4fS & p ) : x( p.x ), y( p.y ), z( p.z ), w( p.w )
/*
 PURPOSE: Initialize a 3D homogenous point with given coordinates.
 RECEIVES: 4D Coordinates.
 */
{
}


/*-----------------------------------------------*/
inline GLfloat Point4f::getX() const
/*
 PURPOSE: Retrieve the X-coordinate.
 RETURNS: The X-coordinate.
 */
{
	return x;
}

/*-----------------------------------------------*/
inline GLfloat Point4f::getY() const
/*
 PURPOSE: Retrieve the Y-coordinate.
 RETURNS: The Y-coordinate.
 */
{
	return y;
}

/*-----------------------------------------------*/
inline GLfloat Point4f::getZ() const
/*
 PURPOSE: Retrieve the Z-coordinate.
 RETURNS: The Z-coordinate.
 */
{
	return z;
}

/*-----------------------------------------------*/
inline GLfloat Point4f::getW() const
/*
 PURPOSE: Retrieve the W-coordinate.
 RETURNS: The W-coordinate.
 */
{
	return w;
}

/*-----------------------------------------------*/
inline const GLfloat * Point4f::getRawFloats() const
/*
 PURPOSE: Access the coordinates as an array.
 RETURNS: A pointer to a four-element array of coordinates, { X, Y, Z, W }.
 */
{
	assert( &y == &x + 1 && &z == &y + 1 && &w == &z + 1 );
	return &x;
}

/*-----------------------------------------------*/
void Point4f::setW( GLfloat iw )
/*
 PURPOSE: Modifies the W-coordinate.
 RECEIVES: W coordinate.
 */
{
	w = iw;
}


/*-----------------------------------------------*/
void Point4f::set( GLfloat ix, GLfloat iy, GLfloat iz )
/*
 PURPOSE: Initialize a point with given coordinates.
 RECEIVES: Coordinates.
 */
{
	x = ix;
	y = iy;
	z = iz;
	w = 1.0f;
}

/*-----------------------------------------------*/
void Point4f::set( GLfloat ix, GLfloat iy, GLfloat iz, GLfloat iw )
/*
 PURPOSE: Initialize a point with given coordinates.
 RECEIVES: 4D Coordinates.
 */
{
	x = ix;
	y = iy;
	z = iz;
	w = iw;
}

/*-----------------------------------------------*/
void Point4f::set( const Point4f & ipoint )
/*
 PURPOSE: Copy coordinates from another point.
 RECEIVES: A 4D point.
 */
{
	x = ipoint.x;
	y = ipoint.y;
	z = ipoint.z;
	w = ipoint.w;
}

/*-----------------------------------------------*/
void Point4f::startRunningMin()
/*
 PURPOSE: Initialize the point to track a running maximum.
 */
{
	x = FLT_MAX;
	x *= 2.0f;
	y = z = w = x;
}


/*-----------------------------------------------*/
void Point4f::startRunningMax()
/*
 PURPOSE: Initialize the point to track a running maximum.
 */
{
	x = FLT_MAX;
	x *= -2.0f;
	y = z = w = x;
}


/*-----------------------------------------------*/
void Point4f::runningMin( const Point4f & ipoint )
/*
 PURPOSE: Accumulate a point into a running minimum in this point.
 RECEIVES: A 4D point to merge into this point.
 */
{
	if ( ipoint.x < x )
		x = ipoint.x;
	if ( ipoint.y < y )
		y = ipoint.y;
	if ( ipoint.z < z )
		z = ipoint.z;
	if ( ipoint.w < w )
		w = ipoint.w;
}

/*-----------------------------------------------*/
void Point4f::runningMax( const Point4f & ipoint )
/*
 PURPOSE: Accumulate a point into a running maximum in this point.
 RECEIVES: A 4D point to merge into this point.
 */
{
	if ( ipoint.x > x )
		x = ipoint.x;
	if ( ipoint.y > y )
		y = ipoint.y;
	if ( ipoint.z > z )
		z = ipoint.z;
	if ( ipoint.w > w )
		w = ipoint.w;
}


/*-----------------------------------------------*/
inline Point4f Point4f::operator * ( GLfloat alpha ) const
/*
 PURPOSE: Scale a point by a scalar.
 RECEIVES: Scale factor.
 REMARKS: Side effect is scaling the coordinates of this point by the given scale factor.
 */
{
	return Point4f( x * alpha, y * alpha, z * alpha, w * alpha );
}

/*-----------------------------------------------*/
inline Point4f Point4f::operator + ( const Point4f & b ) const
/*
 PURPOSE: Sums two points.
 RECEIVES: The right-hand side of the addition expression.
 RETURNS: Sum of this point and the addend
 */
{
	return Point4f( x + b.x, y + b.y, z + b.z, w + b.w );
}

/*-----------------------------------------------*/
inline Point4f Point4f::operator -( const Point4f & b ) const
/*
 PURPOSE: Subtracts two points.
 RECEIVES: The right-hand side of the subtraction expression.
 RETURNS: Difference of this point and the right-hand side point
 */
{
	return Point4f( x - b.x, y - b.y, z - b.z, w - b.w );
}

/*-----------------------------------------------*/
void Point4f::cross( const Point4f & ivector )
/*
 PURPOSE: Generate the cross product of this vector and a right-hand-side vector.
 RECEIVES: The right-hand-side vector.
 REMARKS: Side effect is modifying this point to contain the cross product.
 */
{
	/*
	 |  i  j  k |
	 |  x  y  z |
	 | ix iy iz |
	 */
	set( y * ivector.z - z * ivector.y, z * ivector.x - x * ivector.z, x * ivector.y - y * ivector.x );
}

/*-----------------------------------------------*/
void Point4f::combine( const Point4f & a, GLfloat alpha, const Point4f & b, GLfloat beta )
/*
 PURPOSE: Generate a linear combination of two points.
 RECEIVES: Two points and a weight for each point.
 REMARKS: Side effect is modifying this point to contain the combined coordinates.
 */
{
	*this = a * alpha + b * beta;
}

/*-----------------------------------------------*/
void Point4f::lerp( const Point4f & a, const Point4f & b, GLfloat alpha )
/*
 PURPOSE: Linearly interpolate between two points.
 RECEIVES: Two points and a weight for the second point.
 REMARKS: Side effect is modifying this point to contain the interpolated position.
 */
{
	Point4f v = b - a;
	Point4f u = *this - a;
	*this = a * ( 1 - alpha ) + b * alpha;
}

/*-----------------------------------------------*/
GLfloat Point4f::distToLine( const Point4f & a, const Point4f & b ) const
/*
 PURPOSE: Calculate distance from a point to a line.
 RECEIVES: Two distinct points on the line.
 RETURNS: The distance of this point to two line formed by the two given points.
 */
{
	Point4f v = b - a;
	Point4f u = *this - a;
	u.cross( v );
	return static_cast( sqrt( u.square_length() / v.square_length() ) );
}

/*-----------------------------------------------*/
GLfloat Point4f::length() const
/*
 PURPOSE: Calculate length of vector from origin to the point.
 RETURNS: The distance of this point to the origin.
 */
{
	return static_cast( sqrt( square_length() ) );
}

/*-----------------------------------------------*/
GLfloat Point4f::square_length() const
/*
 PURPOSE: Calculate squared length of (3D) vector from origin to the point.
 RETURNS: The squared distance of this point to the origin.
 */
{
	return x * x + y * y + z * z;
}

/*-----------------------------------------------*/
void Point4f::project()
/*
 PURPOSE: Project the homogenous 3D point.
 REMARKS: Side effect is to modify the point. Point must not be at infinity.
 */
{
	assert( w != 0.0f ); 
	x /= w;
	y /= w;
	z /= w;
	w = 1.0f;
}


/*
 CLASS DEFINITIONS
 */

/*-----------------------------------------------*/
class SurfaceData
/*
 PURPOSE: This class encapsulates the data describing a surface of
 revolution formed by rotating a curve around the Z-axis.
 */
{
private:
	// Stores either the control points or interpolated points of the curve.
	Point4fArray vertices;
	// Stores the knot vector.
	FloatArray knots;
	// A char defined by the CurveType enumeration indicating the kind of curve.
	char type;
	// The maximum sweep angle of the surface of revolution about the Z-axis.
	float angle;
	
	// The minimum extent point of the axis-aligned bounding box of the vertices.
	Point4f min;
	// The maximum extent point of the axis-aligned bounding box of the vertices.
	Point4f max;
	
	// For parsing input files.
	FILE * stream;
	
public:
	SurfaceData();
	
	int readData( const char * file );
	void draw( GLenum fill_mode, bool draw_control_cage );
	
	void getExtent( Point4f & bboxMin, Point4f & bboxMax ) const;
	
	enum CurveType
	{
		CT_Rational_Bezier  = 'R',
		CT_Catmull_Rom = 'C',
		CT_Bessel_Overhauser = 'B',
	};
	
private:
	int ParseError();
	
	char ProcessChar( const char * token );
	float ProcessFloat( const char * token );
	uint ProcessUInt( const char * token );
};

/*-----------------------------------------------*/
SurfaceData::SurfaceData() : type( CT_Rational_Bezier ), angle( 0.0f ), stream( NULL )
/*
 PURPOSE: Construct a SurfaceData object read for use with readData.
 REMARKS: The curve is "empty" after the constructor. Initialize it with readData.
 */
{
}

/*-----------------------------------------------*/
int SurfaceData::readData( const char * file )
/*
 PURPOSE: Parse a profile curve file.
 RECEIVES: A NULL-terminated filename indicating the curve data file.
 RETURNS: Non-zero indicates success; zero indicates failure and the curve 
 may be in an indeterminate state.
 REMARKS: The function reads curve data into the object from the given file.
 The file format is as given in the homework assignment:
 
 X:angle_in_radians:num_points
 point_1
 point_2
 
 X: one of {R, C, B} with interpretation defined by the CurveType enumeration.
 angle_in_radians: a floating-point number.
 
 num_points: a positive integer indicating the number of points on lines 
 following the header.
 
 point_1, point_2, etc: comma-delimited pairs or triples indicating the XZ shape 
 of the profile curve, but prepended with a knot value for the Bessel-Overhauser 
 spline type.
 */
{
	stream = fopen( file, "r" );
	if ( !stream )
	{
		printf( "Missing file %s\n", file );
		return 0;
	}
	
	uint num_points = 0;
	char text[ 256 ];
	if ( !fgets( text, 256, stream ) )
		return ParseError();
	
	type = ProcessChar( strtok( text, ": \t" ) );
	angle = ProcessFloat( strtok( NULL, ": \t" ) );
	num_points = ProcessUInt( strtok( NULL, ": \t" ) );
	
	min.startRunningMin();
	max.startRunningMax();
	
	Point4f data;
	for ( num_points; num_points && fgets( text, 256, stream ) != NULL; --num_points )
	{
		float a = 0.0f, b = 0.0f, c = 0.0f;
		a = ProcessFloat( strtok( text, ", \t" ) );
		b = ProcessFloat( strtok( NULL, ", \t" ) );
		
		if ( type == CT_Rational_Bezier )
		{
			c = ProcessFloat( strtok( NULL, ", \t" ) );
			data.set( a, 0.0f, b, c );
		}
		else
			if ( type == CT_Catmull_Rom )
			{
				data.set( a, 0.0f, b, 1.0f );
			}
			else
			{
				knots.push_back( a );
				c = ProcessFloat( strtok( NULL, ", \t" ) );
				data.set( b, 0.0f, c, 1.0f );
			}
		
		min.runningMin( data );
		max.runningMax( data );
		vertices.push_back( data );
	}
	
	if ( num_points != 0 )
		return ParseError();
	
	fclose( stream );
	return 1;
}

/*-----------------------------------------------*/
int SurfaceData::ParseError()
/*
 PURPOSE: Cancel parsing and return an error code of zero.
 REMARKS: The curve is in an indeterminate state after a parse error, except
 the input stream is guaranteed to be closed.
 */
{
	fclose( stream );
	stream = NULL;
	
	printf( "Parse error!\n" );
	return 0;
}

/*-----------------------------------------------*/
char SurfaceData::ProcessChar( const char * token )
/*
 PURPOSE: Process a single-character token.
 RECEIVES: An input string token.
 RETURNS: The character value, or NULL in case of a parse error such as an 
 invalid (non-character) token.
 REMARKS: The function issues a parse error if the input is not a string
 containing a single character.
 */
{
	if ( strlen( token ) != 1 )
	{
		ParseError();
		return 0;
	}
	
	return token[ 0 ];
}

/*-----------------------------------------------*/
float SurfaceData::ProcessFloat( const char * token )
/*
 PURPOSE: Process a floating-point number token.
 RECEIVES: An input string token containing a floating-point number.
 RETURNS: The floating-point number value.
 REMARKS: The function does not perform any validation of the token, and
 silently fails if the number is malformed.
 */
{
	return atof( token );
}

/*-----------------------------------------------*/
uint SurfaceData::ProcessUInt( const char * token )
/*
 PURPOSE: Process an unsigned integer number token.
 RECEIVES: An input string token containing a non-negative number.
 RETURNS: The integer number value.
 REMARKS: The function does not perform any validation of the token, and
 silently fails if the number is malformed.
 */
{
	uint temp = 0;
	sscanf( token, "%u", &temp );
	return temp;
}

/*-----------------------------------------------*/
void SurfaceData::draw( GLenum fill_mode, bool draw_control_cage )
/*
 PURPOSE: Displays a surface of revolution given the profile curve defined in 
 the SurfaceData object, over the sweep angle limit given by the angle data member.
 RECEIVES: fill_mode: An OpenGL evaluator fill mode.
 draw_control_cage: true indicates draw the surface and the control polygon 
 mesh; false indicates draw only the surface.
 REMARKS: Side effect is to modify OpenGL state to initialize surface tessellation
 and to send the surface to OpenGL.
 */
{
	switch ( type )
	{
		case CT_Rational_Bezier:
			drawCubicRationalBezierPatchSOR( &vertices[ 0 ], NULL,
											vertices.size(), angle, fill_mode, draw_control_cage );
			break;
			
		case CT_Catmull_Rom:
			drawCubicRationalCatmullRomPatchSOR( &vertices[ 0 ], 
												vertices.size(), angle, fill_mode, draw_control_cage );
			break;
			
		case CT_Bessel_Overhauser:
			drawCubicRationalBesselOverhauserPatchSOR( &vertices[ 0 ],
													  knots.size() ? &knots[ 0 ] : NULL, vertices.size(), angle, fill_mode,
													  draw_control_cage );
			break;
	}
}

/*-----------------------------------------------*/
void SurfaceData::getExtent( Point4f & bboxMin, Point4f & bboxMax ) const
/*
 PURPOSE: Retrieves extent of the profile curve. 
 RECEIVES: bboxMin: The method fills this with the axis-aligned bounding box minimum.
 bboxMax: The method fills this with the axis-aligned bounding box maximum.
 */
{
	bboxMin = min;
	bboxMax = max;
}


/*
 GLOBAL VARIABLES
 */

// This data encapsulates the surface data for use by the GLUT callback functions.
SurfaceData surface;

/*
 FUNCTIONS
 */


/*-----------------------------------------------*/
void deCasteljeauCubicRationalBezier( const Point4f & p0, const Point4f & p1,
									 const Point4f & p2, const Point4f & p3, Point4f & r0, Point4f & r1, Point4f & r2,
									 Point4f & s0, Point4f & s1, Point4f & t0, float u )
/*
 PURPOSE: Construct the de Casteljeau interpolation points.
 RECEIVES: p0, p1, p2, p3: The four Bezier curve control points.
 r0, r1, r2: The method constructs the first stage interpolation points in these points.
 s0, s1: The method constructs the second stage interpolation points in these points.
 t0: The method constructs the final stage point on the curve in this point.
 RETURNS: The function modifies the input ri, si, and t0.
 */
{
	r0.lerp( p0, p1, u );
	r1.lerp( p1, p2, u );
	r2.lerp( p2, p3, u );
	
	s0.lerp( r0, r1, u );
	s1.lerp( r1, r2, u );
	
	t0.lerp( s0, s1, u );
}

/*-----------------------------------------------*/
void subdivideCubicRationalBezierFirstHalf( const Point4f & p0, Point4f & p1,
										   Point4f & p2, Point4f & p3, float u )
/*
 PURPOSE: Subdivide a cubic rational curve at the given parameter and return the
 control points for the first half of the subdivided curve.
 RECEIVES: The four Bezier curve control points.
 RETURNS: The function modifies the last three of the input control points to contain
 the subdivided curve. The initial point remains the same.
 */
{
	Point4f r0, r1, r2, s0, s1, t0;
	deCasteljeauCubicRationalBezier( p0, p1, p2, p3, r0, r1, r2, s0, s1, t0, u );
	
	// take the first half
	p1 = r0;
	p2 = s0;
	p3 = t0;
}

/*-----------------------------------------------*/
float findBezierParameterRS( Point4f p0, Point4f p1,
							Point4f p2, Point4f p3, const Point4f & root, float distance_epsilon )
/*
 PURPOSE: Search a 2D cubic rational Bezier curve for a root.
 RECEIVES: p0, p1, p2, p3: The four Bezier curve control points
 root: the root point
 distance_epsilon: the root accuracy distance tolerance
 RETURNS: The parameter value of the root.
 */
{
	Point4f r0, r1, r2, s0, s1, t0, t0_p;
	float u0 = 0.0f, u1 = 1.0f;
	bool bRootFound = false;
	
	distance_epsilon = distance_epsilon * distance_epsilon;
	
	do
	{
		deCasteljeauCubicRationalBezier( p0, p1, p2, p3, r0, r1, r2, s0, s1, t0, 0.5f );
		
		// project the midpoint to compare with the root point
		t0_p.set( t0 );
		t0_p.project();
		
		if ( !( bRootFound = ( t0_p - root ).square_length() <= distance_epsilon ) )
		{
			// decide which half of the curve to subdivide into
			if ( root.getX() > t0_p.getX() )
			{
				// take the first half
				p1 = r0;
				p2 = s0;
				p3 = t0;
				u1 = ( u0 + u1 ) * 0.5f;
			}
			else
			{
				// take the second half
				p0 = t0;
				p1 = s1;
				p2 = r2;
				u0 = ( u0 + u1 ) * 0.5f;
			}
		}
	}
	while ( !bRootFound );
	
	return ( u0 + u1 ) * 0.5f;
}


/*-----------------------------------------------*/
float findBezierCircleParameter( const Point4f & root, float distance_epsilon )
/*
 PURPOSE: Search a 2D cubic rational Bezier circle curve for a root.
 RECEIVES: root: The root point
 distance_epsilon: The root accuracy distance tolerance
 RETURNS: The parameter value of the root.
 */
{
	// Note the point order here gives a hemi-ciricle curve that sweeps from -X towards +Y
	// then to the +X axis. This is opposite the sweep of positive
	// angles in Cartesian coordinates, but doesn't affect the root finding
	return findBezierParameterRS(
								 cubic_xy_py_circ[ 0 ], cubic_xy_py_circ[ 1 ], cubic_xy_py_circ[ 2 ], cubic_xy_py_circ[ 3 ],
								 root, distance_epsilon );
}

/*-----------------------------------------------*/
void initializeCubicSolidOfRevPatch( const CubicRational3DBezierCurveControlPoints circle_arc,
									const Point4f profile[], CubicRational3DBezierPatchControlPoints patch )
/*
 PURPOSE: Given a rational cubic Bezier profile curve to rotate about the
 Z-axis, and a rational cubic Bezier circle arc section, generate a rational
 cubic Bezier patch for the correspondig part of the surface.
 RECEIVES: circle_arc: an arc of the circle (on the Z = 0 plane) defined by a rational cubic Bezier.
 profile: a cubic Bezier profile curve on the Y = 0 plane.
 patch: the function stores the patch control points into this 2D array.
 */
{
	for ( int y = 0; y < CUBIC_RATIONAL_PATCH_VERTS_V; ++y )
	{
		for ( int v = 0; v < CUBIC_RATIONAL_PATCH_VERTS_U; ++v )
		{
			// scale and translate the rational unit circle curve to form a 
			// rational circle around the axis, interpolating the desired
			// homogeneous profile control point. The X- and Y-axis have the
			// same scale factor, at least for a circular solid of revolution.
			float xy_s = profile[ y ].getX();
			float z_s = profile[ y ].getX();
			float w_s = profile[ y ].getW();
			// Only Z-axis has translation
			float z_t = profile[ y ].getZ();
			
			// Perform a homogenous affine transform on the vertices
			patch[ v ][ y ].set(
								circle_arc[ v ].getX() * xy_s,
								circle_arc[ v ].getY() * xy_s,
								circle_arc[ v ].getZ() * z_s + circle_arc[ v ].getW() * z_t,
								circle_arc[ v ].getW() * w_s );
		}
	}
}


/*-----------------------------------------------*/
void drawCubicBezierControlCage( const CubicRational3DBezierPatchControlPoints patch_verts )
/*
 PURPOSE: Displays the wireframe of the control polygon mesh or cage of the patch.
 RECEIVES: The patch vertices.
 REMARKS: The function sends line data to OpenGL.
 */
{
	glColor3fv( GREEN_COLOR );
	
	glBegin( GL_LINES );
	
	for ( int i = 0; i < CUBIC_RATIONAL_PATCH_VERTS_U; ++i )
		for ( int j = 0; j < CUBIC_RATIONAL_PATCH_VERTS_V - 1; ++j )
		{
			glVertex4fv( patch_verts[j][i].getRawFloats() );
			glVertex4fv( patch_verts[j+1][i].getRawFloats() );
		}
	
	for ( int j = 0; j < CUBIC_RATIONAL_PATCH_VERTS_V; ++j )
		for ( int i = 0; i < CUBIC_RATIONAL_PATCH_VERTS_U - 1; ++i )
		{
			glVertex4fv( patch_verts[j][i].getRawFloats() );
			glVertex4fv( patch_verts[j][i+1].getRawFloats() );
		}
	
	glEnd();
	
	glColor3fv( YELLOW_COLOR );
	glPointSize( 3.0f );
	
	glBegin( GL_POINTS );
	
	for ( int i = 0; i < CUBIC_RATIONAL_PATCH_VERTS_U; ++i )
		for ( int j = 0; j < CUBIC_RATIONAL_PATCH_VERTS_V; ++j )
		{
			glVertex4fv( patch_verts[j][i].getRawFloats() );
		}
	
	glEnd();
	
	glColor3fv( RED_COLOR );
	glPointSize( 1.0f );
}

/*-----------------------------------------------*/
bool initializeCircleArcs( float max_angle, Point4f cross_section_1[],
						  Point4f cross_section_2[] )
/*
 PURPOSE: Initialize up-to two rational cubic Bezier curves to cover the sweep
 of angle from zero to max_angle radians around the Z-axis.
 RECEIVES: max_angle: The curves rotate around the Z-axis from zero radians 
 to this angle in radians, so 2*Pi radians specifies a complete surface.
 cross_section_1: The curve covering zero to up-to Pi radians.
 cross_section_2: The curve covering from Pi radians to 2*Pi radians,
 if max_angle > Pi.
 RETURNS: True indicates the second curve is necessary to complete the entire arc
 of the circle; false indicates only the first curve is required.
 */
{
	bool bIncludeSecondHalf = true;
	memcpy( cross_section_1, cubic_xy_py_circ, sizeof( cubic_xy_py_circ ) );
	memcpy( cross_section_2, cubic_xy_ny_circ, sizeof( cubic_xy_ny_circ ) );
	
	if ( max_angle && max_angle < 2 * _M_PI )
	{
		if ( max_angle > _M_PI )
		{
			max_angle -= _M_PI;
			float max_angle_u = findBezierCircleParameter( Point4f( cos( max_angle ), sin( max_angle ), 0.0f ), 1.0e-5f );
			
			subdivideCubicRationalBezierFirstHalf( cross_section_2[ 0 ], cross_section_2[ 1 ], 
												  cross_section_2[ 2 ], cross_section_2[ 3 ], max_angle_u );
		}
		else
		{
			bIncludeSecondHalf = false;
			
			if ( max_angle < _M_PI )
			{
				float max_angle_u = findBezierCircleParameter( Point4f( cos( max_angle ), sin( max_angle ), 0.0f ), 1.0e-5f );
				subdivideCubicRationalBezierFirstHalf( cross_section_1[ 0 ], cross_section_1[ 1 ], 
													  cross_section_1[ 2 ], cross_section_1[ 3 ], max_angle_u );
			}
		}
	}
	return bIncludeSecondHalf;
}

/*-----------------------------------------------*/
void drawCubicRationalBezierPatchSegment( const Point4f cross_section[],
										 const Point4f profile[], GLfloat u0, GLfloat u1, GLenum fill_mode,
										 bool draw_control_cage )
/*
 PURPOSE: Draw a solid-of-revolution surface using piecewise cubic rational 
 Bezier patches. Optionally displays control point frame.
 RECEIVES: cross_section: A cubic rational Bezier circle arc around the Z-axis.
 profile: The rational cubic Bezier segment rotating around the circle.
 u0: The initial knot value.
 u1: The ending knot value.
 fill_mode: One of GL_LINE or GL_FILL, indicating wireframe or filled
 patch.
 draw_control_cage: true indicates draw the cage, false indicates draw
 only the surface patch.
 REMARKS: Requires GL_MAP2_VERTEX_4 enabled.
 */
{
	CubicRational3DBezierPatchControlPoints patch_verts;
	
	initializeCubicSolidOfRevPatch( cross_section, profile, patch_verts );
	glMapGrid2f( SURFACE_TESSELLATION, u0, u1, SURFACE_TESSELLATION, 0.0f, 1.0f );
	glMap2f( GL_MAP2_VERTEX_4,
			u0, u1, HOMOG_3D_COORDS, CUBIC_ORDER,
			0.0f, 1.0f, HOMOG_3D_COORDS * CUBIC_RATIONAL_PATCH_VERTS_U, CUBIC_ORDER,
			patch_verts[ 0 ][ 0 ].getRawFloats() );
	glEvalMesh2( fill_mode, 0, SURFACE_TESSELLATION, 0, SURFACE_TESSELLATION );
	
	if ( draw_control_cage )
		drawCubicBezierControlCage( patch_verts );
}

/*-----------------------------------------------*/
void drawCubicRationalBezierPatchSOR( const Point4f profile_rcbez[], const GLfloat knot_vector[],
									 uint num_points, float max_angle, GLenum fill_mode, bool draw_control_cage )
/*
 PURPOSE: Draw a solid-of-revolution surface using piecewise cubic rational 
 Bezier patches. Optionally displays control point frame.
 RECEIVES: profile_rcbez: A piecewise cubic rational Bezier curve profile.
 num_points: The number of control points in the curve. Neighboring 
 segments share a common vertex, so there must be ( N * 3 + 1 ) 
 points for N piecewise segments in the profile curve.
 max_angle: The surface rotates around the Z-axis from zero radians 
 to this angle in radians, so 2*Pi radians specifies a complete surface.
 fill_mode: One of GL_LINE or GL_FILL, indicating wireframe or filled
 patch.
 draw_control_cage: true indicates draw the cage, false indicates draw
 only the surface patch.
 REMARKS: Requires GL_MAP2_VERTEX_4 enabled.
 */
{
	CubicRational3DBezierCurveControlPoints cross_section_1;
	CubicRational3DBezierCurveControlPoints cross_section_2;
	
	bool bIncludeSecondHalf = initializeCircleArcs( max_angle, cross_section_1, cross_section_2 );
	
	GLfloat u0 = 0.0f, u1 = 1.0f;
	for ( int knot = 0, profile_point = 0; profile_point < num_points - 1;
		 ++knot, profile_point += PROFILE_CUBIC_CONTROL_POINT_STEP )
	{
		if ( knot_vector )
		{
			u0 = knot_vector[ knot ];
			u1 = knot_vector[ knot + 1 ];
		}
		
		drawCubicRationalBezierPatchSegment( cross_section_1, profile_rcbez + profile_point,
											u0, u1, fill_mode, draw_control_cage );
		if ( bIncludeSecondHalf )
			drawCubicRationalBezierPatchSegment( cross_section_2, profile_rcbez + profile_point,
												u0, u1, fill_mode, draw_control_cage );
	}
}

/*-----------------------------------------------*/
void drawCubicRationalCatmullRomPatchSOR( const Point4f profile_points[],
										 uint num_points, float max_angle, GLenum fill_mode, bool draw_control_cage )
/*
 PURPOSE: Draw a solid-of-revolution surface using piecewise cubic rational 
 Bezier patches. Optionally displays control point frame.
 RECEIVES: profile_points: The function interpolates these points with a 
 Catmull-Rom spline, and uses this spline as the profile curve.
 num_points: The number of interpolated points in the curve.
 max_angle: The surface rotates around the Z-axis from zero radians 
 to this angle in radians, so 2*Pi radians specifies a complete surface.
 fill_mode: One of GL_LINE or GL_FILL, indicating wireframe or filled
 patch.
 draw_control_cage: true indicates draw the cage, false indicates draw
 only the surface patch.
 REMARKS: Requires GL_MAP2_VERTEX_4 enabled.
 */
{
	CubicRational3DBezierCurveControlPoints profile_bez_verts;
	
	for ( int profile_point = 1; profile_point < num_points - 2; ++profile_point )
	{
		// p_i
		profile_bez_verts[ 0 ].set( profile_points[ profile_point ] );
		
		Point4f l;
		l = ( profile_points[ profile_point + 1 ] - profile_points[ profile_point - 1 ] ) * ( 1.0f / 2 );
		
		// p^+_i
		profile_bez_verts[ 1 ].set( profile_points[ profile_point ] + l * ( 1.0f / 3 ) );
		profile_bez_verts[ 1 ].setW( 1.0f );
		
		// p^-_i+1
		l = ( profile_points[ profile_point + 2 ] - profile_points[ profile_point ] ) * ( 1.0f / 2 );
		profile_bez_verts[ 2 ].set( profile_points[ profile_point + 1 ] - l * ( 1.0f / 3 ) );
		profile_bez_verts[ 2 ].setW( 1.0f );
		
		// p_i+1
		profile_bez_verts[ 3 ].set( profile_points[ profile_point + 1 ] );
		
		GLfloat knot_vector[ 2 ] = { profile_point, profile_point + 1 };
		drawCubicRationalBezierPatchSOR( profile_bez_verts, knot_vector, 4, max_angle,
										fill_mode, draw_control_cage );
	}
}

/*-----------------------------------------------*/
static void calculateBesselOverhauserHalfKnotDerivative( const Point4f profile_points[], const GLfloat knot_vector[], 
														int i, Point4f & v_i_p_1o2 )
/*
 PURPOSE: Calculate v_i+1/2 for a Bessel-Overhauser spline.
 RECEIVES: profile_points: The points, p_i, the spline interpolates.
 knot_vector: The knot array, u_i.
 i: The function calculates the velocity at the half-knot point of this interval.
 v_i_p_1o2: The function returns the calculated v_i+1/2 in this point.
 RETURNS: The function modifies v_i_p_1o2, but does not otherwise return a value.
 */
{
	v_i_p_1o2.set( ( profile_points[ i + 1 ] - profile_points[ i ] ) * ( 1.0f / ( knot_vector[ i + 1 ] - knot_vector[ i ] ) ) );
}

/*-----------------------------------------------*/
static void calculateBesselOverhauserDerivative( const Point4f profile_points[], const GLfloat knot_vector[], 
												int i, Point4f & v_i, const Point4f & v_im1o2, const Point4f & v_ip1o2 )
/*
 PURPOSE: Calculate v_i for a Bessel-Overhauser spline.
 RECEIVES: profile_points: The points, p_i, the spline interpolates.
 knot_vector: The knot array, u_i.
 i: The function calculates the approximate velocity at this knot.
 v_i: The function returns the approimated velocity in this point.
 v_im1o2: The function uses the prior half knot velocity, v_i-1/2 in 
 the formula approximating the velocity at the knot.
 v_ip1o2: The function uses the next half knot velocity, v_i+1/2 in 
 the formula approximating the velocity at the knot.
 RETURNS: The function modifies v_i_p_1o2, but does not otherwise return a value.
 */
{
	v_i = v_im1o2 * ( knot_vector[ i + 1 ] - knot_vector[ i ] ) + 
	v_ip1o2 * ( knot_vector[ i ] - knot_vector[ i - 1 ] );
	v_i = v_i * ( 1.0f / ( knot_vector[ i + 1 ] - knot_vector[ i - 1 ] ) );
}

/*-----------------------------------------------*/
void drawCubicRationalBesselOverhauserPatchSOR( const Point4f profile_points[],
											   const GLfloat knot_vector[], uint num_points, float max_angle,
											   GLenum fill_mode, bool draw_control_cage )
/*
 PURPOSE: Draw a solid-of-revolution surface using piecewise cubic rational 
 Bezier patches. Optionally displays control point frame.
 RECEIVES: profile_points: The function interpolates these points with a 
 Catmull-Rom spline, and uses this spline as the profile curve.
 knot_vector: The knot vector for the curve, or NULL to automatically
 use the chord-length parameterization.
 num_points: The number of interpolated points in the curve.
 max_angle: The surface rotates around the Z-axis from zero radians 
 to this angle in radians, so 2*Pi radians specifies a complete surface.
 fill_mode: One of GL_LINE or GL_FILL, indicating wireframe or filled
 patch.
 draw_control_cage: true indicates draw the cage, false indicates draw
 only the surface patch.
 REMARKS: Requires GL_MAP2_VERTEX_4 enabled.
 */
{
	CubicRational3DBezierCurveControlPoints profile_bez_verts;
	if ( !knot_vector )
	{
		GLfloat * chord_length_knots = (GLfloat*)alloca( sizeof( GLfloat ) * ( num_points + 1 ) );
		chord_length_knots[ 0 ] = 0.0f;
		for ( int profile_point = 1; profile_point < num_points; ++profile_point )
			chord_length_knots[ profile_point ] = ( profile_points[ profile_point ] - profile_points[ profile_point - 1 ] ).length() + chord_length_knots[ profile_point - 1 ];
		knot_vector = chord_length_knots;
	}
	
	Point4f v_im1o2, v_ip1o2;
	Point4f v_i, v_ip1;
	calculateBesselOverhauserHalfKnotDerivative( profile_points, knot_vector, 0, v_im1o2 );
	calculateBesselOverhauserHalfKnotDerivative( profile_points, knot_vector, 1, v_ip1o2 );
	calculateBesselOverhauserDerivative( profile_points, knot_vector, 1, v_i, v_im1o2, v_ip1o2 );
	
	for ( int profile_point = 1; profile_point < num_points - 2; ++profile_point )
	{
		// p_i
		profile_bez_verts[ 0 ].set( profile_points[ profile_point ] );
		
		
		// p^+_i
		profile_bez_verts[ 1 ].set( profile_points[ profile_point ] + v_i * ( ( knot_vector[ profile_point ] - knot_vector[ profile_point - 1 ] ) / 3 ) );
		profile_bez_verts[ 1 ].setW( 1.0f );
		
		
		v_im1o2.set( v_ip1o2 );
		calculateBesselOverhauserHalfKnotDerivative( profile_points, knot_vector, profile_point + 1, v_ip1o2 );
		calculateBesselOverhauserDerivative( profile_points, knot_vector, profile_point + 1, v_ip1, v_im1o2, v_ip1o2 );
		v_i.set( v_ip1 );
		
		// p^-_i+1
		profile_bez_verts[ 2 ].set( profile_points[ profile_point + 1 ] - v_ip1 * ( ( knot_vector[ profile_point + 1 ] - knot_vector[ profile_point ] ) / 3 ) );
		profile_bez_verts[ 2 ].setW( 1.0f );
		
		// p_i+1
		profile_bez_verts[ 3 ].set( profile_points[ profile_point + 1 ] );
		
		GLfloat temp_knot[ 2 ] = { 0.0f, 5.0f };
		drawCubicRationalBezierPatchSOR( profile_bez_verts, 
										//knot_vector + profile_point,
										temp_knot,
										4, max_angle,
										fill_mode, draw_control_cage );
	}
}

/*-----------------------------------------------*/
void drawBezierPatchSOR(void)
/*
 PURPOSE: Draw a solid-of-revolution surface using piecewise cubic rational Bezier patches.
 */
{
	glColor3fv( RED_COLOR );
	
	glEnable( GL_MAP2_VERTEX_4 );
	glEnable( GL_AUTO_NORMAL );
	
	surface.draw( GL_LINE, true );
	
	
	glDisable( GL_MAP2_VERTEX_4 );
	glDisable( GL_AUTO_NORMAL );
}


/*
 GLUT AND DRIVER CODE
 */

/*-----------------------------------------------*/
void init(void)
/*
 PURPOSE: Used to initializes our OpenGL window
 */
{
	glClearColor(1.0, 1.0, 1.0, 0.0);
	
	surface.readData( FILE_NAME );
}

/*-----------------------------------------------*/
void drawFn(void)
/*
 PURPOSE: GLUT display callback function for this program.
 Draw an ampersand using piecewise cubic Bezier curves.
 */
{
	glClear(GL_COLOR_BUFFER_BIT);
	
	drawBezierPatchSOR();
	
	glFlush();
}

/*-----------------------------------------------*/
void winReshapeFn(int newWidth, int newHeight)
/*
 PURPOSE: Resizes/redisplays the contents of the current window
 RECEIVES: newWidth, newHeight -- the new dimensions (ignored)
 REMARKS: Selects an orthographic projection surrounding the ampersand.
 */
{   
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	
	Point4f min, max, center;
	surface.getExtent( min, max );
	//center.lerp( min, max, 0.5f );
	center.set( 0.0f, 0.0f, 0.0f );
	max = max - min;
	float majorAxis = max.getX();
	if ( majorAxis < max.getZ() )
		majorAxis = max.getZ();
	majorAxis *= 1.1666667f;
	min = center - Point4f( majorAxis, majorAxis, majorAxis );
	max = center + Point4f( majorAxis, majorAxis, majorAxis );
	
	glOrtho(
			min.getX(), max.getX(),
			min.getZ(), max.getZ(),
			-100.0f, 100.0f );
	
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	
	glRotatef( -60, 1.0f, 0.0f, 0.0f );
	
	glClear(GL_COLOR_BUFFER_BIT);
}

int main(int argc, char** argv)
{
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
	
	const int WINDOW_INITIAL_POSITION_X = 50;
	const int WINDOW_INITIAL_POSITION_Y = WINDOW_INITIAL_POSITION_X;
	const int WINDOW_INITIAL_SIZE_X = 500;
	const int WINDOW_INITIAL_SIZE_Y = WINDOW_INITIAL_SIZE_X;
	glutInitWindowPosition(WINDOW_INITIAL_POSITION_X, WINDOW_INITIAL_POSITION_Y);
	glutInitWindowSize(WINDOW_INITIAL_SIZE_X, WINDOW_INITIAL_SIZE_Y);
	glutCreateWindow("Recursive Subdivision of Bezier Curve");
	
	init();
	glutDisplayFunc(drawFn);
	glutReshapeFunc(winReshapeFn);
	glutMainLoop();
	return 0;
}