b-spline学习-系数计算及程序实践

内容

• b-spline的学习网址
• 学习理解
• 现成的代码
• 动手翻译系数计算函数
• 用opencv图形验证效果

b-spline的学习网址

对b-spline的介绍网址：
http://www.qiujiawei.com/b-spline-1/

根据上面的网址，找到下面的讲义，这个讲义，内容已经说的很清楚了：
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/

基本上按照顺序对b-spline的内容读了一次，但其实不够的，最好还是对前因后果都读一下会比较前面。

学习理解

b-spline是一种可以更灵活的表示曲线的数学模型，是对bezier曲线的继承发展。
主要的概念有：knots，control point，degree。

knots

b-spline的定义域是[0,1]之间，而knots是将定义域进行划分的，分成多段子定义域。

degree

degree这个我理解的不深，不好解析。

basic function

b-spline由一系列的basic function组成的，degree的不同，basic function的定义域也不同。

仅仅通过knots，degree就可以计算出basic function，也是系数。

递归形式如下：

control points

control points影响曲线的shape，basic function，得到定义域内某个值计算出来的point，密集的插值就可以得到一条曲线了。
具体的插值公式为：

所以，当knots确定，degree确定，basic function也就确定了，即系数函数与control point无关。

knots，degree，control point的数量必须满足计算关系：
假设knots的数量为m，degree的值为p，control point的数量为n，则：
m=n+p+1

现有代码

在github上有js版的代码，自己用浏览器控制台执行验证过，确实是可以的，其verify是通过的。
js版b-spline曲线计算代码

动手翻译系数计算函数

虽然已经有了现有代码，可以翻译为c++后直接使用，但为了更好的掌握这个知识，还是自己动手实现一遍系数计算的过程。
参考这里的伪代码：
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html

重要的是对着图去看和理解代码：

b_spline.h

#ifndef __B_SPLINE__#define __B_SPLINE__#include <vector>/** * @brief calc_b_spline_coefficient * @param max_control_point_cnt_idx * control point的最大下标，从0开始，所以，control point的数量为 max_control_point_cnt_idx+1 * @param p * degree * @param max_knots_idx * knots的最大下标，所以，knots的数量为 max_knots_idx+1 * @param knots * 具体的konts数量 * @param u * @param coeffis * 计算返回的系数 ，长度为 max_control_point_cnt_idx+1 */void calc_b_spline_coefficient(int max_control_point_cnt_idx,int p,int max_knots_idx,float* knots,float u,std::vector<float>& coeffis);#endif

b_spline.cpp

#include <iostream>#include <assert.h>#include "math/b_spline.h"/** * @brief calc_b_spline_coefficient * @param max_control_point_cnt_idx * control point的最大下标，从0开始，所以，control point的数量为 max_control_point_cnt_idx+1 * @param p * degree * @param max_knots_idx * knots的最大下标，所以，knots的数量为 max_knots_idx+1 * @param knots * 具体的konts数量 * @param u * @param coeffis * 计算返回的系数 ，长度为 max_control_point_cnt_idx+1 */void calc_b_spline_coefficient(int max_control_point_cnt_idx,int p,int max_knots_idx,float* knots,float u,std::vector<float>& coeffis){//    std::cout<<"max_knots_idx="<<max_knots_idx<<",max_control_point_cnt_idx="<<max_control_point_cnt_idx<<",p="<<p<<std::endl;    assert(max_knots_idx==max_control_point_cnt_idx+p+1 && "m=n+1+p");    for(int i=0;i<=max_control_point_cnt_idx;i++){        coeffis[i]=0;    }    int m=max_knots_idx+1;    float um=knots[m-1];    float u0=knots[0];    if(u==u0){        coeffis[0]=1.0;        return;    }    if(u==um){        coeffis[max_control_point_cnt_idx]=1.0;        return;    }    assert(u>=u0&&u<=um && "u is invalid!");    //find u in which segment    int k=-1;    for(int i=0;i<m-1;i++){        if(u>=knots[i] && u<knots[i+1]){            k=i;            break;        }    }//    std::cout<<"control_point_cnt="<<max_control_point_cnt_idx+1<<",m="<<m<<",u0="<<u0<<",um="<<um<<",u="<<u<<",k="<<k<<std::endl;    coeffis[k]=1.0;    for(int d=1;d<=p;d++){//        std::cout<<"coeffis[0]="<<coeffis[0]<<//                ",coeffis[1]="<<coeffis[1]<<",k="<<k<<",d="<<d<<",k+d="<<k+d<<",k-d="<<k-d<<",k+1="<<k+1<<",k-d+1="<<k-d+1<<"\n";        if(k-d>=0){            coeffis[k-d]=(knots[k+1]-u)/(knots[k+1]-knots[k-d+1])*coeffis[k-d+1];// right (south-west corner) term only//            std::cout<<"1 d="<<d<<",k-d="<<k-d<<",coeffis["<<k-d<<"]="<<coeffis[k-d]<<std::endl;        }        for(int i=k-d+1;i<=k-1;i++){            coeffis[i]=(u-knots[i])/(knots[i+d]-knots[i])*coeffis[i]+(knots[i+d+1]-u)/(knots[i+d+1]-knots[i+1])*coeffis[i+1];//            std::cout<<"2 d="<<d<<",coeffis["<<i<<"]="<<coeffis[i]<<std::endl;        }        coeffis[k]=(u-knots[k])/(knots[k+d]-knots[k])*coeffis[k];//  left (north-west corner) term only//        std::cout<<"3 d="<<d<<",coeffis["<<k<<"]="<<coeffis[k]<<std::endl;    }}

简单测试代码：

#include <iostream>#include "math/b_spline.h"using namespace std;int basic_test_b_spline(){    //参考http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-ex-1.html    float knots_arr[]={ 0, 0.25, 0.5, 0.75, 1};    int knots_size=sizeof(knots_arr) / sizeof(knots_arr[0]);    int max_knots_idx=knots_size-1;    std::vector<float> knots(knots_arr,knots_arr+knots_size);    int p=1;    float u=0.6;    for(int i=0;i<knots.size();i++){        std::cout<<"i="<<i<<" knots="<<knots[i]<<std::endl;    }    int max_control_point_cnt_idx=max_knots_idx-p-1;    vector<float> coeffis(max_control_point_cnt_idx+1);    calc_b_spline_coefficient(max_control_point_cnt_idx,p,max_knots_idx,knots_arr,u,coeffis);    std::cout<<"coefficients of u="<<u<<" p="<<p<<" are as follows:\n";    for(int i=0;i<=max_control_point_cnt_idx;i++){        std::cout<<coeffis[i]<<" ";    }    std::cout<<"\n";}int main(){    basic_test_b_spline();}

测试的数据用的是：
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-ex-1.html

用opencv图形验证效果

用了基本的计算系数的函数，就可以进行插值了，用opencv的图形来进行操作与显示。

思路

1 设定好degree，以及knots数组
2 创建一个500*500的空白图
3 鼠标在空白图上点击若干个控制点，数量满足要求：m=n+1+p
如果选多了，则只取前面的点
4 选够了控制点后，双击鼠标左键，进行插值，并将结果显示到界面中；
如果没有选够，则提示信息。

代码

测试代码如下：

#include <iostream>#include "math/b_spline.h"#include <opencv2/opencv.hpp>using namespace std;using namespace cv;int basic_test_b_spline(){    //参考http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-ex-1.html    float knots_arr[]={ 0, 0.25, 0.5, 0.75, 1};    int knots_size=sizeof(knots_arr) / sizeof(knots_arr[0]);    int max_knots_idx=knots_size-1;    std::vector<float> knots(knots_arr,knots_arr+knots_size);    int p=1;    float u=0.6;    for(int i=0;i<knots.size();i++){        std::cout<<"i="<<i<<" knots="<<knots[i]<<std::endl;    }    int max_control_point_cnt_idx=max_knots_idx-p-1;    vector<float> coeffis(max_control_point_cnt_idx+1);    calc_b_spline_coefficient(max_control_point_cnt_idx,p,max_knots_idx,knots_arr,u,coeffis);    std::cout<<"coefficients of u="<<u<<" p="<<p<<" are as follows:\n";    for(int i=0;i<=max_control_point_cnt_idx;i++){        std::cout<<coeffis[i]<<" ";    }    std::cout<<"\n";}//为了让得到的曲线，经过首尾控制点，对u0,和um进行p次重复，即u0和um共有p+1个。float knots_arr[]={ 0,0,0,0.2, 0.25, 0.5, 0.6,0.75,0.8, 1,1,1};int knots_size=sizeof(knots_arr) / sizeof(knots_arr[0]);int max_knots_idx=knots_size-1;int p=2;int max_control_point_cnt_idx=max_knots_idx  - p -1;vector<float> coeffis(max_control_point_cnt_idx+1);vector<Point> control_pts;Mat pict;//mouse event callbackvoid mouse_track_event(int event, int x, int y, int flags, void *param ){    //左键点击选则控制点    if (event==CV_EVENT_LBUTTONDOWN) {        control_pts.push_back(Point(x,y));        cout<<"click,need control point ="<<max_control_point_cnt_idx+1<<",now has:"<<control_pts.size()<<endl;        //控制点画圆        circle(pict,Point(x,y),3,Scalar(255),2);    }else if (event==CV_EVENT_RBUTTONDOWN){        //右键单击清空图片        pict-=pict;        //清空        control_pts.resize(0);    }else if (event==CV_EVENT_LBUTTONDBLCLK){        //左键单击进行拟合曲线        //插值，画图        if(control_pts.size()>=knots_size-p-1){            //可以进行了            //控制点画圆            pict-=pict;            for(int i=0;i<=max_control_point_cnt_idx;i++){                circle(pict,control_pts[i],5,Scalar(255),2);            }            for(int i=max_control_point_cnt_idx+1;i<control_pts.size();i++){                circle(pict,control_pts[i],3,Scalar(0,255,0),2);            }            //            for(float u=0;u<1;u=u+0.01){                calc_b_spline_coefficient(max_control_point_cnt_idx,p,max_knots_idx,knots_arr,u,coeffis);                //计算                float u_x=0;                float u_y=0;                for(int i=0;i<=max_control_point_cnt_idx;i++){                    u_x+=coeffis[i]*control_pts[i].x;                    u_y+=coeffis[i]*control_pts[i].y;                }                cout<<"u="<<u<<",u_x="<<u_x<<",u_y="<<u_y<<endl;                circle(pict,Point((int)u_x,(int)u_y),1,Scalar(0,0,255));            }        }else{            cout<<"not enough control point~!"<<endl;        }    }    imshow("b-spline",pict);}void test_show_b_spline(){    namedWindow("b-spline");    pict=Mat::zeros(500,500,CV_8UC3);    setMouseCallback("b-spline",mouse_track_event);    waitKey(0);}int main(){//    basic_test_b_spline();    test_show_b_spline();}

效果图：