匀速贝塞尔曲线运动（续）

转：http://www.thecodeway.com/blog/?p=749

自从上次我写那篇讨论贝塞尔曲线的文章之后，一直有一些朋友问我关于如何实现更高次贝塞尔曲线匀速运动的问题，比如三次曲线。

    实现匀速运动的关键的一个问题就是如何积分求得曲线长度，其实在实际工程中，很多情况是无法直接求得一个函数积分后的表达式的，这种情况就要用数值积分(Numerical integration)的算法求得近似解。我用Mathmatica模拟了一下这个过程。

    另外我用C++实现了一份代码，其中使用辛普森积分法求得曲线长度，有兴趣的朋友可以参考。更高次的曲线，稍作修改即可实现。

#include <stdio.h>

#include <math.h>

#include <windows.h>



//四个控制点

POINT P0={50,50}, P1={300,600}, P2={600,600}, P3={800,200};



//曲线总长度

double total_length = 0.0;



//曲线分割的份数

const int STEP = 70;



//用于保存绘制点数据的数组

POINT pixels[STEP];





//-------------------------------------------------------------------------------------

//x坐标方程

double beze_x(double t)

{

	double it = 1-t;

	return it*it*it*P0.x + 3*it*it*t*P1.x + 3*it*t*t*P2.x + t*t*t*P3.x;

}



//-------------------------------------------------------------------------------------

//y坐标方程

double beze_y(double t)

{

	double it = 1-t;

	return it*it*it*P0.y + 3*it*it*t*P1.y + 3*it*t*t*P2.y + t*t*t*P3.y;

}



//-------------------------------------------------------------------------------------

//x坐标速度方程

double beze_speed_x(double t)

{

	double it = 1-t;

	return -3*P0.x*it*it + 3*P1.x*it*it - 6*P1.x*it*t + 6*P2.x*it*t - 3*P2.x*t*t + 3*P3.x*t*t;

}



//-------------------------------------------------------------------------------------

//y坐标速度方程

double beze_speed_y(double t)

{

	double it = 1-t;

	return -3*P0.y*it*it + 3*P1.y*it*it - 6*P1.y*it*t + 6*P2.y*it*t - 3*P2.y*t*t + 3*P3.y*t*t;

}



//-------------------------------------------------------------------------------------

//速度方程

double beze_speed(double t)

{

	double sx = beze_speed_x(t);

	double sy = beze_speed_y(t);



	return sqrt(sx*sx+sy*sy);

}



//-------------------------------------------------------------------------------------

//长度方程,使用Simpson积分算法

double beze_length(double t)

{

	//在总长度范围内，使用simpson算法的分割数

	#define TOTAL_SIMPSON_STEP	(10000)

	

	//分割份数

	int stepCounts = (int)(TOTAL_SIMPSON_STEP*t);

	if(stepCounts & 1) stepCounts++;	//偶数

	if(stepCounts==0) return 0.0;



	int halfCounts = stepCounts/2;

	double sum1=0.0, sum2=0.0;

	double dStep = t/stepCounts;



	for(int i=0; i<halfCounts; i++)

	{

		sum1 += beze_speed((2*i+1)*dStep);

	}

	for(int i=1; i<halfCounts; i++)

	{

		sum2 += beze_speed((2*i)*dStep);

	}



	return (beze_speed(0.0)+beze_speed(1.0)+2*sum2+4*sum1)*dStep/3.0;

}



//-------------------------------------------------------------------------------------

//根据t推导出匀速运动自变量t'的方程(使用牛顿切线法)

double beze_even(double t)

{

	double len = t*total_length; //如果按照匀速增长,此时对应的曲线长度

	double t1=t, t2;

	do

	{

		t2 = t1 - (beze_length(t1)-len)/beze_speed(t1);

		if(abs(t1-t2)<0.0000001) break;

		t1=t2;

	}while(true);



	return t2;

}



//-------------------------------------------------------------------------------------

LRESULT CALLBACK _WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)

{

	switch (message) 

	{

	case WM_TIMER:

		{

			static unsigned int nIndex = 0;

			if(nIndex>=0 && nIndex<=STEP)

			{

				double t = (double)nIndex/STEP;

				//求得匀速运动对应的t值

				t = beze_even(t);



				//根据贝塞尔曲线函数，求得取得此时的x,y坐标

				double x = beze_x(t);

				double y = beze_y(t);



				//取整

				pixels[nIndex].x = (int)(x+0.5);

				pixels[nIndex].y = (int)(y+0.5);



				nIndex++;

				InvalidateRect(hWnd, 0, 0);

			}

			else

			{

				KillTimer(hWnd, 101);

			}

		}

		break;

	case WM_PAINT:

		{

			PAINTSTRUCT ps;

			HDC hdc = BeginPaint(hWnd, &ps);

			::MoveToEx(hdc, P0.x, P0.y, 0);

			LineTo(hdc, P1.x, P1.y);

			LineTo(hdc, P2.x, P2.y);

			LineTo(hdc, P3.x, P3.y);



			for(int i=0; i<STEP; i++)

			{

				const POINT &pt = pixels[i];

				if(pt.x==0 && pt.y==0) break;



				::MoveToEx(hdc, pt.x-2, pt.y, 0);

				::LineTo(hdc, pt.x+3, pt.y);

				::MoveToEx(hdc, pt.x, pt.y-2, 0);

				::LineTo(hdc, pt.x, pt.y+3);

			}

			EndPaint(hWnd, &ps);

		}

		break;

	case WM_DESTROY:

		PostQuitMessage(0);

		break;

	default:

		return DefWindowProc(hWnd, message, wParam, lParam);

	}

	return 0;

}



//-------------------------------------------------------------------------------------

int APIENTRY WinMain(HINSTANCE hInstance,

                     HINSTANCE hPrevInstance,

                     LPTSTR    lpCmdLine,

                     int       nCmdShow)

{

	//计算总长度

	total_length = beze_length(1.0);



	//注册窗口类

	WNDCLASSEX wcex;

	ZeroMemory(&wcex, sizeof(WNDCLASSEX));



	wcex.cbSize = sizeof(WNDCLASSEX); 

	wcex.style			= CS_HREDRAW | CS_VREDRAW;

	wcex.lpfnWndProc	= (WNDPROC)_WndProc;

	wcex.hInstance		= hInstance;

	wcex.hCursor		= LoadCursor(NULL, IDC_ARROW);

	wcex.hbrBackground	= (HBRUSH)(COLOR_WINDOW+1);

	wcex.lpszClassName	= "BezierClass";

	RegisterClassEx(&wcex);



	//创建窗口

	HWND hWnd = CreateWindow("BezierClass", "BezierDemo", WS_OVERLAPPEDWINDOW,

      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL);

	ShowWindow(hWnd, nCmdShow);

	UpdateWindow(hWnd);





	//清空绘制点数据

	ZeroMemory(&pixels, sizeof(pixels));



	//设定定时刷新计时器

	SetTimer(hWnd, 101, 10, 0);



	//消息循环

	MSG msg;

	while(GetMessage(&msg, NULL, 0, 0)) 

	{

		TranslateMessage(&msg);

		DispatchMessage(&msg);

	}



	return (int) msg.wParam;

}