粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

① ikawidget2のデータを使ってSplatoon2での編成を評価 ~1.データ取得と集計編~

Abstract

発売から2年半経っていますがSplatoon2を未だにやってます。ここまで長いことやったゲームは初めてでつまり神ゲーってことです。

今回はSplatoon2の編成強い弱い問題を最近遊んだ機械学習ライブラリで解き明かせないかという企画になります。(なお著者はCS専門ではないです。)

  • 自分の戦績データから,「勝敗に関連しやすそうな要素」を抽出する
  • 編成などから「勝敗予想」モデルを作成する

というのを目標に据えてやっていきます。

記事

①ではデータの読み出しとデータ構造の解析に焦点を当てます。とりあえずデータのとり方など技術的な話をするので興味ない方は飛ばして結構です。

Splatoon 2 (スプラトゥーン2) - Switch

Splatoon 2 (スプラトゥーン2) - Switch

  • 発売日: 2017/07/21
  • メディア: Video Game

背景(モチベーション)

Splatoon2では敵味方4対4に分かれて試合を行うのですが,その時のブキの編成が結果を非常に左右しやすいという印象を受けます。

例えば,チャージャー,ローラー,ブラスターなどを2枚引くと塗りや打開が厳しいというのは誰しもが思っていることでしょう。

感覚的に厳しい編成かどうかはわかるのですがうまいことこれをデータで示せないか,というのがモチベーションです。

自分の戦績はikawidget2というスマホアプリから取得したものを用います。

参考情報ですが,2020/4/28時のデータはこの様になってます。

  • ウデマエ:エリア・ホコ・ヤグラ XP2000~2200,アサリ S+ 2000前後
  • 持ちブキ:チャージャー・スピナー以外全般
  • 戦績データ数:エリア: 826戦ヤグラ:845戦ホコ:432戦アサリ:339戦

です。

ikawidget2の戦績データを抽出する

このブログを参考にまずはデータをPythonで読み出します

  1. バックアップを行う → 本体に ~~.ikaxという名前のファイルが保存されるのでそれをPC上にダウンロード。
  2. ~~.ikaxの拡張子を.zipに変更して中身を開く → stats.realmというファイルが出てくる。
  3. Realm StudioをDLしてstats.realmを開き,json形式で保存する。

json形式の戦績データを読む

以下のようにして読み出します。

import json

with open("../data/stats.json","r",encoding="utf-8") as f:
    result_json=json.load(f)

試合結果はresult_json["Results"]にリストの形式で入っています。

試合の情報全てがこのResultsに入っており,

  • ガチパワー,ステージ等の情報
  • プレーヤーのギア・ブキや戦績情報
  • 勝敗

が辞書型で記録されています。

試しにkeysを見てみると,

result_json.keys()
>> dict_keys(['Boss', 'BossCount', 'Brand', 'CoopResult', 'CoopSchedule', 'CoopStage', 'CoopWeapon', 'EventType', 'Fes', 'Game', 'Gear', 'Grade', 'Player', 'PlayerType', 'Result', 'Skill', 'SkillLog', 'Skills', 'Special', 'Stage', 'SubWeapon', 'WaterLevel', 'WaveDetail', 'Weapon', 'Worker'])

となります。主に使う戦績データはResultに入ってます。

ステージ毎の勝率計算

参考ブログでは,戦績からfor文を用いてガチエリアの戦績を集計しています。 試しに苦手なガチヤグラ'game': `を集計してみました。


補足:

  • ルール情報は'game'にある。ヤグラはtower_controlgachi,エリアはsplat_zonesgachietc...
  • ブキ・ギア・ステージ等は番号で管理されており,IDと名前の対応はWeapon・Gear・Stageのkeys下に保管されている
yagura_result={}

for item in result_json["Result"]:
    if item["udemae"]<10 and item["game"]!="tower_controlgachi": #S以上の戦績に限定&ガチヤグラ以外のルールを排除
        continue
    for i in result_json["Stage"]:
        if i["ID"]==item["stage"]:
            stage_name=i["name"]

    # 勝ちと負けをステージごとにカウント
    if not stage_name in yagura_result.keys():
        yagura_result[stage_name]={"win":0,"lose":0}
    if item["win"]:
        yagura_result[stage_name]["win"]+=1
    else:
        yagura_result[stage_name]["lose"]+=1

# 勝率高い順にソート        
for key in sorted(yagura_result.keys(),key = lambda x:yagura_result[x]["win"]/(yagura_result[x]["win"] +yagura_result[x]["lose"]),reverse=True):
    print(key,yagura_result[key]["win"]/(yagura_result[key]["win"] +yagura_result[key]["lose"]))

結果は

デボン海洋博物館 0.6153846153846154
(中略)
モズク農園 0.41346153846153844

となりました。この時点での情報はikawidgetでも見れますね。

プレーヤー情報の取得

ここから勝ち負けを評価するためにプレーヤーの情報を抽出します。

プレーヤー情報はResults以下の myMembers(味方),otherMember(敵),player(自分自身)の3つのkeyに保存されています。使いそうなのは

  • kii / allKill: キル数 / アシスト込のキル
  • assist: アシスト数
  • death: デス数
  • paintpoint: 塗りポイント
  • special: スペシャル回数
  • weapon: ブキ

くらいでしょうか。

偽相関が出そうですが,以下の情報も使えそうです。

  • elapsed time:試合時間
  • gachiEstimatePower / gachiEstimateXPower :ガチパワー / Xパワー (Player以下のudemaeIsXの変数フラグにより切り替え)

味方ブキ毎の勝率計算

例としてステージ毎の集計と似たような感じで味方のブキ毎の勝率を計算することができます。

#エリアで勝ったときの味方のブキを集める
win_weapon_ally = {}

# プレイヤー情報の入った辞書データから勝敗を抽出する関数
def addwin(win_weapon,result_json,item,win):
    w_name = "null"
    for i in result_json["Weapon"]:
        if i["ID"]==item:
            w_name=i["name"]
    if not w_name in win_weapon.keys():
        win_weapon[w_name]={"win":0,"lose":0}
    if win:
        win_weapon[w_name]["win"] +=1
    else:
        win_weapon[w_name]["lose"] +=1

# 試合結果から結果を抜き出すループ
for item in result_json["Result"]:
    #if item["udemae"]<10 and item["game"]!="tower_controlgachi": #S以上の戦績に限定&ガチヤグラ以外のルールを排除
    if item["udemae"]<10 and item["game"]!="splat_zonesgachi": #S以上の戦績に限定&ガチヤグラ以外のルールを排除
        continue
    for ally in item["myMembers"]:
        if item["win"]:
            addwin(win_weapon_ally,result_json,ally["weapon"],1)
        else:
            addwin(win_weapon_ally,result_json,ally["weapon"],0)

結果がこれ。(長文注意)

出会ったことの少ないマイナーなブキが上位に来てますが50戦以上データがあるなかで勝率が良かったのはなんとボールドマーカー7(約62%)でした。 最下位はソイチューバー2種。かわいそう...

H3リールガンチェリー 1.0 2
ヒーロースピナー レプリカ 0.8888888888888888 9
14式竹筒銃・丙 0.8181818181818182 11
スプラローラーコラボ 0.7619047619047619 21
スパイガジェット 0.7368421052631579 19
スパイガジェットベッチュー 0.7368421052631579 19
スクリュースロッシャーネオ 0.7 10
プロモデラーRG 0.6666666666666666 12
スプラスコープベッチュー 0.6666666666666666 21
ノーチラス79 0.6666666666666666 15
ノヴァブラスターネオ 0.6666666666666666 3
ロングブラスターカスタム 0.6666666666666666 3
スプラマニューバーベッチュー 0.6410256410256411 39
ホクサイ 0.6363636363636364 22
ホットブラスターカスタム 0.6285714285714286 35
L3リールガンベッチュー 0.625 8
リッター4Kカスタム 0.625 8
ジェットスイーパーカスタム 0.6222222222222222 45
ボールドマーカー7 0.6197183098591549 71
ホクサイベッチュー 0.6176470588235294 34
カーボンローラー 0.6153846153846154 13
プライムシューター 0.6111111111111112 36
スプラマニューバー 0.6086956521739131 23
プロモデラーPG 0.6071428571428571 28
.52ガロン 0.6 10
ヴァリアブルローラー 0.6 5
パーマネント・パブロ 0.6 10
クラッシュブラスターネオ 0.5897435897435898 78
ヒッセン・ヒュー 0.5882352941176471 51
ジェットスイーパー 0.5882352941176471 68
クーゲルシュライバー 0.5882352941176471 17
ケルビン525デコ 0.5869565217391305 46
ヒッセン 0.5853658536585366 41
ヒーローローラー レプリカ 0.5842696629213483 89
ヒーローチャージャー レプリカ 0.5833333333333334 60
ヒーローブラシ レプリカ 0.5833333333333334 12
バケットスロッシャーソーダ 0.5833333333333334 24
ダイナモローラーベッチュー 0.5789473684210527 133
ボールドマーカーネオ 0.5740740740740741 54
ノヴァブラスター 0.5737704918032787 61
スプラスコープコラボ 0.5714285714285714 28
Rブラスターエリートデコ 0.5714285714285714 7
デュアルスイーパー 0.5645161290322581 62
L3リールガン 0.5625 32
スパッタリークリア 0.5609756097560976 123
スプラローラー 0.5598290598290598 234
.96ガロンデコ 0.5555555555555556 36
.96ガロン 0.5555555555555556 27
バレルスピナーデコ 0.5555555555555556 27
ケルビン525ベッチュー 0.5555555555555556 9
ロングブラスター 0.55 60
ボトルガイザーフォイル 0.55 40
プライムシューターコラボ 0.5454545454545454 44
ヒーローマニューバー レプリカ 0.5454545454545454 11
ハイドラントカスタム 0.5443786982248521 169
オーバーフロッシャーデコ 0.5384615384615384 39
H3リールガンD 0.5384615384615384 26
スプラスピナーベッチュー 0.5384615384615384 13
リッター4K 0.5368421052631579 95
オクタシューター レプリカ 0.5344827586206896 58
もみじシューター 0.5333333333333333 75
デュアルスイーパーカスタム 0.5329512893982808 349
オーバーフロッシャー 0.5307692307692308 130
スプラスピナー 0.5277777777777778 36
パラシェルターソレーラ 0.5263157894736842 19
4Kスコープ 0.5263157894736842 95
14式竹筒銃・甲 0.5254237288135594 59
スクイックリンα 0.5217391304347826 23
バレルスピナーリミックス 0.5190839694656488 131
スプラシューターコラボ 0.5181818181818182 220
N-ZAP89 0.5163398692810458 153
プライムシューターベッチュー 0.5150684931506849 365
ボールドマーカー 0.5142857142857142 140
クアッドホッパーブラック 0.5135135135135135 74
ホクサイ・ヒュー 0.5128205128205128 39
N-ZAP85 0.5114285714285715 350
スプラシューターベッチュー 0.5032258064516129 155
わかばシューター 0.5023041474654378 217
ラピッドブラスターベッチュー 0.5 70
ノーチラス47 0.5 24
ハイドラント 0.5 20
スパイガジェットソレーラ 0.5 6
ヴァリアブルローラーフォイル 0.5 48
バケットスロッシャー 0.5 8
スクリュースロッシャー 0.5 6
ヒーロースロッシャー レプリカ 0.5 2
ラピッドブラスターデコ 0.5 14
ケルビン525 0.5 6
スプラスピナーコラボ 0.5 10
シャープマーカーネオ 0.49382716049382713 81
パブロ 0.49019607843137253 51
キャンピングシェルターカーモ 0.4864864864864865 37
L3リールガンD 0.4827586206896552 58
スプラチャージャー 0.4803921568627451 102
クーゲルシュライバー・ヒュー 0.47706422018348627 109
エクスプロッシャーカスタム 0.47297297297297297 74
スパッタリー・ヒュー 0.4722222222222222 72
スプラシューター 0.4714285714285714 70
ラピッドブラスター 0.46153846153846156 13
Rブラスターエリート 0.46153846153846156 13
スプラマニューバーコラボ 0.4578313253012048 83
スプラチャージャーコラボ 0.45454545454545453 44
クラッシュブラスター 0.45454545454545453 22
ヒーローブラスター レプリカ 0.45454545454545453 11
シャープマーカー 0.45454545454545453 33
スクリュースロッシャーベッチュー 0.45098039215686275 102
ロングブラスターネクロ 0.45098039215686275 51
.52ガロンベッチュー 0.4482758620689655 29
おちばシューター 0.4444444444444444 63
スプラローラーベッチュー 0.44 25
カーボンローラーデコ 0.44 50
ダイナモローラーテスラ 0.4318181818181818 44
.52ガロンデコ 0.42857142857142855 7
スクイックリンβ 0.42857142857142855 21
スプラチャージャーベッチュー 0.42857142857142855 35
クアッドホッパーホワイト 0.42857142857142855 7
ダイナモローラー 0.42857142857142855 14
パブロ・ヒュー 0.425 40
ホットブラスター 0.4166666666666667 24
スプラスコープ 0.4084507042253521 71
H3リールガン 0.4 5
4Kスコープカスタム 0.4 5
スパッタリー 0.39285714285714285 28
ヒーローシューター レプリカ 0.38461538461538464 13
プロモデラーMG 0.3793103448275862 29
パラシェルター 0.375 24
スクイックリンγ 0.375 16
N-ZAP83 0.36 25
ノヴァブラスターベッチュー 0.34782608695652173 46
エクスプロッシャー 0.3333333333333333 36
バケットスロッシャーデコ 0.3333333333333333 21
ボトルガイザー 0.3333333333333333 3
14式竹筒銃・乙 0.3333333333333333 9
バレルスピナー 0.29411764705882354 17
キャンピングシェルターソレーラ 0.2857142857142857 7
ヒーローシェルター レプリカ 0.2 10
キャンピングシェルター 0.14285714285714285 7
ソイチューバー 0.14285714285714285 7
ソイチューバーカスタム 0.0 2

この結果は私の持ちブキや立ち回りとの相性という話なので一般的ではないですが,傾向を見る上では参考になりそうです。

より大規模にはika.statsというサイトがあります。

結構ちゃんと分析されてるようなので私の出る幕があるかどうかは微妙です。まぁしばらくは車輪を再発明していきましょう。

github.com

解析に使うデータを整形する

今回の生データはjsonで記述されていますが,これは今後の解析を考える上ではあまり使いやすいものとは言えません。

Kaggleで遊んだときのようなcsvデータのような形式で試合結果を保存するのが良いでしょう。

ossyaritoori.hatenablog.com

行列形式で1行に試合結果を入れるとすると,ひとまず列に必要そうな要素としては

  • ステージ
  • 自分のブキ
  • 自分のキル・デス
  • 自チームの合計キル・デス , 相手チームの合計キル・デス
  • 自チームの合計SP回数, 相手チームの合計SP回数
  • 自チームの合計塗りポイント,相手チームの合計塗りポイント(以後”相手~”の部分省略)
  • 自チームの武器ごとの編成数, 相手~
  • 勝ち負け

あたりでしょうか。 Resultsからこれを抽出する関数を書きましょう。

使用するのはPandasのDataFrameとします。

結果からスコアを抜き出すテスト

試しにいくつかの要素をピックアップしてpandasのdataframeに変換してみましょう。ブキの分類はちょっと考察がいるので後回しにすることにします。

import pandas as pd
import numpy as np

# 一時保存用
dfarray = []

# For ループを回してデータをArray状にする
for item in result_json["Result"]:
    # 1. game information
    game = item["game"]
    # 2. gachi power if X or not
    gachipower = item["gachiEstimateXPower"] if item["udemaeIsX"] else item["gachiEstimatePower"] 
    # 3. extract elapsed time
    battletime = item["elapsedTime"]
    # 4. stage
    stage = item["stage"]
    # 5, Mypower
    mypower = item["xPower"]
    # 6. My data
    mykill = item["player"]["kill"]
    mydeath = item["player"]["death"]    
    mysp = item["player"]["special"]
    myweapon = item["player"]["weapon"]
    mypaintpt = item["player"]["paintPoint"]
    # 7. Our data
    ourkill = mykill
    ourdeath = mydeath
    oursp = mysp
    ourweapon = [myweapon]
    ourpaintpt = mypaintpt
    for ally in item["myMembers"]:
        ourkill += ally["kill"]
        ourdeath += ally["death"]    
        oursp += ally["special"]
        ourweapon.append(ally["weapon"])
        ourpaintpt += ally["paintPoint"]
    # 8. Enemy data
    theirkill = 0
    theirdeath = 0
    theirsp = 0
    theirweapon = []
    theirpaintpt = 0
    for ene in item["otherMembers"]:
        theirkill += ene["kill"]
        theirdeath += ene["death"]    
        theirsp += ene["special"]
        theirweapon.append(ene["weapon"])
        theirpaintpt += ene["paintPoint"]
    # 9. Results
    win = 1 if item["win"] else 0
    
    dfarray.append([game,gachipower,battletime,stage,mypower,mykill,mydeath,mysp,mypaintpt,ourkill,ourdeath,oursp,ourpaintpt,theirkill,theirdeath,theirsp,theirpaintpt,win])

column_name = ["game","gachipower","battletime","stage","mypower","mykill","mydeath","mysp","mypaintpt","ourkill","ourdeath","oursp","ourpaintpt","theirkill","theirdeath","theirsp","theirpaintpt","win"]

# テスト用のデータ格納
testdata=pd.DataFrame(data=dfarray,columns=column_name)

f:id:ossyaritoori:20200429233406p:plain
抜き出したのはこんな感じ。昔のデータがあるのでS+時代のデータから始まってるのはご愛嬌。

テスト解析:sklearnの決定木でガチエリアの集計

以前やった決定木から試していきましょう。(Localの環境構築を忘れていたのでsklearnから始めます。)

ガチエリアのX帯の戦績に着目してデータを抜き出します。

areadata = testdata[testdata["game"]=="splat_zonesgachi" ]
areaXdata = areadata[areadata["mypower"]>0]
areaXdata = areaXdata.drop(["game","stage","battletime"],axis=1)

f:id:ossyaritoori:20200429235659p:plain
整形後のガチエリアX帯のデータ

さて,以前のコードを流用して簡単に分析ごっこをしてみましょう。

KaggleのTitanicデータに対してsklearnの決定木を試してみる - 粗大メモ置き場

まずはモデルを作成。

import sklearn.tree as tree
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

clf = tree.DecisionTreeClassifier(max_depth=10)

そして交差検証で学習です。

# 訓練と結果に分ける
X_train = areaXdata.drop(["win"],axis=1)
y_train = areaXdata["win"]


# 5分割交差検証を指定し、インスタンス化
from sklearn.model_selection import KFold
kf = KFold(n_splits=5, shuffle=True, random_state=0)

# スコアとモデルを格納するリスト
score_list = []
models = []

# 各分割ごとに評価
for fold_, (train_index, valid_index) in enumerate(kf.split(X_train, y_train)):    
    print(f'fold{fold_ + 1} start')
    train_x = X_train.iloc[train_index]
    valid_x = X_train.iloc[valid_index]
    train_y = y_train.iloc[train_index]
    valid_y = y_train.iloc[valid_index]
    
    ## 分割データで学習・予測・評価
    clf = tree.DecisionTreeClassifier(max_depth=10)
    model = clf.fit(train_x, train_y)
    
    # データを用いて予測,記録
    predicted = model.predict(valid_x)
    score_list.append(accuracy_score(predicted,valid_y))
    models.append(model)
print(score_list, '平均score', round(np.mean(score_list), 3))

結果は以下の通りで,85%の精度で結果を推定できました。

[0.9099099099099099, 0.8018018018018018, 0.8198198198198198, 0.8918918918918919, 0.8090909090909091] 平均score 0.847

ここでこの木における重要度,すなわちどの要素が結果を大きく左右したかについて見てみます。

import matplotlib.pyplot as plt

importances = np.zeros(model.feature_importances_.size)

for model in models:
    importances += model.feature_importances_

fig, ax = plt.subplots()
plt.grid()
ax.bar(train_x.columns,importances)
fig.autofmt_xdate() # make space for and rotate the x-axis tick labels
plt.show()

f:id:ossyaritoori:20200430024808p:plain
エリアの勝敗への各変数の影響力:敵味方のキルデスの他に自身のデスなども比較的重要と出ました。


データ数と解析手法がおもちゃなのでこのまま信用することはできませんが,言われてみればちょっと納得できるような結果になったのではないでしょうか。

これを拡張して編成を評価できるようにしていきたいですね。

決定木による検証2:他のルールでのキル・デスの関係性

同様の解析をルール毎にかけてみました。

ここから,ホコ・ヤグラ・アサリはオブジェクトの動きに依存するためエリアのようにリザルトだけを見ても勝敗を予想しづらいという仮説が建てられます。

これらのルールそれぞれでの各変数の影響力を比較してみましょう。xticksという機能を使って折れ線グラフのx軸を各変数名に変えています。

fig, ax = plt.subplots()
plt.grid()
xrange = range(len(area_importances))
ax.plot(xrange,area_importances,marker = 'o',label="area")
ax.plot(xrange,yagura_importances, marker='x',label="yagura")
ax.plot(xrange,hoko_importances,marker = 's',label="hoko")
ax.plot(xrange,asari_importances, marker='*',label="asari")
plt.legend()
plt.xticks(xrange,train_x.columns)
fig.autofmt_xdate() # make space for and rotate the x-axis tick labels
plt.show()

f:id:ossyaritoori:20200430030646p:plain
各ルール毎の変数の重要度の比較:エリアでは味方のキルが重要であることがわかります。また,ヤグラでは敵のキルと試合結果に大きく相関があります。


深さ10,5回交差検定の決定木による解析結果まとめ

  • エリアでは味方のキル数が試合結果に大きく影響する
  • ヤグラでは敵のデス数が試合結果に大きく影響する
  • アサリ以外では味方のキル・デスと敵のキルデスの占める割合が大きい
  • アサリは他に比べて塗りポイントの重要性が高い?

アサリは特に複雑性の高いルールですのでちょっとなんとも言えないですが,言われてみればそこそこ的を射ているのではないでしょうか。

まとめ

長くなりましたが,第一回としては

  • ikawidget2のデータをrealm→json→pandas(csv)形式に変換
  • データ内容,解析の手法・方向性やテスト

について書きました。

  • 戦績のリザルトからある程度勝敗について予測できるのは直感通り
  • どのスコアがどの程度影響するか可視化できるのは面白い

戦績データをためてて解析してみてほしいという人は一報ください。

補足:LightGBMを用いた際のリザルト重要度比較

今度はKaggleでよく使われるLightGBMの決定木を用いて重要度を比較してみましょう。()内はsklearnからの増分です。

  • ガチエリア:平均86.8%の的中率(+2%)
  • ガチヤグラ:平均81.5%の的中率(+5%)
  • ガチホコ:平均75.9%の的中率(+5%)
  • ガチアサリ:平均70.7%の的中率(+7%)

このときの各変数の重要度の比較は,以下のようになりました。

f:id:ossyaritoori:20200501173809p:plain
LightGBMで学習したアンサンブル木での重要度:味方のキルデスの重要度が高いのはそのままですが,敵のキルデスよりもこちらの塗りや自分のデスなどに比重が多くなる場合もあります。

LightGBMの決定木による解析結果まとめ

  • 全般的に味方のキル・デス数が試合結果に大きく影響する
  • SPの回数はあまり勝敗を決しない
  • ヤグラ以外ではデスの方がキルよりも試合への影響が大きい
  • アサリは他に比べて塗りポイントの重要性が高い

TODO

ブキの分類

流石にブキの種類が多いため,そのままブキ名をデータに加えても有意義な解析ができないと思います。 したがって種類ごとに分類します。

  • 公式の分類通りに分類
  • 射程毎に分類
  • 塗りの強さ毎に分類

最も方針が立てやすいのは公式の分類でしょう。ika statsでも似たような表があります。

f:id:ossyaritoori:20200429230610p:plain
ika stats様から引用

塗りの評価はika statsさんが出してくれているのでそれを流用するのも良いかもしれません。

Kaggle-titanicコンペのデータを用いてRandomForest,SVM,LightGBMを試す

前回のつづきです。

(注)本記事はコンペ目的ではなく,ライブラリの使用感を確かめているところで個人メモの範疇です。

ossyaritoori.hatenablog.com

前処理やデータなど

KaggleのTitanicを使っています。 前処理は前記事を参照してください。

RandomForest

例のごとくscikit learnのライブラリを用いていきます。

scikit-learn.org

前回の木構造を複数学習して多数決をとるような手法です。

Boostingなどのキーワードは調べておいても良いでしょう。

コード

デフォルトでは100本の木でアンサンブルしようとするので10本程度に収まるようにさせました。深さ4を上限としています。

import sklearn.tree as tree
from sklearn.ensemble import RandomForestClassifier

# RandomForest 10個のモデルの重ね合わせ
clf = RandomForestClassifier(max_depth=4,random_state=0,n_estimators=10)

ここから交差検証に入ります。5分割くらいでも良かったかも。

# 3分割交差検証を指定し、インスタンス化
from sklearn.model_selection import KFold
kf = KFold(n_splits=3, shuffle=True, random_state=0)

# スコアとモデルを格納するリスト
score_list = []
models = []

# 各分割ごとに評価
for fold_, (train_index, valid_index) in enumerate(kf.split(X_train, y_train)):    
    print(f'fold{fold_ + 1} start')
    train_x = X_train.iloc[train_index]
    valid_x = X_train.iloc[valid_index]
    train_y = y_train.iloc[train_index]
    valid_y = y_train.iloc[valid_index]
    
    ## 分割データで学習・予測・評価
    model = clf.fit(train_x, train_y)
    
    # データを用いて予測,記録
    predicted = model.predict(valid_x)
    score_list.append(accuracy_score(predicted,valid_y))
    models.append(model)
print(score_list, '平均score', round(np.mean(score_list), 3))

このように交差検証法で多数決をとることでデータの偏りの影響を低減できます。(ランダムフォレストの多数決とは別の多数決,ややこしいですね。)

もちろん図示も可能ですが,複数の木の多数決構造をとっているので可視化をする時はそのうちの一つに絞ります。 model.estimators_にリストで入ってるので適当な番号のを覗いてみましょう。

図は前回を参照にしてください。

重要度は全部の木について計算してくれているのかmodelから直接とることができます。

# bar plot
fig, ax = plt.subplots()
plt.grid()
ax.bar(train_x.columns,model.feature_importances_)
fig.autofmt_xdate() # make space for and rotate the x-axis tick labels
plt.show()

f:id:ossyaritoori:20200421010909p:plain
RandomForestで見た変数重要度。決定木と比べて Fareに重点が置かれています。

何本アンサンブルすればいいの?

木を増やせば汎化性能は上がりそうですが,計算量の問題もあります。 モデルの複雑性が増すと過学習しないかなという懸念もあります。

結構議論されているようですが,いくつか見て以下の意見を参考にしました。

How to determine the number of trees to be generated in Random Forest algorithm?

f:id:ossyaritoori:20200421011153p:plain
要約:原則木が多いほど性能が上がるが,上がり幅が小さくなっていく。個々のケースでは木が少ない場合が良い結果を出す可能性はあるかもね。10,30,100あたりで決めたら?

10本 RF 交差検証:0.808 本番:0.77990
30本 RF 交差検証:0.805 本番:0.76555
100本 RF 交差検証:0.805 本番:0.78468

確かに,若干は上がるようです。この辺は試行錯誤する必要があるので後で説明するgrid searchをすると良いでしょう。

SVMサポートベクターマシン

木構造とは異なりますが,SVMのような分類手法も有名です。

簡単に言うとデータ群を最適に分割する超平面を求める手法です。

www.slideshare.net

  • カーネルトリック:平面で分割とは言ってもデータの変数をかさ増しすることで非線形な分割も可能になります。例えば元の変数がx,yなら,x^2, xy ,y^2も変数に加えた空間で分割すれば元の空間では直線ではなくなります。ポイントは内積計算が楽になるようにこのかさ増しを工夫することです。

sklearnではデフォルトでRBFカーネルというものが使われています。

参考:

[python 機械学習初心者向け] scikit-learnでSVMを簡単に実装する - Qiita

カーネル近似(クラス分類)【Pythonとscikit-learnで機械学習:第2回】

コード

from sklearn.svm import SVC

clf = SVC(kernel='rbf', C=1, gamma=0.1)


# 3分割交差検証を指定し、インスタンス化
from sklearn.model_selection import KFold
kf = KFold(n_splits=3, shuffle=True, random_state=0)

# スコアとモデルを格納するリスト
score_list = []
models = []

# 各分割ごとに評価
for fold_, (train_index, valid_index) in enumerate(kf.split(X_train, y_train)):    
    print(f'fold{fold_ + 1} start')
    train_x = X_train.iloc[train_index]
    valid_x = X_train.iloc[valid_index]
    train_y = y_train.iloc[train_index]
    valid_y = y_train.iloc[valid_index]
    
    ## 分割データで学習・予測・評価
    model = clf.fit(train_x, train_y)
    
    # データを用いて予測,記録
    predicted = model.predict(valid_x)
    score_list.append(accuracy_score(predicted,valid_y))
    models.append(model)
print(score_list, '平均score', round(np.mean(score_list), 3))

結果は[0.7104377104377104, 0.6868686868686869, 0.6666666666666666] 平均score 0.688と少し低めに出てしまいました。

ハイパーパラメータ調整の指針

以下のサイトから引用します。

scikit-learn.org

upura.hatenablog.com

  • γ:「一つの訓練データが与える影響の範囲」を意味します。小さいほど「遠く」、大きいほど「近く」まで影響します。
  • C:分類ミスに対するペナルティと分離平面からの距離のトレードオフです。大きいほど誤りに厳しく,データすれすれに線をひきます。したがって小さいほど大雑把でシンプルな分類になりやすいようです。

実際のシチュエーションでは複数のハイパーパラメータを試して最も良さそうな組み合わせをチョイスする必要があります。 for文でも書けますが,sklearnにはgrid_searchという関数があります。

バージョンによって関数名が変わるのですが,Ver1.0では以下のように書きます。

from sklearn.model_selection import GridSearchCV

# グリッドサーチを用いたオートチューニング
parameters = [{'kernel':['rbf'], 'C':np.logspace(-4, 4, 9), 'gamma':np.logspace(-4, 4, 9)}] # RBF kernel でグリッドサーチ

# 交差検定5回でSVCを用いてやる
clf = GridSearchCV(cv=5,estimator=SVC(), param_grid=parameters, n_jobs = -1)
out = clf.fit(X_train, y_train)

print(out.best_estimator_)

このbest estimatorが最も成績の良かった学習器です。

SVC(C=10000.0, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma=0.0001, kernel='rbf',
    max_iter=-1, probability=False, random_state=None, shrinking=True,
    tol=0.001, verbose=False)

このときの精度は0.8417508417508418と前回より随分高くなりました。

LightGBM

LightGBMは決定木のアルゴリズムを用いた機械学習フレームワークです。以下のサイトが勉強になります。

www.codexa.net

コード

使い方はsklearnのものと非常によく似ています。

import lightgbm as lgb

# LightGBMの分類器をインスタンス化
gbm = lgb.LGBMClassifier(objective='binary')

# trainとvalidを指定し学習
gbm.fit(train_x, train_y, eval_set = [(valid_x, valid_y)],
        early_stopping_rounds=20,  # 20回連続でlossが下がらなかったら終了
        verbose=10  # 10round毎に、lossを表示
) ;

一応推論をしてみると,score 81.48 %と出ます。

# valid_xについて推論
oof = gbm.predict(valid_x, num_iteration=gbm.best_iteration_)  # oofはout of fold
print('score', round(accuracy_score(valid_y, oof)*100,2), '%')  # 正解率の表示

図示

こちらはデフォルトでそこそこ見やすい図が出ます。

lgb.create_tree_digraph(gbm)

f:id:ossyaritoori:20200428003643p:plain
一部分のみ。

重要度は以下のように簡単にplotできます。

lgb.plot_importance(gbm, figsize=(12, 6))

f:id:ossyaritoori:20200428004244p:plain
LightGBMでの変数重要性:船賃と年齢が圧倒的で,性差は低めに見積もられていることがわかります。

提出に用いた交差検証 LightGBM

こちらも交差検証を用いて複数のモデルの多数決をとります。

# 5分割交差検証を指定し、インスタンス化
from sklearn.model_selection import KFold
kf5 = KFold(n_splits=5, shuffle=True, random_state=0)

# スコアとモデルを格納するリスト
score_list = []
models = []
test_predicts = []

# 各分割ごとに評価
for fold_, (train_index, valid_index) in enumerate(kf5.split(X_train, y_train)):    
    print(f'fold{fold_ + 1} start')
    train_x = X_train.iloc[train_index]
    valid_x = X_train.iloc[valid_index]
    train_y = y_train.iloc[train_index]
    valid_y = y_train.iloc[valid_index]
    
    ## 分割データで学習・予測・評価
    gbm.fit(train_x, train_y, eval_set = [(valid_x, valid_y)],
            early_stopping_rounds=20,  # 20回連続でlossが下がらなかったら終了
            verbose=10  # 10round毎に、lossを表示
    ) ;
    # データを用いて予測,記録
    test_predicts.append(gbm.predict(test, num_iteration=gbm.best_iteration_) )
    predicted = gbm.predict(valid_x, num_iteration=gbm.best_iteration_) 
    score_list.append(accuracy_score(predicted,valid_y))
    
print(score_list, '平均score', round(np.mean(score_list), 3))
[0.8491620111731844, 0.8089887640449438, 0.8314606741573034, 0.8033707865168539, 0.8089887640449438] 平均score 0.82

コンペ結果

気休めですが,LightBGMが僅差でランダムフォレストを上回っていました。

f:id:ossyaritoori:20200428004543p:plain
そこまで各手法で劇的に差が出るわけではありませんね。

他にも手法があるので,概略と動かし方をわかっていれば今後の解析の時に役立ちそうですね。

scikit-learn.org

[https://scikit-learn.org/stable/images/sphx_glr_plot_classifier_comparison_001.png:image=https://scikit-learn.org/stable/images/sphx_glr_plot_classifier_comparison_001.png]

KaggleのTitanicデータに対してsklearnの決定木を試してみる

(注)コンペ目的でない素人の備忘録です。参考になるかは不明ですがコメントは歓迎します。

Kaggleとは

www.kaggle.com

KaggleはHPにあるとおり機械学習/データサイエンスに携わる人や企業などのコミュニティです。

正直始めたてなのでざっくりとしかわからないですが,Competitionという、企業や政府などがコンペ形式で課題を提示し、賞金と引き換えに最も制度の高い分析モデルを買い取るといったシステムが有名です。

勉強用としてのKaggle(Titanic)

機械学習は教科書で深層学習やパーセプトロンなどをざっくり読んだだけなのでコンペに参加することはしばらくないでしょうが,勉強用にはいくつか使えそうな点があります。

例えば,練習用としてTitanicというタイタニック号事件のデータが提供されています。 構造もシンプルで初学者にとっても分かりやすいです。

また,人によってはnotebookをわかりやすくコメント付きで投稿してくれていたりします。 www.kaggle.com

他にも,ブラウザ上からアクセスするnotebook形式ではデフォルトでsklearnやlightBGMなどのライブラリがインストールされてる状態で始まるので環境構築の心配が現状なさそうです。

なお,真面目に入門するなら以下の解説が懇切丁寧でよかったです。(英語)

www.kaggle.com

私個人の目的としてはスプラトゥーン2の編成格差問題を自前のデータから解析してみようかなと思っています。(別記事にて記載予定)

メモ:Notebookの作成

初めてKaggleに登録して開くと何をしていいのか本当によくわからなくなりますね。 コンペを探してそれに参加するというのが本来でしょうが,今回はKaggleの環境とデータで練習をするという目的でNotebookを作ります。なお,PythonやPandasなどの知識はどこか別のところで勉強してください。

コンペから探して作成

上のSearchでtitanicを検索すればtitanicのページに移動します。

f:id:ossyaritoori:20200419162920p:plain
Titanicでサーチした結果のnotebookタブ

右上の「New Notebook」から始めます。

なお,notebookのタブには他の人のnotebookがありさんこうになります。

新規作成→Notebookにあとからデータを追加

f:id:ossyaritoori:20200420011520p:plain
Competition Dataから追加します。



sklearnの分類木ライブラリを用いた分類

決定木の説明は教科書等にに投げます。(Qiitaなら以下のリンクなんかが読みやすいのでは?)

[入門]初心者の初心者による初心者のための決定木分析 - Qiita

今回はsklearnというライブラリを用いて決定木のうち,分類木 (classification tree)を用いて学習します。

import pandas as pd
import numpy as np
import warnings
warnings.filterwarnings('ignore')

train = pd.read_csv('../input/titanic/train.csv')
test = pd.read_csv('../input/titanic/test.csv')
sample_submission = pd.read_csv('../input/titanic/gender_submission.csv')

trainは訓練用,testは評価用のデータです。

f:id:ossyaritoori:20200419180806p:plain
trainの中身はこんな感じ。氏名とか公開されてるんですね。

データの前処理

まずはデータの前処理をします。 機械学習のデータ前処理備忘録 - Qiita

主に,以下のようなことをします。

  • 訓練に必要なデータの抜き出し,加工
  • 空白データの補間,削除

データ形式の変換と不要なデータの削除

# データの変換
##maleを0に、femaleを1に変換
train["Sex"] = train["Sex"].map({"male":0,"female":1})
test["Sex"] = test["Sex"].map({"male":0,"female":1})
# EmbarkedのOne-Hotエンコーディング
train = pd.get_dummies(train, columns=['Embarked'])
test = pd.get_dummies(test, columns=['Embarked'])

## 不要な列の削除
train.drop(['PassengerId', 'Name', 'Cabin', 'Ticket'], axis=1, inplace=True)
test.drop(['PassengerId', 'Name', 'Cabin', 'Ticket'], axis=1, inplace=True)
  • One-Hotエンコーディング:[A,B,C]のような3つ以上の値のある要素をAの有無としてBoolean値に変換する。

f:id:ossyaritoori:20200419184512p:plain
必要そうなデータを抜き出したあとのテストデータ。

欠損データの確認と対処

ところで,データにはところどころ欠損(NaN)があり,学習に悪影響を及ぼすので取り除いたり穴埋めをしたいです。 pandasの機能でどの量に欠損があるのか簡単に割り出せます。 例えばNanを確認して取り除く場合は以下のように書くことができます。

# NaN の存在確認 と 除去
print(train.isnull().sum())

train2 = train.dropna()
test2 = test.dropna()

一応結果を示すと「年齢」の項が抜けているのがわかりますね。

Survived        0
Pclass          0
Sex             0
Age           177
SibSp           0
Parch           0
Fare            0
Embarked_C      0
Embarked_Q      0
Embarked_S      0
dtype: int64

他にも中央値や平均値,ランダム値などで埋める手法があるようです。 参考:平均値でNanを埋めるタイプの方

学習

さて,いよいよ学習のフェーズです。まずは与えられたデータを入出力に分けます。ここでの出力(=予測したい値)は生存したかどうかなのでSurvivedの列を抽出します。

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


# devide train input and output
X_train = train2.drop(['Survived'], axis=1)  # X_trainはtrainのSurvived列以外
y_train = train2['Survived']  # y_trainはtrainのSurvived列

# X_trainとY_trainをtrainとvalidに分割
train_x, valid_x, train_y, valid_y = train_test_split(X_train, y_train, test_size=0.20, random_state=0)

最後の実行では訓練用のデータを更に分けて4:1に分けて訓練の評価をします。 訓練用データ内でのフィッティング結果と実際の評価用データ内でのフィッティング結果を分けたいがためです。(過学習とかで検索してください。)

訓練用データを2つに分けるのをホールドアウト法,3つ以上に分けるのを交差検証法と言います。(引用

これで訓練用の入力と出力を分けられたのでsklearnの分類木ライブラリの学習で学習します。

import sklearn.tree as tree

# 分類木だからClassifier,Regressorもある
clf = tree.DecisionTreeClassifier(max_depth=4)

# データを用いて学習
model = clf.fit(train_x, train_y)

# データを用いて予測
predicted = model.predict(valid_x)
print(accuracy_score(predicted,valid_y))

なんと,これだけです。 結果は0.7902097902097902になりました。約8割程度の評価です。これは訓練用データに8割程度の正答率で合致しているという解釈になります。

補足:交差検証法による評価

先程の結果はデータのどの部分を訓練用に使うかにより左右されそうです。

したがってよりきちんとやる場合は,交差検証法で訓練データと評価用データをうまくローテーションして平均値などでモデルを評価します。

# 3分割交差検証を指定し、インスタンス化
from sklearn.model_selection import KFold
kf = KFold(n_splits=3, shuffle=True, random_state=0)

# スコアとモデルを格納するリスト
score_list = []
models = []

# 各分割ごとに評価
for fold_, (train_index, valid_index) in enumerate(kf.split(X_train, y_train)):    
    print(f'fold{fold_ + 1} start')
    train_x = X_train.iloc[train_index]
    valid_x = X_train.iloc[valid_index]
    train_y = y_train.iloc[train_index]
    valid_y = y_train.iloc[valid_index]
    
    ## 分割データで学習・予測・評価
    clf = tree.DecisionTreeClassifier(max_depth=4)
    model = clf.fit(train_x, train_y)
    
    # データを用いて予測,記録
    predicted = model.predict(valid_x)
    score_list.append(accuracy_score(predicted,valid_y))
    models.append(model)
print(score_list, '平均score', round(np.mean(score_list), 3))

結果は以下の通り,データによってばらつきがあるということがわかりますね。

[0.8067226890756303, 0.7647058823529411, 0.819327731092437] 平均score 0.797

混合行列による予測結果の視覚化

予測結果の成否をよりわかりやすくするために混合行列を用いた評価ができます。

from sklearn.metrics import confusion_matrix

#混同行列の作成
cmatrix = confusion_matrix(valid_y,predicted)

#pandasで表の形に
df = pd.DataFrame(cmatrix,index=["actual_died","actual_survived"],columns=["pred_died","pred_survived"])

print(df)
                 pred_died  pred_survived
actual_died            131             14
actual_survived         29             64

ざっくりというと,死んだという予測は18.1%で間違っており,生きていたという予測は17.9%で間違っていたということになります。 誤差のようなものですが,生きていたという予測のほうがかろうじて得意と言えそうです。

木構造の図示

学習した木を評価する手法にはいくつかあります。まずは木の構造を見たいと思うので木を図示します。

# 図示その1
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(30, 30))  # whatever size you want
tree.plot_tree(model, ax=ax)
plt.show()

ちょっと見づらい感じになります。

f:id:ossyaritoori:20200419234211p:plain

より良く見せたいならば,graphvizdtreevizなどを使うと良いと思います。

木構造のかっこいい視覚化

Kernelに以下のモジュールをインストールします。

!pip install dtreeviz pydotplus
# graphvizによる視覚化
import pydotplus as pdp

file_name = "./tree_visualization.png"
dot_data = tree.export_graphviz(model, # 決定木オブジェクトを一つ指定する
                                out_file=None, # ファイルは介さずにGraphvizにdot言語データを渡すのでNone
                                filled=True, # Trueにすると、分岐の際にどちらのノードに多く分類されたのか色で示してくれる
                                rounded=True, # Trueにすると、ノードの角を丸く描画する。
                                feature_names=train_x.columns, # これを指定しないとチャート上で特徴量の名前が表示されない
                                class_names=['Dead','Sruvived'], # これを指定しないとチャート上で分類名が表示されない
                                special_characters=True # 特殊文字を扱えるようにする
                                )
graph = pdp.graph_from_dot_data(dot_data)
graph.write_png(file_name)

# dtreevizによる視覚化
from dtreeviz.trees import dtreeviz

viz = dtreeviz(
    model,
    train_x, 
    train_y,
    target_name='alive',
    feature_names=train_x.columns,
    class_names=['Dead','Sruvived']
) 

#viz.view() # 外窓で画像を開くのでKaggleのKernelではエラーがでる?
viz.save(”dtreeviz.svg”)

f:id:ossyaritoori:20200421003827p:plain
graphvizによる可視化

f:id:ossyaritoori:20200421004012p:plain
dtreevizによる視覚化:どの様な分類をしたのかかなり分かりやすいですね。

木構造の解析:変数の重要性

では,解析の結果どの変数が生存に大きな影響を及ぼしたのでしょうか。 関数を叩くと一発ですが,以下のワードくらいは抑えておくとします。

  • ジニ不純度: ノードによる分類のバラケ方(≒エントロピーみたいな)を数値化したもの→半々で最大。「1 - 各分類木のノードの分類割合の2乗和 」の値。CARTで使われている。→ sklearnでのデフォルト。
  • 変数重要度:”ある特徴量で分割することでどれくらいジニ不純度を下げられるか。”を表す量。分割後のジニ不純度のサンプル数の分の和が,分割前のジニ不純度のサンプル数和よりどの程度小さくなったかを表す。

詳細な説明はこちらの記事に投げます。


treeオブジェクトでは変数重要度は.feature_importances_というメンバに入っています。

# 重要度を表示
print(dict(zip(train_x.columns, model.feature_importances_)))

# bar plot
fig, ax = plt.subplots()
plt.grid()
ax.bar(train_x.columns,model.feature_importances_)
fig.autofmt_xdate() # make space for and rotate the x-axis tick labels
plt.show()

f:id:ossyaritoori:20200420003916p:plain
1に性別,2に乗客の等位,3に年齢となっていますね。

確かタイタニックでは女性と子供が優先的に救命ボートに乗せられたはずなのでこの洞察は合っていると言えるでしょう。

ただし重要度は,如何にクラスを上手く抽出できたかを表す無方向性の値なので,年齢が幼い方が良いのか女性の方が良いのかはわかりません。どうすればよいかちょっと調べてみます。(TODO)

テスト提出

とりあえず提出用のデータを作ります。

# テスト  これを提出

## test の中身
print(test.isnull().sum())
## 中央値で埋める
test = test.fillna(test.median())

## 予測結果
test_predicted = model.predict(test)

## 提出用データ
sample_submission['Survived'] = test_predicted
sample_submission.to_csv('submission.csv',index=False)

このcsvをダウンロードしてSubmit Predictionの該当箇所にUploadします。

f:id:ossyaritoori:20200420014616p:plain
Upload中の様子

本当はNotebookでSubmitしたいんですが,Submitできそうなボタンがないような?(@2020年4月20日)なんでやろ?

f:id:ossyaritoori:20200420015626p:plain
提出した結果は74%くらい。訓練時より低く出ました。そりゃそうか。

スコアボードを見ると全く同じスコアの人々がいるので多分テストコードを走らせただけの人々がそれなりにいるのでしょう。 なお,正答率100%の人々もかなりいます。だって答えわかってるもんね。

最後に一応My Notebookを載せます。大してきちんと書いてないですけどね。

www.kaggle.com

補足・TODO

気分がアガり次第書きます。スプラトゥーンの記事は書いたまま塩漬けになってるので今度ちゃんと書き起こします…

TODOs:

今回は最も単純な分類木を用いましたが,

  • スプラトゥーンの記事との連携
  • ランダムフォレストやlightBGMなどの実装と評価
  • SVM等の他の手法に関する考察

補足:今後の参考

素晴らしい記事をありがとうございます。後ほど参考にさせて頂く予定です。

maruo51.com

qiita.com

pythondatascience.plavox.info

ボツデータ:欠損値をどのように埋めるか問題

欠損値をどのように埋めるかですが,3番めに重要と思われるAgeが抜けておりました。

絶対適当に埋めるべきではない気がするのですが,中央値,平均値,それとも削除して無視するのどれが良いのか比較しようと試みました。

結論から言うと結果は全く変わらず。なんか違いが出ると思ったんですがね…

# 欠損値除去 VS 補填 (なお,比較のためここでやっているが,本来はデータの前処理の段階でやるべき)

## train_yのうちtrain_xがNanを含む行をdropする
train_y_dropna = train_y.drop(train_x[train_x['Age'].isnull()].index.tolist())
train_x_dropna = train_x.dropna()

## 埋める方 shallow copy
train_x_fillmed = train_x
train_x_fill0 = train_x

train_x_fillmed["Age"] = train_x_fillmed["Age"].fillna(train_x["Age"].median())
train_x_fill0["Age"] = train_x_fill0["Age"].fillna(0)

# それぞれFit
model_dropna = clf.fit(train_x_dropna, train_y_dropna)
model_fillmed = clf.fit(train_x_fillmed, train_y)
model_fill0 = clf.fit(train_x_fill0, train_y)

# 予測

## validの欠損値補間はMedianで
valid_x["Age"] = valid_x["Age"].fillna(valid_x["Age"].median())

pred_dropna = model_dropna.predict(valid_x)
pred_fillmed = model_fillmed.predict(valid_x)
pred_fill0 = model_fill0.predict(valid_x)

## 評価
print(accuracy_score(pred_dropna,valid_y),
      accuracy_score(pred_fillmed,valid_y),
      accuracy_score(pred_fill0,valid_y))

追記:ちゃんと変化しました。Medianで適当に埋めるほうが少し精度が下がりそうです。

[0.8148148148148148, 0.7811447811447811, 0.7542087542087542] 平均score 0.783

f:id:ossyaritoori:20200420102911p:plain
Medianで適当にNanを埋めて学習した時の重要度。他の変数にも重みが振られていることがわかります。

柏キャンパス-江戸川台駅 のグルメ攻略ガイド<Maniac>

どうも。本日柏とお別れしましたが,今後は終身栄誉柏民です。

さるネイティブカシワリアンの書いた記事に触発(挑発)されたので, 柏キャンパス-江戸川台のエリアにフォーカスしたグルメ攻略ガイドのManiac版を公開します。

TX民?ララポに行きなさい。

柏キャンパスのグルメがわかる記事

渋谷・本郷と暮らしてきた人々からすれば柏の葉は不毛の地と思われがちですが,よく見ると柏キャンパス周辺にはたくさん食事処があります。

以前簡潔に書いた柏キャンパス民のためのグルメガイド記事から約2年が経ちました。

ossyaritoori.hatenablog.com

その後,優秀なネイティブカシワリアン?氏に以下の記事にて補足(と煽り)をいただきました。

note.com

ということで,本記事の趣旨は上記の記事の店も含めよりManiacな 表通り(通称グルメストリート)にないお店を紹介することです。

f:id:ossyaritoori:20200412021626p:plain
ということで一覧図。

和食

突然ですが,柏キャンパスには和食があまりありません。 チェーンを除くと基本中華の勢力が強いんですね。

鳥一理彦などは必修ですので他の和食店を紹介します。 といっても実質2店舗くらいしかないんですけど。

わっ嘉

インドカレー屋の向かいの路地裏あたりにあるお店。昼はランチ,夜は飲み屋らしいです。

  • ○:小皿の料理が美味しい
  • ○(✖):店主がとてもフレンドリーでアットホームな雰囲気(苦手な人はいると思う)
  • ✖:学生にはちと高い
  • △:ランチしかやってない/日曜空いてない

”グルメ”としては結構完成度高いので自分へのご褒美としてどうぞ。夜も安く飲ましてくれるそうです。

f:id:ossyaritoori:20200412005946p:plain
鯖の塩焼き定食の図

f:id:ossyaritoori:20200412005640p:plain
昼のメニュー。他にも釜飯など美味しそうなメニューがたくさんあります。

12品の彩り定食が魅力の「わっ嘉」は、車いす料理人が江戸川台にオープンした和風ダイニング | リビングかしわWeb

美こま・お食事処

すごく評判が良さげなので行きたかったのですが,ついぞ行けなかったので名前だけの紹介にします。 レビューの感じもすごく良さそうなので期待を込めて入れています。

だんらん食堂

どうやらチェーンの類のようですがあまり知られていないのでここに記します。 自分でおかずを選んでご飯,豚汁と組み合わせるタイプの店です。

f:id:ossyaritoori:20200412011223p:plain

  • ○:比較的手軽な値段(600円程度~)で健康的な定食がを自分で選んで食べれる。
  • ○:250円くらいでご飯&豚汁を無限おかわりできる。うまい。
  • ✖:江戸川台駅から西へチャリで10分。超遠い。

洋食

江戸川台駅の北のエリアには実はフレンチのビストロがあったりします。うっかり柏に来た友達などを接待してやる時に使うと良いでしょう。

フレンチ(コンツェルト・ブラッスリータケ)

駅から北に歩いて5分程度の位置に並んで2店舗があります。2店舗の違いはそんなに通ってないのでなんとも言えないですが個人的にはコンツェルトの方のメニューが印象に残っています。

f:id:ossyaritoori:20200412012107p:plain
こんな感じのちゃんとしたコース料理が食べれます。

  • ○:駅から近い
  • ○:1500円くらい出せばちゃんとしたコースが食べれる
  • △:昼遅めになるとお得なコースがなくなってたりするので注意
  • あくまで接待用??

イタリアン

イタリアンはいくつかあるようですが,行ったことのある店のみにします。

ボスコロッソ

駅の広場のすぐ裏にあるカフェっぽいイタリアンです。 1000円前後でランチが食べられます。

  • ○:駅から近い
  • ○:イタリアンなので基本はずれない
  • △: 値段からして学生が通う店ではない気がする

ひるふぁ~む

みどり台郵便局の所にあるピザの美味しいお店。 ランチは1000円強くらいです。

  • ○:ピザがうまい。サラダもうまい。おそらくパスタもうまいんだろう。
  • ○:健康的な食事をした!という気分になる
  • △:通うには高いなぁ…

イタリアンにはコスパサイゼリアがあるのでなかなかよっしゃ食べに行こうとはならないのが残念なところです。

中華

中華は激戦区です。イチオシは華記食府ですが,個人的には小満も好きです。 ホントはもっといろいろあるんですが紹介したい店だけにします。あと,ラーメンも同列です。

ラーメン食いたければ柏駅に行け。

ほんわ華屋

ヨークマートの近くにある品の良い中華屋という感じです。

f:id:ossyaritoori:20200412013925p:plain
品の良いとは言ったが夏に食ったジャージャー麺の写真しかなかった…

  • ○:胃に優しいタイプの中華
  • △:1000円前後でドリンク付きランチ。やはりちょっと高い

親御さんなど連れてきてもいいでしょう。

ラーメンイヨマンテ

おじいちゃんが経営しているラーメン屋さん。ベーシックだけど結構好きな味です。

f:id:ossyaritoori:20200412014405p:plain
Captured From Google Map

  • ○:値段の割に量が多め。お腹いっぱいになれる。
  • △:場所を知っている人は殆ど居ないだろう。近くに住んでいるならどうぞ。

さくらミニ

駅前北側にあるラーメン屋さん。意外と気づかないような位置にある。 まぜそばが人気のようだが,個人的には味噌ラーメンもイケる。

f:id:ossyaritoori:20200412014848p:plain
Captured from Google Map

f:id:ossyaritoori:20200412014949p:plain
味噌ラーメンは太麺と濃いスープが合う。

  • ○:結構美味しい
  • ○:駅に近い
  • ✖:柏駅などに強い競合が多い。駅前のこの店に来れるなら柏駅に行っちゃう。

その他

その他の雑多なグルメ情報です。

ケーキ・パン屋

パン屋として必修なのはサフランですが,付近のちょっと路地を入った所にミレーというケーキ屋さんがあります。

f:id:ossyaritoori:20200412015518p:plain
可愛げのある店内。Captured from Google Map

ケーキのデザインも良く,値段も高くないので気になる教授にケーキを買う時に是非(?)

配達系

知っているか。柏キャンパスにUber eatsは来ない(2020年4月現在)。

配達系はピザか寿司,かわりものでパエリア程度しかないのでおうちで何か注文するには絶望的に適していない。ご注意。

スーパー

新人はわかりやすい駅前の”かどや”やヨークマートに行きがちですが,より食品の鮮度・コスパで優れているのはフーズマーケットセレクション西原店です。 位置はほぼ初石ですが,果物系の鮮度や安さ,暴力的な弁当のコスパ(~300円でそこそこ)を考えると行くべきはここ! 頑張って自転車を漕ぎましょう。

なお,近くの西原小学校では毎週日曜日の昼に有志でバドミントンをやっています。興味ある方はそちらもどうぞ。

16号沿い

キャンパスの裏から高校の脇の道を通って国道16号線に出れます。キャンパスから最も近い牛丼屋である”すき家”他よくあるチェーンやラーメン屋もあるので頑張って足を伸ばしてみましょう。 もつ煮太郎は必修な。

おわりに

柏キャンパス-江戸川台付近にはもっともっとお店があるのですが,ある程度書いて満足したのでこの辺で終わりにします。 みなさんもぜひお店を開拓して見てください。

【雑談】宝石をじゃらじゃらしたいのでAmazonUSで100カラット注文してみた話

この記事は私が疲れている時に勢いで書いたものになります。

読まれる方もどうか勢いをつけて読んでください。

追記:後に冷静になってから貴石達が届きました。石に詳しい人に是非助言をいただきたいです。

宝石じゃらじゃら欲

追記:この章は茶番です。結局AmazonUSで買いました。

宝石じゃらじゃらしたい。これは石に興味のある子供なら誰しも思ったことがあるはずである。

というかキラキラ輝く宝石を手に収めて眺めていたい欲は多分結構普遍的に人類の持っている欲だと思う。

https://iroishi-kaitori.jp/img/og-image-min.jpg

昔は道端で透明な石英を拾って喜んでいたちびっこだった私はその後もビー玉を集めたりしたものだった。ラムネの瓶から抜いてくるやつがちょっとマットな質感があって好きだった。 鉱物に関する本も結構読んだ気がする(ただし小学生レベル)。

小さい頃はそれで我慢できたかもしれない。が,ちょっと大人になった今なんだか急に宝石をじゃらじゃらしたい欲が湧いてきた。

アクリルでもガラスでもない。「本物」をただじゃらじゃらしたいのだ。 はい,というわけで探してきました。安価にじゃらじゃらする方法を。


え?安価に??

https://imasoku.com/wp-content/uploads/2019/08/%E3%82%B9%E3%83%BC%E3%83%91%E3%83%BC%E3%81%8F%E3%81%84%E3%81%97%E3%82%93%E5%9D%8A.jpg

サブスクリプションサービス

デザインアトリエカケラというサブスクリプションサービスがあることを最近になって知った。

コンセプトも夢があって,加工ミスや品質の問題により値段があまり付けられないような宝石を定期的にランダムに送り届けてくれるというものだ。

community.camp-fire.jp

  • 本来価値の高いものを訳ありということで安価に
  • ある程度ランダム性をもって届けられる

という点において「訳あり高級カニ」系統に劣らないワクワク感を演出してくれる。


このランダムというのがおみくじ(ガチャ)感あってとても好きだが,私は今すぐじゃらじゃらしたいので今回は保留とした。

なお,こちらのショップはバラ売りも行っている。定期的に眺めたい感じはある。

www.atelierkakera.biz

ネットショッピング

こういう掘り出し物的なのはネットで買うのが一番である。どうやら「天然石」「ルース」「訳あり」などのワードでサーチすると良いらしい。

専門店で買う

普通に「天然石」「ルース」「訳あり」などのワードでサーチすると次のようなページに辿り着く。

nature-guidance.jp

まずページがずるい。きれいな形の宝石がじゃらっと並んでるだけでIQが下がってしまう。 しかも一個あたり結構安いのでホイホイ買えてしまうではないか。

https://nature-guidance.jp/pic-labo/rondonbluetopazobal641.jpg

ただし実際のサイズより写真ではかなり大きく見えていることに注意。

また,ぶっちゃけどの宝石をどの割合でとかを考えるのも面倒だ。 私は「じゃらじゃら」が注文したいんであって最適な石の選定には興味がなかった。たくさんじゃらるには結構値も張るので今回は保留とした。

ヤフオク・メルカリ・楽天etc

上記の専門店のような業者は一部フリマにも進出している。 友人曰く交渉次第で良いものを手に入れられるらしいが,私は宝石をじゃらじゃらしたいチンパンなので今回はパス。

しかし,メルカリって無限に時間潰せちゃうよね。

Amazonで買う

今や偽レビューや謎業者によって信頼性がガタ落ちしてはいるが腐ってもネットショッピング大手。Amazonにはこの世の半分くらいのものは揃っている。


たまたまそれっぽいのを見つけたので筆者は最終的にAmazonで買うことにした。

日本のAmazon

自分が探した中だと以下のようなランダムミックスが面白そうと思った。他にもいくつかあったがアクリルが混じっているなどの悪いレビューがあったりするので気をつけてみるべし。

5000円は少し買うのに勇気がいる値段だがサブスクリプション毎月待たずとも一気にじゃらじゃらできる量が手に入る。 具体的にレビューの写真を見ていこう。

f:id:ossyaritoori:20200124231124p:plain
結構コンパクトに送られてくる。あと保証書てきなのもあるらしい。
f:id:ossyaritoori:20200124230933p:plain
じゃらじゃら~(語彙力)

料金は5000円くらい,納期は2週間といったところか。

f:id:ossyaritoori:20200124230349p:plain
日本のAmazonで買った場合

すでにありなのだが海外から輸入になるので仲介業者の存在により高くなっている場合がある。

英語で同様のショップを調べた所アメリカのAmazonにも同じ商品を見つけた。

アメリカのAmazon

海外のAmazonからものを買う手法は

  1. AmazonGlobalというサービスを使う
  2. 普通にアメリカのAmazonにアカウント登録して日本に送らせる

があるが,目当てのものはAmazonGlobalで探せなかったので仕方なくアメリカのAmazonのアカウントを作成してみた。

mvno.xsrv.jp

これの良いところは仲介の店舗を介さないのでより安く買える点である。

https://www.amazon.com/gp/product/B006C4VWX2/ref=ppx_od_dt_b_asin_title_s00?ie=UTF8&psc=1

f:id:ossyaritoori:20200124230544p:plain
自分が買った時の時価はこんなかんじ。2000円近く安いだけでなく届く日にちも早くなっているのがわかる。ウケる。

というか納期も日本より早いのでアカウント作る手間を惜しまなければアメリカのAmazonで買ったほうが良さげである。

なおコメントにもあるのだが,これはロットによって当たり外れが大きく,外れた人は低評価をつけている。 曰く,

  • 欠けや傷が存在している
  • 色や大きさに偏りがある
  • 思ったより少ない

ということらしい。 しかしまぁ「本来価値のあるものバルクで扱ってみたい」という本趣旨を考えるにそこまで問題にはならなそうではある。

新年まだやっていないおみくじのつもりで買ってみることにする。ガチャっとな。

https://www.fate-go.jp/manga_fgo/images/comic03/comic_point.png

追記:届いたものについて

茶番終了です。 こちらが届いたものになります。アメシスト多めですが各色揃っているのがわかります。

f:id:ossyaritoori:20200203215339p:plain
石が漏れてるのは袋が破れたから

内訳について

鑑定に関しては素人ですが多い順に並べていきます。

アメシスト紫水晶):27個

一番多かったです。 手前から奥に向けて濃くなるように並べました。

f:id:ossyaritoori:20200203215604p:plain
おそらく全部紫水晶

ガーネット?:15個

次に多かったのが暗赤色の宝石です。多分ガーネットと思います。

f:id:ossyaritoori:20200203215829p:plain
ピントがあってませんね。加工傷が一番目立っていました。

水晶?:11個

石英(クォーツ)のうち特に透明度の高いものをクリスタルと言います。

他にも,トパーズなんかは透明のものもあるようです。 【ダイヤ】無色透明の宝石一覧【珍しい石】 - NAVER まとめ

以下の写真は少し自信がないので色ごとに分けています。

f:id:ossyaritoori:20200203220757p:plain
左と右では明らかに色合いが違います。シトリンの一種なのですかね?

でもバルクで売ってる以上全部クォーツだと思います。

スモーキークォーツ(煙水晶):5個

水晶が放射線を浴びたり中にアルミニウムイオンなどが含まれている時にできると言われています。

f:id:ossyaritoori:20200203222743p:plain
一番濃い色のやつはガーネットと紛らわしかったです。

シトリン(黄水晶):3つ + 何か1つ

紫水晶を加熱することで得られるのがシトリンらしいです。こちらは鉄イオンが関係しているとか。

f:id:ossyaritoori:20200203224924p:plain
下の小さいのは明らかに色が違いますし別の鉱物に見えます。

残りの1つはなんでしょう?ペリドット,イエローベリル,トパーズあたりでしょうか?

翡翠?:2つ

濁り具合が翡翠を想起したので翡翠としましたがよくわかりません。 大粒で入っているのでそこまで貴重な石ではなさそうですが…

f:id:ossyaritoori:20200203223432p:plain
個人的にはなんかグミっぽくて美味しそうと感じました()

ということで袋に詰め替えて終了です。 じゃらじゃらは数回して満足してしまったので鑑定の方に時間を費やしました。


やはり福袋的な感じがして楽しいのですが正解が非常に気になる…石に詳しい人に聞きたいですね。 化学系の研究室なら紫外線ライトとかがあるんでしょうが…。

f:id:ossyaritoori:20200203224034p:plain
計64個。26とは良い数ですね。

布教コーナー

宝石の国は良いぞ。

宝石の国(1) (アフタヌーンコミックス)

宝石の国(1) (アフタヌーンコミックス)

RICOH THETA_SCで多重露光・合成をしてPixel4みたいにクリアな星空を撮りたい

この記事はたまたま集まった技術&もの作りが好きな人たち Advent Calendar 2019の記事です。

はじめに

それはゴミのような論文を泣く泣く提出した日の夜でした。デスマを終えて見上げた空は信じられないくらい澄んでいてそれはもう最高に晴れやかな気分でした。

『こら撮らな!』と思いたって家からTHETAを引っ張り出してきて撮ったのが以下の写真です。

f:id:ossyaritoori:20191207235116p:plain
オリオンをなぞる

何度かパシャパシャ撮ってから帰りましたが,やっぱり後でビューワーで視るとまだ少し物足りない感じがします。

そこで,Google Pixel4でも使われている多重露光を行って合成してより精細な写真を撮るという手法を試してみたくなりました。

以下がPixel4で撮った写真の一例らしいです。カッコよすぎんか。私も撮りたい。 https://image.itmedia.co.jp/news/articles/1910/17/ki_1609376_pixel4camera01.jpg

戦略

暗い夜空をきれいに取るためには露光時間を長くすればよいのですが,長すぎると星々が動いてしまいます。

Google Pixel4の戦略としては短めの露光を複数回撮って,それらを重ね合わせてノイズを除去しているようです。

japanese.engadget.com

ということで私も以下のような流れをとることにしました。

  • 360°カメラで星空を複数回撮影
  • 複数の画像を合成してノイズを落とす
  • 星空の回転に合わせて画像を補正

諸々の処理のプログラムはPythonで書きます。

Step0: 星空を撮る

撮影に使っているのはRICOHの全天球カメラのTHETA_SCです。

RICOH 360度カメラ RICOH THETA SC (ブルー) 全天球カメラ 910743

RICOH 360度カメラ RICOH THETA SC (ブルー) 全天球カメラ 910743

  • 発売日: 2016/10/28
  • メディア: エレクトロニクス

扱いやすく個人的には360度カメラの入門にはもってこいと思います。 撮像素子はデジカメでは標準的な 1/2.3型,レンズのF値2.0とまぁまぁの明るさなので15秒程度露光すれば都心の夜空は十分に綺麗に撮れます。

都市郊外で空を撮る感覚だと 露光20-30秒 ,ISO:100-160程度が向いていそうです。

当日は小型の三脚に固定して撮影しました。

全天球カメラで撮った画像はEquirectalgular (正距円筒図法)という形式で保存されます。

f:id:ossyaritoori:20191208001928p:plain
スマホのパノラマで一周して撮ったときも同じような画像になりますね。


本稿ではこの画像の中心を原点とします。

横方向:水平方向へは一周で2π,縦方向は緯度方向なので一周でπのスケールがあります。 したがってざっくり画像座標→正距円筒座標の変換を書くと以下のようになります。

hei,wid,_ = image.shape
cx,cy = wid/2,hei/2
Equilecpoints = (np.array(polarpoints)-np.array([cx,cy]))/np.array([wid/2/pi,-hei/pi])

Step1:多重露光を合成する

ということで多重露光をした画像を合成する段階に入ります。 本来なら露光やISOを統一するべきなのですが,合成を思いついたのが撮影会後なので設定がバラバラな画像になってしまいました。

何も考えずにピクセル輝度の平均をとると以下のようにぼやけた画像になってしまいます。

f:id:ossyaritoori:20191210093417p:plain
建物と空の境界が曖昧に白っぽくなってしまう

そこでOpenCVにあるHDR合成(ハイダイナミックレンジ合成)のライブラリを使用してみることにしました。

HDR合成とは簡単に言うと露光条件の異なる画像を組み合わせて白飛び黒飛びを軽減した画像を撮像することです。

丁寧にまとめてくれている方がこちら(文献[1]): [CV]High-Dynamic-Range(HDR) Imagingについて - Qiita

今回は最も楽そうなMergeMertensという手法を用いました。 Pythonで書くなら以下のような感じです。

import cv2
merge=cv2.createMergeMertens()
proc1 =merge.process([img1,img2,img3])

f:id:ossyaritoori:20191210112502p:plain
MergeMertensを用いた複数画像の合成結果


先ほどのようなモヤが生じる問題を回避できましたが,よく見ると星の動きが軌跡となっています。これはこれで乙ですがくっきりとした星空が見たいですね。

実際オリオン座の周辺に着目すると,星空が動いているのがよくわかります。

f:id:ossyaritoori:20191208001147g:plain
短い間にも空は動いている

露光を繰り返した3分強の間にも星々はよく動いていることがわかります。

Step 2 星空の回転を補正する

先程のGIF画像の様に星空は一日で天球を一回転します。回転中心は北極星Polaris)近辺の極北にあたります。

OpenCVの既存関数による位置ずれ補正

OpenCVにもAlignMTBというクラスがあり,calculateShiftという関数で位置ずれを補正できるそうですがプログラムの中身を見た感じこれは微小並進量のみの対応です。

また,一般の平面画像であれば射影変換によって回転ズレを補正できます。 ossyaritoori.hatenablog.com

しかし,Equirectangular画像においては回転操作は非線形変換になります。

Equirectangular画像における三次元回転

回転操作がどの様に影響するか,次の座標変換から導出することができます。

f:id:ossyaritoori:20191210112100p:plain
3D直交座標と球面極座標の関係。文献[2]から引用

逆変換は以下の通りです。


\begin{align}
x = cos \phi cos \theta, \\
y = cos \phi sin \theta, \\
z = sin \phi
\end{align}

星空はrが変化しない無限遠にあるので,変換後のEquirectangular画像の \theta’,\phi’は天球上での三次元回転行列を用いて


\begin{pmatrix}
 x' \\ y' \\ z'
\end{pmatrix}
= R\begin{pmatrix}
 x\\ y \\  z
\end{pmatrix},  \  \theta' = tan^{-1}\frac{x’}{y'} ,\ \phi'=sin^{-1}z'

と計算できます。

実際変形する際には (\theta’,\phi’)=(i,j)に対応する元画像の座標を逆変換から求め,ピクセル値をBilinear補間で求めます。

ロドリゲスの回転行列の計算

回転行列はOpenCVでサポートされているロドリゲスの回転行列を用います。


星空は原点と極北≒北極星を通る回転軸の反時計回りに回転します。1秒あたり1/240度なので微小のように思われますが,先程のGIFのように複数回露光をすると動きとして効いてきます。

変化を補正するために手動で極北の座標を決定した後,画像のExif情報から撮像の間隔を取得して回転変換をかけます。

# polarw :Equirectangularでの極北のtheta,phi
nx = np.cos(polarw[1])*np.cos(polarw[0])
ny = np.cos(polarw[1])*np.sin(polarw[0])
nz = np.sin(polarw[1])

# omega = time_diff/3600/24*2*pi
Rrod,_ = cv2.Rodrigues(omega*np.array([nx,ny,nz]))

このRrodは先程の式における回転行列Rの転置になります。

変換&重ね合わせの結果

ということで重ね合わせの結果です。

f:id:ossyaritoori:20191210161026p:plain
星空が合うようにマッチングしたので建物などがぶれて見える。

重ね合わせが成功して満足してしまいましたが,見えないものが見えるわけではないのでもっと暗い所に行って多重露光合成の威力を試したくなりました。 まずは休暇が欲しい…

STEP3 :自動で極北を抽出(途中)

SIFTやAKAZE特徴点でマッチングした画像は案の定,空と建物の境界にハマります。 星々をマッチングできれば回転を検出できますが,やはり空と建物・地面の分離を事前にやる必要がありそうです。

f:id:ossyaritoori:20191210212145p:plain
案の定有効な特徴点は空との境界に寄るため星空間の対応は難しそう。

夜空部分の抽出

先程のマッチングでは信頼度の低いマッチング結果を除去していました。具体的には距離情報で縛りをかけています。

良いマッチングを得るためにユークリッド距離の比を元に選別するコードがよく書かれています。

for m,n in self.matchings: 
   if m.distance < 0.75 * n.distance:
         #good.append()        

この選別部分をなくした場合のマッチング点(AKAZE)はいくつか空の星々を捉えています。

f:id:ossyaritoori:20191213110200p:plain
青:画像1のAKAZE特徴点,オレンジ:画像2のAKAZE特徴点(マッチング実施後の対応点)

したがって,最初の選別有りのマッチングにて建物と空の協会となるY座標を抽出し,(今回はY=619)Y座標がその閾値以下,すなわちカメラから見てより上方にある箇所を空の領域と仮定しました。 その結果絞り込まれたのが以下の17点の対応です。

f:id:ossyaritoori:20191213112316p:plain
Y座標を元に夜空の部分の特徴点のみを抽出した

3次元座標での回転推定

Equirectangular平面内ではお互いの変換は非線形になってしまうので,一度求めた特徴点の対応を三次元球面極座標写像してから対応する回転行列を求めます。

三次元回転の最適化といえば金谷先生の教科書を参照するのが良いでしょう。

3次元回転: パラメータ計算とリー代数による最適化

3次元回転: パラメータ計算とリー代数による最適化

また,一応英語のテキストが無料で転がっていたりします。

https://igl.ethz.ch/projects/ARAP/svd_rot.pdf

 \sum Rx_i-y_iを最小化する回転Rは対応点x_iとy_iを列方向に並べた行列X,Yを用いて表した

 S=XY^\top特異値分解して得た U\Sigma V^\top=Sから,

 \hat{R} = V diag(1,1,..,det(VU^\top)) U^\top

のように計算できます。なお,ここで推定したRはデータによっては厳密な回転行列にならない場合があります。それを補正するために同様の操作を繰り返します。 こちらは金谷先生の教科書に記述があるのでそちらを確認してください。

numpyでは以下のように書きます。np.linalgのsvdがVを転置の状態で出すことには注意です。

def estimate3dR(X,Y):
    S=np.dot(X,Y.T)
    U,W,V=np.linalg.svd(S)
    Rhat_ = np.dot(V.T,U.T)
    Rhat = np.linalg.multi_dot([V.T,np.diag([1,1,np.linalg.det(Rhat_)]),U.T])
    # additional Step For Optimization
    U,W,V=np.linalg.svd(Rhat)
    Rhat_ = np.dot(V.T,U.T)
    Rhat = np.linalg.multi_dot([V.T,np.diag([1,1,np.linalg.det(Rhat_)]),U.T])
    return Rhat

ランダムサンプルによるロバスト推定

ところで上記の手法でそのまま回転を計算すると,外れ値の影響を受けてしまいます。

この外れ値を除去するためにランダムサンプルによる手法が数多く提案されており,有名なRANSACなどはopenCVでの推定アルゴリズムにも採用されています。

二次文献ですが以下のブログがわかりやすいと思います。

www.sanko-shoko.net

RANSACの完全な実装はなかなかに面倒で,以下のブログを参考にしました。 クラスの継承を有効に使うコードを初めて書いた気がします。これはこれで長いので別記事で。

一目で分かるRANSAC - Qiita

その結果,

  • 手動設定の回転行列
array([[ 0.99989064, -0.00811519, -0.01236295],
       [ 0.00808042,  0.99996327, -0.00285925],
       [ 0.0123857 ,  0.00275904,  0.99991949]])
  • RANSACで推定した回転行列
array([[ 0.99994968, -0.00893819, -0.00455381],
       [ 0.00895058,  0.99995627,  0.00270865],
       [ 0.0045294 , -0.00274927,  0.99998596]])

はい。そこそこまともな値が抽出できてきました。 観測にノイズが混じっているのでこれは厳密には回転行列ではないので,補正が必要になります。(TODO:もう1つSubsectionが必要)

とりあえず,これをこのサイトにて分解して回転軸と回転角(rad)を抽出してみると,

[ -0.2625027, -0.4368636, 0.8603735 ], 0.0103961

これは,約143秒の露光という判定です…うーん,一分ほど短い? 北極星の位置も微妙にずれた検出をしているのでマッチングの精度が悪かったかRANSACの調整ミスですね。(TODO)

また,回転中心の位置推定もやはり手動には敵いません。これはなかなか手ごわいです。

f:id:ossyaritoori:20191214002406p:plain
白:本来の北極星,黄色:自動推定した空の回転中心


とまぁこれ以上詳細に書くとめちゃ長くなるので別記事にします。 とりあえず,

  • 星空の領域抽出
  • 星空の運動の推定

がネックとなるため最初持っていたイメージよりかは簡単には自動化できないという感覚を得ました。 手動で北極星さえ指定すればExif情報で基本的に推定できるのでこの頑張りがそもそも必要かというのも少し怪しげですが気力があれば続きも記事にします。

TODO

  • 空の領域とそうでない領域の分離:雑にDone
  • 回転中心の自動抽出(特徴点マッチングを用いる):雑にDone
  • 三次元回転推定まわりを整形:多分やる
  • 公開ようにプログラムを整形:やる?

おまけプログラム's

PILでgif画像を保存する関数

from PIL import Image
def saveasGif(images,filename='out.gif',duration=500,loop=0):
    imlist = []
    for img in images:
        img_array = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        imlist.append(Image.fromarray(img_array))
        
    imlist[0].save(filename, save_all=True, append_images=imlist[1:], duration=duration,loop=loop)

Exifデータから撮影間隔を取得する関数

from PIL.ExifTags import TAGS
from datetime import datetime

# Exif データを取得
def get_exif_of_image(file):
    """Get EXIF of an image if exists.
    指定した画像のEXIFデータを取り出す関数
    @return exif_table Exif データを格納した辞書
    """
    im = Image.open(file)

    try:
        exif = im._getexif()
    except AttributeError:
        return {}

    exif_table = {}
    for tag_id, value in exif.items():
        tag = TAGS.get(tag_id, tag_id)
        exif_table[tag] = value
    return exif_table

def get_diff_time(file1,file2):
    """Calculate capture interval of images from  EXIF
    指定した2つの画像ファイル間の撮影時間の間隔を計算する
    """
    FMT = '%H:%M:%S'
    exif1=get_exif_of_image(file1)['DateTime']
    exif2=get_exif_of_image(file2)['DateTime']
    dtime = datetime.strptime(exif2[-8:], FMT) - datetime.strptime(exif1[-8:], FMT)
    return min(3600*24-dtime.seconds,dtime.seconds)

多くなりすぎたのでStep3に関しては別記事に移行する気がします。

参考文献等

[1] HDRなどに関して:[CV]High-Dynamic-Range(HDR) Imagingについて - Qiita

[2] 球面極座標に関する図表などの引用:360度パノラマ画像から平面画像への変換手法 - LASTMILE WORKS / DYNAMO TECH - R&D Project - Medium

[3] パノラマ画像から矩形画像を切り取るプログラム:GitHub - NitishMutha/equirectangular-toolbox: Handy tool for equirectangular images

他余裕でき次第追記

VSCodeでコマンドをショートカットから実行する設定方法 -docmuteで個別コンパイルができない問題の解決-

TexはずっとTexstudioで書いてきましたが,docmuteで分割コンパイルするときに困ったのでVScodeも併用します。

VScodeでコマンドをショートカットから実行する

ほぼ以下の記事の通りです。

qiita.com

keybindings.jsonを開く

まずは,設定を書き込むkeybindings.jsonを開きます。

コマンドパレット(ctrl+shift+p)を開き'key'と入力し「基本設定:キーボードショートカットを開く(JSON)」を選びます。

設定の書き方

開いたjsonファイルは[]で囲まれた記述が連なっていると思います。 下のような表記を中括弧[]の中に追記します。

    {
        "key": "F5",
        "command": "workbench.action.terminal.sendSequence",
        "args":{
            "text": "echo ${file}\n"
        }
    }

これはworkbench.action.terminal.sendSequenceというコマンドを介してtext: に連なるターミナルのコマンドを実行しています。

ここでいうターミナルはVSCodeの下にあるターミナルのことで,実行結果や経過を確認することができます。

f:id:ossyaritoori:20191117175325p:plain
Windowsの場合はPowershellを起動している。UbuntuとかならBashかな?

引数について

textの引数には以下のURLにあるような変数が使用できます。

例えば上の例の${file}は現在VSCode上で編集しているファイルのフルパスとファイル名です。

Visual Studio Code Variables Reference

自分の設定

後半で解説する問題を解決するために自分が書いたファイルは以下のとおりです。

    {
        "key": "F5",
        "command": "workbench.action.terminal.sendSequence",
        "args":{
            "text": "cd ${fileDirname}\n platex ${fileBasenameNoExtension}.tex \n dvipdfmx ${fileBasenameNoExtension}.dvi \n SumatraPDF ${fileBasenameNoExtension}.pdf \n"
        }
    }

「自分の今開いているファイルのディレクトリに移動」「platexとdvipdfmxでPDF作成」「SumatraPDFで表示」という順になっています。

docmuteパッケージの個別コンパイルができない問題

docmuteとは端的に言うとTexファイルの分割コンパイルのためのパッケージです。

ossyaritoori.hatenablog.com

以下のような構造で,mainがcahpter1,2をincludeないしinputしているときにmainでの全体コンパイルとchapter1単体のコンパイルをサポートしてくれます。

root/
 ├ main.tex
 ├ chapter/
   ├ chapter1.tex
   └ chapter2.tex

確認している問題

TexStudioのオートコンパイルの機能を使っている前提ですが, 他のファイルをコンパイルする前にmainをコンパイルした場合,chapter1のファイルを開いてコンパイルしても間違って全体がコンパイルされてしまう現象が起きてしまいます。


おそらく,main.texコンパイルで生成されたファイルの依存関係を参照してしまうためと思われ,ターミナルで必要なコマンドを直打ちするか,一度中間生成ファイルをすべて取り払うことで想定した動作を実現できます。

VSCodeを使った解決

TexStudioで編集しているファイルをVSCode同時に開いて全体コンパイルと個別コンパイルでエディタを分けて対応する様にしてやろうというのが今回の解決策です。

この解決策は最適ではきっとないでしょうが,VSCodeとは今後長いこと付き合っていきそうなのでできそうなことはVSCodeにやってもらおうかなと思っています。