维诺图 (Voronoi Diagram)

#include"Heder.h" bool key1 = true; bool key2 = false; bool key3 = true; bool key4 = false; double width = WIDTH; double height = HEIGHT; Mesh* pDelaunayMesh; Mesh* pVoronoiMesh; int delaunayVertexId = 0; int delaunayFaceId = 0; int voronoiVertexId = 0; int voronoiFaceId = 0; //用来存储DelaunayMesh的所有点 vector<CVertex*> allVertex; //用来存储ConvexHull算法产生的包围圈的点 vector<CVertex*> vecConvexHull; //用来存储凸包内部的点 vector<CVertex*> toBeAdd; //用来存储face对应的外接圆圆心 map<CFace*, MyPoint2> mpFaceCenter; //随机数初始化Vertex void initVertex(); //插入排序，从小到大，以x坐标大小排序，若x坐标相同，则按y坐标大小排序 void insertSort(); //向量积，大于0表示从oa到ob为逆时针旋转 GLfloat cross(CVertex* o, CVertex* a, CVertex* b); //Andrew's Monotone Chain void ConvexHull(); //计算圆心，输入为三个点 MyPoint2 getCenter(vector<CVertex*> vec); //计算圆心，输入为一个面 MyPoint2 getTriCenter(CFace* pFace); //将Face和相应的圆心Vertex加入mpFaceCenter中 void addMpFaceCenter(CFace* pFace); //将Face和相应的圆心Vertex从mpFaceCenter中删除 void delMpFaceCenter(CFace* pFace); //Euclidean Distance的平方 double eDistance(double x1, double y1, double x2, double y2); //判断一个点是否在一个三角面的外接圆中 bool isCircleContain(CFace* pFace, CVertex* pVertex); //根据ConvexHull生成的点集初始化Delaunay图 void initDelaunayMesh(); //判断是否为DelaunayTriangle bool isDelaunayTriangle(); //Delaunay Triangulation void DelaunayTriangulation(); //VoronoiDiagram算法 void VoronoiDiagram(); //帮助 void help(); //绘制点 void drawVertex(Mesh* pMesh); //绘制边 void drawEdge(Mesh* pMesh); //绘制面 void drawFace(Mesh* pMesh); //OpenGL图像绘制 void renderScene(); //窗口改变 void changeSize(int w, int h); //响应鼠标事件 void mouseClick(int button, int state, int x, int y); //响应键盘事件 void keyBoardFunc(unsigned char key, int x, int y); //初始化OpenGL void init_openGL(int argc, char *argv[]); //随机数初始化Vertex void initVertex() { CVertex* tmp = NULL; while (delaunayVertexId < VERTEXNUM) { tmp = (*pDelaunayMesh).createVertex(delaunayVertexId++); tmp->point().coord()[0] = (rand() % (int)(2 * ACCURACY) - ACCURACY) / ACCURACY; tmp->point().coord()[1] = (rand() % (int)(2 * ACCURACY) - ACCURACY) / ACCURACY; allVertex.push_back(tmp); } } //插入排序，从小到大，以x坐标大小排序，若x坐标相同，则按y坐标大小排序 void insertSort() { int i, j; CVertex* tmp = NULL; for (i = 0; i < allVertex.size(); i++) { tmp = allVertex[i]; j = i - 1; while (j >= 0 && allVertex[j]->point().coord()[0] >= tmp->point().coord()[0]) { if (allVertex[j]->point().coord()[0] == tmp->point().coord()[0] && allVertex[j]->point().coord()[1] <= tmp->point().coord()[1]) { break; } allVertex[j + 1] = allVertex[j]; j--; } allVertex[j + 1] = tmp; } } //向量积，大于0表示从oa到ob为逆时针旋转 GLfloat cross(CVertex* o, CVertex* a, CVertex* b) { return (a->point().coord()[0] - o->point().coord()[0]) *(b->point().coord()[1] - o->point().coord()[1]) - (a->point().coord()[1] - o->point().coord()[1]) *(b->point().coord()[0] - o->point().coord()[0]); } //Andrew's Monotone Chain void ConvexHull() { //将Vertex的顶点id排序 insertSort(); //记录包围点的个数 int m = 0; int size = allVertex.size(); //从左到右，得出下半部分包围点 for (int i = 0; i < size; i++) { while (m >= 2 && cross(vecConvexHull[m - 2], vecConvexHull[m - 1], allVertex[i]) < 0) { vecConvexHull.pop_back(); m--; } vecConvexHull.push_back(allVertex[i]); m++; } //从右到左，得出上半部分包围点，多存储了一次起点id for (int i = size - 2; i >= 0; i--) { while (cross(vecConvexHull[m - 2], vecConvexHull[m - 1], allVertex[i]) < 0) { vecConvexHull.pop_back(); m--; } vecConvexHull.push_back(allVertex[i]); m++; } vecConvexHull.pop_back(); } //计算圆心，输入为三个点 MyPoint2 getCenter(vector<CVertex*> vec) { //求出外接圆的圆心 double x0 = vec[0]->point().coord()[0]; double y0 = vec[0]->point().coord()[1]; double x1 = vec[1]->point().coord()[0]; double y1 = vec[1]->point().coord()[1]; double x2 = vec[2]->point().coord()[0]; double y2 = vec[2]->point().coord()[1]; double a0 = (x0 + x1) / 2; double b0 = (y0 + y1) / 2; double a1 = (x2 + x1) / 2; double b1 = (y2 + y1) / 2; //横纵坐标 double x; double y; if (y0 != y1 && y1 != y2) { double k1 = -(x0 - x1) / (y0 - y1); double c1 = b0 - k1 * a0; double k2 = -(x1 - x2) / (y1 - y2); double c2 = b1 - k2 * a1; x = (c2 - c1) / (k1 - k2); y = k1 * x + c1; } else if (y0 == y1) { x = a0; double k2 = -(x1 - x2) / (y1 - y2); double c2 = b1 - k2 * a1; y = k2 * x + c2; } else { x = a1; double k1 = -(x0 - x1) / (y0 - y1); double c1 = b0 - k1 * a0; y = k1 * x + c1; } MyPoint2 point2(x, y); return point2; } //计算圆心，输入为一个面 MyPoint2 getTriCenter(CFace* pFace) { vector<CVertex*> vec; for (fvIter iter(pFace); !iter.end(); iter++) { vec.push_back(*iter); } return getCenter(vec); } //将Face和相应的圆心Vertex加入mpFaceCenter中 void addMpFaceCenter(CFace* pFace) { MyPoint2 point2 = getTriCenter(pFace); mpFaceCenter.insert(pair<CFace*, MyPoint2>(pFace, point2)); } //将Face和相应的圆心Vertex从mpFaceCenter中删除 void delMpFaceCenter(CFace* pFace) { mpFaceCenter.erase(pFace); } //Euclidean Distance的平方 double eDistance(double x1, double y1, double x2, double y2) { return (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2); } //判断一个点是否在一个三角面的外接圆中 bool isCircleContain(CFace* pFace, CVertex* pVertex) { MyPoint2 point2 = mpFaceCenter.at(pFace); fvIter iter(pFace); double r = eDistance(point2.x, point2.y, (*iter)->point().coord()[0], (*iter)->point().coord()[1]); double distance = eDistance(point2.x, point2.y, pVertex->point().coord()[0], pVertex->point().coord()[1]); return r > distance; } //根据ConvexHull生成的点集初始化Delaunay图 void initDelaunayMesh() { vector<CVertex*> vecTmp = vecConvexHull; bool isContainOther = false; while (vecTmp.size() > 3) { int size = vecTmp.size(); for (int i = 0; i < size; i++) { isContainOther = false; vector<CVertex*> fv; fv.push_back(vecTmp[i]); fv.push_back(vecTmp[(i + 1) % size]); fv.push_back(vecTmp[(i + 2) % size]); MyPoint2 center = getCenter(fv); double r = eDistance(center.x, center.y, fv[0]->point().coord()[0], fv[0]->point().coord()[1]); for (int j = 0; j < vecConvexHull.size(); j++) { CVertex* pVertexTmp = vecConvexHull[j]; if (pVertexTmp != fv[0] && pVertexTmp != fv[1] && pVertexTmp != fv[2] && r > eDistance(center.x, center.y, pVertexTmp->point().coord()[0], pVertexTmp->point().coord()[1])) { isContainOther = true; break; } } if (!isContainOther) { mpFaceCenter.insert(pair<CFace*, MyPoint2>((*pDelaunayMesh).createFace(fv, delaunayFaceId++), center)); vecTmp.erase(vecTmp.begin() + (i + 1) % size); break; } } } addMpFaceCenter((*pDelaunayMesh).createFace(vecTmp, delaunayFaceId++)); } //判断是否为DelaunayTriangle bool isDelaunayTriangle() { for (mfIter iter(pDelaunayMesh); !iter.end(); iter++) { fvIter iter2(*iter); CVertex* v1 = *iter2; iter2++; CVertex* v2 = *iter2; iter2++; CVertex* v3 = *iter2; for (int i = 0; i < allVertex.size(); i++) { if (v1 != allVertex[i] && v2 != allVertex[i] && v3 != allVertex[i] && isCircleContain(*iter, allVertex[i])) { return false; } } } return true; } //Delaunay Triangulation void DelaunayTriangulation() { //初始化toBeAdd for (int i = 0; i < allVertex.size(); i++) { int num = count(vecConvexHull.begin(), vecConvexHull.end(), allVertex[i]); if (num == 0) { toBeAdd.push_back(allVertex[i]); } } //用包围点和里面的一点创建三角面 init