OpenGL - Hermite算法多点画光滑曲线

本来以为矩阵求方程是最难的，没想到写个高斯消元一下就解决了，十分钟就写完了解方程的部分。然后花了将近四个小时查bug（QAQ）

说一下算法思路：
根据上一篇博客（传送门：点击打开链接），我们可以根据两点用Hermite算法绘制三次曲线。但是考虑多点问题时，光滑连接就是主要问题了，如果我们能求出中间点的切矢，那么就可以两两点绘制了。
所以我们的主要问题就是求出中间点的切矢。

我们知道，中间点和其左右点的光滑连线上，法向量应该是相等的，所以我们利用这一点列出两个方程，然后把P′′midP_{mid}^{''} 消去，最后得到一个方程：（一撇表示切矢，两撇表示法矢）

P′i−1+4P′i+P′i+1=3(Pi+1−Pi−1)

P_{i-1}^{'}+4P_{i}^{'}+P_{i+1}^{'}=3(P_{i+1}-P_{i-1})

然后我们利用这个公式，可以得到矩阵：

⎡⎣⎢⎢⎢⎢⎢⎢411411....114⎤⎦⎥⎥⎥⎥⎥⎥∗⎡⎣⎢⎢⎢⎢⎢⎢⎢P′1P′2..P′n−2⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢3(P2−P0)−P′03(P3−P1)..3(Pn−1−Pn−3−)P′n−1⎤⎦⎥⎥⎥⎥⎥⎥⎥

\begin{bmatrix} 4 & 1 & & & \\ 1 & 4 & 1 & & \\ & 1 & . & . & \\ & & . & . & 1\\ & & & 1 & 4 \end{bmatrix} * \begin{bmatrix} P_{1}^{'}\\ P_{2}^{'}\\ .\\ .\\ P_{n-2}^{'} \end{bmatrix} = \begin{bmatrix} 3(P_{2}-P_{0})-P_{0}^{'}\\ 3(P_{3}-P_{1})\\ .\\ .\\ 3(P_{n-1}-P_{n-3}-)P_{n-1}^{'} \end{bmatrix}

（LateX公式真好用！）其中有n-2个方程，但是有n个未知数，所以方程没法解。那么我们还要考虑考虑端点的切线，老师上课讲的是把两个端点看做自由点（即端点切矢为0⃗ \vec{0}），然后就是n-2个方程求解n-2个未知数了。

跑一发高斯消元（保证了一组解的条件，简直是天堂啊）。然后就得到了每个点的坐标和切矢。两两点用Hermite算法就OK啦。

程序使用说明：
绘图默认是自由点，也就是控制点和端点重合，使用时可以用左键把控制点拖出来。然后右键移动型值点，查看结果就好啦。
另外，如果按照标准的切矢计算的话，图形变化不是很明显，这时候可以认为的把端点的切矢放大4倍（实验中4倍是最棒的结果）。

运行结果：

代码如下：

#include<Windows.h>
#include<GL/glut.h>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<vector>
#include<cstring>
using namespace std;#define MatrixMaxSize 50
#define ACCURACY 200
#define CLR(a,b) memset(a,b,sizeof(a))bool mouseLeftIsDown = false;
bool mouseRightIsDown = false;
int cnt;        //点的数量struct Point
{double x, y;double tX, tY;      //该点处的切线(tangent)(相对矢量)Point(){}Point(double tx, double ty){x = tx;y = ty;}
};
struct Matrix
{double a[MatrixMaxSize][MatrixMaxSize];int w, h;void init(int th, int tw, double b[])       //初始化系数矩阵并形成增广矩阵{h = th;w = tw;CLR(a, 0);      //非全局定义不能初始化，所以记得初始化为0for (int i = 0; i < h; i++){a[i][i] = 4;if (i + 1 < h)a[i + 1][i] = 1;if (i + 1 < w - 1)a[i][i + 1] = 1;}for (int i = 0; i < h; i++){a[i][w - 1] = b[i];}}void Guass()        //高斯消元（默认只有一组解）{for (int i = 0; i < h-1; i++){//把系数最大的换到当前行，减小误差int t = i;for (int j = t+1; j < h; j++){if (a[j][i] > a[i][i])t = j;}if (t != i){for (int j = i; j < w; j++)swap(a[i][j], a[t][j]);}//然后消元for (int j = i + 1; j < h; j++){double mul = a[j][i] / a[i][i];     //消元因子for (int k = i; k < w; k++)     //变为行阶梯型{a[j][k] -= a[i][k] * mul;}}}/*puts("====");for (int i = 0; i < h; i++){for (int j = 0; j < w; j++){printf("%.2lf ", a[i][j]);}puts("");}puts("====");*///最后从下到上变为行最简型（a[i][w-1]为解）for (int i = h - 1; i >= 0; i--){for (int j = i + 1; j < h; j++)     //减去已知解a[i][w - 1] -= a[i][j] * a[j][w - 1];a[i][w - 1] /= a[i][i];     //最后除以系数}}
};
vector<Point> p;        //点集
int caculateSquareDistance(double x1,double y1,double x2,double y2)
{return (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
}
void mouse(int button, int state, int x, int y)       //监听鼠标动作
{if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN){mouseLeftIsDown = true;}if (button == GLUT_LEFT_BUTTON && state == GLUT_UP){mouseLeftIsDown = false;}if (button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN){mouseRightIsDown = true;}if (button == GLUT_RIGHT_BUTTON && state == GLUT_UP){mouseRightIsDown = false;}
}
void motion(int x, int y)       //移动点
{if (mouseLeftIsDown)        //左键移动控制点{if (caculateSquareDistance(x, y,p[0].x+p[0].tX,p[0].y+p[0].tY) < 400)     //防止鼠标移动过快点位无法及时读取，经测试，400为一个较适合的值{p[0].tX = x-p[0].x;p[0].tY = y-p[0].y;}else if (caculateSquareDistance(x, y, p[cnt - 1].x + p[cnt - 1].tX, p[cnt - 1].y+p[cnt - 1].tY) < 400){p[cnt - 1].tX = x-p[cnt - 1].x;p[cnt - 1].tY = y-p[cnt - 1].y;}}if (mouseRightIsDown)      //右键移动型值点{int tp = -1;double MinDis = caculateSquareDistance(x, y, p[0].x, p[0].y);for (int i = 0; i < cnt; i++){double t = caculateSquareDistance(x, y, p[i].x, p[i].y);if (t < 400){if (tp == -1){tp = i;MinDis = t;}else if (t < MinDis){tp = i;MinDis = t;}}}if (tp != -1){p[tp].x = x;p[tp].y = y;}}glutPostRedisplay();        //重新构图
}void getTangent()       //求出所有点的切线
{Matrix m;double arr[MatrixMaxSize];int n = p.size();       //总点数//先求x坐标for (int i = 0; i < n - 2; i++)arr[i] = 3.0*(p[i + 2].x - p[i].x);arr[0] -= p[0].tX;arr[n - 3] -= p[n - 1].tX;m.init(n - 2, n - 1, arr);m.Guass();for (int i = 1; i < n - 1; i++)p[i].tX = m.a[i - 1][m.w - 1];//再求y坐标for (int i = 0; i < n - 2; i++)arr[i] = 3.0*(p[i + 2].y - p[i].y);arr[0] -= p[0].tY;arr[n - 3] -= p[n - 1].tY;m.init(n - 2, n - 1, arr);m.Guass();for (int i = 1; i < n - 1; i++)p[i].tY = m.a[i - 1][m.w - 1];
}void initWindow(int &argc, char *argv[], int width, int height, char *title)    //初始化并显示到屏幕中央
{glutInit(&argc, argv);glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);glutInitWindowPosition((GetSystemMetrics(SM_CXSCREEN) - width) >> 1, (GetSystemMetrics(SM_CYSCREEN) - height) >> 1);       //指定窗口位置glutInitWindowSize(width, height);       //指定窗口大小glutCreateWindow(title);glClearColor(1, 1, 1, 0.0);glShadeModel(GL_FLAT);
}
void Reshape(int w, int h)      //两个参数：窗口被移动后大小
{glViewport(0, 0, w, h);glMatrixMode(GL_PROJECTION);glLoadIdentity();gluOrtho2D(0, w, h, 0);glMatrixMode(GL_MODELVIEW);glLoadIdentity();
}void Hermite(Point a,Point b)
{//求出相对于控制点的向量（切线）double delTa = 1.0 / ACCURACY;glBegin(GL_LINE_STRIP);for (int i = 1; i < ACCURACY; i++){double t = i * delTa;double t1 = 2 * pow(t, 3) - 3 * pow(t, 2) + 1;double t2 = -2 * pow(t, 3) + 3 * pow(t, 2);double t3 = pow(t, 3) - 2 * pow(t, 2) + t;double t4 = pow(t, 3) - pow(t, 2);glVertex2d(a.x*t1 + b.x*t2 + a.tX*t3 + b.tX*t4, a.y*t1 + b.y*t2 + a.tY*t3 + b.tY*t4);}glEnd();
}void myDisplay()
{glClear(GL_COLOR_BUFFER_BIT);       //清除。GL_COLOR_BUFFER_BIT表示清除颜色if (p.empty())      //如果是第一次定义{Point tp;for (int i = 0; i < cnt; i++){tp.x = 100 + i*800.0 / cnt;tp.y = 250;p.push_back(tp);}p[0].tX = 0;p[0].tY = 0;p[cnt - 1].tX = 0;p[cnt - 1].tY = 0;}glPointSize(10.0f);glColor3f(0, 0, 1);glBegin(GL_POINTS);     //画点（蓝色）for (int i = 0; i < cnt; i++)glVertex2d(p[i].x, p[i].y);glVertex2d(p[0].x + p[0].tX, p[0].y + p[0].tY);glVertex2d(p[cnt - 1].x + p[cnt - 1].tX, p[cnt - 1].y + p[cnt - 1].tY);glEnd();glLineWidth(3);glColor3f(1, 0, 0);glBegin(GL_LINES);      //两控制点和两端点连线glVertex2d(p[0].x, p[0].y);  glVertex2d(p[0].x+p[0].tX, p[0].y+p[0].tY);glVertex2d(p[cnt - 1].x, p[cnt - 1].y);  glVertex2d(p[cnt-1].x + p[cnt-1].tX, p[cnt-1].y + p[cnt-1].tY);glEnd();getTangent();       //得出所有切线for (int i = 0; i < cnt - 1; i++)Hermite(p[i], p[i+1]);/*for (int i = 0; i < cnt; i++)printf("%d %lf %lf %lf %lf\n",i, p[i].x, p[i].y, p[i].tX,p[i].tY);*/glFlush();glutSwapBuffers();
}int main(int argc, char *argv[])
{initWindow(argc, argv, 1000, 500, "Hermite");puts("\n\t使用Hermite算法，多点绘制光滑曲线。");puts("\t左键移动控制点，右键移动型值点(此处默认两端为自由端，也就是控制点和端点重合)");printf("\n\t请输入点的数量(2-20): ");while (~scanf("%d", &cnt)){if (cnt >= 2 && cnt <= 20){glutDisplayFunc(myDisplay);glutReshapeFunc(Reshape);glutMouseFunc(mouse);glutMotionFunc(motion);break;}else{puts("\t您的输入有误！请重新输入！");printf("\n\t请输入点的数量: ");}}glutMainLoop();return 0;
}

OpenGL - Hermite算法多点画光滑曲线相关推荐

java曲线平滑算法_JFreeChart简单实现光滑曲线绘制
用JFreeChart绘制光滑曲线,利用最小二乘法数学原理计算,供大家参考,具体内容如下绘制图形: 代码: FittingCurve.java package org.jevy; import ja ...
OpenGL实现Hermite算法绘制三次曲线
首先是推导:节省功夫我就直接贴照片了. 程序加了鼠标的监听器,可以移动控制点和型值点. 注意:图片中矩阵第二行第二列应该为3,当时笔误程序效果: 代码如下: #include<gl/glut. ...
WPF将点列连接成光滑曲线——贝塞尔曲线
背景最近在写一个游戏场景编辑器,虽然很水,但是还是遇到了不少问题.连接离散个点列成为光滑曲线就是一个问题.主要是为了通过关键点产生2D的赛道场景.总之马路不可能是直线相连的,当然需要曲线光滑相连.现 ...
折线图转成光滑曲线并求拐点
目录 1. 折线转成光滑曲线 2. 求拐点 3. paper 3.1 基于图像轮廓曲线的拐点检测算法摘要 0. 引言 1. APTD 拐点检测算法 1.1 算法原理 1.2 拐点判别函数的推导 3. ...
OpenGL Cubic Bezier三次贝塞尔曲线修补实例
OpenGL Cubic Bezier三次贝塞尔曲线修补先上图,再解答. 正常显示按下C键按下W键按下X键完整主要的源代码源代码剖析先上图,再解答. 正常显示按下C键按下W键
光滑曲线_对第一/二型曲线/曲面积分的小总结
公式第一型曲线积分(Line Integrals): 第二型曲线积分(Line Integrals of Vector Fields): 第一型曲面积分(Surface Integrals): 第二 ...
【基础】光滑曲线什么意思？以及 n次方差、n次方和公式、二项式定理（和的n次方）
目录一.光滑曲线二.n次方的差.和 1. 记忆 2. n次方差为什么有两个公式? 三.二项式定理一.光滑曲线同济七版注释:当曲线上每一点都具有切线,且切线随切点的移动而连续转动,这样的曲线是光 ...
光滑曲线_光滑曲线可求长定理证明
同济大学<高等数学>(第七版)中对于光滑曲线求长以及可求长定理没有做出过多的讨论,只是以可求长为前提推导了长度公式. 若想严谨地证明存在性, 首先需要严格的定义. (1) 设是一光滑的曲 ...
光滑曲线_极简微积分——函数的曲线的描绘
我们曾为导数是什么,导数如何计算付出了很多的努力去搞明白它到底是什么一回事.导数在物理学,经济学中都发挥了重要的作用,下面将讲述一些导数在函数图像分析方面的简单应用,以证明我们学过的导数不只是理论上的 ...

OpenGL - Hermite算法多点画光滑曲线

OpenGL - Hermite算法多点画光滑曲线相关推荐

最新文章

热门文章