Unraveling Bezier Splines 

本文深入浅出地介绍了Bezier曲线的基础数学原理及其编码实现方法,并探讨了如何利用Bezier曲线创建曲面。

Unraveling Beziér Splines
by Justin Reynen

原文网址:http://www.gamedev.net/reference/articles/article888.asp

简介

毫无疑问,您已经听说过这些用于处理曲线和曲面的新技术,二次Bezier曲线,三次Bezier曲线,非均匀有理B样条(NURBS),并且非常迫切的想要把所有这些东西搞清楚,甚至可能尝试将这些技术中的一种加入到您正在制作的游戏或演示中。

这篇文章是帮您揭开这些谜团的最好的开始。

在本文中,我会向你解释Bezier曲线的数学原理,编写程序画出一条二次Bezier样条(最简单的样条),然后你就可以自己实现任意你所需要次数的Bezier曲线了。

Bezier曲线的数学原理

好了,我们将要作的第一件事是温习一下Bernsteins基本方程。在下面我会给出一个公式,别着急,你看,她很简单!

1 = t + (1 - t)

T为任意值等式都成立,但是我们所关注的T值是在01之间。

好,现在我们已经了解Bernsteins基本方程了,下面将要开始生成曲线的工作。我们选择二次Bezier曲线作为切入点,因为它是最容易实现的。。

好了,现在已经确定我们想要一条二次Bezier曲线,我们需要从万能的基本方程派生出我们需要的方程,现在就开始吧!

对于二次Bezier曲线,我们需要将Bernsteins基本方程的两段同时取平方;对于三次Bezier曲线,我们就去立方,以此类推。下面是两边取平方后的Bernsteins基本方程

1^2 = (t + (1 - t))^2

 

1 = t^2 + 2*t*(1-t) + (1-t)*(1-t)

可能你注意到有很多项我没有展开,因为这样更有助于理解和编码。

好的,现在我们获得了从基本方程衍生出来的3个函数。你可能会问他们在哪儿,就在你的面前!我们的函数就是等式右侧的那3个多项式。我们将这3个函数取名为Bx(t)

#define B1(t)          (t*t)

#define B2(t)          (2*t*(1-t))

#define B3(t)          ((1-t)*(1-t))

你应该注意到,对于二次Bezier曲线有三个函数,对于三次Bezier曲线有四个参数,以此类推。在本文的演示程序中我已经计算了三次和四次曲线的方程,请查看"Functions.txt"

好了,现在我们已经搞清楚Bezier曲线的简单数学意义,接下来让我们谈一下编码吧。

编码实现Bezier曲线

好的,我们已经准备好了基本函数,现在是充分利用他们的时候了。

我认为现在你应该意识到与这3个函数相对应,有3个控制点,我们将这些点起名叫Cx。你可以向这样来存储这些点:

typedef struct sCPoint

{

        int x;

        int y;

}  C_POINT;

注:我们将要生成的是2DBezier曲线,如果你想要生成3DBezier曲线的话可以在结构中加入一个z元素。

C_POINT controlP[3];

好了,现在我们将要接触在屏幕上绘制Bezier曲线的真正代码了。目前,我们的曲线只存在于想象之中,而且我们只有数学方程。我们需要做的就是利用数学方程,沿着曲线将点画在屏幕上。

这就给我们带来了最后一个问题“采样率”。采样率就是我们在曲线上画点时,每绘制一个点需要沿着曲线移动的距离。采样率是个小数,你可以简单的把它理解为百分比。例如采样率0.5的含义就是:从曲线的一端开始,在这里绘制一个点;然后沿曲线移动50%,然后绘制一个点,然后再移动50%(到达曲线的末端),绘制最后一个点。这样就会在曲线的左端、中间和右端共绘制了3个点。我们需要一个计数器来记录已经绘制的曲线的百分比,然后我们每次对计数器累加,直到计数器到达1为止。定制的采样率越小,我们绘制的点就越多,曲线看起来就越精细。

//Note: it is assumed that each controlP has already been filled out with

//an x and y coordinate.

 

double count;  //used as our counter

double detailBias; //how many points should we put on our curve.

 

double x,y; //used as accumulators to make our code easier to read

 

 

detailBias = 1 / 50; //we'll put 51 points on out curve (0.02 detail bias)

 

do

{

x = controlP[0].x*B1(count) + controlP[1].x*B2(count) + controlP.x[2]*B2(count);

y = controlP[0]*B1.y(count) + controlP[1].y*B2(count) + controlP.y[2]*B2(count);

 

        PutPixel(x,y,RGB(255,255,255));

 

        count += detailBias;

 

}while( count <= 1);

请注意我们是如何获得xy的值。我们用控制点与相应的基本方程相乘,然后将他们加在一起。我们对xy坐标进行了这样的操作,如果我们是工作在3维空间,那么我们也要对z坐标进行一样的运算。

Bezier曲线创建曲面

我现在就要揭开如何绘制Bezier曲面的神秘面纱。这是个很大的课题,完全可以把它单独拿出来写一篇文章,而且它也超出了本文的研究范畴。我将只介绍二次Bezier曲面,并且没有提供演示代码。我将会介绍如何生成曲面上的点,你可以很容易的将其改造成3次或者4次曲面。

上面这幅图是最基本的差值曲面,它由2Bezier曲线组成。它有9个控制点,你可以想象这幅图的x轴沿789方向,y轴沿741方向。为了让事情简单,我把采样率定为0.25,这样我们在每条线上都会产生5个点。在整个面上将会生成25个点。我们可以简单的把这些点以行和列分组,这样我们的数据结构就会像这样:

C_POINT points[5][5];

我们要做的工作就是对每一条曲线进行插值,产生新的控制点,然后再由新生成的控制点继续插值,这样就生成了我们需要的曲面。

In Conclusion

我希望这篇文章可以帮助你了解Bezier曲线的用法并且能够将他们应用到你的游戏或者演示当中。在本文结束之前,提出一些关于Bezier曲线的优点和缺点:

Pro's:

在层次细节模型当中,曲面是一个非常好的用于减少绘制多边形数量的工具。你可以很容易的将这些面减少到只有4个多边形(对于2次方程),而且你可以很容易的做到这一点(通过改变采样率)。它也可以很容易的使表面变得平滑。另外,如果你计划添加摄像机移动功能的话,Bezier曲线对于生成平滑的移动路线,动画移动来说是至关重要的。

Con's:

曲面对于大多数引擎来说最主要的缺点就是要花费太长的时间将曲面变成多边形,你在载入时做这些事情也不会有太大改变。但是如果你每一帧都进行计算的话,你可能会注意到它存在巨大的性能问题。这就是大多数引擎开发者对曲面发牢骚的最主要原因。

(happykevins)

 


``` function [T21, err90, err95, err99] = causality_est(xx1, xx2, np) % % function [T21, err90, err95, err99] = causality_est(x1, x2, np) % % Estimate T21, the information transfer from series X2 to series X1 % dt is taken to be 1. % % On input: % X1, X2: the series (n by 1 colum vectors) % np: integer >=1, time advance in performing Euler forward % differencing, e.g., 1, 2. Unless the series are generated % with a highly chaotic deterministic system, np=1 should be % used. % % On output: % T21: info flow from X2 to X1 (Note: Not X1 -> X2!) % err90: standard error at 90% significance level % err95: standard error at 95% significance level % err99: standard error at 99% significance level % % Citations: % X. San Liang, 2014: Unraveling the cause-effect relation between time series. Phys. Rev. E 90, 052150. % X. San Liang, 2015: Normalizing the causality between time series. Phys. Rev. E 92, 022126. dt = 1; [nm, one] = size(xx1); dx1(:,1) = (xx1(1+np:nm, 1) - xx1(1:nm-np, 1)) / (np*dt); x1(:,1) = xx1(1:nm-np, 1); dx2(:,1) = (xx2(1+np:nm, 1) - xx2(1:nm-np, 1)) / (np*dt); x2(:,1) = xx2(1:nm-np, 1); clear xx1 xx2; N = nm-np; C = cov(x1, x2); dC(1,1) = sum((x1 - mean(x1)) .* (dx1 - mean(dx1))); dC(1,2) = sum((x1 - mean(x1)) .* (dx2 - mean(dx2))); dC(2,1) = sum((x2 - mean(x2)) .* (dx1 - mean(dx1))); dC(2,2) = sum((x2 - mean(x2)) .* (dx2 - mean(dx2))); dC = dC / (N-1); % C_infty = cov(x1, x2); C_infty = C; detc = det(C); a11 = C(2,2) * dC(1,1) - C(1,2) * dC(2,1); a12 = -C(1,2) * dC(1,1) + C(1,1) * dC(2,1); % a21 = -C(1,2) * dC(1,2) + C(1,1) * dC(2,2); % a22 = C(2,2) * dC(1,2) - C(1,2) * dC(2,2); a11 = a11 / detc; a12 = a12 / detc; % a21 = a21 / detc; % a22 = a22 / detc; f1 = mean(dx1) - a11 * mean(x1) - a12 * mean(x2); % f2 = mean(dx2) - a21 * mean(x1) - a22 * mean(x2); R1 = dx1 - (f1 + a11*x1 + a12*x2); % R2 = dx2 - (f2 + a21*x1 + a22*x2); Q1 = sum(R1 .* R1); % Q2 = sum(R2 .* R2); b1 = sqrt(Q1 * dt / N); % b2 = sqrt(Q2 * dt / N); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % covariance matrix of the estimation of (f1, a11, a12, b1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NI(1,1) = N * dt / b1/b1; NI(2,2) = dt/b1/b1 * sum(x1 .* x1); NI(3,3) = dt/b1/b1 * sum(x2 .* x2); NI(4,4) = 3*dt/b1^4 * sum(R1 .* R1) - N/b1/b1; NI(1,2) = dt/b1/b1 * sum(x1); NI(1,3) = dt/b1/b1 * sum(x2); NI(1,4) = 2*dt/b1^3 * sum(R1); NI(2,3) = dt/b1/b1 * sum(x1 .* x2); NI(2,4) = 2*dt/b1^3 * sum(R1 .* x1); NI(3,4) = 2*dt/b1^3 * sum(R1 .* x2); NI(2,1) = NI(1,2); NI(3,1) = NI(1,3); NI(3,2) = NI(2,3); NI(4,1) = NI(1,4); NI(4,2) = NI(2,4); NI(4,3) = NI(3,4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% invNI = inv(NI); var_a12 = invNI(3,3); % % approx. variance of a12, corr. to variance of T21 % % % Information transfer: T21 = C12/C11 * a12 % T12 = C21/C22 * a21 % T21 = C_infty(1,2)/C_infty(1,1) * (-C(2,1)*dC(1,1) + C(1,1)*dC(2,1)) / detc; % T12 = C_infty(2,1)/C_infty(2,2) * (-C(1,2)*dC(2,2) + C(2,2)*dC(1,2)) / detc; % var_T21 = (C_infty(1,2)/C_infty(1,1))^2 * var_a12; % % From the standard normal distribution table, % significance level alpha=95%, z=1.96 % 99%, z=2.56 % 90%, z=1.65 % z99 = 2.56; z95 = 1.96; z90 = 1.65; err90 = sqrt(var_T21) * z90; err95 = sqrt(var_T21) * z95; err99 = sqrt(var_T21) * z99;```解释这段代码
最新发布
03-27
<think>好的,我现在需要解释用户提供的MATLAB代码。首先,我需要理解这段代码的功能和结构。代码名为causality_est,估计两个时间序列X2到X1的信息传递,并输出T21以及不同显著性水平下的标准误差。 首先,我会通读代码,了解整体流程。函数输入是xx1、xx2和np,输出是T21和三个误差项。根据注释,该方法是基于Liang在2014和2015年的论文,涉及信息流(因果)的估计。 接下来,我会逐步分析代码的各个部分: 1. **输入处理**:检查输入参数,xx1和xx2是列向量,np是用于前向差分的时间步长。dt设为1,说明时间间隔为1单位。 2. **数据预处理**:计算dx1和dx2,即使用np步长的前向差分来近似导数。例如,dx1是xx1的未来np步减去当前值,除以np*dt。同时,x1和x2被截断为原始数据的前nm-np个点,以便与导数对齐。 3. **协方差矩阵计算**:计算x1和x2的协方差矩阵C,然后计算dC,即协方差的变化率。这里通过元素相乘求和后除以N-1得到。 4. **参数估计**:根据协方差矩阵和其变化率,计算系数a11和a12,这些可能对应于动态系统中的参数。然后计算残差R1,以及其平方和Q1,进而得到b1,即噪声项的标准差估计。 5. **信息矩阵和方差估计**:构建NI矩阵(可能是Fisher信息矩阵的负矩阵),然后求逆得到参数估计的协方差矩阵。这里关注的是a12的方差var_a12,进而计算T21的方差var_T21。 6. **信息传递计算**:根据公式,T21由C_infty的元素和a12等计算得出。这里的公式可能基于Liang的论文中的推导,将信息流表示为协方差与参数的乘积。 7. **标准误差计算**:使用正态分布的分位数(z值)乘以标准差,得到不同置信水平下的误差范围。 在理解过程中,我需要验证每一步的数学背景。例如,协方差矩阵C和其变化率dC的作用,参数a11和a12的推导,以及信息流T21的具体公式来源。同时,需要确认NI矩阵的构造是否符合最大似然估计或最小二乘法的信息矩阵计算方式。 可能存在的疑问点包括:为什么使用前向差分而不是其他数值微分方法?协方差变化率dC的计算是否有特定理论依据?信息流T21的公式如何从论文中的理论转化为代码中的表达式?此外,NI矩阵的具体形式是否与参数估计的协方差有关? 还需注意,代码中的注释提到C_infty被赋值为C,但原注释中C_infty = cov(x1, x2),而实际代码中是C_infty = C,可能因为前面已经计算过C,这里直接使用。需要确认是否正确。 最后,总结各部分如何组合起来实现信息流的估计,并解释输出结果的含义,即T21的大小及其在不同置信水平下的误差范围,用于判断因果关系的显著性。</think>以下是对该MATLAB代码的逐步解释: ### 1. **函数功能** 该函数用于估计两个时间序列$X_2 \to X_1$的因果信息流$T_{21}$,并计算不同显著性水平下的标准误差。基于Liang提出的信息流理论(参考文献见代码注释)。 --- ### 2. **输入参数** - `xx1`, `xx2`: 时间序列数据(列向量) - `np`: 前向差分的时间步长(用于数值微分) --- ### 3. **输出参数** - `T21`: $X_2$对$X_1$的因果信息流强度 - `err90`, `err95`, `err99`: 90%/95%/99%显著性水平下的标准误差阈值 --- ### 4. **关键步骤解析** #### 4.1 **数据预处理** - 计算导数近似: ```matlab dx1 = (xx1(np+1:end) - xx1(1:end-np)) / (np*dt) % X1的导数 x1 = xx1(1:end-np) % 截断后的X1 ``` 同理处理`dx2`和`x2`,通过前向差分近似导数。 #### 4.2 **协方差计算** - 协方差矩阵: ```matlab C = cov(x1, x2) % x1和x2的协方差矩阵 ``` - 协方差变化率`dC`: ```matlab dC(1,1) = sum((x1 - mean(x1)) .* (dx1 - mean(dx1))) / (N-1) % dC11/dt ``` 其他元素类似,反映协方差随时间的变化。 #### 4.3 **参数估计** - 动态系统参数: ```matlab a11 = [C(2,2)*dC(1,1) - C(1,2)*dC(2,1)] / det(C) a12 = [-C(1,2)*dC(1,1) + C(1,1)*dC(2,1)] / det(C) ``` 这些参数描述$X_1$和$X_2$对$X_1$变化的贡献。 #### 4.4 **残差与噪声估计** - 残差计算: ```matlab R1 = dx1 - (f1 + a11*x1 + a12*x2) % 实际导数与模型预测的差 ``` - 噪声标准差: ```matlab b1 = sqrt(sum(R1.^2) * dt / N) % 白噪声强度估计 ``` #### 4.5 **信息矩阵与方差** - 构建Fisher信息矩阵`NI`: ```matlab NI(2,2) = dt/(b1^2) * sum(x1.^2) % 参数a11的Fisher信息 NI(3,3) = dt/(b1^2) * sum(x2.^2) % 参数a12的Fisher信息 ``` - 求逆得协方差矩阵`invNI`,提取`a12`的方差`var_a12`。 #### 4.6 **信息流计算** - 最终信息流: ```matlab T21 = (C_infty(1,2)/C_infty(1,1)) * (-C(2,1)*dC(1,1) + C(1,1)*dC(2,1)) / det(C) ``` 该公式源于Liang理论,量化$X_2$对$X_1$的因果影响。 #### 4.7 **标准误差** - 通过正态分布分位数计算显著性阈值: ```matlab err95 = sqrt(var_T21) * 1.96 % 95%置信水平的误差阈值 ``` --- ### 5. **使用建议** - **数据要求**:输入序列需平稳,避免剧烈波动 - **参数`np`选择**:通常设为1,除非系统高度混沌 - **结果解读**:若`T21 > err95`,则认为$X_2 \to X_1$因果关系在95%水平显著 --- ### 6. **参考文献** - Liang, X. S. (2014). *Phys. Rev. E* 90, 052150 - Liang, X. S. (2015). *Phys. Rev. E* 92, 022126 该代码实现了一个基于协方差分析的因果推断方法,适用于气候、经济等复杂系统的因果分析。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值