根據地圖數據程序化城市建模 Houdini - Procedural Modeling - Procedural City

花了大概一個禮拜左右的時間把城市建模的插件做出來了，還有一些瑕疵，但是已經達到預期想要的效果了。說在最前面的是，做這個東西完全沒有任何商業目的，只是希望在houdini里面掌握更多程序化的思維。

先看效果：

 

整體思路分為：

1；抓取地圖數據

2：分析地圖數據，並根據坐標位置來生成最基本的建築和道路的幾何形狀

3：根據道路和現有建築布局，用點填充其他空余空間

4：程序化生成建築樓，並將上一步的點提換成建築

 

下面將每一步的具體方法和思路講解一下：

1：抓取地圖數據：

現在市面上能夠免費得到地圖所有元素坐標和類別信息的機構我只找到了openstreetmap.org這一家，他們這個項目的目的也是為了能夠向世界所有人創造並提供免費的地理數據。我們可以直接在他的地圖網頁上框選一定的地圖空間然后導出相應經緯度的地圖坐標數據。非常重要的一點事他們目前只提供osm格式文件，其實就是xml的一種，如果想用編輯器直接查看里面的內容，可以吧后后綴名.osm改為.xml就可以打開了。

地圖相關資源：

地圖鏈接：https://www.openstreetmap.org

關於地圖到處的xml數據結構：http://wiki.openstreetmap.org/wiki/OSM_XML

關於xml幾個重要標記的解釋：http://wiki.openstreetmap.org/wiki/Elements

地圖中的元素一覽表：http://wiki.openstreetmap.org/wiki/Map_Features

 

這里簡單講解一下.osm格式的結構:

<osm>
	<bounds minlat="***" minlon="***" maxlat="***" maxlon="***"/>
	<node id="***" ******** lat="40.4444343" lon="-79.9486641"/>
	<node/>
	....
	<node/>
	<way id="27574365" ....>
  		<nd ref="1705867277"/>
  		<nd ref="302440336"/>
  		<tag k="highway" v="service"/>
 	</way>
 	....
 	<relation/>
 	...
 </osm>

 

bounds定義了整個數據的上下左右邊界經緯度，Node相當於是地圖中記錄下來了的每一個點的數據，他可以是道路上的也可以是建築上的。way是把所有點集結起來的一種方法，通過對他的tag做標簽能夠確定單個way下面所有的點是某一條街道還是建築或者其他物體。relation是node和way之間的所有關系，我們這里沒有用上所以不做介紹了。

 

2：分析地圖數據，並根據坐標位置來生成最基本的建築和道路的幾何形狀

地圖數據分析到這就基本上差不多了，下面我們需要一種機制能夠將地圖的xml數據傳入到houdini里面。houdini自身的程序語言Vex實現不了這個，但是好在houdini已經非常好的整合了Python進來，通過Python那問題就容易多了。

首先科普一下python查詢xml的參考：

1：http://www.cnblogs.com/fnng/p/3581433.html

2：http://www.cnblogs.com/xuxm2007/archive/2011/01/16/1936610.html

 

針對於houdini里面建立python節點並抓取數據我放了源代碼，這里我只講一講方法：

1：根據bounds做好邊界限定（源數據有很多超出邊界的雜點），並建立好縮放。

2：讀取每個node標簽，根據每一個node上面的經緯坐標增加一個對應的幾何點。並把所有有效node的id存入一個數組中。

3：讀取每個way標簽，遍歷每個way下面的所有ref中提到的id，判定該id是否為有效id，如果是調出該id對應的點，在點上增加vertex並用vertex生成多邊形線。

4：根據way里面對該對象描述的不同，確定是路還是建築，是封閉多邊形還是開放多邊形。

5：如果是路就是開放多邊形，並增加寬度給每個點；如果是建築則是封閉多邊形，這個多邊形則是之后擠出建築的底座形狀。

源碼：

# This code is called when instances of this SOP cook.
geo = hou.pwd().geometry()

# Add code to modify the contents of geo.

def mapCreate():
    import xml.dom.minidom

    mapScale = hou.pwd().parm("map_scale").eval()

    mapFile = hou.pwd().parm("map_data").evalAsString()
    map = xml.dom.minidom.parse(mapFile)

    #test whether the file path is legal
    if not (mapFile.endswith((".osm",".xml"))):
        print "This node just except xml or osm file."
        return

    mapIDs = []
    mapPoints = []

    streetWithScale = 0.0001 * hou.pwd().parm("street_width_scale").eval()

    polyAttrib  = geo.addAttrib(hou.attribType.Prim, "type", 0)
    streetWidth = geo.addAttrib(hou.attribType.Prim, "width", 0.0)


    #define the boundary of this xml file
    def getBoundary(mapData, mapScale):
        boundary = mapData.getElementsByTagName("bounds")[0]
        minLat = float(boundary.getAttribute("minlat"))
        maxLat = float(boundary.getAttribute("maxlat"))
        minLon = float(boundary.getAttribute("minlon"))
        maxLon = float(boundary.getAttribute("maxlon"))

        #return [scale, minlon, minlat, maxlon, maxlat]
        return [(mapScale / max((maxLon - minLon), (maxLat - minLat))), minLon, minLat, maxLon, maxLat]

    #draw points
    def drawNodes(mapNode):
        lontitude = float(mapNode.getAttribute("lon"))
        latitude  = float(mapNode.getAttribute("lat"))

        scaledLontitude = boundary[0] * (lontitude - 0.5 *(boundary[1] + boundary[3]))
        scaledLatitude  = boundary[0] * (latitude  - 0.5 *(boundary[2] + boundary[4]))
        if (boundary[1] <= lontitude <= boundary[3] and boundary[2] <= latitude <= boundary[4]):
            point = geo.createPoint()
            point.setPosition((scaledLontitude,0,-scaledLatitude))

            pointID = int(mapNode.getAttribute("id"))
            #this two lists are prepared for create ways
            mapIDs.append(pointID)
            mapPoints.append(point)

    #create ways no matter is closed or open loop
    def createWays(mapWays):
        
        #loop for each way tag
        for mapWay in mapWays:
            nodes = mapWay.getElementsByTagName("nd")
            tags  = mapWay.getElementsByTagName("tag")

            #put all points of one way into the wayPoint
            wayPoints = []
            #flag the define this "way" is closed
            isClosed = 0

            #flag for different elements
            isHighway   = 0
            isBuilding  = 0
            #these way is the "v" value of "highway"
            isFirstWay  = 0
            isSecondWay = 0
            isThirdWay  = 0
            isLastWay   = 0


            #define the way should be closed or open, and get different way types
            #the feature list of every elements in this xml map could be found in this link:
            #http://wiki.openstreetmap.org/wiki/Map_Features
            for tag in tags:
                keyType = str(tag.getAttribute("k"))
                if keyType == "building" or keyType == "craft":
                    isBuilding = 1
                    isClosed = 1
                elif keyType == "highway":
                    isHighway = 1
                    valueType = str(tag.getAttribute("v"))
                    if valueType == "motoway" or valueType == "trunk":
                        isFirstWay = 1
                    elif valueType == "primary" or valueType == "secondary" or valueType == "tertiary":
                        isSecondWay = 1
                    elif valueType == "unclassified" or valueType == "residential" or valueType == "service":
                        isThirdWay = 1
                    else:
                        isLastWay =1


            if isBuilding == 1 or isHighway ==1:
                #find every nodes in each way
                for node in nodes:
                    ref = int(node.getAttribute("ref"))
    
                    #test wether this point could be found in the point list we created
                    #if we can find it, get the index of this point in this list,
                    #we will use this index to get the actual point
                    try:
                        index = mapIDs.index(ref)
                    except:
                        index = -1
    
                    if(index != -1):
                        point = mapPoints[index]
                        wayPoints.append(point)

                        
                #create the polygon using the way point list
                poly = geo.createPolygon()
                
                for point in wayPoints:
                    poly.addVertex(point)
                poly.setIsClosed(isClosed)

                #type : 0 - building ; 1 - first_way ; 2 - second_way ; 3 - last_way
                if isBuilding :
                    poly.setAttribValue(polyAttrib, 0)
                elif isHighway :
                    if isFirstWay:
                        poly.setAttribValue(polyAttrib, 1)
                        streetWidth_temp = streetWithScale * boundary[0] * 1
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    elif isSecondWay: 
                        poly.setAttribValue(polyAttrib, 2) 
                        streetWidth_temp = streetWithScale * boundary[0] * 0.6
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    elif isThirdWay: 
                        poly.setAttribValue(polyAttrib, 3) 
                        streetWidth_temp = streetWithScale * boundary[0] * 0.3
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    else: 
                        poly.setAttribValue(polyAttrib, 4)
                        streetWidth_temp = streetWithScale * boundary[0] * 0.2
                        poly.setAttribValue(streetWidth, streetWidth_temp)


    boundary = getBoundary(map, mapScale)

    mapNodes = map.getElementsByTagName("node")
    mapWays  = map.getElementsByTagName("way")

    for mapNode in mapNodes:
        drawNodes(mapNode)

    createWays(mapWays)

mapCreate()

 

3：根據道路和現有建築布局，用點填充其他空余空間

將地圖數據提取到houdini之后，發現還有很多道路之間本來有建築的地方是空的，主要是因為該地圖的數據對於建築還是遠遠不夠，尤其是美國以外的國家，基本上除了大條的路其他的信息都非常不全。免費的質量就是這樣沒辦法。

接下來就是在Houdini里面的處理了。

整體結構：

根據線生成街道：

生成原理街道的建築點：

刪除距離小於一定閾值的點：

根據每個點相互間的距離以及其到街道的距離判定每個建築的最大寬度：

4：程序化生產成本每個不一樣的建築（把整個已經打成了一個包）：

這個自己做的節點生成的建築其實還是有很多問題的，沒有做更深入的控制了，先就做成這樣了。

基本上思路就是這樣了，中間遇到不不少困難，第一個就是在houdini里面Python寫起來一點都沒感覺，主要是因為用Python調houdini實例的方法還不熟悉，剛接觸難免會這樣，第二個是整個過程中每一個都步奏都是自己接觸的少的，每一個環節在開始做之前需要不停的查找文獻，尤其是在siggraph上發表的一些學術文章，現在算是真正體會了別人technical director平時都是在干嘛了，尋找發現方法，用藝術家的角度將方法實現，並確保別人也能夠使用。我在這方面還是差太遠了，總是覺得怎么做都還是不夠，這種不夠不是說沒做好還是怎樣，只是像一同水，不論怎么加水都還是看不見滿，也許明年去base學習后會找到答案吧。

這是回國實習前最后做的一個東西了，下一步就是回國適應工作環境，開始實戰。