任意多边形三维屋顶自动生成算法

任意多变性三维屋顶自动生成算法

        已经在CSDN猫了一年了，通过这里找到许多资源，特别是Android方面的。但是一直没有写过博客，很感激在一些牛人博客里找到自己问题的答案，所以从今天起我也将效仿大家，写一些自己做过的东西，感觉比较有价值技术探讨吧，初次写会有表达不清的地方，欢迎批评讨论。

       第一篇，写一下今年10月份做的一个三维屋顶算法。这个算法是为三维GIS根据地图数据自动构建建筑物的预研项目做基础的。当时自己遇到许多难题，最近要写一些算法文档，顺便就贴出来，相互讨论学习一下。这个算法中涉及几何库CGAL的编译使用，希望能给在第一次编译这个库的朋友一点提示。此外，这个算法虽然是以CGAL为基础的，但是我对于CGAL库的使用仍不熟练。CGAL是一个非常给力的库，是做几何图形这一块的不可或缺的优质资源。下面介绍一下任意多边形三维屋顶自动生成算法实现的过程。

一、任务要求：

1.      根据提供的多边形序列顶点和输入的参数类型自动生成坡面屋顶。

2.      输入的多变形内允许存在洞，洞的个数为任意。

3.      将生成的多边形屋顶以OSG模型格式输出。

4.      算法要求稳定、高效，能接受任意测试。

二、算法设计

1.设计思路

根据任务需要，选取CGAL开源库中的Straight_skeleton_2任意多边形骨架算法为底层骨架生成函数，通过对得到的骨架边进行最短路径生成，得打屋顶面的多变形拓扑关系。通过设定屋顶坡面的角度，经过得到骨架点的高度值。将生成的骨架点坐标和拓扑关系渲染到OSG场景中生成OSG格式模型。

2.过程步骤

   （1）编译CGAL开源库。

   （2）从Straight_skeleton_2算法中得多边形到骨架线的拓扑关系。

   （3）通过最短路径算法对骨架边进行遍历检查生成屋顶多边形面。

   （4）根据骨架点相对临近骨架线的距离设定骨架点的高度。

   （5）根据设定的屋面高度对屋顶进行截切，生成对应的屋顶平面。

   （6）将生成的屋顶拓扑关系输入到OSG场景中，渲染得到OSG模型格式。

三、算法实现过程

1.CGAL开源库的编译

  CGAL是比较经典的开源几何算法库，算法经典，稳定高效。其编译有一定难度，需要按照相应的步骤进行配置才能成功。其难点主要有三：1需要依赖BOOST库，QT环境，LibQGLviewer 。2.相应的环境搭建比较困难，设置问题直接导致编译失败。3由于BOOST库编译过程产生的版本问题需要修改相应的库文件名称，否则编译失败。（概括一点，就是麻烦，不过总能编出来！认真+细心+毅力=CGAL编译成功，加油！）网上的例子给出很多，但是不同的环境总是出现不同的情况，下面给出自己的编译过程……。

（1）       准备工作

本机操作系统：Win7+VS2010

【boost_1_56_0】

【CGAL-4.2】

【cmake-2.8.0-win32-x86】

【qt-opensource-windows-x86-vs2010-4.8.6】

【libQGLViewer】

       工欲善其事必先利其器，虽然有网友说需要部分功能可以不装某些模块就行，但是对于没有开源库编译经验的朋友还是老老实实做好以上充分准备，上面那些东西最好一样都不要少。否则，你将会懂得……我用的是这个版本，亲测成功，但不是唯一方案，可以参看其他网友，版本的下载链接请找度妈。

备注：有网友说BOOST库最好用在线安装版，但是好像那个网页失效了，就使用boost_1_56_0吧。还有CGAL，安装时需要联网，有时网速慢提示某些文件下载失败，我的解决办法是将其地址拷贝到浏览器中下载，然后放到相应的文件夹下，下面会具体提到。

（2）       第一批安装配置

为难于其易，先从软柿子开始。依次安装QT、libQGLViewe、Cmake，解压BOOST到一个文件夹里，这里我的文件夹是D:\Instruct。

为防止路径变量设置出错，下面给出我的路径安装环境变量配置，这个相当重要！如有不知道环境变量怎么设置，先到墙角面壁思过一会儿。

PATH：D:\CGAL-4.2;D:\CGAL-4.2\auxiliary\gmp\lib;

D:\Instruct\cmake\cmake-2.8.10-win32-x86\bin;

C:\Qt\4.8.6\bin;C:\Qt\bin;

D:\Instruct\boost_1_56_0\boost_1_56_0;

      C:\CGAL\lib;

 

BOOST_INCLUD：D:\Instruct\boost_1_56_0\boost_1_56_0

BOOST_LIBRAR：D:\Instruct\boost_1_56_0\boost_1_56_0\stage\lib

BOOST_ROOT：D:\Instruct\boost_1_56_0\boost_1_56_0

(3)BOOST库自编译

这个路径D:\Instruct\boost_1_56_0\boost_1_56_0\stage\lib，在编译之前是不存在的。需要首先运行bootstrap.bat生成，然后运行生成的bjam.exe；具体就不讲了，参见其他网友过程。本部分需要等待20分钟。

(4)安装CGAL，并设置环境变量路径。

   安装时候需要联网，如果网速慢，有文件下载失败，就根据提示用浏览器下载。然后放到相应的文件夹目录下。

(5)开始Cmake

                                   

图1 Cmake路径设置

根据上图1设置相应的路径。然后点击configure，会让你选择编译器，我用的是VS2010，选择2010这个过程大概两分钟，如果没有错，就点击generate，会在你选择的输出路径下有一个CGAL.sln。问题来了，很多情况下会提示各种异常，原因就在于路径设置有误，正确编译会显示下面三行内容，如果路径都正确，还有错误，就根据问题错误提示修改文件名，检查下面三条信息中是否配置正确。

D:/CGAL-4.2/auxiliary/gmp/lib/libgmp-10.lib

MPFR include: D:/CGAL-4.2/auxiliary/gmp/include

MPFR libraries:D:/CGAL-4.2/auxiliary/gmp/lib/libmpfr-4.lib

（6）打开CGAL.sln，在Debug和Release模式下都运行一遍，显示的是成功5，失败0。输出路径下的lib文件夹里会有CGAL-vc100-mt.lib、CGAL-vc100-mt-gd.lib等lib生成。

这就算编译完成了。

如果您一次编译成功，就太TMD幸运了。如果一天编出来说明您聪明，如果两天说明您能干，如果三天才编出来，说明您有毅力！我就很有毅力，呵呵！当初configure成功眼泪都快出来了！

2. 通过Straight_skeleton_2骨架生成算法得到骨架边拓扑关系

参考：

http://www.cgal.org/Manual/3.2/doc_html/cgal_manual/Straight_skeleton_2/Chapter_main.html

如下图2，

黑色的圈代表输入的顶点，蓝色的圈表示生成的骨架点，其中红色的线条为骨架边。

 

图2

通过Straight_skeleton_2 例子给出的骨架输出函数(代码如下)可以输出生成的骨架点和骨架边。

template<class K>

void print_straight_skeleton(CGAL::Straight_skeleton_2<K>const&ss )           （1）

{

  typedefCGAL::Straight_skeleton_2<K> Ss ;

 

  typedef typename Ss::Vertex_const_handle     Vertex_const_handle ;

  typedef typename Ss::Halfedge_const_handle   Halfedge_const_handle ;

  typedef typename Ss::Halfedge_const_iteratorHalfedge_const_iterator ;

 

  Halfedge_const_handle null_halfedge ;//半边控制handle

  Vertex_const_handle   null_vertex ;//顶点控制handle

 

  std::cout << "Straightskeleton with " << ss.size_of_vertices() //顶点数

            << "vertices, " << ss.size_of_halfedges()//半边的条数

            << "halfedges and " << ss.size_of_faces()//骨架面数

            << "faces" << std::endl ;

           

  for (Halfedge_const_iterator i = ss.halfedges_begin(); i != ss.halfedges_end(); ++i)

  {

    print_point(i->opposite()->vertex()->point());//输出骨架边起点坐标

    std::cout << "->";

    print_point(i->vertex()->point());//输出骨架边终点坐标

std::cout << "" << ( i->is_bisector() ?"bisector": "contour" ) << std::endl;

//"bisector" :"contour" 原始的边：生成的轮廓边

  }

}

例子中只给出了骨架边的信息输出，没有骨架面的输出，我试图找到骨架面的控制函数，

typedef typenameSs::Face_const_handle    Face_const_handle  ;

typedef typenameSs::Face_const_iterator    Face_const_iterator  ;

但是没有智能提示其中骨架面的输出方法。找了很久，找到了关于face的一些函数，但是没有操作成功，随后就放弃了直接得到骨架面拓扑关系想法，因为可以根据找到原始边和骨架边的关系，采用最短路径算法找到每条原始边对应的骨架面。

3.最短路径算法得到屋顶多边形

如图3所示，顶点1—8为原始多边形顶点，9—13为生成的骨架点，其连接关系如图中所示，可以由Straight_skeleton_2中的函数（1）得到。采用最短路径算法，将原始链接边的权重设为无穷大（黑色边：1000000），将生成的骨架边的权重设为1（红色边），然后依次计算原始边1->2,2->3，…，8->1的最短路径。即可得到每条原始边对应的骨架面。

最短路径算法参考：DijkstraShortestPathAlg。

 

图3

4、骨架点的高程值计算

由Straight_skeleton_2计算出的骨架点坐标为平面坐标，想得到三维屋顶，只要给生成的骨架点赋一个高度值。但是这个高度值是一个主观意识的设定，也许你想要每个点的高度不一样，这样可能更加真实。因此这里给出一个通过设定屋顶的坡面角度自动计算生成屋顶点高度的方法。如图3在一个骨架面中的五边形（1 —11—12—13—2）中，可以根据顶点11、12、13到边（1—2）的距离乘以屋面角度的正切值得到他们的高度值。当然这是也只是一种规定，假定原始点的高度值相等。你也可以设定为其他规则。

5设定屋顶平面高度得到屋顶切截面

图4(a)

如图4(a),是生成的全坡面屋顶，但有时候我们的屋顶不都是坡面，而是有坡面的房檐部分，有平的屋顶，像一个棱台那样。 所以根据以上的算法，我们得到屋顶点（即骨架点）的高度，然后用一个一定高度的平面来与屋顶相切，得到一个切截面。当然这个切截面有时候有多个，这需要我们设定一定的方法将他们找出来，如图4(b)所示，红色部分即为切截面。

s

图4(b)

6渲染成OSG模型格式

通过前几个步骤可以得到屋顶面的拓扑关系了，随后就是读入OSG渲染的过程了，当然OPENGL足可以胜任，只要有模型的拓扑关系用什么引擎渲染，就根据自己的喜好了。想要最后的屋顶模型为OSG格式，当然要有编译OSG库了，这个过程参见其他网友博客。

保存为OSG模型用这个函数 ：osgDB::writeNodeFile(*(root.get()),"RzXpdRoof.osg");一行代码，即可将绘制场景保存为OSG模型。

四、算法主要代码实现

（1）Straight_skeleton_2骨架生成算法得到骨架边拓扑关系

参照CGAL中example中示例代码：http://doc.cgal.org/latest/Straight_skeleton_2/index.html

在函数（1）中可以添加如下代码查看顶点编号：

std::cout <<i->opposite()->vertex()->id()<<" "<< i->vertex()->id();

（2）屋顶生成及切截面寻找

详见源码print.h中函数：ROOFgetxpd_straight_skeleton( CGAL::Straight_skeleton_2<K> const& ss ,PLOYGON merofploy ,float Angle,float maxhight )

 

（3）       OSG渲染参看工程C:\Users\zyp\Desktop\1111111111中main.cpp文件

（4）       最后整合工程为MFC截面，进行测试，详见工程C:\Users\zyp\Desktop\SSZRidegLine

五、算法测试分析

实现这一整套流程，用时近三周，现在思路清晰，觉得很简单，但当时摸着石头过河非常艰难。难点：1.CGAL编译；2.发现骨架边的拓扑关系；3.找到骨架面（我是用最短路径算法得到的。至今仍不知道怎么操作骨架面，对CGAL熟悉的朋友请您赐教）；4.找出切截面。

特别是找到切截面，我开始的算法又复杂，又有Bug，这个地方调试就卡了近一两周，真实痛苦，有80%时候得到的结果正确，其他就是三角片连错。开始怀疑是我渲染模式不正确，后来大量测试，才意识到找切截面算法身不严密。意识到问题出在那里后，一天就搞定了。

 

图5不含有洞的任意多边形

图6含有洞的情况

在SSZRidegLine工程中定义了ROOF结构体用来存储生成的模型拓扑关系。

工程C:\Users\zyp\Desktop\SSZRidegLine\cinploys.txt文件中定义了多边形输入读取格式（多边形最外多边形排序为顺时针，洞为逆时针）

PloygonNum:  6 10 3 12 5 7 3

6：表示总共有六个多边形//

ploygonouter:表示最外边的多边形//

ploygonhole://表示多边形中的洞

10 3 12 5 7 3分别依次表示6个多边形的顶点数。

ploygonouter:

134.000 97.0000.000 ///坐标（x，y）

273.000 68.0000.000

561.000 62.0000.000

530.000 298.0000.000

455.000 381.0000.000

363.000 130.0000.000

279.000 260.0000.000

148.000 206.0000.000

127.000 320.0000.000

91.000 117.0000.000

ploygonhole:

227.000 116.0000.000

197.000 149.0000.000

264.000 187.0000.000

ploygonhole:

481.000 119.0000.000

445.000 151.0000.000

453.000 193.0000.000

487.000 192.0000.000

510.000 183.0000.000

515.000 151.0000.000

517.000 114.0000.000

505.000 91.0000.000

449.000 81.0000.000

391.000 104.0000.000

406.000 124.0000.000

440.000 118.0000.000

ploygonhole:

464.000 253.0000.000

455.000 273.0000.000

469.000 293.0000.000

486.000 289.0000.000

501.000 264.0000.000

ploygonhole:

153.000 185.0000.000

204.000 209.0000.000

220.000 195.0000.000

137.000 148.0000.000

150.000 126.0000.000

118.000 121.0000.000

115.000 153.0000.000

ploygonhole:

336.000 99.0000.000

295.000 117.0000.000

317.000 133.0000.000

 

输出模型文件格式

C:\Users\zyp\Desktop\SSZRidegLine\Debug\faceoutxpd.txt

RoofObject: 1 3 11 14

///

1:表示该屋顶有削平顶，该位置次数为0则表示严格坡面屋顶

3：表示该屋顶有由三个多边形生成（即有两个洞）

11：表示最外边的多边形有11个顶点

14：表示共有14个屋顶面

///

PosId: 35 4 3 4

0 907.000 140.000 0.000///序号0；坐标（x,y,z）

1 916.000 459.000 0.000

2 116.000 420.000 0.000

3 228.000 100.000 0.000

4 382.000 197.000 0.000

5 332.000 259.000 0.000

6 507.000 332.000 0.000

7 814.000 257.000 0.000

8 709.000 287.000 0.000

9 705.000 390.000 0.000

10 811.000 390.000 0.000

11 675.253 418.615 28.592

12 624.759 414.220 30.522

13 842.971 422.701 32.675

14 380.184 143.822 34.769

15 861.972 194.873 46.521

16 862.564 404.457 51.835

17 855.390 190.053 52.959

18 230.129 284.846 63.023

19 640.432 380.605 64.833

20 285.195 185.860 82.282

21 618.737 216.125 92.871

22 250.119 329.674 96.673

23 589.450 242.684 121.084

24 882.637 163.628 25.000

25 890.228 432.694 25.000

26 150.684 396.641 25.000

27 245.378 126.087 25.000

28 380.694 158.763 25.000

29 291.590 269.253 25.000

30 603.455 399.345 25.000

31 839.780 223.614 25.000

32 684.702 267.921 25.000

33 678.990 415.020 25.000

34 835.462 415.020 25.000

Face: 4//屋顶面有4个顶点

3 27 24 0屋顶面的连接关系；记录顶点序号

Face: 4

0 24 25 1

Face: 4

1 25 26 2

Face: 4

2 26 27 3

Face: 4

6 30 28 4

Face: 4

4 28 29 5

Face: 4

5 29 30 6

Face: 4

10 34 31 7

Face: 4

7 31 32 8

Face: 4

8 32 33 9

Face: 4

9 33 34 10

Face: 4

27 24 25 26 27

Face: 3

30 28 29 30

Face: 4

34 31 32 33 34

算法介绍结束了！

感兴趣请留言评论，提供代码！