Morry's Blog

A wizzard is never late.

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++;
    }
}

最后的结果截图如下:

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)

,其中

在Mac环境下跑汇编

今天汇编作业做到第七章,就想在Mac下跑自己的asm程序,看到了一个很好的教程

虽然对自己没有用,但是对汇编又有了一次了解,一些嵌入式设备系统如ios之类的是主要基于arm体系结构的操作系统,而pc机上的大多是基于intel体系结构的操作系统。查到了几种编译汇编代码的方式:

  1. 用nasm直接编译.asm文件:

    nasm -f macho hello.asm

    生成可执行文件:

    ld -o hello -e mystart hello.o

    运行

    ./hello

    检查退出

    echo $?

  2. 用gcc在.c文件中嵌套asm代码:

    gcc -Wall -m32 -O3 -fasm-blocks convert.c -o convert

    运行

    ./convert

  3. 本来想直接运行书里的代码,但是发现书里的代码都是用masm语法编程的,所谓masm就是windows平台下的汇编编译器(对应linux下是nasm),所以在自己的电脑上就坑爹的跑不起来了,但是查到还是可以用一些如DOSbox的工具编译的asm文件然后和c文件连起来的。下了DOSbox,还要配置实在有些麻烦就想着还是明天去实验室做大程的时候顺便写吧,懒得折腾了。

  4. 生成list文件

nasm -f elf 6-2.asm -l 6-2.lst

顺带出来的.o里面有一堆机器码,比原来的程序长很多 .lst文件是.asm文件翻译成机器码的结果就是一对一的翻译

反汇编 ndisasm 6-2.o

Install Opencv on Mac

下载好opencv之后

  1. 在文件夹下新建一个release或build的文件夹;

  2. cmake .

 make

Python中运行:

在该build文件夹下

nano .bash_profile

把python的路径加下去

export PYTHONPATH=/usr/local/lib/python2.7/site-packages/:$PYTHONPATH

Xcode中运行,可能碰到的一些问题:

  • usr/local/include和usr/local/lib的问题就不说了,编译之前要把头文件和动态库都连进去否则会产生dyld: Library not loaded的错误
  • dyld: Library not loaded: cv2.so,有放到cv2.so檔案的, 則是要記得到Build Phases(project target下)內,把選項改為optional
  • Undefined symbols for architecture x86_64: C++库和opencv库编译不一致,将XCode的C++ Standard Library改为GNU C++ standard library即可正常编译显示图片。Project Build settings –> APPLE LLVM compiler 4.2 – Language –> C++ Standard Library

在这里赞一下这个网站,里面有很多例程参考,每个都有不同语言的版本

Steganography

今天做DAM的作业,做到图片水印的时候,想起来当初小调同学把言页之庭的种子通过图片发给我。看到下面这个新闻真是觉得碉堡了!!

下面分享一下自己用python实现的隐写术,后面如果有时间会再更深入研究一下解码的问题!!

基本思路是:

1.每张图片都有RGB三个通道,每个像素点的值为3个8位的R/G/B值。

2.把每个RGB值的最后两位&11111100(0xFC)移掉,因为最后两位对图片的影响不大,肉眼基本不可见

3.把要隐藏的图片的RGB值都除以85(255/3=85),这样正好可以通过2bits(00-11,0-3)来表示图片

4.把得到的值填到之前的2个空bit中

5.等解码的时候只要把之前在原图中藏进去的2位提取出来(&0x3)再*85就可以得到原图信息

还有一种方法是把信息写道Alpha通道里并按RBG格式编码,解码时用L格式。

但是图片会失真一点,因为毕竟相比8位的存储,2位的存储信息会少很多,但是却已经可以看到80%的图片信息了。所以不得不再赞一下二八法则,20%的编码信息里存储了80%的图片信息啊。

import Image,ImageChops

def waterMark(originFile,markFile):
    origin = Image.open(originFile)
    mark = Image.open(markFile)
    size = mark.size #record the original size of the marked image
    mark = mark.resize(origin.size) #set the pics in same size

    origin.load();
    mark.load()
    source1=origin.split() #split the pic in 3 channel
    source2=mark.split()

    image=[(),(),()]
    waterMark=[(),(),()]
    for x in [0,1,2]: #loop in 3 channel, RGB
        image[x]=source1[x].point(lambda i:i & 0xFC) #remove the last 2 bits of the original image, 0xFC=11111100
        waterMark[x]=source2[x].point(lambda i:i/85) #reduce the RGB value of each channel, 255/3=85

    mark=Image.merge("RGB",waterMark)
    origin=Image.merge("RGB",image)

    result=ImageChops.add(mark,origin)

    return result,size

def deCode(waterMark,size):
    waterMark.load()
    source = waterMark.split()
    originSize=size 

    mark=[(),(),()]
    for x in [0,1,2]:
        mark[x]=source[x].point(lambda i:(i & 0x3)*85) #get back the RGB value of the marked image

    result=Image.merge("RGB",mark)
    result = result.resize(originSize)

    return result

if __name__ == '__main__':
    img1,size = waterMark("avatar2.jpg","green.jpg")
    img2 = deCode(img1,size)

    img1.show() #original image with watermark
    img2.show() #embeded image