opengl:凸包算法

准备工作

判断点在有向线段的左侧

可以通过叉积判断，如下为k在有向线段ab的左侧代码描述：

double multiply(Point a, Point b, Point k)
{
    double x1 = b.x-a.x;
    double y1 = b.y-a.y;
    double x2 = k.x-a.x;
    double y2 = k.y-a.y;
    return x1*y2-x2*y1;
}

bool toLeft(Point a, Point b, Point k)
{
    return multiply(a,b,k)>0;
}

判断点在三角形的内部

给三角形abc定义一定的次序，按照一般习惯，假设abc是逆时针的，则判断k是否在三角形内部，只需要判断k是否在有向线段ab,bc,ac的左侧：

bool inTriangle(Point a, Point b, Point c, Point k)
{
    bool abLeft = toLeft(a,b,k);
    bool bcLeft = toLeft(b,c,k);
    bool caLeft = toLeft(c,a,k);
    return (abLeft==bcLeft)&&(bcLeft==caLeft);
}

几种典型的算法

极点算法

凸包上的顶点称为极点，极点有一个特性，总可以找到过极点的一条直线使得其他所有的顶点，在这个直线的一侧。所以极点不可能在某一个顶点三角形的内部，则可以在初始化时，标示所有的顶点为极点，然后遍历所有的顶点组成的三角形，排除掉三角形内部的顶点，则剩下的顶点则为凸包的极点。该算法时间复杂度为O(N^4)，算法描述如下：

极边算法

凸包上的边称为极边，所有的顶点都在极边的一侧，所以可以遍历所有的边，检查它是否为极边，算法时间复杂度为O(N^3),算法描述如下：

GiftWrapping算法

两个相邻的极边之间有一个共同的极点，所以一条极边的尾端也是另一条极边的顶端。如果已知一个极点，则可以寻找以该极点作为顶端的极边的尾端极点。方法是任取一个点作为候选点，如果下一个点在已知点与候选点组成的有向线段的右端，则把这个点作为候选点，这样不断的更新。因为最下的点肯定是一个极点，所以可把最下点作为初始点。算法的复杂度为O(N*W)(W是凸包的边数)，算法描述如下：

Graham Scan算法

算法需要借助一次排序，和两个栈：

下图描述了整个流程：

opengl实现

geometry.h文件：

#include <GL/glut.h>
#include <vector>
#include <algorithm>
#include <stack>
using namespace std;

class Point
{
public:
    Point(){};
    Point(double a,double b,double c):x(a),y(b),z(c),extreme(true){};
public:
    double x;
    double y;
    double z; //平面凸包，此项为0
    bool extreme;   //EE算法中用到，标识该点是否为极点
    pair<double,double> ref; //Graham Scan算法中，排序的参考点
    bool operator < (const Point &a);

};

class Edge
{
public :
    Edge(Point a,Point b):s(a),e(b){};
    Point s,e;
};

enum  PLOTMODE
{
    EP=0,
    EE,
    GW,
    GS
};

double multiply(Point a, Point b, Point c);

bool toLeft(Point a, Point b, Point c);

bool inTriangle(Point a, Point b, Point c, Point k);

 void EPAlgorithm(vector<Point> &list);

 vector<Edge> EEAlgorithm(vector<Point> list);

 vector<Edge> GWAlgorithm(vector<Point> list);

 stack<Point> GSAlgorithm(vector<Point>list);

geometry.cpp文件：

#include "Geometry.h"

bool Point::operator<(const Point &a)
{
    Point p = Point(ref.first,ref.second,0);
    return toLeft(p,*this,a);
}

double multiply(Point a, Point b, Point k)
{
    double x1 = b.x-a.x;
    double y1 = b.y-a.y;
    double x2 = k.x-a.x;
    double y2 = k.y-a.y;
    return x1*y2-x2*y1;
}

bool toLeft(Point a, Point b, Point k)
{
    return multiply(a,b,k)>0;
}

bool inTriangle(Point a, Point b, Point c, Point k)
{
    bool abLeft = toLeft(a,b,k);
    bool bcLeft = toLeft(b,c,k);
    bool caLeft = toLeft(c,a,k);
    return (abLeft==bcLeft)&&(bcLeft==caLeft);
}

//极点算法
void EPAlgorithm(vector<Point> &list)
{
    //三重循环遍历所有三角形
    for(int i=0;i<list.size();i++)
    {
        for(int j=i+1;j<list.size();j++)
        {
            for(int k=j+1;k<list.size();k++)
            {
                for(int m=0;m<list.size();m++)
                {
                    if(!list[m].extreme)
                        continue;
                    if(m==i||m==j||m==k)
                        continue;
                    if(inTriangle(list[i],list[j],list[k],list[m]))
                        list[m].extreme = false;
                }
            }
        }
    }
}

//极边算法
vector<Edge> EEAlgorithm(vector<Point> list)
{
    vector<Edge> res;
    //二重循环遍历所有边
    for(int i=0;i<list.size();i++)
    {
        for(int j=i+1;j<list.size();j++)
        {
            bool left = true, right = true;
            for(int k=0;k<list.size();k++)
            {
                if(k!=i&&k!=j)
                    (toLeft(list[i],list[j],list[k])>0?left:right) = false;
            }
            if(left|right)
                res.push_back(Edge(list[i],list[j]));
        }
    }
    return res;
}

//GiftWrapping算法
vector<Edge> GWAlgorithm(vector<Point> list)
{
    vector<Point> listCopy = list;
    vector<Edge> res;
    if(listCopy.size()<=2)
        return res;
    int ltl = 0;
    //找出lowest-then-leftest的点
    for(int i=1;i<listCopy.size();i++)
    {
        if(listCopy[i].y<listCopy[ltl].y||(listCopy[i].y==listCopy[ltl].y&&listCopy[i].x<listCopy[ltl].x))
            ltl = i;
    }
    int p = ltl;
    //找出下一条极边
    while(1)
    {
        int q = -1;
        for(int i=0;i<listCopy.size();i++)
        {
            if(i!=p&&(q<0||!toLeft(listCopy[p],listCopy[q],listCopy[i])))
                q = i;
        }
        res.push_back(Edge(listCopy[p],listCopy[q]));
        if(q==ltl)
            break;
        p = q;
    }
    return res;
}

Point getPoint(stack<Point>s, int num)
{
    stack<Point> temp = s;
    for(int i=0;i<num;i++)
    {
        temp.pop();
    }
    return temp.top();
}

stack<Point> GSAlgorithm(vector<Point>list)
{
    vector<Point> listCopy = list;
    stack<Point> S,T;
    if(listCopy.size()<3)
        return S;
    int ltl = 0;
    //找出lowest-then-leftest的点
    for(int i=1;i<listCopy.size();i++)
    {
        if(listCopy[i].y<listCopy[ltl].y||(listCopy[i].y==listCopy[ltl].y&&listCopy[i].x<listCopy[ltl].x))
            ltl = i;
    }
    //给所有定点附加ref属性
    for(int i=0;i<listCopy.size();i++)
    {
        listCopy[i].ref = pair<double,double>(listCopy[ltl].x,listCopy[ltl].y);
    }
    S.push(listCopy[ltl]);
    listCopy.erase(listCopy.begin()+ltl);
    //对定点进行排序
    sort(listCopy.begin(),listCopy.end());
    //构造初始的S和T
    S.push(listCopy[0]);
    for(int i = listCopy.size()-1;i>=1;i--)
    {
        T.push(listCopy[i]);
    }
    while(!T.empty())
    {
        while(!toLeft(getPoint(S,1),getPoint(S,0),getPoint(T,0)))
        {
            S.pop();
        }
        S.push(getPoint(T,0));
        T.pop();
    }
    return S;
}

convexHull.cpp文件：

#include <iostream>
#include <vector>
#include <GL/glut.h>
#include "Geometry.h"

using namespace std;

GLsizei width = 600,height = 600;
vector<Point> list;

PLOTMODE mode ;

void init()
{
    glClearColor(0.0f,0.0f,0.0,1.0f);
    glViewport(0,0,width,height);
    glMatrixMode(GL_PROJECTION);
    gluOrtho2D(0,width,0,height);
}

//选择菜单
void selectMenu(GLint option)
{
    switch (option)
    {
    case 1:
        list.clear();
        glutPostRedisplay();
        break;
    case 2:
        mode = EP;
        glutPostRedisplay();
        break;
    case 3:
        mode = EE;
        glutPostRedisplay();
        break;
    case 4:
        mode = GW;
        glutPostRedisplay();
        break;
    case 5:
        mode = GS;
        glutPostRedisplay();
        break;
    default:
        break;
    }
}

void plotPoints(vector<Point> list)
{
    glBegin(GL_POINTS);
    for(int i=0;i<list.size();i++)
    {
        glVertex2i(list[i].x,list[i].y);
    }
    glEnd();
}

void plotEP(vector<Point>list)
{
    EPAlgorithm(list);
    glPointSize(3.0);
    glBegin(GL_POINTS);
    for(int i=0;i<list.size();i++)
    {
        if(list[i].extreme)
            glVertex2d(list[i].x,list[i].y);
    }
    glEnd();
    glPointSize(1.0);
}

void plotEE(vector<Point>list)
{
    vector<Edge> edge = EEAlgorithm(list);
    glBegin(GL_LINES);
    for(int i=0;i<edge.size();i++)
    {
        glVertex2d(edge[i].s.x,edge[i].s.y);
        glVertex2d(edge[i].e.x,edge[i].e.y);
    }
    glEnd();
}

void plotGW(vector<Point>list)
{
    vector<Edge> edge = GWAlgorithm(list);
    glBegin(GL_LINES);
    for(int i=0;i<edge.size();i++)
    {
        glVertex2d(edge[i].s.x,edge[i].s.y);
        glVertex2d(edge[i].e.x,edge[i].e.y);
    }
    glEnd();
}

void plotGS(vector<Point> list)
{
    stack<Point> s = GSAlgorithm(list);
    glBegin(GL_LINE_LOOP);
    {
        while (!s.empty())
        {
            glVertex2d(s.top().x,s.top().y);
            s.pop();
        }
    }
    glEnd();
}

void displayFunc()
{
    glClear(GL_COLOR_BUFFER_BIT);
    glColor3f(1.0,0.0,0.0);
    glPointSize(1.0);
    plotPoints(list);
    cout<<mode<<endl;
    switch (mode)
    {
    case EE:
        plotEE(list);
        break;
    case EP:
        plotEP(list);
        break;
    case GW:
        plotGW(list);
        break;
    case GS:
        plotGS(list);
        break;
    default:
        break;
    }
    glFlush();
}

void reshapeFunc(GLint newWidth,GLint newHeight)
{
    glViewport(0,0,newWidth,newHeight);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0.0,GLdouble(newWidth),0.0,GLdouble(newHeight));
    width = newWidth;
    height = newHeight;
}

void mouseFunc(GLint button, GLint action, GLint x,GLint y)
{
    if(button==GLUT_LEFT_BUTTON&&action==GLUT_DOWN)
    {
        list.push_back(Point(x,height-y,0));
        glutPostRedisplay();
    }
}

int main(int argc,char* argv[])
{
    glutInit(&argc,argv);
    glutInitDisplayMode(GLUT_SINGLE|GLUT_RGB);
    glutInitWindowPosition(100,100);
    glutInitWindowSize(width,height);
    glutCreateWindow("Convex Hull");

    init();
    glutCreateMenu(selectMenu);
    glutAddMenuEntry("清除点",1);
    glutAddMenuEntry("极点算法",2);
    glutAddMenuEntry("极边算法",3);
    glutAddMenuEntry("GiftWrapping算法",4);
    glutAddMenuEntry("Graham Scan算法",5);
    glutAttachMenu(GLUT_RIGHT_BUTTON);
    glutDisplayFunc(displayFunc);
    glutReshapeFunc(reshapeFunc);
    glutMouseFunc(mouseFunc);

    glutMainLoop();
}
时间： 2024-10-05 04:04:48

opengl:凸包算法的相关文章

计算几何-凸包算法 Python实现与Matlab动画演示

凸包算法是计算几何中的最经典问题之一了.给定一个点集,计算其凸包.凸包是什么就不罗嗦了 本文给出了<计算几何——算法与应用>中一书所列凸包算法的Python实现和Matlab实现,并给出了一个Matlab动画演示程序. 啊,实现谁都会实现啦╮(╯▽╰)╭,但是演示就不一定那么好做了. 算法CONVEXHULL(P)  输入:平面点集P  输出:由CH(P)的所有顶点沿顺时针方向组成的一个列表 1.   根据x-坐标,对所有点进行排序,得到序列p1, …, pn 2.   在Lupper中加入p

Graham Scan凸包算法

获得凸包的算法可以算是计算几何中最基础的算法之一了.寻找凸包的算法有很多种,Graham Scan算法是一种十分简单高效的二维凸包算法,能够在O(nlogn)的时间内找到凸包. 首先介绍一下二维向量的叉积(这里和真正的叉积还是不同的):对于二维向量a=(x1,y2)和b=(x2,y2),a×b定义为x1*y2-y1*x2.而它的几何意义就是|a||b|sin<a,b>.如果a与b夹角小于180度(逆时针),那么这个值就是正值,大于180度就是负值.需要注意的是,左乘和右乘是不同的.如图所示:

openlayer的凸包算法实现

最近在要实现一个openlayer的凸多边形,也遇到了不小的坑,就记录一下 1.具体的需求: 通过在界面点击,获取点击是的坐标点,来绘制一个凸多边形. 2.思考过程: 1)首先,我们得先获取点击事件发生时,触发的点的坐标 map.events.register('click', map, function (e) { var pixel = new OpenLayers.Pixel(e.xy.x,e.xy.y); var lonlat = map.getLonLatFromPixel(pixel

线段余弦角+凸包算法

/// /// 根据余弦定理求两个线段夹角 /// /// 端点 /// start点 /// end点 /// double Angle(PointF o, PointF s, PointF e) { double cosfi = 0, fi = 0, norm = 0; double dsx = s.X - o.X; double dsy = s.Y - o.Y; double dex = e.X - o.X; double dey = e.Y - o.Y; cosfi = dsx * de

凸包算法

先理解下凸包 说凸包首先要说凸性的定义,简单点说就是平面邻域中任意两点所在的线段上的点都在该邻域中,则该邻域具有凸性.简单推敲一下,就可以发现如果邻域中存在一阶导数不连续的点一定无法被某点集线性表示出来.再往下的内容属于数学分析了,对我们的算法设计帮助不大,暂时先不管. 一般的计算几何问题都是处理的离散点集形成的平面域,所以我们感兴趣的是怎样找一个包含这个点集的面积最小的凸多边形,这就是凸包.作为常识也应该知道凸包上的顶点必然是该点集的子集,所以根据此性质我们就可以设计高效算法. 下面将介绍三种

平面上圆的凸包算法

平面上圆的凸包算法 我们之前探讨过这个有趣的问题: 平面上有若干圆,求包含这些圆的所有凸集的交. 根据之前讨论的结果,直接按圆心排序做扫描线的方法是错误的.我们需要考虑圆上的每一个点是否可以作为凸包上的一部分. 然而圆上的点的数目是无限多的.我们需要借助离散化的思想:因为我们发现凸包一定是由若干圆弧和线段构成的.而其中线段一定是切线段.由此,我们只需要将每两两圆的切点取出来求一个凸包.显然,在圆凸包上出现的切线段一定会在切点凸包中出现:而切点凸包中其余的线段则一定是弧弦. 但是,这个算法需要枚举

凸包算法的应用——数一数图形中共有多少三角形

一.问题引入 网络上经常会遇到判断图形个数的题目,如下例: 如果我们要把图中所有三角形一个一个选出来,在已知每个交点的前提下,该如何用代码判断我们选的图形是否是三角形呢.如下图,如何把图3筛选出来呢? 这里需要用到两步: 1.得到所选图形(阴影部分)所包含的所有小图形的顶点集合,求集合的凸包,根据凸包顶点个数判定凸包围成的图形是否是三角形,若顶点个数不为3则不是三角形,如图(1). .2.若凸包围成的图形是三角形,判断凸包的面积与所选图形(所有选中的小图形面积之和)是否相等,若相等则所选图形是三

OpenGL深度缓存算法

在OpenGL渲染过程中,深度测试对于多个物体之间的显示至关重要,如果不做适当处理,显示的结果将会与预期截然不容.所以深度缓存算法(depth buffer method)是一个比较常用的判别对象表面可见性的空间算法.它在投影面上的每一个像素位置比较场景中所有面的深度.该算法对场景各个对象表面单独进行处理,且在表面的逐点进行.该算法通常应用于只包含多边形面的场景,因为这些场景适合于很快地计算出深度值且算法易于实现,当然,该算法也可以应用于非平面的对象表面.由于通常沿着观察系统的z轴来计算对各对象

凸包算法-GrahamScan+暴力+分治

RT.求平面上点集的凸包. 1. GrahamScan算法,<算法导论>上的例子,先找到y最小的点O,以O建立极坐标,其它点按极角排序后再从头开始扫描(配合stack实现). 2.BruteForce算法,依赖定理:如果一个点在平面上某三个点组成的三角形内,那么这个点不可能是凸包上的点. 所以暴力的思路是平面上的点每4个进行枚举,并判断是否满足定理,若满足,则删除这个点继续找:一直找到没有满足定理的点为止.. 3.DivideAndConquer思路:有很多种,这里我实现的只是严格意义上最后一