ぷろぐ((>ω<))

ぷろぐらみんぐ関係のメモ

Raspberry Pi 4B の GPU で OpenCL (VC4CL)

RaspberryPiの映像出力用のGPUOpenCLから利用できると知ったので、興味本位でトライ。
結果から言うと、行列演算サンプルを実行したところ、GPU処理よりCPU処理の方が若干早かった…。

材料

    • 64bitのつもりだったが32bitのRaspbian入れてた
$ getconf LONG_BIT
32

$ lsb_release -a
No LSB modules are available.
Distributor ID: Raspbian
Description:    Raspbian GNU/Linux 10 (buster)
Release:        10
Codename:       buster

参考サイト

  1. Raspberry PiのGPUをOpenCLで使う方法
  2. Raspberry PiでOpenCLを触ってみた
    • これらのサイトでは Raspberry Pi Zero W や Raspberry Pi 2B で試している
    • また、debパッケージをインストールしているが、私の環境では例のおまじない sudo apt --fix-broken install を実行しても全然ダメだった
    • OpenCLのテスト実行のリポジトリ紹介あり。本記事でも試す。

ビルド&インストール

というわけで、手元で必要物をビルドする。
具体的には下記サイトの"VC4CL installation"のセクションの通り実行する。

確認

$ sudo clinfo

Number of platforms                               1
  Platform Name                                   OpenCL for the Raspberry Pi VideoCore IV GPU
  Platform Vendor                                 doe300
  Platform Version                                OpenCL 1.2 VC4CL 0.4
  Platform Profile                                EMBEDDED_PROFILE
  Platform Extensions                             cl_khr_il_program cl_khr_spir cl_altera_device_temperature cl_altera_live_object_tracking cl_khr_icd cl_vc4cl_performance_counters
  Platform Extensions function suffix             VC4CL

  Platform Name                                   OpenCL for the Raspberry Pi VideoCore IV GPU
Number of devices                                 1
  Device Name                                     VideoCore IV GPU
  Device Vendor                                   Broadcom
  Device Vendor ID                                0xa5c
  Device Version                                  OpenCL 1.2 VC4CL 0.4
  Driver Version                                  0.4
  Device OpenCL C Version                         OpenCL C 1.2
  Device Type                                     GPU
  Device Profile                                  EMBEDDED_PROFILE
  Device Available                                Yes
  Compiler Available                              Yes
  Linker Available                                Yes
  Max compute units                               1
  Max clock frequency                             500MHz
  Device Partition                                (core)
    Max number of sub-devices                     0
    Supported partition types                     None
    Supported affinity domains                    (n/a)
  Max work item dimensions                        3
  Max work item sizes                             12x12x12
  Max work group size                             12
  Preferred work group size multiple              1
  Preferred / native vector sizes
    char                                                16 / 16
    short                                               16 / 16
    int                                                 16 / 16
    long                                                 0 / 0
    half                                                 0 / 0        (n/a)
    float                                               16 / 16
    double                                               0 / 0        (n/a)
(以下略)

CPU vs GPU 速度比較

vc4clを触るときはsudoが必要

$ git clone git://github.com/HandsOnOpenCL/Exercises-Solutions.git
$ cd Exercises-Solutions/Solutions/Exercise06/C
$ make
$ sudo ./mult
Using OpenCL device: VideoCore IV GPU

===== Sequential, matrix mult (dot prod), order 1024 on host CPU ======
 25.80 seconds at 83.3 MFLOPS

===== OpenCL, matrix mult, C(i,j) per work item, order 1024 ======
 30.00 seconds at 71.6 MFLOPS

 Errors in multiplication: 246368803225600.000000
terminate called without an active exception
中止
  • CPUの方が若干早かった。
  • 途中で落ちてしまったが、原因究明未。

更に調べていると

Raspberry Pi 4 のGPUについて言及されている記事を発見。

RPi4になってGPUが変わったようですが (VideoCore VI, VC6)、本記事で入れたライブラリはRPi3以前向け (VideoCore IV, VC4) ということでした。
VC6向けのライブラリが待たれる…。

Google Colab 上で Keras 入門(SegNet-Basic を実装・学習・推論)

4年前ぐらいに Deep Learning 入門をしてからしばらく触っていなかったが、コロナ外出自粛で時間が余っていることもあって久しぶりにトライ。
Google Colab でやる手順をメモ。試したのはSegNetの軽量版のSegNet-Basic。

  1. SegNet: 画像セグメンテーションニューラルネットワーク - Qiita
  2. GitHub - 0bserver07/Keras-SegNet-Basic: SegNet-Basic with Keras

(1)でSegNet-Basicの情報源として(2)にリンクを貼っているが、(2)のHow-to が書きかけだったので、(2)のページを参考にしながら下の手順でKerasで実装したところ、うまくいった。

ちなみに、何故SegNetではなくてSegNet-Basicかというと、実はSegNetも試してみたが、Google Colabではリソースの制約上クラッシュして仕方がなかったので諦めただけ。

環境確認

以下の記事は、下記の環境で実行した。

import tensorflow as tf
import keras

print("tf.__version__ is", tf.__version__)
print("tf.keras.__version__ is ", tf.keras.__version__)
print("keras.__version__ is ", keras.__version__)
Using TensorFlow backend.
tf.__version__ is 2.2.0-rc4
tf.keras.__version__ is  2.3.0-tf
keras.__version__ is  2.3.1

下準備

ノートブックを新規作成

  • 「ファイル」 > 「ノートブックを新規作成」
  • ここではSegNetBasic.ipynbとした

Google Drive マウント

  1. Google DriveGoogle Colab にマウント
    • 左端の「ファイル」から操作できる
    • ↓こういう状態になるはず

f:id:presan:20200508144101p:plain
Google Drive をマウント

  1. 下記で公開されているSegNetリポジトリ内のデータ(CamVidフォルダ以下)をGoogle Drive にコピー。
    • 今回は drive/My Drive/Colab Notebooks/にSegNetというフォルダを作り、その中にCamVidをコピーした。
  2. SegNetフォルダ下に次のフォルダを作成
    • data
    • models
    • weights
    • valid  (おまけ)

モデル準備

SegNetフォルダ下に下記をmodel.pyとして保存。

import keras.models as models
from keras.models import Model
from keras.layers import Input
from keras.layers.core import Layer, Dense, Dropout, Activation, Flatten, Reshape, Permute
from keras.layers.convolutional import Convolution2D, MaxPooling2D, UpSampling2D, ZeroPadding2D
from keras.callbacks import ModelCheckpoint

def buildSegnetBasicModel(input_shape, n_labels, kernel=3, pool_size=(2, 2), pad=1, output_mode="softmax"):
    # encoder
    inputs = Input(shape=input_shape)

    conv_1 = ZeroPadding2D(padding=(pad,pad))(inputs)
    conv_1 = Convolution2D(64, (kernel, kernel), padding="valid")(conv_1)
    conv_1 = BatchNormalization()(conv_1)
    conv_1 = Activation("relu")(conv_1)

    pool_1 = MaxPooling2D(pool_size)(conv_1)

    conv_2 = ZeroPadding2D(padding=(pad,pad))(pool_1)
    conv_2 = Convolution2D(128, (kernel, kernel), padding="valid")(conv_2)
    conv_2 = BatchNormalization()(conv_2)
    conv_2 = Activation("relu")(conv_2)

    pool_2 = MaxPooling2D(pool_size)(conv_2)

    conv_3 = ZeroPadding2D(padding=(pad,pad))(pool_2)
    conv_3 = Convolution2D(256, (kernel, kernel), padding="valid")(conv_3)
    conv_3 = BatchNormalization()(conv_3)
    conv_3 = Activation("relu")(conv_3)
    
    pool_3 = MaxPooling2D(pool_size)(conv_3)

    conv_4 = ZeroPadding2D(padding=(pad,pad))(pool_3)
    conv_4 = Convolution2D(512, (kernel, kernel), padding="valid")(conv_4)
    conv_4 = BatchNormalization()(conv_4)
    conv_4 = Activation("relu")(conv_4)

    print("Build SegNet-Basic enceder done..")

    # decoder
    conv_5 = ZeroPadding2D(padding=(pad,pad))(conv_4)
    conv_5 = Convolution2D(512, (kernel, kernel), padding="valid")(conv_5)
    conv_5 = BatchNormalization()(conv_5)

    unpool_1 = UpSampling2D(pool_size)(conv_5)

    conv_6 = ZeroPadding2D(padding=(pad,pad))(unpool_1)
    conv_6 = Convolution2D(256, (kernel, kernel), padding="valid")(conv_6)
    conv_6 = BatchNormalization()(conv_6)
    
    unpool_2 = UpSampling2D(pool_size)(conv_6)

    conv_7 = ZeroPadding2D(padding=(pad,pad))(unpool_2)
    conv_7 = Convolution2D(128, (kernel, kernel), padding="valid")(conv_7)
    conv_7 = BatchNormalization()(conv_7)

    unpool_3 = UpSampling2D(pool_size)(conv_7)

    conv_8 = ZeroPadding2D(padding=(pad,pad))(unpool_3)
    conv_8 = Convolution2D(64, (kernel, kernel), padding="valid")(conv_8)
    conv_8 = BatchNormalization()(conv_8)

    conv_9 = Convolution2D(n_labels, (1, 1), padding="valid")(conv_8)
    conv_9 = Reshape(
        (input_shape[0] * input_shape[1], n_labels),
        input_shape=(input_shape[0], input_shape[1], n_labels),
    )(conv_9)

    outputs = Activation(output_mode)(conv_9)
    
    print("Build SegNet-Basic decoder done..")

    model = Model(inputs=inputs, outputs=outputs, name="SegNetBasic")

    return model

データ準備

次を実行して、画像データをnumpyバイナリ形式(.npy)に変換。SegNet/data内にnpyファイルが6個、合計4.5GB程度生成されているはず。

import cv2
import numpy as np

from helper import *

import os
import gc

# Copy the data to this dir here in the SegNet project /CamVid from here:
# https://github.com/alexgkendall/SegNet-Tutorial
RootPath = 'drive/My Drive/Colab Notebooks/SegNet'
DataPath = 'drive/My Drive/Colab Notebooks/SegNet/CamVid/'
data_shape = 360*480

def normalized(rgb):
    return rgb / 255.

def one_hot_it(labels):
    w, h = labels.shape[:2]
    x = np.zeros([w,h,12], dtype=np.uint8)
    for i in range(w):
        for j in range(h):
            x[i,j,labels[i][j]]=1
    return x

def load_data(mode):
    data = []
    label = []
    with open(DataPath + mode +'.txt') as f:
        txt = f.readlines()
        txt = [line.split(' ') for line in txt]
    for i in range(len(txt)):
        datapath = RootPath + txt[i][0][7:]
        print('(', i, '/', len(txt), ') Loading data: ', datapath)
        img = cv2.imread(datapath)
        data.append(np.rollaxis(normalized(img),2))

        labelpath = RootPath + txt[i][1][7:][:-1]
        print('(', i, '/', len(txt), ') Loading label: ', labelpath)
        img = cv2.imread(labelpath)
        label.append(one_hot_it(img[:,:,0]))
    return np.array(data), np.array(label)



train_data, train_label = load_data("train")
train_label = np.reshape(train_label,(367,data_shape,12))
np.save(RootPath + "/data/train_data", train_data)
np.save(RootPath + "/data/train_label", train_label)
del train_data
del train_label
gc.collect()

test_data, test_label = load_data("test")
test_label = np.reshape(test_label,(233,data_shape,12))
np.save(RootPath + "/data/test_data", test_data)
np.save(RootPath + "/data/test_label", test_label)
del test_data
del test_label
gc.collect()

val_data, val_label = load_data("val")
val_label = np.reshape(val_label,(101,data_shape,12))
np.save(RootPath + "/data/val_data", val_data)
np.save(RootPath + "/data/val_label", val_label)
del val_data
del val_label
gc.collect()

学習

モデルを作成

次を実行

# Reference:
# https://qiita.com/cyberailab/items/d11862852eccc17585e8
# https://github.com/0bserver07/Keras-SegNet-Basic

import keras.models as models
from keras.models import Model
from keras.layers import Input
from keras.layers.core import Layer, Dense, Dropout, Activation, Flatten, Reshape, Permute
from keras.layers.convolutional import Convolution2D, MaxPooling2D, UpSampling2D, ZeroPadding2D
from keras.layers.normalization import BatchNormalization
from keras.callbacks import ModelCheckpoint

import cv2
import numpy as np
import matplotlib.pyplot as plt
import time

# Import SegNet/SegNetBasic modules
import sys
RootPath = '/content/drive/My Drive/Colab Notebooks/SegNet/'
sys.path.append(RootPath)
from model import buildSegnetModel, buildSegnetBasicModel

# Start time
start_time = time.time()

# Fix seed for reproducibility
np.random.seed(0)

# Parameters
class_weighting= [0.2595, 0.1826, 4.5640, 0.1417, 0.5051, 0.3826, 9.6446, 1.8418, 6.6823, 6.2478, 3.0, 7.3614]
input_shape = (360, 480, 3)
n_labels =  12
kernel = 3
pool_size = 2
pad = 1
output_mode = 'softmax'
data_shape = input_shape[0] * input_shape[1]

# load the model:
print("Building model...")
model_segnet = buildSegnetBasicModel(input_shape, n_labels, kernel, pool_size, pad, output_mode)
print("done.")

print("Compiling model...")
model_segnet.compile(loss="categorical_crossentropy", optimizer='adadelta', metrics=["accuracy"])
print("done.")

# Visualize model
model_segnet.summary()

# Calculate erapsed time
model_load_time = time.time()
print('Erapsed time: ', (model_load_time - start_time), '[s]')

npyデータを読み込み

次を実行。5~10分程度時間がかかる。

# load the data
print("Loading data...")
print("- Loading train_data...")
train_data = np.load(RootPath + 'data/train_data.npy').transpose((0, 2, 3, 1)) # NCHW to NHWC
print("  -> shape", train_data.shape)
print("- Loading train_label...")
train_label = np.load(RootPath + 'data/train_label.npy')
print("  -> shape", train_label.shape)
print("- Loading test_data...")
test_data = np.load(RootPath + 'data/test_data.npy').transpose((0, 2, 3, 1)) # NCHW to NHWC
print("  -> shape", test_data.shape)
print("- Loading test_label...")
test_label = np.load(RootPath + 'data/test_label.npy')
print("  -> shape", test_label.shape)
print("done.")

# Calculate erapsed time
data_load_time = time.time()
print('Erapsed time: ', (data_load_time - model_load_time), '[s]')

学習開始

ランタイム選択

  • 「ランタイム」 > 「ランタイムのタイプを変更」 で GPU or TPU を選ぶ

学習開始

次を実行。

  • 途中、標準割り当てメモリ(12.5[GB])ではメモリ不足でクラッシュするかも。
    • その場合は2倍の25[GB]をGoogleから有難く割り当ててもらって再実行。
    • 25[GB]あれば大丈夫
  • そのときの割り当てリソースにもよるが、ある日のGPU選択時は1epoch 110[s]程度だった。
    • つまり100epochの学習に2[h]程度
    • batch_sizeを2倍の12にすると1epoch 30[s]で100epoch成功するときもあれば、ResourceExhaustedErrorでクラッシュすることもあった
  • 過去最高の結果になるたび、SegNet/weightsフォルダに重みを上書き保存。(checkpoint)
# Parameter
nb_epoch = 100
batch_size = 6

# checkpoint
print("Deifining callbacks...")
filepath = RootPath + "weights/segnet_weights.best.hdf5"
checkpoint = ModelCheckpoint(filepath, monitor='val_accuracy', verbose=1, save_best_only=True, mode='max')
callbacks_list = [checkpoint]
print("done.")

# Calculate erapsed time
def_cb_time = time.time()
print('Erapsed time: ', (def_cb_time - data_load_time), '[s]')
print('-----------------')

# Fit the model
print("Fitting model...")
hist = model_segnet.fit(train_data, train_label, callbacks=callbacks_list, batch_size=batch_size, epochs=nb_epoch,
                    verbose=2, class_weight=class_weighting, validation_data=(test_data, test_label), shuffle=True) # validation_split=0.33
print("done.")

# Calculate erapsed time
fit_model_time = time.time()
print('Erapsed time: ', (fit_model_time - def_cb_time), '[s]')
print('-----------------')

モデル・重み保存

  • 次を実行して、推論時に使うためにモデルと100epoch目の重みをファイルに保存。SegNetフォルダ下のmodels, weightsフォルダにそれぞれ保存。
  • ただ、重みは100epoch目よりも、過去最高(checkpoint)の結果を使うべき
    • checkpointと100epoch目の重みファイルのサイズがかなり違うのは何故だろう?
# This save the trained model weights to this file with number of epochs
print("Saving model and weights...")
model_segnet.save(RootPath + 'models/segnet_model.hdf5')
model_segnet.save_weights(RootPath + 'weights/segnet_weight_{}.hdf5'.format(nb_epoch))
print("done.")

# Calculate erapsed time
fit_model_time = time.time()
print('Erapsed time: ', (fit_model_time - data_load_time), '[s]')
print('-----------------')

** loss/accuracyの可視化
次を実行してloss/accuracyの推移をグラフ描画。

>|python|
# Visualize
epochs = range(1, len(hist.history['accuracy']) + 1)

plt.plot(epochs, hist.history['loss'], label='Training loss', ls='-')
plt.plot(epochs, hist.history['val_loss'], label='Validation loss')
plt.title('Training and Validation loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()
plt.show()

plt.plot(epochs, hist.history['accuracy'],  label='Training acc')
plt.plot(epochs, hist.history['val_accuracy'], label="Validation acc")
plt.title('Training and Validation accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.show()
学習結果

validationの精度がmaxで約80%ぐらいになった

f:id:presan:20200508141014p:plain
Loss
f:id:presan:20200508141032p:plain
Accuracy

推論

モデルと重みを読み込み

次を実行

from keras.models import load_model
from google.colab import files
from PIL import Image
import cv2
import numpy as np
import time

RootPath = '/content/drive/My Drive/Colab Notebooks/SegNet/'

# Parameter
input_shape = (360, 480, 3)

# Start time
start_time = time.time()

# Load model and weights
print('Loading model and weights...')
model_segnet = load_model(RootPath + 'models/segnet_model.hdf5')
model_segnet.load_weights(RootPath + 'weights/segnet_weights.best.hdf5')
print('done')

# Calculate erapsed time
model_load_time = time.time()
print('Erapsed time: ', (model_load_time - start_time), '[s]')

検証用データを読み込み

次を実行

val_data = np.load(RootPath + '/data/val_data.npy').transpose((0, 2, 3, 1))  # NCHW to NHWC
val_label = np.load(RootPath + '/data/val_label.npy')
print('done')

評価

次を実行

batch_size = 12

# estimate accuracy on whole dataset using loaded weights
scores = model_segnet.evaluate(val_data, val_label, verbose=0, batch_size=batch_size)
print("%s: %.2f%%" % (model_segnet.metrics_names[1], scores[1]*100))
評価結果
accuracy: 87.45%

可視化

次を実行。
検証用データ全部にやっては大変なので、最大10個まで可視化するようにした。

import matplotlib.pyplot as plt

Sky = [128,128,128]
Building = [128,0,0]
Pole = [192,192,128]
Road_marking = [255,69,0]
Road = [128,64,128]
Pavement = [60,40,222]
Tree = [128,128,0]
SignSymbol = [192,128,128]
Fence = [64,64,128]
Car = [64,0,128]
Pedestrian = [64,64,0]
Bicyclist = [0,128,192]
Unlabelled = [0,0,0]

label_colours = np.array([Sky, Building, Pole, Road, Pavement,
                          Tree, SignSymbol, Fence, Car, Pedestrian, Bicyclist, Unlabelled])

def visualize(temp, plot=True):
    r = temp.copy()
    g = temp.copy()
    b = temp.copy()
    for l in range(0,11):
        r[temp==l]=label_colours[l,0]
        g[temp==l]=label_colours[l,1]
        b[temp==l]=label_colours[l,2]

    rgb = np.zeros((temp.shape[0], temp.shape[1], 3))
    rgb[:,:,0] = (r/255.0)#[:,:,0]
    rgb[:,:,1] = (g/255.0)#[:,:,1]
    rgb[:,:,2] = (b/255.0)#[:,:,2]
    if plot:
        plt.imshow(rgb)
    else:
        return rgb

img_size = (input_shape[0], input_shape[1])

# Start time
start_time = time.time()

# Predict
print('Predicting...')
output = model_segnet.predict(val_data)
print('done')

# Calculate erapsed time
end_time = time.time()
print('Erapsed time: ', (end_time - start_time), '[s]')
print('-----------------')

# Visualize
count = min([10, len(output)])
for i in range(count):
  pred_class = np.argmax(output[i], axis=1).reshape(img_size)
  img_ret = visualize(pred_class, False)
  plt.figure(i * 2)
  plt.imshow(val_data[i])
  plt.figure(i * 2 + 1)
  plt.imshow(img_ret)
plt.show()
可視化結果

それっぽいsegmentation結果になった

f:id:presan:20200508142426p:plain
元画像
f:id:presan:20200508142440p:plain
SegNetBasic結果

(おまけ) 自前の画像で可視化

  • SegNet/validフォルダに480x360のRGB画像をtest0~8.pngという名前で保存
  • 次を実行
    • さっきまでSegNet公開リポジトリの検証用データを読み込んでいた変数val_dataに自前の画像を読み込む
val_data = []
for i in range(9):
  # Load image
  img_path = RootPath + 'valid/test{}.png'.format(i)
  print('Loading ', img_path)
  img = cv2.imread(img_path)
  if img is None:
    print('Failed to load ', img_path)
    continue
  #img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
  print('done')
  val_data.append(img)

val_data = np.array(val_data)
可視化結果

うーん。。

f:id:presan:20200508143222p:plain
元画像
f:id:presan:20200508143234p:plain
SegNetBasic結果

マウスカーソルを見失ったときに所定の場所に強制移動させる

マルチディスプレイしているとマウスカーソルを見失うことも多いが、そのたびにマウスをバタバタと左右に動かして画面のどこにあるか探す時間が勿体ない。

そこで、マウスカーソルが画面上の所定の位置に来るようなbatを作成した。これをpathが通るところにおいて、Winキー→bat名入力で即座にマウスカーソルを取り戻す。

@powershell -NoProfile -ExecutionPolicy Unrestricted "$s=[scriptblock]::create((gc \"%~f0\"|?{$_.readcount -gt 1})-join\"`n\");&$s" %*&goto:eof

add-type -AssemblyName System.Windows.Forms
[System.Windows.Forms.Cursor]::Position = new-object System.Drawing.Point(960, 540)

cv::cornerSubPix (サブピクセル推定) アルゴリズムの詳細

OpenCVにはcv::cornerSubPixという関数があるが、これの原理に関する説明が非常に簡素でわかりづらい。

特徴検出 — opencv 2.2 documentation


上記ページに書いてある3つの式の意味についてまとめる。

  1. \displaystyle{\epsilon_{i}=DI^{T}_{p_{i}}\cdot(q-p_{i})}
  2. \displaystyle{\sum_{i} (DI_{p_{i}}\cdot DI^{T}_{p_{i}})-\sum_{i} (DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot p_{i})}
  3. \displaystyle{q=G^{-1}\cdot b}

式1について

上記OpenCVの説明ページで「各点pに向かう画像勾配と直交する」にあたる部分。これについては下記記事がわかりやすかった。

Fundamentals of Features and Corners: Subpixel Corners: Increasing accuracy - AI Shack

式2について

式1から式2の間がちょっとぶっ飛びすぎているし、いろいろと間違っていると思う。

  • \epsilon_{i}を0とすることで」(英語版では"A system of equations may be set up with \epsilon_{i} set to zero")は違うと思う
  • 突然qが式2から消えてなくなっているが間違いだと思う

式1~式2の間には次の論理が省略されていると思う。一言で言えば「二乗誤差の最小化」である。

詳細

探索窓内のすべての点について式1で定義される\epsilon_{i}の二乗の総和をEとする。

\displaystyle{
E=\sum_{i} \epsilon_{i}^{2}
}

これに式1を代入すると、

\displaystyle{
E=\sum_{i} \left(DI^{T}_{p_{i}}\cdot(q-p_{i})\right)^{2}
}

この総和Eが最小値をとるようなqを求めればよい。それはつまり、qについての微分dE/dqがゼロになるときである。よって、

\displaystyle{
\frac{dE}{dq}=2\sum_{i} \left(DI^{T}_{p_{i}}\cdot(q-p_{i})\right)DI^{T}_{p_{i}}=0
}

(a\cdot b)a=(aa^T)bと書き換えられるので、

\displaystyle{
\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot(q-p_{i})=0
}

これを展開すると

\displaystyle{
\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot q - \sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot p_{i}=0
} (★)

この左辺が式2である。繰り返しになるが、qが書かれてしかるべきだと思う。

式3について

突然1次勾配G、2次勾配bという言葉が出てくるが、次の論理が省略されている。

詳細

上述の式(★)より

\displaystyle{
\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot q = \sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot p_{i}
}

q\displaystyle{\sum_{i}}には依存しないので、

\displaystyle{
\left(\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\right) \cdot q = \sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot p_{i}
}

ここで、

\displaystyle{
G=\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}
}
\displaystyle{
b=\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}\cdot p_{i}
}

と置くと、上述の式は次のように書き換えらえる。

\displaystyle{
G\cdot q=b
}

ゆえに、今求めたいq

\displaystyle{
q=G^{-1}\cdot b
}

となり、式3となる。

おまけ

ところで、G=\sum_{i} DI_{p_{i}}\cdot DI^{T}_{p_{i}}は計算すると近似的に次のようになるが、これを structure tensor matrix と呼ぶ。

DI_{p_{i}}=(I_{x_{p_{i}}}, I_{y_{p_{i}}})^{T}とすると、

\displaystyle{
G=\sum_{i} \begin{pmatrix}
I_{x_{p_{i}}}I_{x_{p_{i}}} & I_{x_{p_{i}}}I_{y_{p_{i}}} \\
I_{x_{p_{i}}}I_{y_{p_{i}}} & I_{y_{p_{i}}}I_{y_{p_{i}}} \\
\end{pmatrix} \approx 
\sum_{i} \begin{pmatrix}
I_{xx_{p_{i}}} & I_{x_{p_{i}}}I_{y_{p_{i}}} \\
I_{x_{p_{i}}}I_{y_{p_{i}}} & I_{yy_{p_{i}}} \\
\end{pmatrix} 
}

C, V, H キーが打てない問題の解消法 (When C, V, H keys don't work...)

TL;DR (English)

  • It might depend on each environment, but in my case, "Right Windows" key is regarded as pressed down.
    • Of course my keyboard doesn't have such a key... So I can't press it.
  • To solve this problem, simulate "R Win" key (VK_RWIN) with some tool.
    • You can create it by yourself using Win32API keybd_event()
  • My guess is that this problem might occur due to the malfunction of keyboards.

TL;DR (Japanese)

  • この問題の発生原因は人それぞれかもしれないが、私の場合は「右Windowsキー」が押しっぱなし状態と判定されていることが原因だった。
    • ちなみに使っているキーボードには「左Windowsキー」しかなく、当然「右Windowsキー」を押した覚えはない。
      • たぶん今時「右Windowsキー」のついたキーボードの方が珍しいはず
    • なおWindowsに標準でインストールされているソフトウェアキーボードも「左Windowsキー」しかないので、解消できない。
  • キーボード入力をシミュレートするツールで「右Windowsキー」を話したことにして対処。
  • (恐らくだが)原因はキーボードの誤動作。無線キーボードを使っているのだが、伝送する信号を間違えたのかな…あくまで想像。

詳細

結構長い間、下記の問題に悩んでいた:

  • テキスト入力欄にC, V, Hキーが打てない
  • C, V, H, Aキーを押すと変な動作をする
  • etc

ネットで調べても同様の症状の人は世界中でいるようである。
検索欄に「c v h」と入力したら他の人も検索しているみたい。

f:id:presan:20190622222311p:plain
c v h suggestion

実際いくつかの質問サイトでQ&Aがあるが、人それぞれ原因は違うようでいろんな対策が載っている。

ちなみに私は検索して出てくるやつでは一向に改善しなかった。
これまでは、しばらく我慢してたらどういうわけか直ったりしていた。
(いつか直るまでどうにも我慢できないときは、一度ログオフしたり再起動したりしていた)

ちょっと真面目に原因を探ろうと思い、何か特殊キーが押された状態になっているのではないかと仮説を立ててキー入力状態を見てみたら、当たりだった。

Direct Sparse Odometry (DSO) の要点まとめ

過去のエントリではWindowsへのインストール方法をまとめたが、今回は論文から理解したことをまとめる。
数学的な細かな話は最小限に抑え、私なりに要点だけまとめる。
# 間違ってるところがあるかもしれないので随時修正していくつもり。

従来手法

従来のSFM/SLAMはざっくり次の3ステップ:

  1. 複数画像から対応する特徴点を抽出
  2. 各対応特徴点に対して三角測量を行い3D再構成
  3. 3D再構成した点やカメラ位置姿勢に対してバンドル調整 (Bundle Adjustment : BA)
    • 3D点の画像上への再投影点と元の2D特徴点のx,yのズレ(再投影誤差)という空間的誤差 (geometric error) が最小化されるように、カメラ位置姿勢や3D点座標の微小修正を繰り返す

DSOの手法

これに対して、DSOは入力画像群に対して最初から直接バンドル調整(BA)を実施しているイメージ。特徴点は求めない。ざっくり次の3ステップ:

  1. 一方の画像の(勾配がある程度大きい)とある座標の画素値(A)はもう一方の画像のとある座標の画素値(B)に対応するはずという計算式を立てる
  2. その計算式から、画素値(A)と(B)の差の計算式を立てる
  3. その画素値誤差 (photometric error) が最小化されるようにカメラ位置姿勢や3D点(*1)の微小修正を繰り返す

ちなみに、処理対象とする画像:キーフレーム(KF)や3D点群を随時不要になったら削除するなど上手に管理するところにもテクニックがありそうだ。

※画像/画質処理屋さんのセンスのままで考えると理解できないかも。DSOは数学の塊で成立しているイメージ。
前回エントリを例にとると、私のイメージでは(あくまで個人の感想):

  • マッチング方式 … 画像処理屋さん的センス
  • オプティカルフロー方式 … 数学屋さん的センス

前提条件

Webカメラスマホカメラのような「人の目に美しく映るような画像」ではなく、「マシンビジョンにとって最適な画像(やカメラ)」であること。具体的には:

  • グローバルシャッター
    • 論文では一言、ローリングシャッターの影響をいくらか緩和する手法もあるようなことも書かれてはいるが…
  • キャリブレーションその1 (geometric calibration) (※言わゆる内部パラメタ)
  • キャリブレーションその2 (photometric calibration)
    • 露光時間
    • レンズ減衰・口径食
    • 放射照度の逆関数

DSOの嬉しいところ

  • 従来の特徴点方式SFM/SLAMは所謂「コーナー」を3D化するに留まるが、DSOは処理方式の違いから「エッジ」も再構成しやすい
  • つまり従来方式より多くの3D点を生成しやすい


今後より技術的に詳しく見ていこうと思う。

画像間の特徴点対応付け~マッチング方式とオプティカルフロー方式~

最終出力として2画像間の特徴点の対応関係を求めたいという場合、その方式には大きく2種類がある。
この観点からまとめられている記事をあまり見た事がないので投稿。

マッチング方式

「特徴点 対応付け」や「特徴点 マッチング」と検索すると引っかかるのは大体この方式のこと。
特徴点は英語では feature (point) と言ったり keypoint と言ったりする。
ブロック図にすると下記のようになる:

マッチング方式のブロック図
BlockDiagram

つまり、処理手順としては

  1. 双方の画像で「特徴点」を見つける
  2. 双方の画像で見つけた特徴点について「特徴量(特徴記述子と言ったりもする)」を求める
  3. 双方の特徴量を比較して一番近しい(似ている)ペアを見つける

という流れ。

特徴点抽出手法や特徴量記述手法にはいろんなものが提案されている。

  • 2つセットで提案されているものとして、例えばSIFT, SURF, ORB, KAZE, AKAZE, etc...
  • 特徴点抽出(検出)手法単体として提案されているものとして、例えばHarrisコーナー, Shi-Tomasiコーナー(=GFTT), FAST, etc...
  • 大抵の特徴点抽出手法や特徴量記述手法はグレースケール画像に対して実施するが、RGB画像をそのまま使用するものもあるようだ。例えば、SIFTに色情報を追加する論文がいくつか見つかる。

処理が理解しやすいように、ブロック図にイメージ図も付け足すとこんな感じ:

マッチング方式のイメージ図
BlockDiagram

ちなみに近しい特徴量が複数ある場合はマッチングに失敗することもある。例えばこんな感じ:

マッチング失敗の例
マッチング失敗の例

このような誤マッチングの回避策として、例えばOpenMVGでは下記のようなチェックがなされている:

誤マッチング回避方法のコンセプト
誤マッチングの回避策

つまり特徴が突出しているわけではなければ、マッチング不成立としよう、という考え方。
他にも、双方向にマッチングしてみて、結果が1対1対応しなければ、棄却する、などという手法もある。
ただ、どんどん計算量が増えていってしまうところが悩みどころ。

オプティカルフロー方式(疎)

ここではLucas-kanade法を扱うこととする。
また、説明の簡単化のためにPyramidの件は省略する。
ブロック図にすると下記のようになる:

疎なオプティカルフロー方式のブロック図
疎なオプティカルフロー方式のブロック図

つまり、処理手順としては

  1. 片方の画像で「特徴点」を見つける
  2. 片方の画像について微分画像」を求める
  3. 片方の画像の特徴点と微分画像ともう片方の画像から、特徴点の移動先を推定する(追跡する)

という流れ。

特徴点を追跡するので英語ではtrackingと言ったりする。matchingしているわけではない。
どうして「こんな処理でもう一方の画像上で対応点が見つかるのか?」という疑問に対する雑な回答としては、
「数学的に計算するとこれでできる」となる。下記でわかりやすく解説されている:
OpenCVでとらえる画像の躍動、Optical Flow - Qiita
オプティカルフロー推定の原理・特徴・計算式 | アルゴリズム雑記


処理が理解しやすいように、ブロック図にイメージ図も付け足すとこんな感じ:

疎なオプティカルフロー方式のイメージ図
疎なオプティカルフロー方式のイメージ図

オプティカルフロー方式でも誤った点対応(誤トラッキング)が発生する可能性がある。
具体的には「同じような模様の繰り返し」である。