#include <stdio.h>
#include <string.h>
#include <math.h>
#include <dist.h>

Point2
beval2(Point2 b[4], double t)
{
	double s, c[4];
	Point2 p;

	s = 1.0-t;

	c[0] = s*s*s;
	c[1] = 3*s*s*t;
	c[2] = 3*s*t*t;
	c[3] = t*t*t;

	p.x = c[0]*b[0].x + c[1]*b[1].x + c[2]*b[2].x
		+ c[3]*b[3].x;

	p.y = c[0]*b[0].y + c[1]*b[1].y + c[2]*b[2].y
		+ c[3]*b[3].y;

	return p;
}

double
bldev2(Point2 b[4])
{
	double d1, d2, r;
	Point2 b2[4];

	/* determine deviation from a straight line */

	b2[0] = point2(0,0);
	b2[1] = psub2(b[1], b[0]);
	b2[2] = psub2(b[2], b[0]);
	b2[3] = psub2(b[3], b[0]);

	/* special case */
	if (b2[3].x == 0.0 && b2[3].y == 0.0) {
		d1 = sqrt(b2[1].x*b2[1].x + b2[1].y*b2[1].y);
		d2 = sqrt(b2[2].x*b2[2].x + b2[2].y*b2[2].y);

		if (d1 < d2)
			r = d2;
		else
			r = d1;

		return r;
	}
	
	d1 = b2[1].x*b2[3].y - b2[1].y*b2[3].x;
	d2 = b2[2].x*b2[3].y - b2[2].y*b2[3].x;

	if (d1 < 0)
		d1 = -d1;
	if (d2 < 0)
		d2 = -d2;

	if (d2 < d1)
		r = d1;
	else
		r = d2;

	r /= sqrt(b2[3].x*b2[3].x + b2[3].y*b2[3].y);

	return r;
}

void
bsubdiv2(Point2 b[4], Point2 b0[4], Point2 b1[4])
{
	Point2 t;

	/* page 221 of bartels, beatty, and barsky */

	b0[0] = b[0];
	b0[1].x = 0.5*(b[0].x+b[1].x);
	b0[1].y = 0.5*(b[0].y+b[1].y);
	t.x = 0.5*(b[1].x+b[2].x);
	t.y = 0.5*(b[1].y+b[2].y);
	b0[2].x = 0.5*(b0[1].x+t.x);
	b0[2].y = 0.5*(b0[1].y+t.y);
	b1[3] = b[3];
	b1[2].x = 0.5*(b[2].x+b[3].x);
	b1[2].y = 0.5*(b[2].y+b[3].y);
	b1[1].x = 0.5*(t.x+b1[2].x);
	b1[1].y = 0.5*(t.y+b1[2].y);
	b0[3].x = b1[0].x = 0.5*(b0[2].x+b1[1].x);
	b0[3].y = b1[0].y = 0.5*(b0[2].y+b1[1].y);
}