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

其他 2018-06-20 09:44:56 阅读次数: 2

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

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

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

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

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

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

（LateX公式真好用！）其中有n-2个方程，但是有n个未知数，所以方程没法解。那么我们还要考虑考虑端点的切线，老师上课讲的是把两个端点看做自由点（即端点切矢为 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);      //非全局定义不能初始化，所以记得初始化为0
        for (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;
}

猜你喜欢

转载自blog.csdn.net/wyg1997/article/details/78251887
今日推荐
周排行
Leetcode简单题61~80
解决zookeeper磁盘IO高的问题
多线程相关方法详解
Maven-setting.xml文件详解
Maven 项目的 classpath 理解
渊亭科技大数据笔试题
配置JVM内存分配
计算机网络个人学习笔记 （三）网络层 ：第三部分 连载
js中两个等号(==)和三个等号(===)的区别
用C程序自动打开电脑上的程序
每日归档
更多
2024-09-18(0)
2024-09-17(0)
2024-09-16(0)
2024-09-15(0)
2024-09-14(0)
2024-09-13(0)
2024-09-12(0)
2024-09-11(0)
2024-09-10(0)
2024-09-09(0)

代码天地 © 2018 codetd.com

境外服务器   最新电视剧2022