butterfly01 弧长函数就是一个曲线积分的问题 C2连续：有相同的公共点且有相同的公共点且公共点处有相同的一阶导数和二阶导数（切线斜率的变化率），

弧长函数就是一个曲线积分的问题


C2连续：有相同的公共点且有相同的公共点且公共点处有相同的一阶导数和二阶导数（切线斜率的变化率），


http://www.cnblogs.com/sparkmorry/


2014年1月30日

曲线参数化的Javascript实现（理论篇）

　　在关键帧动画的制作过程中，动画师在k物体运动的过程中，一般要确定2个参数：
　　　　1)运动轨迹（表示物体运动的路径）；
　　　　2)速度曲线（表示物体随时间的速度变化）。
　　对于运动轨迹通常选用一定的样条曲线，通过动画师给出关键点位置，通过曲线插值产生，如之前介绍的一种Cardinal样条曲线。速度曲线就是速度根据时间变化的曲线，速度曲线在一段时间上的积分累加就是运动过的路程（不是位移）。对于简单的运动类型如匀速直线运动（速度曲线是一个常数），匀加速直线运动（速度曲线是一条加速度为斜率的直线方程）都可以直接获得路程随时间变化的方程，从而得到每个时间点在曲线上通过的路径长度。因此就需要在样条曲线上根据运动过的路径的长度来确定物体在曲线上的位置（由坐标(u, P(u))表示），这个过程就是曲线的参数化。
　　如果速度曲线为V(t)，对应的路程曲线为s=A(t)，对于给定的一段路程s，需要在Q(u)上找到对应的u点，设样条曲线函数为Q(u)，上面的弧长函数为s=A(u)，，那么参数化的过程就是把弧长函数表达成反函数的形式，然后代人Q(u)，就可以得到一个根据弧长算出曲线上函数值的方程。
　　弧长函数s=A(u)是路程在u上的积分方程，所以没有解析解（所谓解析解就是通过一个求根方程直接表达成自变量u根据方程的根，如二次函数的求根方程），因此采取数值求解的方法。求解的过程分为两步：
　　1)计算弧长函数s＝A(u)
　　2)对给定的弧长s，通过二分查找法在A(u)中计算u的值

1.弧长函数s＝A(u)的计算
弧长函数就是一个曲线积分的问题，就二维图像来说，如图：

在方程Q(u)～Q(u+du)上面每一段，对s进行积分，则
 其中，

因为上述积分函数并非解析可积，一般采用数值积分技术求解，简单的采用扩展是Simpson方法（在以后篇中介绍），则

这里的每个f(u)就是我们的样条函数Q(u)，一般高阶系数就直接忽略了，对于一段积分s就可以用u值和展开式计算s的值。
2.二分查找
因为弧长函数是严格的单调递增函数（路径长度的累加），那么如图


     (1)如果S3<Sa, 则u3<ua (红点)，
          搜索区间更新成[u3,u2].
     (2)如果S3>Sa, 则u3>ua (蓝点)，
         搜索间更新成[u1,u3].
      重复上述递归过程，直到||A(u1)-A(u2)||<指定精度ε，就把u1当作匹配的u值
　　这样我们就得到了每个s值对应的u值。另外，因为对于每一段曲线的弧长函数不同，连接起来的曲线可以用一个弧长表来纪录，如纪录下第一段曲线的总弧长，在第二段曲线上就在第一段曲线的基础上加上去，查找的方法也可以减小计算的时间复杂度，方便二分查找。

posted @ 2014-01-30 00:11 苏打草莓 阅读(13) 评论(0) 编辑

 

2014年1月29日

Cardinal样条曲线的Javascript实现（代码篇）

由上一篇文章得到了Cardinal曲线的矩阵表达式，下面就这个矩阵表达式就可以来对曲线进行插值了。
这里选用了JS来实现，完全是因为之前交作业的时候还不知道怎么在Xcode里建完整的C++OpenGL的项目，所以就用js在浏览器这个最大的跨平台UI上写。。
先看之前得到的Cardinal插值的样条曲线的矩阵，对每一段曲线PkPk+1来说，都可以通过拟合一个被参数化为P(u)函数（0<= u <=1）的图形的形式：
，其中
同样，M*P得到的是系数矩阵[a b c d]的转置，s是引入的表示每个连接点处曲线尖锐程度的变量－张量
所以我们现在的期望是
输入：n个控制点的xy轴坐标数组x[],y[]，张量tension，每两个控制点之间插入的密度grain
输出：插值好的样条曲线上控制点的一个数组
1.构造函数 function CSpline(x, y, grain, tension, n)
拥有的成员变量和成员函数：

function CSpline(x, y, grain, tension, n) {   
    var jd = new Array(100);  //用户输入的控制点位置数组，临时变量存储
    this.n0=n;   //控制点个数
    this.m = new Array(16);  //Mc矩阵
    this.Spline = new Array(1024);  //插值好的控制点数组
    CSpline.prototype.GetCardinalMatrix = function(a1) {}    //通过给定的tension值生成Mc矩阵
    CSpline.prototype.Matrix = function(a, b, c, d, u){}   //计算P(u)的值，xy分量均相同
    CSpline.prototype.CubicSpline = function(np, knots, grain, tension) {}   //插值函数

    //临时函数，建立曲线控制点数组jd，并调用CubicSpline函数计算插值曲线
    function CSpline(x, y, grain, tension, n)      
    { 
        for(var i=1; i<=np; i++){
        　　jd[i] = new CPT(x[i-1],y[i-1]);   //内部点
        }
        jd[0] = new CPT(x[0],y[0]);  //补上隐含的整条曲线起始点
        jd[np+1] = new CPT(x[np-1],y[np-1]);    //补上隐含的整条曲线终止点
        np=np+2;   
        this.CubicSpline(np, jd, grain, tension);
    }
}      

这里成员变量通过构造函数分别在实例中存储，而成员函数则在对应的原型函数中存储。
2.计算Mc矩阵函数GetCardinalMatrix = function(a1)
根据输入的平滑度 tension 值计算Mc矩阵

CSpline.prototype.GetCardinalMatrix = function(a1) {
    this.m[0]=-a1; this.m[1]=2.0-a1; this.m[2]=a1-2.; this.m[3]=a1;                 
    this.m[4]=2.*a1; this.m[5]=a1-3.; this.m[8]=-a1; this.m[9]=0.;
    this.m[12]=0.; this.m[13]=1.; this.m[6]=3.-2*a1; this.m[7]=-a1;
    this.m[10]=a1; this.m[11]=0.; this.m[14]=0.; this.m[15]=0.;
} 

这里m是一个4*4的矩阵，a1表示tension值
3.计算P(u)的值
输入4个点Pk-1,Pk,Pk+1,Pk+2的值(x分量或者y分量，以及参数化好的u值)，返回计算得到的P(u)的值

CSpline.prototype.Matrix = function(p0, p1, p2, p3, u){ 
    //求解系数
    var a, b, c, d;
    a=this.m[0]*p0+this.m[1]*p1+this.m[2]*p2+this.m[3]*p3;    
    b=this.m[4]*p0+this.m[5]*p1+this.m[6]*p2+this.m[7]*p3; 
    c=this.m[8]*p0+this.m[9]*p1+this.m[10]*p2+this.m[11]*p3; 
    d=this.m[12]*p0+this.m[13]*p1+this.m[14]*p2+this.m[15]*p3; 
    return(d+u*(c+u*(b+u*a))); //au^3+bu^2+cu+d
}

4.主要绘制函数CSpline.prototype.CubicSpline = function(np, knots, grain, tension)

CSpline.prototype.CubicSpline = function(np, knots, grain, tension) {
    alpha = new Array(50);
    var k0, kml, k1, k2;
        //获取Mc矩阵
        this.GetCardinalMatrix(tension);
        //对每两个关键点之间的插值点进行参数化到0～1之间
    for(var i=0; i<grain; i++){
        alpha[i]=i*1.0/grain;
    }
        //从最开始的四个点开始，给第一段曲线插值
    kml = 0;
    k0 = 1;
    k1 = 2;
    k2 = 3;
    this.s = 0; //纪录总共插值后的点数
        //两次循环第一次对输入的控制点遍历，第二次对每两个控制点之间插值，分别计算xy分量上得出的插值后的函数值，k值分别＋1
    for(var i =1; i<this.n0; i++){
        for(var j=0; j<grain; j++){
            cpx = this.Matrix(knots[kml].x, knots[k0].x, knots[k1].x, knots[k2].x, alpha[j]);
            cpy = this.Matrix(knots[kml].y, knots[k0].y, knots[k1].y, knots[k2].y, alpha[j]);
            this.Spline[this.s] = new CPT(cpx, cpy);
            this.s++;
        }
        kml++; k0++; k1++; k2++;
    }
}

最后的结果截图如下：

完整的代码正在上传github。。

posted @ 2014-01-29 11:08 苏打草莓 阅读(9) 评论(0) 编辑

 

Cardinal样条曲线的Javascript实现（理论篇）

　　首先，要对样条曲线进行插值的原因是：希望通过给定的关键帧点生成一条希望的直线或者曲线。
1.直线插值
　　生成一条直线，给定直线首尾的关键点P0,P1，就能确定这条直线的特性，比如y=kx+b中的斜率k和y轴偏移值b。通过线性（P0,P1线性相关）插值（线性的给中间插上一定数量的点使看起来连续）的方式就可以得到我们要的线段。

　　　　　　　　　　图1.1
2.曲线插值
　　但是对于曲线来说比较难确定，我们要对于给定的参数生成唯一的一条曲线并且可以进行方便的调整。这里要确定每一小段曲线，我们需要4个参数，首尾点的位置(参数化后的P(u))及他们的方向（在代数上表现为一阶导数P'(u)的值），并进行平滑的过渡，通过插值来实现对中间点的计算。再将每个小段的曲线进行连接。

　　　　　图2.1
　　曲线的生成分为两个部分：
　　　　1)对每给定的两个关键帧点之间进行插值得到一段曲线（如图2.1）
　　　　2)连接两两关键帧点的曲线段（如图2.2对曲线段进行连接）

　　　　　　　　图2.2
2.1 三次多项式插值
　　先就两个关键帧点之间的插值做计算。我的理解主要的思想是通过一个参数化的函数曲线去拟合，因此要通过给定的条件（如图2.1中的P(0),P(1),P'(0),P'(1)）去确定这个函数的几个系数。然后通过在中间等间距的插点计算根据函数得到的值确定点的位置。这里选择参数化的三次多项式(）去拟合这条曲线，通过插值来求得中间点的位置。这里，选择三次多项式的原因是：灵活性和计算速度的折中方案。
　　1)与更高次多项式相比，三次样条只需较少的计算与存储空间，并且较为稳定；（只需要一个4*4的基函数矩阵）
　　2)与更低次多项式相比，三次样条在模拟任何曲线形状时显得更灵活。（能满足大多数情况下的曲率调整，比如二次抛物线太对称，一次就是直线）
2.2 曲线连续性
　　为了保证分段参数曲线从一段到另一端平滑过渡，可以在连接点处要求各种连续性（可以参考微积分中的函数连续性与导数的关系）。
　　C0连续：有相同的公共点，如图(a)
　　C1连续：有相同的公共点且公共点处有相同的一阶导数（切线斜率），如图(b)
　　C2连续：有相同的公共点且有相同的公共点且公共点处有相同的一阶导数和二阶导数（切线斜率的变化率），如图(c)

　　　　　　　　　　图2.3
2.3 自然样条曲线
　　一种简单的插值方法是自然样条插值，就是对如图2.2的一条有n+1个控制点（P0～Pn）n个分段的曲线求三次多项式系数这样就有4n个系数要确定（需要列出4n个方程才能求解出来），通过曲线内部的n－1个点每相邻曲线段公共点的一阶及二阶导数相等且都通过该点，能确定4n－4 [4*(n-1)]个方程，再加上两个“隐含”控制点P0和Pn的二阶导数为0得到4n个方程来求解。缺点是一个点发生了变化其他点都跟着受到影响，这里不做主要展开了。　　　　　　　　　　　　
2.4 Hermite插值
　　为了解决自然样条曲线不能局部控制的缺点，Hermite插值对每个控制点的切线进行了限制为D（由用户给出），这样切线就变成了D*Pk（D斜率＋经过该点），使每个曲线段仅依赖于端点位置的约束。如对每小段曲线的边界条件表示为：
     －－－－－－－公式1
这里P是向量，对于xy分量来说分别都符合这个边界条件。那么如何保证曲线连续的呢？即对于一个连接点Pk-1来说，这段曲线相应的终点为Pk，在下一段曲线里Pk则作为这段曲线的起始点，保证Pk这个点的函数值和导数值相等就可以。现在点的位置都是Pk，导数值都是DPk，二阶导数值都是D。
可以写出向量方程：

其中P和abcd都是向量，可以再xy分量上（2维情况下）分别表达成这种形式，如，y分量也是类似，组合成的表达，可以获得等价的矩阵形式表达：

以及导数的矩阵表达式：


通过公式1，用u＝0和u＝1替代上面的u就可以得到以下方程：

如Pk=P(0)=[0 0 0 1]*[a b c d]'，这里每一个曲线段上的点都被参数化到0～1之间，所以起始点就是P(0),终止点就是P(1)，利用线代方法（参考线性方程求解）通过4个参数对多项式系数求解

最后边界条件的表达式如下

其中M*P部分就是通过对Hermite函数推导得到的系数矩阵[a b c d]'，Pk和Pk+1由用户给定，DPk和DPk+1由给定D值计算获得，根据不同的控制方式可以生成有不同曲线，Cardinal曲线就是Hermite曲线的一种变形，通过对切线斜率进行控制来控制曲线的平滑度。
2.5 Cardinal样条曲线
 　　因为在Hermite曲线中斜率D还是需要用户给出，Cardinal样条曲线可以由相邻两控制点坐标来计算斜率从而使曲线完全由控制点Pk来决定，但是在这里引入了一个t作为在每个公共点上尖锐程度的控制值。
Cardinal曲线的边界条件方程：

这样控制点Pk和Pk+1处的斜率和弦PkPk+2和PkPk+1成正比（理解为表示比较接近的控制点对曲线曲率的影响）

其中t（称为张量）控制Cardinal样条和输入控制点之间的松紧程度，对曲线的影响程度可以通过下图看出来

通过跟Hermite类似的系数求解方式，将边界条件代入可以得到矩阵表达式：
P(u)

，其中

posted @ 2014-01-29 00:46 苏打草莓 阅读(21) 评论(0) 编辑

 

2014年1月20日