Project Plateauと3D(Blender)、数理最適化を組み合わせて安心・安全な街を設計する

0．初めに

消費電力が最適化されたまち

0.1 改版履歴

v 1.0 4/29 リリース
v 1.1 4/30 誤記修正、画像追加、表現修正、タイトル変更 など
v 1.2 5/7 複数の光源の強さを反映した結果画像に置き換え
　　　　　　*プログラムは複雑化してわかりにくいので、更新していません。
v 1.3 5/8 VR編 可視化を追加

0.2 この記事を作成した背景

先日、素晴らしいページ（サービス）がリリースされました。
それは「Project Plateau」！

多くの人がUE4やBlenderに取り込みレンダリングなどされて、Twitterに投稿され話題となりました。
半面、活用となるとなかなか難しいといった感じですかね。

Aboutから引用

PLATEAU は、国土交通省が進める 3D都市モデル整備・活用・オープンデータ化 のリーディングプロジェクトである。都市活動のプラットフォームデータとして 3D都市モデルを整備し、 そのユースケースを創出。さらにこれをオープンデータとして公開することで、誰もが自由に都市のデータを引き出し、活用できるようになる。

ACRiもそうですが、大きな組織を動かして成果を出す方々に心からの敬意をはらい、最大限の感謝を申し上げずにはいられません！
そんな大変なこと、私は絶対やりたくない、そもそもできない…

そして、最大の恩返しは、活用することだと思いましたので、この記事を作りました。

0.3 本記事の全体的な流れ

何ができるかなぁと想像して、私の得意分野を生かして遊ぶことにしました。
ここには全体の流れを書いておきますので、興味がある章だけ見るでもよいでしょう。

この記事の全体の処理フローです。

1. Project Platauのデータを取り込み
2. ライティングをシミュレーションする（デジタルツイン）
3. 組合せ最適化（整数計画）でもっともよい配置を割り出す
4. 可視化

この記事で説明することは、大手企業などの資料風に言うと「DX(大規模都市データ）、デジタルツインを活用し、もっとも効率よく、安全、安心な”まち”を実現します！」です。

具体的には「（費用的に、電力的に）最も効率よく暗い部分をなくしたまち」を実現します。
徐々にLEDに切り替わっているとはいえ、日本の都会では数mごとに設置されている街灯。
これ、田舎やアジアの都市部に行ったことがある人なら共感いただけると思いますが、実は非常にありがたいことなんです！
どこに行ってもこれだけ明るいというのは非常に、安全にも役立っているのです。

*私はいい年のオッサンですが、出張で訪れた安全なまちでも日が落ちると薄暗く、不安を感じました。

Image from フリー素材ぱくたそ（www.pakutaso.com）

街を明るくする街灯の設置を最高の効率で実施できたら大きなコストダウンになり、もっと使いたいところにお金を回せますよね！
しかし、設置する場所によっては建物の陰になったりして道を十分に照らすことができず、街灯の明るさや数だけでは単純には決めることができません。
たとえば、下図、緑は建物です。建物のすぐ横に街灯を設置すると、さえぎられた部分が影になり、近くでも非常に暗い部分ができてしまいます。

ということで、3D街情報（Project Plateau）をもとにライティングをシミュレーション（デジタルツイン）する、街灯の配置を組み合わせ最適化（整数計画）で解決します。

最初に何が出来上がるか、結果を上げておきます。
ゴールデンウィークの間、プログラムをむにゃむにゃして、50-150WのAreaLightを使うようにしました。

結果の図（緑は建物）

1. Project Platauのデータを取り込み

1.1. データ入手元

なぜか、データですがなぜか「Project Plateau」ページから直接見つけることができません。「Project Plateau fbx download」などのキーワードで見つけることができます。

どれか一つのLINKからデータをダウンロードしてくるとよいと思います。
私は適当に選んだデータが「533956.zip」　足立区でした。

1.2.FBXデータの展開

先ほどダウンロードしたFBXのzipファイルですが、展開するといくつかの圧縮データが格納されています。

これも展開しましょう。

LOD1.zipは足立区をさらに細かく区切った区画ごとのデータ（FBX)がを含んでいます。
なお、LODは、Level of Detailを意味しており後ろの数字が細かくなるほど、より詳細な建物の形状が格納されています。

dem.zipは一つのFBXを含みます。これは地面のデータになります。

足立区は非常にシンプルな構造で以下の通りでした。データによっては他にbirdやLOD2などを含む場合があります。
*詳細はマニュアルをご参照ください。

533956.zip/
　├ dem/  <地面モデル>
　└ LOD1/ <LOD1モデル>

1.3. FBXの結合

さて、FBXを何らかのツールに取り込めばいいのですが、Objectがバラバラなためか、非常に非常に重いファイルになっています。
困ったときは先人の知恵を借りるということで、以下のページを参考にマージします。

Project PLATEAU の LOD1 ファイルを結合する

LOD1をすべて結合しても重くなるため、私は「53395612_bldg_6677.fbx」「533956_dem_6677.fbx」の2つだけを別フォルダにコピーし、結合しました。

1.4. Blenderでの取り込み

Blenderでの取り込みは非常に簡単です。「File」→ 「Import」→「FBX(.fbx)」を選択して、先ほど結合したファイルを取り込みます。
Blenderの説明はメインではないので細かい話は飛ばしますが、注意点をいくつか。

• 原点から遠い部分に取り込まれる場合があります。あわてずZoomOutしましょう。また、ViewのClipEndなどでクリッピングされているかもしれませんので、どうしても見えない場合はこのサイズを大きくしましょう。
• デフォルトでは1/100のようです。位置を移動してから100倍するなど現実のサイズに合わせるとよいでしょう。

Blenderで取り込んだら、不要な地面を削除し位置を移動するなど処理しやすいようにしましょう。
グイっと移動して、背景にHDRIを設定したのがこちら。
かなり、やり遂げた気分になります（笑）。

2. ライティングをシミュレーションする（デジタルツイン）

さて、先ほど取り込んだFBX適当な住宅が並んでいる部分を探します。

適当に空いている空間があって、家が立ち並んでいるような場所を探します。
私は以下のような部分を選びました。

2.1. カメラ、ライト、レンダラの設定

2.1.1. カメラの設定

まずは、真上からレンダリングできるようにカメラを設定しましょう。今回初めて知りましたが、Previewだけでなく、カメラでも「OrthoGraphic」が選べます。デフォルトの「perspective」（透視図法かな？）だと後続の分析で難しくなるので、「OrthoGraphic」を使いましょう。

真上から、おおよそ100メートル四方が入るようにカメラ設定をしました。

2.1.2. 街灯の設定

さて、実は歩道に求められる明るさや、街灯というのは規格が決まっているようです。
参考：
https://www.mlit.go.jp/road/sign/data/chap10.pdf

そのあたり、本番では詰める必要がありますが、今回のお遊びで確認するだけであればそこまでは不要でしょう。
とりあえず、250Wの点光源を街灯と見立てましょう。

これを、街灯の高さ4.5mに設置して、適当な場所に配置したものが以下の図です。

環境光を消し街灯だけにします。
実際にはこの画像を使用します。

雰囲気出ましたね。

2.1.3. レンダラの設定

反射光も考慮して、大好きなCyclesレンダラを用いていましたが、途中からEeveeに変更しました。このレベルと内容では速さが正義です。(1枚当たり数分から1秒に画像作成が短縮されました）
後ほど気が付きましたが、Cyclesじゃないと目的が達成できないので、現在の3.0αで実装されているCyclesXでバンバン画像を生成しました。

2.2. 街灯を移動する。

さて、あとはこの点光源をあちこちに配置したときの画像を生成します。
Blenderが素晴らしいところは、Pythonでこれらの制御ができるところです。
私は以下のようなPythonで、5m刻み計400パターン生成しました。

LightMovement

import bpy
import os
import numpy as np

filepath = bpy.data.filepath
directory = os.path.dirname(filepath)
lamp_object = [o for o in bpy.data.objects if o.name == 'Light']

for i in range(-50,50,5):
    for j in range(-50,50,5):
        lamp_object[0].location = (i,j,4.5)
        bpy.context.scene.render.filepath = directory+f"{50+i:02}_{50+j:02}.png"
        bpy.ops.render.render(write_still = True)
        bpy.ops.render.render()

一枚一枚レンダリングするので、多少は時間はかかりますが待っておけばいいだけなので、焦らずほっておくと以下のような画像を得ることができます。

3．組合せ最適化（整数計画）でもっともよい配置を割り出す

BlenderのPythonで一気に処理してしまってもいいのですが、後続の処理はVSCodeで実装しています（すいません、なれの問題です）。

3.1. 最適化の考え方

最適な街灯の配置が得る戦略を以下のように考えました。
1. 100m四方を細かく分割(街灯の移動距離のさらに半分で分割しました)
2. 各街灯がある場合、ない場合を表す変数xを用意する。xの添え字は街灯の位置とする。
3. 消費電力的、投資的に費用は街灯は少ないほうが良い。
4. ある分割した1つのエリアを見たとき、光を提供する街灯のxを合計する。光が届くということはそれが1以上となる

ここまで準備すると定式化ができます。

3.1.1. 定式化

数学っぽい表記にしていますが、表現しているのは直前で書いたことと同じです。

目的関数
目的関数は、なるべく街灯が少なくする（効率よくすること）こと。

$$ minimize \sum_{i \in LIGHTAREA} x_i$$

*ただし、ここでLIGHTAREAは街灯を試しに設置したすべての場所

制約式
全ての分割したエリアにおいて、いずれかの街灯の光が届いていること。

$$\sum_{k \in LIGHTS_j} x_k >=1 , (\forall j\in DIVIDEDAREA) $$

*ただし、ここでDIVIDEDAREAは分割したすべてのエリア、$LIGHTS_j$はエリアjに光を届けることができる街灯の位置の集合

3.2. 地面はどこだ？ 地面を決定する処理

地面だけを対象としたいので、ここだけはSunLightを追加して、建物と地面を分離します。
後々マスクとして使用するため、2値化してます。

from PIL import Image
import matplotlib.pyplot as plt
import numpy as np
import glob

img_ground = Image.open('data/ground.png')
#地面は二値化して反転
img_ground = img_ground.point(lambda x:255 if x < 10 else 0).convert("1")
plt.imshow(np.array(img_ground),cmap="binary")

3.3. レンダリング画像の読み込み

レンダリングした画像、合計400枚を読み込みます。後ほど連続で処理できるようにPillowのImage型に変換して、街灯の位置を添え字とした辞書型に変換します。

# 画像データを順次読みこむ
files = glob.glob("./data/250w/*.png")
sim_images = {}
for reading_file in files:
    file_name = reading_file.split("/")[3]
    if (file_name == "ground.png"):
        continue
    #ファイル名が00_00.pngであること前提
    x = int(int(file_name[0:2]))
    y = int(int(file_name[3:5]))
    #グレースケールにして保持
    sim_images[(x,y)] = Image.open(reading_file).convert("L")
#Lightを設置した位置パターン
PATTERN_X = list(set(np.array(list(sim_images.keys()))[:,0]))
PATTERN_Y = list(set(np.array(list(sim_images.keys()))[:,1]))
LIGHT_X_NUM = len(PATTERN_X)
LIGHT_Y_NUM = len(PATTERN_Y)


#ライトの移動距離の1/2のサイズで画像を分割する
SPLIT_X_NUM = len(PATTERN_X) * 2
SPLIT_X = img_ground.size[0] / SPLIT_X_NUM
SPLIT_Y_NUM =  len(PATTERN_Y) * 2
SPLIT_Y = img_ground.size[1] / SPLIT_Y_NUM

3.4. 画像ごとに各エリアに光が届いているかの配列を作成する

渡された画像を決められたサイズで分割し、それぞれのエリアに光が届いていれば１，届いていなければ0を格納する2次元配列を返す関数を定義します。

#地面のうち、分割されたエリアの特定の割が明るければ光が届いているという扱いにする

#明るさがBRIGHT_THRESHOLD以上であれば明るいとカウント
BRIGHT_THRESHOLD = 5
#面積のCOUNT_RATIO以上明るいと判定されると、その区画はOKとする。
COUNT_RATIO = 0.7
#分割エリアの3割を地面の面積が超えない場合、対象外とする
IGNORE_RATIO = 0.3


def check_brightness_by_area(img):
    #戻り値変数を確保
    ret = np.zeros((SPLIT_Y_NUM,SPLIT_X_NUM),dtype=int)

    for y in range(SPLIT_Y_NUM):
        for x in range(SPLIT_X_NUM):
            cropped_area = (SPLIT_X*x,SPLIT_Y*y,SPLIT_X*(x+1)-1,SPLIT_Y*(y+1)-1)

            mask_img = img_ground.crop(cropped_area)
            #地面（2値画像のTrue）のPixel数
            road_points = np.sum(np.array(mask_img))
            #地面が存在しないエリアはスキップ
            if road_points < 0.3*SPLIT_X*SPLIT_Y:
                continue

            #レンダリング画像を分割したもの
            cropped_rendered_img = img.crop(cropped_area)

            #レンダリング画像から地面だけを抜き出す
            target = np.array(cropped_rendered_img) * np.array(mask_img)

            #明るさ40以上をカウント
            brighted_area = np.sum(target > BRIGHT_THRESHOLD)
            #print(brighted_area,road_points,brighted_area/road_points)
            if (brighted_area/road_points > COUNT_RATIO):
                #print(y,x,brighted_area,road_points)
                ret[y][x] = 1
    return ret

3.5. 最適化準備

このレベルであれば、定式化した通り整数計画で定義、解くことができる。
Pulpの流儀に従って問題定義、変数定義、目的関数定義を定義する。ただし、どうあがいても光が届かない場所が出る可能性があるので、そのような場所を排除するために、エリアに一度でも光が届いたかを確保するMAPも作成する。

from pulp import LpProblem,lpSum,LpStatus,LpBinary,LpMinimize,LpVariable

#最も少ないライトで街を照らすという設定にする
p = LpProblem(name = "SafetyCity",sense=LpMinimize)

#(y,x)のライトを使うかどうかを、0,1で表現する変数
x = LpVariable.dict(name="x",indexs=(PATTERN_Y,PATTERN_X),lowBound=0,upBound=1,cat=LpBinary)

#街灯の数最小化（目的関数）
p += lpSum(x[i] for i in sim_images.keys())

#全ての街灯の位置において、細分化されたエリアの明るさMatrixを作成する。
#同時に、細分化エリアに光が届く可能性があるかを別でマッピングする（どうあがいても光が届かない場所があると問題として解けないため）

#最適化のターゲットエリアを決める
target_area = np.zeros(shape=(SPLIT_Y_NUM,SPLIT_X_NUM),dtype=int)

#各ライトの位置で明るくなるエリア
light_map = {}
for light_location in sim_images.keys():
    light_map[light_location] = check_brightness_by_area(sim_images[light_location])

    #これまでのLightの位置 、今回のLightの位置の or　で１になるエリアをターゲットとする
    target_area = target_area | light_map[light_location]

#全ての位置が明るくなる制約式
for area_y in range(SPLIT_Y_NUM):
    for area_x in range(SPLIT_X_NUM):

        #もし、最適化の対象エリアでなければ条件式に入れない
        if target_area[area_y][area_x] == 0:
            continue

        #エリア(y、x) を照らすライトを足し合わせる
        tmp_sum = 0
        for light_location in sim_images.keys():
            if light_map[light_location][area_y][area_x] == 1:
                tmp_sum += x[light_location]

        #いずれかのライトで明るくなれば良い
        p += tmp_sum >= 1

3.6. 最適化

あとは最適化を解く(solve関数を実行する)だけ。

p.solve()
cnt = 0
best_locations = []
for light_location in sim_images.keys():
    if x[light_location].value() == 1:
        best_locations.append(light_location)
        cnt = cnt + 1
    #print(light_location,x[light_location].value()) 
print(cnt)
print(best_locations)

最後、かなり手抜きだが、最適化した街灯の位置を出力する。

出力結果

54
[(30, 65), (5, 75), (0, 15), (20, 90), (15, 65), (0, 60), (70, 5), (25, 25), (70, 0), (70, 80), (10, 50), (85, 70), (70, 30), (15, 20), (60, 85), (40, 50), (60, 50), (85, 40), (55, 10), (90, 50), (80, 10), (20, 15), (50, 90), (90, 25), (0, 40), (20, 5), (20, 80), (45, 70), (10, 95), (50, 40), (40, 20), (85, 45), (80, 95), (90, 10), (30, 45), (50, 55), (10, 80), (5, 0), (10, 30), (55, 50), (60, 65), (60, 25), (30, 5), (95, 70), (30, 90), (95, 5), (95, 85), (40, 30), (75, 55), (15, 40), (75, 60), (5, 20), (65, 95), (40, 10)]

4. 可視化

4.1 Blenderでの可視化

街灯の位置（数字）は分かったが、どうなっているかちゃんと結果を見たいところ。
Blenderに戻って出力されたListをコピペして、PointLightを追加するPythonを作成します。
（すいません、このあたりやっつけです）

街灯を追加する処理

import bpy

light_locations = [(30, 65), (5, 75), (0, 15), (20, 90), (15, 65), (0, 60), (70, 5), (25, 25), (70, 0), (70, 80), (10, 50), (85, 70), (70, 30), (15, 20), (60, 85), (40, 50), (60, 50), (85, 40), (55, 10), (90, 50), (80, 10), (20, 15), (50, 90), (90, 25), (0, 40), (20, 5), (20, 80), (45, 70), (10, 95), (50, 40), (40, 20), (85, 45), (80, 95), (90, 10), (30, 45), (50, 55), (10, 80), (5, 0), (10, 30), (55, 50), (60, 65), (60, 25), (30, 5), (95, 70), (30, 90), (95, 5), (95, 85), (40, 30), (75, 55), (15, 40), (75, 60), (5, 20), (65, 95), (40, 10)]

for light_location in light_locations:
    x = light_location[0]
    y = light_location[1]

    bpy.ops.object.light_add(type="POINT", radius=0.4, align='WORLD', location=(x-50, y-50, 4.5), scale=(1, 1, 1))
    bpy.context.object.data.energy = 250

「うーん、そこに必要かなぁ」というところもありますが、ばらけているようには見えます。

さて、レンダリングしましょう。！

お、いい感じです。かなり細かい路地まで光が届いていますね。
最初の図ですが、建物を緑にして道だけわかりやすくします。

ループ回数を間違って街灯が届かなかった右上は光が届いていない路地ができてしまいました。実装を確認いただけるとわかりますが、意味がなさそうな地面があまりに少ないエリアなども切り捨てているため、細い路地(?)は暗い部分がありそうです。

もうすこし街をあかるくしないといけないのでは？ というときは、BRIGHT_THRESHOLDの値を上げてください。街灯は増えますが全体的に明るくなります。

4.2 VR可視化

では本当に人から見るとどのように見えるのでしょうか？

こういう時に便利なものがVRです。
詳細をかなり飛ばしますが、UE4のブループリント内で動的にライトを作成するのは今回はお勧めしません。
非常に重いうえに、このように固定のライトはライティングをビルドしたほうが見た目がきれいになるからです。
*今回のLOD1レベルでは何でもよいかもしれませんが。

もうひと頑張りPythonでライトをWorldに追加しましょう。そのために、まずPython Script関連を有効にします。

次にPythonScriptでRectLightを追加します。
*下記の結果ではLightの強さも含めて3つの配列に位置が分かれていますが、UE4でIntensityの設定方法がわからなかったため、デフォルトのIntensityで追加しています。

EditorScript

import unreal

light_locations = []
light_locations.append([[85.0, 35.0], [70.0, 40.0], [85.0, 40.0], [95.0, 60.0], [10.0, 30.0], [5.0, 20.0], [70.0, 80.0], [10.0, 80.0], [0.0, 25.0], [50.0, 95.0], [75.0, 45.0], [40.0, 45.0], [30.0, 70.0], [35.0, 15.0], [40.0, 85.0], [70.0, 50.0], [95.0, 100.0], [20.0, 75.0], [65.0, 65.0], [15.0, 20.0], [50.0, 30.0], [0.0, 30.0], [85.0, 90.0], [70.0, 5.0], [65.0, 35.0], [75.0, 15.0], [20.0, 65.0], [30.0, 45.0], [35.0, 0.0], [15.0, 95.0], [100.0, 35.0], [70.0, 75.0], [40.0, 30.0], [0.0, 60.0], [15.0, 90.0], [80.0, 90.0], [0.0, 10.0]])
light_locations.append([[30.0, 60.0], [85.0, 70.0], [75.0, 60.0], [40.0, 10.0], [45.0, 70.0], [75.0, 30.0], [15.0, 40.0], [60.0, 5.0], [65.0, 25.0], [55.0, 10.0], [95.0, 10.0], [40.0, 20.0], [60.0, 85.0], [80.0, 10.0], [65.0, 95.0], [85.0, 50.0], [100.0, 70.0], [60.0, 15.0], [25.0, 15.0], [60.0, 50.0], [25.0, 25.0], [55.0, 40.0], [45.0, 55.0], [25.0, 35.0], [95.0, 50.0], [20.0, 5.0], [95.0, 25.0], [0.0, 45.0], [0.0, 0.0]])
light_locations.append([[100.0, 80.0], [25.0, 90.0], [5.0, 0.0], [90.0, 80.0], [60.0, 95.0], [10.0, 65.0], [30.0, 95.0], [65.0, 85.0], [55.0, 70.0]])

cnt = 0
for i in range(len(light_locations)):
    for light_location in light_locations[i]:
        x = light_location[0]
        y = light_location[1]

        sun_location = unreal.Vector(x*100-5000, 5000-y*100, 450.0)
        lightObj = unreal.EditorLevelLibrary.spawn_actor_from_class(unreal.RectLight, sun_location, rotation=[0,0,0])
        lightObj.set_actor_label(f"light{cnt}")
        cnt = cnt + 1

注意点は座標系とScaleがBlenderとUE4では異なることでしょうか。上記Scriptではそのあたりは適当に変換しています。
BlenderからFBX書き出し、UE4への読み込み、その他もろもろ実施するとこんな感じ。

で、VRかぶって人を適当に立たせてみてみると…。
Oculusで見ると数m離れていてなんとなく顔が見えるくらいなので、セーフかな。

5. 今後の発展

実は以下のところまで実施してから公開しようと思ましたが、サクッとやって出したほうがいいなと判断して、このレベルで公開しています。
反応あれば続きを実施するかもしれません。

• ◆済 ワット数の異なる街灯を準備してより現実的な問題にする
• ◆済 近接の複数の街灯の影響を顧慮する（複数の街灯の影響でThresholdを超えるエリアが出るかも）
• ◆済 レンダリング設定をもう少し詰める。　→ Cycles必須
• もっと広い範囲を最適化する。
• ◆済 VRで結果見る（そもそもBlenderのUV展開で時間食ってる）

6. 最後に

使ってこそのデータです。
「Project Plateauの裏には、きっと燃えるような情熱をもち血を吐くような苦労・調整をされた方がいらっしゃるのかな。」と想像します。

それと、ほかにも色々応用が考えられますので、ぜひ一緒にやりたいという、都道府県 市町村の職員の方、または面白そうな活動を指定らしゃる方、声かけてください！！
（個人的にゆかりのある、泉佐野とか熊本とか川崎とか福島とかインドとかフィリピンとかだとさらにうれしいです）

あまりにデータなどが汚かったため、Blender、iPython Notebookは現時点で未公開です。
もし、LGTMが多いか、コメントが多いようでしたら、整理して公開しようと思います。