【转】三次贝塞尔曲线绘制算法
原文：http://www.cnblogs.com/flash3d/archive/2012/01/30/2332176.html
源码：https://files.cnblogs.com/flash3d/bezier.rar
====================================================
这学期学图形学，就把自己的一些粗浅的理解发上去让大家拍砖。前些天做三次贝塞尔曲线绘制的上机练习，正好将从直线扫描算法中启发得来的n次多项式批量计算用上了，自认为优化得还可以。原来写的版本是C++，为了便于网上观看特改写成AS3，对这方面内行的朋友，希望还多多指点！
在讲三次贝塞尔曲线之前，先给出n次贝塞尔曲线的一般式：
[R]( t ) = ( +..i=0,n ) ( [R]i * Bi,n( t ) ) , t属于[0,1]
Bi,n( t ) = Cni * pow( 1 - t , n - i ) * pow( t , i )
这里我把符号说明下，中括号[]括起来的表示矩阵。比如[R]就是个矩阵，如果我们讨论的贝塞尔曲线的点在平面内，那么矩阵里面有两个元素，比如我们熟知的x,y，如果是空间的曲线，就是三个元素，元素再多。。那就不知道是几次元了，反正是个数学概念。在计算的时候，将矩阵的元素单独拎出来计算也是可以的，因为公式中并没有出现矩阵和矩阵相乘，也就是说，各元素计算时其实各不相干。
( +..i=0,n ) 这个运算符，这个符号是我自己发明的，其实是累加符号，表示公式后面括号内将i分别用0,1,2,3..n代替，每个代替运算的结果全部加起来。
Cni是二项式系数，其值为 n! / ( i! * ( n - i )! )
pow是求指数函数。
其中 Bi,n( t ) = Cni * pow( 1 - t , n - i ) * pow( t , i ) 为n次二项式展开后的一项。其意义在后面的三次贝塞尔曲线里面讨论更容易理解
那么，现在先讨论3次贝塞尔曲线，我们把公式中的n全部替换成3，得到三次贝塞尔曲线的一般式如下：
[R]( t ) = ( +..i=0,3 ) ( [R]i * Bi,3( t ) ) , t属于[0,1]
Bi,3( t ) = C3i * pow( 1 - t , 3 - i ) * pow( t , i )
直接把Bi,3( t )代入吧，得到式子：
[R]( t ) = ( +..i=0,3 ) ( [R]i * C3i * pow( 1 - t , 3 - i ) * pow( t , i ) ) , t属于[0,1]
对于不熟悉线性代数或者不习惯用矩阵表示数据的朋友，我们就直接把矩阵拆开好了，[R]矩阵中总共就两个元素，分别是x和y，这里就把x和y分别计算，得到关于x和y的参数曲线：
X( t ) = ( +..i=0,3 ) ( Xi * C3i * pow( 1 - t , 3 - i ) * pow( t , i ) ) , t属于[0,1]
Y( t ) = ( +..i=0,3 ) ( Yi * C3i * pow( 1 - t , 3 - i ) * pow( t , i ) ) , t属于[0,1]
如果对参数曲线不熟悉，可以参考直线的参数方程
比如一条直线经过(x1,y1)，(x2,y2)，那么其参数方程可表示为：
X(t)=x1+t*(x2-x1)
Y(t)=y1+t*(y2-y1)
如果是线段那么t在0和1之间，如果是直线，t是实数。
言归正传，下面我们开始处理参数曲线，由于x的参数方程和y的形式差不多，我们就以x的参数方程做例子来分析
X( t ) = ( +..i=0,3 ) ( Xi * C3i * pow( 1 - t , 3 - i ) * pow( t , i ) ) , t属于[0,1]
首先我们把连加拆开，那么i分别等于0,1,2,3，得到一个式子：
X( t ) =
X0 * C30 * pow( 1 - t , 3 - 0 ) * pow( t , 0 )
+
X1 * C31 * pow( 1 - t , 3 - 1 ) * pow( t , 1 )
+
X2 * C32 * pow( 1 - t , 3 - 2 ) * pow( t , 2 )
+
X3 * C33 * pow( 1 - t , 3 - 3 ) * pow( t , 3 )
t属于[0,1]
看这个式子，我们可以发现，C30
* pow( 1 - t , 3 - 0 ) * pow( t , 0 )+C31 * pow( 1 - t , 3 - 1 ) * pow(
t , 1 )+C32 * pow( 1 - t , 3 - 2 ) * pow( t , 2 )+C33 * pow( 1 - t , 3 -
3 ) * pow( t , 3 )=pow( 1 - t + t , 3)=1
那就是说，x0，x1，x2，x3相加前的系数和为1，我们称这些系数为各点的权重。当然，点的权重是随着t的值连续变化的，这也是为什么三次贝塞尔曲线能光滑地通过四个点控制的原因。
下面我们进一步对式子进行化简。
其中C30等于1，C31等于3，C32等于3，C33等于1，把能合并的都合并，能计算值的都计算出来，将式子化简：
X( t ) =
X0 * pow( 1 - t , 3 )
+
X1 * 3 * pow( 1 - t , 2 ) * t
+
X2 * 3 * ( 1 - t ) * pow( t , 2 )
+
X3 * pow( t , 3 )
t属于[0,1]
这样的式子应该很好求了，他是一个最高次为3次的多项式，在《Alchemy的使用和多项式批量计算的优化》中提到，a*pow(x,n)+b*pow(x,n-1)+..这样形式的式子，都能通过多项式批量计算方法来快速求得结果。
对于最高三次的式子来说，我们需要前三个原始数据，以及式子求导三次的常数，算法细节就不讲了，细节参考上面那篇博文。
package
{
import flash.display.Bitmap;
import flash.display.BitmapData;
import flash.display.Shape;
import flash.display.Sprite;
import flash.events.Event;
import flash.events.MouseEvent;
import flash.geom.Point;
import flash.geom.Rectangle;
/**
*贝塞尔曲线绘制算法
* @author boycy815
*/
public class Main extends Sprite
{
/*C30到C33的值总共才4个直接枚举*/
private const C30:int = 1;
private const C31:int = 3;
private const C32:int = 3;
private const C33:int = 1;
/*定义曲线上点的个数（包括两端）*/
private const LENGTH:int = 10000;
/*btC30到btC33分别是R0 R1 R2 R3四个点当t取前三个值时的权重值*/
private var btC30:Vector.<Number> = new Vector.<Number>(3);
private var btC31:Vector.<Number> = new Vector.<Number>(3);
private var btC32:Vector.<Number> = new Vector.<Number>(3);
private var btC33:Vector.<Number> = new Vector.<Number>(3);
/*画点的位图*/
private var data:BitmapData;
private var map:Bitmap;
/*用于连接四个点的四条直线*/
private var lin:Shape;
/*分别是四个点的显示对象*/
private var ct0:Sprite;
private var ct1:Sprite;
private var ct2:Sprite;
private var ct3:Sprite;
public function Main():void
{
if (stage) init();
else addEventListener(Event.ADDED_TO_STAGE, init);
}
private function init(e:Event = null):void
{
removeEventListener(Event.ADDED_TO_STAGE, init);
/*初始化位图和一些显示对象*/
data = new BitmapData(800, 600,false);
map = new Bitmap(data);
this.addChild(map);
lin = new Shape();
this.addChild(lin);
ct0 = new Sprite();
ct0.buttonMode = true;
ct0.graphics.lineStyle(0, 0x00000);
ct0.graphics.beginFill(0xffffff);
ct0.graphics.drawCircle(0, 0, 5);
ct0.graphics.endFill();
ct0.addEventListener(MouseEvent.MOUSE_DOWN, mouseDown);
ct0.addEventListener(MouseEvent.MOUSE_UP, mouseUp);
this.addChild(ct0);
ct1 = new Sprite();
ct1.buttonMode = true;
ct1.graphics.lineStyle(0, 0x00000);
ct1.graphics.beginFill(0xffffff);
ct1.graphics.drawCircle(0, 0, 5);
ct1.graphics.endFill();
ct1.addEventListener(MouseEvent.MOUSE_DOWN, mouseDown);
ct1.addEventListener(MouseEvent.MOUSE_UP, mouseUp);
this.addChild(ct1);
ct2 = new Sprite();
ct2.buttonMode = true;
ct2.graphics.lineStyle(0, 0x00000);
ct2.graphics.beginFill(0xffffff);
ct2.graphics.drawCircle(0, 0, 5);
ct2.graphics.endFill();
ct2.addEventListener(MouseEvent.MOUSE_DOWN, mouseDown);
ct2.addEventListener(MouseEvent.MOUSE_UP, mouseUp);
this.addChild(ct2);
ct3 = new Sprite();
ct3.buttonMode = true;
ct3.graphics.lineStyle(0, 0x00000);
ct3.graphics.beginFill(0xffffff);
ct3.graphics.drawCircle(0, 0, 5);
ct3.graphics.endFill();
ct3.addEventListener(MouseEvent.MOUSE_DOWN, mouseDown);
ct3.addEventListener(MouseEvent.MOUSE_UP, mouseUp);
this.addChild(ct3);
/*初始化计算的一些准备数据*/
btInit();
/*默认进行第一次定位和绘制*/
ct0.x = 100;
ct0.y = 100;
ct1.x = 700;
ct1.y = 100;
ct2.x = 700;
ct2.y = 500;
ct3.x = 100;
ct3.y = 500;
lin.graphics.lineStyle(0, 0x000000);
lin.graphics.moveTo(ct0.x, ct0.y);
lin.graphics.lineTo(ct1.x, ct1.y);
lin.graphics.lineTo(ct2.x, ct2.y);
lin.graphics.lineTo(ct3.x, ct3.y);
drawBezier3(new Point(ct0.x, ct0.y), new Point(ct1.x, ct1.y), new Point(ct2.x, ct2.y), new Point(ct3.x, ct3.y));
}
/*按钮按下便开始可以拖动点*/
private function mouseDown(e:MouseEvent):void
{
e.target.startDrag();
stage.addEventListener(MouseEvent.MOUSE_MOVE, mouseMove);
}
/*鼠标移动便不断绘制曲线*/
private function mouseMove(e:MouseEvent):void
{
/*绘制连接四个点的直线*/
lin.graphics.clear();
lin.graphics.lineStyle(0, 0x000000);
lin.graphics.moveTo(ct0.x, ct0.y);
lin.graphics.lineTo(ct1.x, ct1.y);
lin.graphics.lineTo(ct2.x, ct2.y);
lin.graphics.lineTo(ct3.x, ct3.y);
/*绘制曲线*/
data.fillRect(new Rectangle(0, 0, 800, 600), 0xffffffff);
drawBezier3(new Point(ct0.x, ct0.y), new Point(ct1.x, ct1.y), new Point(ct2.x, ct2.y), new Point(ct3.x, ct3.y));
}
/*释放按钮停止拖动*/
private function mouseUp(e:MouseEvent):void
{
var n:int;
e.target.stopDrag();
stage.removeEventListener(MouseEvent.MOUSE_MOVE, mouseMove);
}
/***********算法部分**********/
/*初始化计算所需数据，主要是btC30到btC33*/
private function btInit():void
{
var add:Number;
/*由于t的值域为[0,1]，要绘制LENGTH个点，则必须计算出每次t的增量*/
add = 1.0 / (LENGTH - 1.0);
/*三次白赛尔曲线Ri的权重公式为：C3i * pow( t - 1 , 3 - i ) * pow( t , i )
*前三次t的值分别为0，add，add*2
*则代入公式就能求出R0 R1 R2 R3的前三次的权重值
* */
btC30[0] = 1;
btC30[1] = C30 * Math.pow(1 - add, 3);
btC30[2] = C30 * Math.pow(1 - add * 2, 3);
btC31[0] = 0;
btC31[1] = C31 * Math.pow(1 - add, 2) * add;
btC31[2] = C31 * Math.pow(1 - add * 2, 2) * add * 2;
btC32[0] = 0;
btC32[1] = C32 * (Math.pow(add, 2) - Math.pow(add, 3));
btC32[2] = C32 * (Math.pow(add * 2, 2) - Math.pow(add * 2, 3));
btC33[0] = 0;
btC33[1] = C33 * Math.pow(add, 3);
btC33[2] = C33 * Math.pow(add * 2, 3);
}
/*利用准备好的初始数据，通过递推方法求出四个点对应的贝塞尔曲线上每个点的坐标，并绘制
*计算中用到多项式批量计算优化策略，上一篇博文中有提及具体的方法和思想
* */
private function drawBezier3(p0:Point,p1:Point,p2:Point,p3:Point):void
{
/*qun_x和qun_x本别为x坐标和y坐标的临时数据数组*/
var qun_x:Vector.<Number> = new Vector.<Number>(3);
var qun_y:Vector.<Number> = new Vector.<Number>(3);
/*dd_x和dd_y分别为x坐标的多项式和y坐标式求导三次后得到的常数*/
var dd_x:Number;
var dd_y:Number;
var n:int;
/*计算前三次种子数据*/
dd_x = 6 * (3 * p1.x - 3 * p2.x - p0.x + p3.x) / Math.pow(LENGTH - 1.0, 3);
qun_x[0] = p0.x * btC30[0] + p1.x * btC31[0] + p2.x * btC32[0] + p3.x * btC33[0];
dd_y = 6 * (3 * p1.y - 3 * p2.y - p0.y + p3.y) / Math.pow(LENGTH - 1.0, 3);
qun_y[0] = p0.y * btC30[0] + p1.y * btC31[0] + p2.y * btC32[0] + p3.y * btC33[0];
data.setPixel(qun_x[0], qun_y[0], 0x000000);
qun_x[1] = p0.x * btC30[1] + p1.x * btC31[1] + p2.x * btC32[1] + p3.x * btC33[1];
qun_x[0] = qun_x[1] - qun_x[0];
qun_y[1] = p0.y * btC30[1] + p1.y * btC31[1] + p2.y * btC32[1] + p3.y * btC33[1];
qun_y[0] = qun_y[1] - qun_y[0];
data.setPixel(qun_x[1], qun_y[1], 0x000000);
qun_x[2] = p0.x * btC30[2] + p1.x * btC31[2] + p2.x * btC32[2] + p3.x * btC33[2];
qun_x[1] = qun_x[2] - qun_x[1];
qun_x[0] = qun_x[1] - qun_x[0];
qun_y[2] = p0.y * btC30[2] + p1.y * btC31[2] + p2.y * btC32[2] + p3.y * btC33[2];
qun_y[1] = qun_y[2] - qun_y[1];
qun_y[0] = qun_y[1] - qun_y[0];
data.setPixel(qun_x[2], qun_y[2], 0x000000);
/*进行递推计算的循环*/
for (n = 3; n < LENGTH; n++)
{
qun_x[0] += dd_x;
qun_x[1] += qun_x[0];
qun_x[2] += qun_x[1];
qun_y[0] += dd_y;
qun_y[1] += qun_y[0];
qun_y[2] += qun_y[1];
data.setPixel(qun_x[2], qun_y[2], 0x000000);
}
}
}
}
