
#import <Foundation/Foundation.h>
#import <AppKit/AppKit.h>
#import <ansi/math.h>

#import "Spline.h"
#import <MiscAppKit/MiscUserPath.h>
#import <MiscAppKit/MiscHitPath.h>
#import <MiscAppKit/MiscControlPoints.h>
#import "BezierORShape.h"


static float distance(float x1, float y1, float x2, float y2);
static void tridiagonalSolve(int dim,const float *a,const float *b, const float *c, const float *r,float *x);
static float *allocUnitParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints);
static float *allocChordalParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints, float ub, float ue);
static float *allocZentripedalParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints, float ub, float ue);
static float *allocFoleyParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints);
static NSPoint *allocDPoints(NSZone *zone,const NSPoint *points,int numberOfPoints,float *delta, NSPoint b1, NSPoint b3LMinus1);
static void cubicBSplineInterpolation(NSZone *zone,const NSPoint *points,
	int numberOfPoints,NSPoint *bezierPoints,ParametrisationType pType,EndPointType eAType,EndPointType eBType);
static NSPoint quadraticFit(NSPoint x0, NSPoint x1, NSPoint x2, float d0, float d1);
static float foleyAngle(NSPoint p0, NSPoint p1, NSPoint p2);

							
@implementation Spline
	
- init
{
    [super init];
	linePath = [[MiscUserPath allocWithZone:[self zone]] init];
	splinePath = [[MiscUserPath allocWithZone:[self zone]] init];
	bezierControlPath = [[MiscUserPath allocWithZone:[self zone]] init];
	points = NSZoneMalloc([self zone],sizeof(NSPoint) * 30);
	bezierPoints = NSZoneMalloc([self zone],sizeof(NSPoint) * 88);
	maxNumberOfPoints = 30;
	numberOfPoints = 1;
	MiscPSWDefineFont("ControlPoints");
	parametrisation = UnitParametrisationType;
	endPointA = endPointB = FirstDerivateType;
    return self;
}

- (void)dealloc
{
	if(linePath)
		[linePath release];
	if(splinePath)
		[splinePath release];
	if(bezierControlPath)
		[bezierControlPath release];
	NSZoneFree([self zone],points);
	NSZoneFree([self zone],bezierPoints);
	return [super dealloc];
}

- (void)empty
{
	[linePath resetFill];
	[splinePath resetFill];
	[bezierControlPath resetFill];
	numberOfPoints = 1;
}

- (void)drawWithPoint:(NSPoint)aPoint :(int)flags
{
	int lengthOps, lengthParams, i;
	float *params, *bbox;
	
	lengthOps = [linePath lengthOps];
	lengthParams = [linePath lengthParams];
	params = [linePath params];
	bbox = [linePath bbox];
	if(numberOfPoints == lengthOps){
		params[lengthParams - 2] = aPoint.x;
		params[lengthParams - 1] = aPoint.y;
		if(aPoint.x < bbox[0])
			bbox[0] = aPoint.x;
   		if(aPoint.y < bbox[1])
			bbox[1] = aPoint.y;
   		if(aPoint.x > bbox[2])
			bbox[2] = aPoint.x;
   		if(aPoint.y > bbox[3])
			bbox[3] = aPoint.y;	
	} else {
		if(lengthOps)
			[linePath lineto:aPoint.x :aPoint.y];
		else
			[linePath moveto:aPoint.x :aPoint.y];
	}
	points[numberOfPoints - 1].x = aPoint.x;
	points[numberOfPoints - 1].y = aPoint.y;
	if(flags & DRAW_LINE){
		PSsetgray(NSLightGray);
		[linePath send:dps_ustroke cached:NO];
	}
	if(numberOfPoints >= 3){
		cubicBSplineInterpolation([self zone],points,numberOfPoints,bezierPoints,
				parametrisation,endPointA,endPointB);
		[splinePath resetFill];
		[bezierControlPath resetFill];
		[splinePath moveto:bezierPoints[0].x :bezierPoints[0].y];
		[bezierControlPath moveto:bezierPoints[0].x :bezierPoints[0].y];
		for(i = 1; i < numberOfPoints;i++){
			[splinePath curveto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y 
				:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y
				:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
			[bezierControlPath lineto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y];
			[bezierControlPath lineto:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y];
			[bezierControlPath lineto:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
		}
		if(flags & DRAW_CURVE){
			PSsetgray(NSBlack);
			[splinePath send:dps_ustroke cached:NO];
		}
		if(flags & DRAW_BEZIER){
			PSsetrgbcolor(1.0,0.0,0.0);
			[bezierControlPath send:dps_ustroke cached:NO];
		}
	}	
}

- (void)drawWithFlags:(int)flags
{
	int i,l;
	char *controlChars;
	char *controlPrimChars;
	float *controlPts;
	
	if(flags & DRAW_LINE){
		PSsetgray(NSLightGray);
		[linePath send:dps_ustroke cached:NO];
	}
	if(numberOfPoints >= 3){
		l = numberOfPoints;
		controlChars = NSZoneMalloc([self zone], l + 1);
		controlPrimChars = NSZoneMalloc([self zone], l + 1);
		controlPts = NSZoneMalloc([self zone],sizeof(float) * l * 2);
		for(i = 0; i < l - 1; i++){
			controlChars[i] = 'a';
			controlPrimChars[i] = 'e';
			controlPts[2 * i] = points[i + 1].x - points[i].x;
			controlPts[2 * i + 1] = points[i + 1].y - points[i].y;
		}
		controlPts[2 * (l - 1)] = 0;
		controlPts[2 * (l - 1) + 1] = 0;
		controlChars[l - 1] = 'a';
		controlPrimChars[l - 1] = 'e';
		controlChars[l] = '\0';
		controlPrimChars[l] = '\0';
		if(flags & DRAW_CURVE){
			PSsetgray(NSBlack);
			[splinePath send:dps_ustroke cached:NO];
		}
		PSselectfont("ControlPoints",8.0);
		MiscPSWDrawControlPoints(points[0].x,points[0].y,controlPts,2 * l,controlChars);
		PSsetgray(NSDarkGray);
		MiscPSWDrawControlPoints(points[0].x,points[0].y,controlPts,2 * l,controlPrimChars);
		NSZoneFree([self zone], controlChars);
		NSZoneFree([self zone],	controlPrimChars);
		NSZoneFree([self zone],	controlPts);
		if(flags & DRAW_BEZIER){
			PSsetrgbcolor(1.0,0.0,0.0);
			[bezierControlPath send:dps_ustroke cached:NO];
		}
	}
}
	
- (void)addPoint:(NSPoint)aPoint
{
	numberOfPoints++;
}

- (int)hitWithPoint:(NSPoint)aPoint
{
	int lengthParams;
	float *params;
	int result = -1, i = 0;
	BOOL found = NO;
	NSRect hitRect;
	NSPoint hotSpot;
	
	lengthParams = [linePath lengthParams];
	params = [linePath params];
	hitRect.size.width = hitRect.size.height = 4;
	hotSpot = aPoint;
	while(!found && (2 * i + 1 < lengthParams)){
		hitRect.origin.x = params[2 * i] - 2;
		hitRect.origin.y = params[2 * i + 1] - 2;
		if(NSMouseInRect(hotSpot , hitRect , NO)){
			found = YES;
			result = i;
		}		
		i++;
	}
	return result;
}

- (void)changePointNr:(int)nr to:(NSPoint)aPoint
{
	int lengthParams,i;
	float *params, *bbox;
	
	lengthParams = [linePath lengthParams];
	if((2 * nr + 1 >= lengthParams) || (nr < 0))
		return;
	params = [linePath params];
	bbox = [linePath bbox];
	params[2 * nr] = aPoint.x;
	params[2 * nr + 1] = aPoint.y;
	if(aPoint.x < bbox[0])
		bbox[0] = aPoint.x;
   	if(aPoint.y < bbox[1])
		bbox[1] = aPoint.y;
   	if(aPoint.x > bbox[2])
		bbox[2] = aPoint.x;
   	if(aPoint.y > bbox[3])
		bbox[3] = aPoint.y;
	if(numberOfPoints >= 3){
		points[nr].x = aPoint.x;
		points[nr].y = aPoint.y;
		cubicBSplineInterpolation([self zone],points,numberOfPoints,bezierPoints,
				parametrisation,endPointA,endPointB);
		[splinePath resetFill];
		[bezierControlPath resetFill];
		[splinePath moveto:bezierPoints[0].x :bezierPoints[0].y];
		[bezierControlPath moveto:bezierPoints[0].x :bezierPoints[0].y];
		for(i = 1; i < numberOfPoints;i++){
			[splinePath curveto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y 
				:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y
				:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
			[bezierControlPath lineto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y];
			[bezierControlPath lineto:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y];
			[bezierControlPath lineto:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
		}
	}	
}

- (ParametrisationType)parametrisationType
{
	return parametrisation;
}

- (void)setParametrisationType:(ParametrisationType)aParametrisationType
{
	int i;
	
	if(aParametrisationType != parametrisation){
		parametrisation = aParametrisationType;
		if(numberOfPoints >= 3){
			cubicBSplineInterpolation([self zone],points,numberOfPoints,bezierPoints,
				parametrisation,endPointA,endPointB);
			[splinePath resetFill];
			[bezierControlPath resetFill];
			[splinePath moveto:bezierPoints[0].x :bezierPoints[0].y];
			[bezierControlPath moveto:bezierPoints[0].x :bezierPoints[0].y];
			for(i = 1; i < numberOfPoints;i++){
				[splinePath curveto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y 
					:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y
					:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y];
				[bezierControlPath lineto:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
			}
		}
	}
}
	
- (EndPointType)endPointAType
{
	return endPointA;
}

- (void)setEndPointAType:(EndPointType)aEndPointType
{
	int i;
	
	if(aEndPointType != endPointA){
		endPointA = aEndPointType;
		if(numberOfPoints >= 3){
			cubicBSplineInterpolation([self zone],points,numberOfPoints,bezierPoints,
				parametrisation,endPointA,endPointB);
			[splinePath resetFill];
			[bezierControlPath resetFill];
			[splinePath moveto:bezierPoints[0].x :bezierPoints[0].y];
			[bezierControlPath moveto:bezierPoints[0].x :bezierPoints[0].y];
			for(i = 1; i < numberOfPoints;i++){
				[splinePath curveto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y 
					:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y
					:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y];
				[bezierControlPath lineto:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
			}
		}
	}
}


- (EndPointType)endPointBType
{
	return endPointB;
}

- (void)setEndPointBType:(EndPointType)aEndPointType
{
	int i;
	
	if(aEndPointType != endPointB){
		endPointB = aEndPointType;
		if(numberOfPoints >= 3){
			cubicBSplineInterpolation([self zone],points,numberOfPoints,bezierPoints,
				parametrisation,endPointA,endPointB);
			[splinePath resetFill];
			[bezierControlPath resetFill];
			[splinePath moveto:bezierPoints[0].x :bezierPoints[0].y];
			[bezierControlPath moveto:bezierPoints[0].x :bezierPoints[0].y];
			for(i = 1; i < numberOfPoints;i++){
				[splinePath curveto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y 
					:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y
					:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 2].x :bezierPoints[3 * i - 2].y];
				[bezierControlPath lineto:bezierPoints[3 * i - 1].x :bezierPoints[3 * i - 1].y];
				[bezierControlPath lineto:bezierPoints[3 * i].x :bezierPoints[3 * i].y];
			}
		}
	}
}


- (NSRect)splineBounds
{
	NSRect b;
	NSRect result;
	
	b = [bezierControlPath bounds];
	result.origin.x = b.origin.x;
	result.origin.y = b.origin.y;
	result.size.width = b.size.width;
	result.size.height = b.size.height;
	return result;
}

- (void)drawSpecial
{
	if(numberOfPoints >= 3){
		PSsetgray(NSBlack);
		[splinePath send:dps_ustroke cached:NO];
	}
}

- (BezierORShape *)rot3DShape
{
         return [BezierORShape shapeWithPoints:bezierPoints :(3 * numberOfPoints - 2)];
}

@end

static float distance(float x1, float y1, float x2, float y2)
{
	float dx,dy;
	
	dx = x1 - x2; dy = y1 - y2;
	return (float)sqrt((double)(dx * dx + dy * dy));
}

static void tridiagonalSolve(int dim,const float *a,const float *b, const float *c, const float *r,float *x)
{
	float aux, *auxArray;
	int i;
	
	auxArray = NSZoneMalloc(NSDefaultMallocZone(),sizeof(float) * dim);	
	aux = 1 / b[0];
	x[0] = r[0] * aux;
	auxArray[0] = 0;
	for(i = 1; i < dim;i++){
		auxArray[i] = c[i - 1] * aux;
		aux = 1 / (b[i] - a[i - 1] * auxArray[i]);
		x[i] = (r[i] - a[i - 1] * x[i - 1]) * aux;
	}
	for(i = dim - 2; i >= 0; i--)
		x[i] -= auxArray[i + 1] * x[i + 1];
	NSZoneFree(NSDefaultMallocZone(),auxArray);
}

static float *allocUnitParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints)
{
	float *resultArray;
	int i;
	
	resultArray = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	for(i = 0; i < numberOfPoints - 1;i++){
		resultArray[i] = 1;
	}
	return resultArray;
}

static float *allocChordalParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints, float ub, float ue)
{
	float *resultArray, sum;
	int i;
	
	resultArray = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	sum = 0;
	for(i = 0; i < numberOfPoints - 1;i++){
		resultArray[i] = distance(points[i].x,points[i].y,points[i+1].x,points[i+1].y);
		sum += resultArray[i];
	} 
	for(i = 0; i < numberOfPoints - 1;i++){
		if(resultArray[i] == 0)
			resultArray[i] = 1;
		resultArray[i] *= (ue - ub) / sum;
	}
	return resultArray;
}

static float *allocZentripedalParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints, float ub, float ue)
{
	float *resultArray, sum;
	int i;
	
	resultArray = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	sum = 0;
	for(i = 0; i < numberOfPoints - 1;i++){
		resultArray[i] = (float)sqrt((double)distance(points[i].x,points[i].y,points[i+1].x,points[i+1].y));
		sum += resultArray[i];
	} 
	for(i = 0; i < numberOfPoints - 1;i++){
		if(resultArray[i] == 0)
			resultArray[i] = 1;
		resultArray[i] *= (ue - ub) / sum;
	}
	return resultArray;
}

static float *allocFoleyParametrisation(NSZone *zone,const NSPoint *points,int numberOfPoints)
{
	float *resultArray, *d, *omega;
	int i;
	
	resultArray = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	d = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	omega = NSZoneMalloc(zone,sizeof(float) * (numberOfPoints - 1));
	d[0] = distance(points[0].x,points[0].y,points[1].x,points[1].y);
	for(i = 1; i < numberOfPoints - 1;i++){
		d[i] = distance(points[i].x,points[i].y,points[i+1].x,points[i+1].y); 
		omega[i] = foleyAngle(points[i - 1],points[i], points[i + 1]);	
	}
	for(i = 1; i < numberOfPoints - 2;i++)
		resultArray[i] = d[i] * (1 + 1.5 * omega[i] * d[i - 1] / (d[i - 1] + d[i]) + 
				1.5 * omega[i + 1] * d[i + 1] / (d[i] + d[i + 1]));
	resultArray[0] = resultArray[1] * d[0] / d[1];
	resultArray[numberOfPoints - 2] = resultArray[numberOfPoints - 3] * d[numberOfPoints - 2] / d[numberOfPoints - 3];
	NSZoneFree(zone,d);
	NSZoneFree(zone,omega);
	return resultArray;
}

static NSPoint *allocDPoints(NSZone *zone,const NSPoint *points,int numberOfPoints,float *delta, NSPoint b1, NSPoint b3LMinus1)
{
	float *dv,*av,*bv,*cv,*rv,sym,alpha,gamma;
	int i,L;
	NSPoint *resultArray;
	
	L = numberOfPoints - 1;
	resultArray = NSZoneMalloc(zone,sizeof(NSPoint) * numberOfPoints);
	resultArray[0] = b1;
	resultArray[L] = b3LMinus1;
	if(numberOfPoints > 3){
		av = NSZoneMalloc(zone,sizeof(float) * (L - 2));
		cv = NSZoneMalloc(zone,sizeof(float) * (L - 2));
	} else {
		av = NULL;
		cv = NULL;
	}
	bv = NSZoneMalloc(zone,sizeof(float) * (L - 1));
	dv = NSZoneMalloc(zone,sizeof(float) * (L - 1));
	rv = NSZoneMalloc(zone,sizeof(float) * (L - 1));
	bv[0] = delta[1] * delta[0] / (delta[0] + delta[1]) + 
			   delta[0] * (delta[1] + delta[2]) / (delta[0] + delta[1] + delta[2]);
	bv[0] /= delta[0] *  delta[0] * delta[1] *  delta[1];
	if(numberOfPoints > 3){
		cv[0] = delta[0] * delta[0] / (delta[0] + delta[1] + delta[2]);
		cv[0] /= delta[0] *  delta[0] * delta[1] *  delta[1];
		for(i = 2; i < L - 1; i++){
			av[i - 2] = delta[i] * delta[i] / (delta[i - 2] + delta[i - 1] + delta[i]);
			bv[i - 1] = delta[i] * (delta[i - 2] + delta[i - 1]) / (delta[i - 2] + delta[i - 1] + delta[i]) + 
				   delta[i - 1] * (delta[i] + delta[i + 1]) / (delta[i - 1] + delta[i] + delta[i + 1]);
			cv[i - 1] = delta[i - 1] * delta[i - 1] / (delta[i - 1] + delta[i] + delta[i + 1]);
			sym = delta[i - 1] *  delta[i - 1] * delta[i] *  delta[i];
			av[i - 2] /= sym;
			bv[i - 1] /= sym;
			cv[i - 1] /= sym;
		}
		av[L - 3] = delta[L - 1] * delta[L - 1] / (delta[L - 3] + delta[L - 2] + delta[L - 1]);
		av[L - 3] /= delta[L - 2] *  delta[L - 2] * delta[L - 1] *  delta[L - 1];
	}
	bv[L - 2] = delta[L - 1] * (delta[L - 3] + delta[L - 2]) / (delta[L - 3] + delta[L - 2] + delta[L - 1]) + 
			   		delta[L - 2] * delta[L - 1] / (delta[L - 2] + delta[L - 1]);
	bv[L - 2] /= delta[L - 2] *  delta[L - 2] * delta[L - 1] *  delta[L - 1];				
	alpha = delta[1] * delta[1] / (delta[0] + delta[1]) / (delta[0] *  delta[0] * delta[1] *  delta[1]);
	gamma = delta[L - 2] * delta[L - 2] / (delta[L - 2] + delta[L - 1]) 
			/ (delta[L - 2] *  delta[L - 2] * delta[L - 1] *  delta[L - 1]);
	for(i = 1; i < L;i++){	
		rv[i - 1] = (delta[i - 1] + delta[i]) * points[i].x;
		sym = delta[i - 1] *  delta[i - 1] * delta[i] *  delta[i];
		rv[i - 1] /= sym;
	}
	rv[0] -= alpha * b1.x;
	rv[L - 2] -= gamma * b3LMinus1.x;
	if(numberOfPoints > 3)
		tridiagonalSolve(L - 1,av,bv,cv,rv,dv);
	else
		tridiagonalSolve(L - 1,&alpha,bv,&gamma,rv,dv);
	for(i = 1; i < L;i++)
		resultArray[i].x = dv[i - 1];
	for(i = 1; i < L;i++){	
		rv[i - 1] = (delta[i - 1] + delta[i]) * points[i].y;
		sym = delta[i - 1] *  delta[i - 1] * delta[i] *  delta[i];
		rv[i - 1] /= sym;
	}
	rv[0] -= alpha * b1.y;
	rv[L - 2] -= gamma * b3LMinus1.y;
	if(numberOfPoints > 3)
		tridiagonalSolve(L - 1,av,bv,cv,rv,dv);
	else
		tridiagonalSolve(L - 1,&alpha,bv,&gamma,rv,dv);
	for(i = 1; i < L;i++)
		resultArray[i].y = dv[i - 1];
	if(numberOfPoints > 3){
		NSZoneFree(zone,av);
		NSZoneFree(zone,cv);
	}
	NSZoneFree(zone,bv);
	NSZoneFree(zone,dv);
	NSZoneFree(zone,rv);
	return resultArray;
}

static void cubicBSplineInterpolation(NSZone *zone,const NSPoint *points,
	int numberOfPoints,NSPoint *bezierPoints,ParametrisationType pType,EndPointType eAType,EndPointType eBType)
{
	int L,i;
	float *delta = NULL;
	NSPoint *dPoints,b1,b3LMinus1;
	
	L = numberOfPoints - 1;
	switch(pType){
		case UnitParametrisationType:
			delta = allocUnitParametrisation(zone,points,numberOfPoints);
			break;
		case ChordalParametrisationType:
			delta = allocChordalParametrisation(zone,points,numberOfPoints,0,numberOfPoints);
			break;
		case ZentripedalParametrisationType:
			delta = allocZentripedalParametrisation(zone,points,numberOfPoints,0,numberOfPoints);
			break;
		case FoleyParametrisationType:	
			delta = allocFoleyParametrisation(zone,points,numberOfPoints);
			break;
	};
	switch(eAType){
		case SecondDerivateType:
			b1 = points[0];
			break;
		case FirstDerivateType:
			b1 = points[0];
			break;
		case QuadraticType:
			b1 = quadraticFit(points[0],points[1],points[2],delta[0],delta[1]);
			break;
	};
	switch(eBType){
		case SecondDerivateType:
			b3LMinus1 = points[L];
			break;
		case FirstDerivateType:
			b3LMinus1 = points[L];
			break;
		case QuadraticType:
			b3LMinus1 = quadraticFit(points[L],points[L - 1],points[L - 2],delta[L - 1],delta[L - 2]);
			break;
	};	
	dPoints = allocDPoints(zone,points,numberOfPoints,delta,b1, b3LMinus1);
	for(i = 1; i < L - 1; i++){
		bezierPoints[3 * i] = points[i]; 
		bezierPoints[3 * i + 1].x = ((delta[i] + delta[i + 1]) * dPoints[i].x + delta[i - 1] * dPoints[i + 1].x) / 
				(delta[i - 1] + delta[i] + delta[i + 1]);
		bezierPoints[3 * i + 1].y = ((delta[i] + delta[i + 1]) * dPoints[i].y + delta[i - 1] * dPoints[i + 1].y) / 
				(delta[i - 1] + delta[i] + delta[i + 1]);		 
		bezierPoints[3 * i + 2].x = (delta[i + 1] * dPoints[i].x + (delta[i - 1] + delta[i]) * dPoints[i + 1].x) /
				(delta[i - 1] + delta[i] + delta[i + 1]);
		bezierPoints[3 * i + 2].y = (delta[i + 1] * dPoints[i].y + (delta[i - 1] + delta[i]) * dPoints[i + 1].y) /
				(delta[i - 1] + delta[i] + delta[i + 1]);		
	}
	bezierPoints[0] = points[0];
	bezierPoints[1] = b1;
	bezierPoints[2].x = (delta[1] * dPoints[0].x + delta[0] * dPoints[1].x)/(delta[0] + delta[1]);
	bezierPoints[2].y = (delta[1] * dPoints[0].y + delta[0] * dPoints[1].y)/(delta[0] + delta[1]);
	bezierPoints[3 * (L - 1)] = points[L - 1];
	bezierPoints[3 * (L - 1) + 1].x = (delta[L - 1] * dPoints[L - 1].x + delta[L -2] * dPoints[L].x) / 
				(delta[L - 2] + delta[L - 1]);
	bezierPoints[3 * (L - 1) + 1].y = (delta[L - 1] * dPoints[L - 1].y + delta[L -2] * dPoints[L].y) / 
				(delta[L - 2] + delta[L - 1]);			
	bezierPoints[3 * (L - 1) + 2] = b3LMinus1;
	bezierPoints[3 * L]	= points[L];
	NSZoneFree(zone,delta);
	NSZoneFree(zone,dPoints);
}

static NSPoint quadraticFit(NSPoint x0, NSPoint x1, NSPoint x2, float d0, float d1)
{
	NSPoint result;
	float c0,c1,c2;
	
	c0 = d1 / (2 * (d0 + d1));
	c1 = (d0 + d1) / (2 * d1);
	c2 = - d0 * d0 / (2 * d1 * (d0 + d1));
	result.x = c0 * x0.x + c1 * x1.x + c2 * x2.x;
	result.y = c0 * x0.y + c1 * x1.y + c2 * x2.y;	
	return result;
}
	
static float foleyAngle(NSPoint p0, NSPoint p1, NSPoint p2)
{
	double a2,b2,c2;
	float dx,dy,result,pi;
	double ausdruck;
	
	pi = 3.14159;
	dx = p1.x - p0.x; dy = p1.y - p0.y;
	a2 = (double)(dx * dx + dy * dy);
	dx = p1.x - p2.x; dy = p1.y - p2.y;
	b2 = (double)(dx * dx + dy * dy);
	dx = p0.x - p2.x; dy = p0.y - p2.y;
	c2 = (double)(dx * dx + dy * dy);	
	ausdruck = (a2 + b2 - c2) / (2 * sqrt(a2 * b2));
	result = (float)acos(ausdruck);
	return MIN(pi - result, pi / 2);
}
	
