おぺんcv

画像処理エンジニアのブログ

CUDA 8がVisual Studio 2015 Update 3に対応した

雑記

今年の6月頃からCUDA 8.0 RC(Release Candidate)は出てたのですが
このたび正式なCUDA 8.0がリリースされたようです
CUDA Toolkit | NVIDIA Developer

Parallel Forall Blog - Nvidiaでも紹介されていて
個人的に嬉しかったのがVisual Studio 2015 Update 2以降に対応したこと(RCの時は未対応だった)

Expanded developer platform support including Microsoft Visual Studio 2015 (updates 2 and 3) and GCC 5.4 (Ubuntu 16.04).


これで最新のVisual StudioでCUDA8を試すことができますね！
Pascal GPUが欲しくなるなあ

ARM NEONの使い方 組み込み関数一覧編

ARM

組み込み関数一覧

こちらのページにNEONの組み込み関数一覧が載ってます
慣れてきたらここから必要な演算を探すと良いでしょう
infocenter.arm.com

命名規則

NEONの組み込み関数は大抵のものが以下の規則で名前が付けられています

<opname>[q]_<type>

<opname>はvaddやvmul等、演算の内容を表します(必ず頭に'v'が付きます)
[q]はオプションフラグで、qが付く場合は128bit型に対する演算を表します
<type>はオペランドの型を表します

例えば、 vadd_s16 は int16x4_t(64bit型)同士の加算を、
vaddq_s16 は int16x8_t(128bit型)同士の加算を表します

次回

次回からは個別の演算に対して解説していきます
ARM NEONの使い方 ロード・ストア編

ARM NEONの使い方 データ型編

ARM

NEONデータ型の概要

NEONのデータ型は以下の規則で名前が付けられています

<type><size>x<number of lanes>_t

例えば、 int16x4_t は符号付き16bit整数を4個保持する型となります
それぞれの要素を「レーン」と呼びます

NEONのデータ型は64bitのものと128bitのものがあります

NEONデータ型一覧

NEONのデータ型一覧を示します

64bit型 128bit型
int8x8_t int8x16_t
int16x4_t int16x8_t
int32x2_t int32x4_t
int64x1_t int64x2_t
uint8x8_t uint8x16_t
uint16x4_t uint16x8_t
uint32x2_t uint32x4_t
uint64x1_t uint64x2_t
float16x4_t float16x8_t
float32x2_t float32x4_t
poly8x8_t poly8x16_t
poly16x4_t poly16x8_t

以上をさらに配列にした型もあり、以下の規則で名前が付けられています

<type><size>x<number of lanes>x<length of array>_t

例えば、 int16x4x2_t は int16x4_t を2つ要素に持つ型です
各要素はvalというメンバに格納されています

struct int16x4x2_t
{
    int16x4_t val[2];
};

これらの配列型は一部のNEON命令で使用されます
使いどころについては今後の記事で解説したいと思います

polyって何?

データ型一覧の中に、見慣れないpolyという型がありました
気になってググったところ、以下の投稿を発見しました

stackoverflow.com

どうやら多項式(polynomial)算術演算というものをするための型らしく
誤り検出符号や暗号化に利用されているようです

参考

ARM Information Center

次回

個別の組み込み関数について説明する前に、全体像に触れておきたいと思います
ARM NEONの使い方 組み込み関数一覧編

ARM NEONの使い方 予告編

ARM

はじめに

ここ最近、NEON(ARMのSIMD命令)を使ったプログラムを初めて書きました
NEONについて調べていて思ったのが、NEONの使い方を初心者向けに
日本語で説明したページがあまり見当たらないということ
そんなわけで、自分なりにNEONの使い方をまとめておこうと思いました

書く予定のもの

とりあえず自分がある処理をNEON化する上で必要だった知識を取り上げます
すべてのNEON命令をカバーしているわけではありませんので悪しからず

次回

まずはNEONのデータ型について解説します
ARM NEONの使い方 データ型編

今使ってるvimプラグイン

雑記

最近の私はIDE(Visual Studio、Eclipse)を中心に開発することが多かったのですが
休日は久々にvimで遊びました
というわけで(?)、現在私が使っているvimプラグインをご紹介します
私は普段c++を使うことが多いので、c++向けの設定も入ってます

目次

プラグイン管理ツール

まずは管理ツールを入れなきゃねってことで
これはNeoBundleで決まりでしょ!
…って思ってたら気になる記事がいくつか

このvim-plug、導入が楽なので気に入った!
NeoBundleも十分楽だけどこっちはもっと楽でした

それでは本編

scrooloose/nerdtree

GitHub - scrooloose/nerdtree: A tree explorer plugin for vim.

言わずと知れたファイルエクスプローラ
私はF2でエクスプローラを開閉するように設定しています

nnoremap <F2> :NERDTreeToggle<CR>

tomasr/molokai

GitHub - tomasr/molokai: Molokai color scheme for Vim

カラースキーマ
Sublime Textで使われてるやつ
molokai_originalの方が柔らかくて好き

colorscheme molokai
let g:molokai_original = 1

tyru/caw.vim

GitHub - tyru/caw.vim: Vim comment plugin: supported operator/non-operator mappings, repeatable by dot-command, 300+ filetypes

自動コメント/コメント解除
ノーマルモード + <C-k>で現在行をコメント
ビジュアルモード + <C-k>で複数行を一括コメント!

nmap <C-k> <plug>(caw:i:toggle)
vmap <C-k> <plug>(caw:i:toggle)

itchyny/lightline.vim

GitHub - itchyny/lightline.vim: A light and configurable statusline/tabline for Vim

ステータスライン
現在の編集モードなどをカラフルに表示してくれる

octol/vim-cpp-enhanced-highlight

GitHub - octol/vim-cpp-enhanced-highlight: Additional Vim syntax highlighting for C++ (including C++11/14)

c++のシンタックスハイライト強化版
デフォルトのハイライトが寂しいと感じたらおススメ

Townk/Vim-autoclose

GitHub - Townk/vim-autoclose: This plugin for Vim enable an auto-close chars feature for you. For instance if you type an '(', ``autoclose`` will automatically insert a ')' and put the cursor between than.

最近のIDEなら必ずあるであろう
対応する括弧の自動補完

junegunn/vim-easy-align

GitHub - junegunn/vim-easy-align: A Vim alignment plugin

選択範囲を指定した文字(<Space>, =, :, etc.)で揃えてくれる
これが決まると超絶気持ちがいい
vimを使ってない人にドヤりたくなる(やらないけどw)プラグイン

justmao945/vim-clang

GitHub - justmao945/vim-clang: Clang completion plugin for vim

こっからは補完系プラグインの紹介
vim-clangはclangを使ったc++の補完プラグイン
標準ライブラリの補完ができるのは強い

Shougo/neocomplete

GitHub - Shougo/neocomplete.vim: Next generation completion framework after neocomplcache

言わずと知れた補完プラグイン
タイプしたその場で補完候補が出るので一気にIDEっぽくなる
ファイルパスを保管してくれるのも地味に便利

Shougo/neosnippet

GitHub - Shougo/neosnippet.vim: neo-snippet plugin contains neocomplcache snippets source

コードスニペット補完
今日入れましたが、なにこれ超便利
挿入モードでp<C-k>を打つと↓が挿入されて感動!

std::cout <<  << std::endl;

Shougo/neosnippet-snippets

GitHub - Shougo/neosnippet-snippets: The standard snippets repository for neosnippet

スニペットの辞書
neosnippetと一緒にいれておく

Free Space Computation

画像処理 車載

以前の記事で紹介したStixelのアルゴリズムに登場する
Free Space Computationを実装しようと頑張ってます

Free Space Computationとは

視差画像をもとにシーンのFree Space(移動可能な領域、路面)を推定する手法です
推定結果のイメージです(赤く塗った領域がFree Space)

f:id:mizunashi:20160820225810j:plain

より具体的には、Free Space Computationは
画像の列ごとに「路面と物体の下端との境界」となる行を見つける問題となります

参考文献

[1]Efficient Representation of Traffic Scenes by Means of Dynamic Stixels
[2]Free Space Computation Using Stochastic Occupancy Grids and Dynamic Programming

考え方

下図に、画像のある列におけるステレオ観測値(=奥行き)の分布を示します

f:id:mizunashi:20160821230953p:plain

緑線が路面、赤線が物体、青線が背景の領域です
この図から、

  1. 「路面と物体の下端との境界」から手前側(カメラ側)は路面の領域
  2. 物体領域の奥行きはほぼ一定(物体は大抵の場合直立しているため)

ということがわかると思います

これをより数学的に書くと次のようになります

  1. v_bを「路面と物体の下端との境界」となる行、Vを画像下端の行とすると、[v | v_b <= v <= V]を満たすvは路面領域内にあり、その視差はd_R(v)で表される
  2. h_vを物体の高さ(画素単位)とすると、[v | v_b - h_v <= v <= v_b]を満たすvは物体領域内にあり、その視差はd_R(v_b)で一定である

ここでd_R(v)はある行vに対応する路面上の視差を返す関数であり
論文では観測した路面の視差をB-Splineでフィッティングすることでd_R(v)を求めています

この考えをもとに、境界の位置で最小となるようなコスト関数を定義します

コスト関数

コスト関数は以下のように、2つのサブコスト \Gamma_v^R と \Gamma_v^O の重み付き和となります
\Gamma_v = \alpha_1 \Gamma_v^R + \alpha_2 \Gamma_v^O

\Gamma_v^Rは路面(Road)のコストで、以下のように表されます
\Gamma_v^R = \sum_{v=v_b}^V | d_R(v) - d_v |

\Gamma_v^Oは物体(Object)のコストで、以下のように表されます
\Gamma_v^O = \sum_{v={v_b - h_v}}^{v_b} | d_R(v_b) - d_v |

これを各画素について計算します
コードで書くとこんな感じです

for (int u = 0; u < disp.cols; ++u)
{
	for (int vb = 0; vb < disp.rows; ++vb)
	{
		float objectCost = 0, roadCost = 0;

		for (int v = vb; v < disp.rows; ++v)
			roadCost += fabsf(disp(v, u) - roadDisp[v]);

		for (int v = max(vb - hv, 0); v < vb; ++v)
			objectCost += fabsf(disp(v, u) - roadDisp[vb]);

		score(vb, u) = alpha1 * roadCost + alpha2 * objectCost;
	}
}

DPによる境界の計算

最終的な出力である「路面と物体の下端との境界」を計算します
これは先ほど求めたスコアの値だけではなく、境界の空間的な連続性も考慮する必要があります
そこで画像の各画素をノードとして、左端と右端を結ぶ最適な経路を求める問題を考えます

これは動的計画法的に以下のように書けます

for (int v = 0; v < score.rows; v++)
{
	dpcost(v, 0) = score(v, 0)
}

for (int u = 1; u < score.cols; u++)
{
	for (int v = 0; v < score.rows; v++)
	{
		float mincost = FLT_MAX;
		int minpath = 0;

		for (int vv = 0; vv < score.rows; vv++)
		{
			int jump = abs(disp(v, u) - disp(vv, u - 1));
			int penalty = std::min(P1 * jump, P1 * P2);
			float cost = score(v, u) + dpcost(vv, u - 1) + penalty;

			if (cost < mincost)
			{
				mincost = cost;
				minpath = vv;
			}
		}

		dpcost(v, u) = mincost;
		dppath(v, u) = minpath;
	}
}

// back track
float mincost = FLT_MAX;
int minv = 0;
for (int v = 0; v < dpcost.rows; v++)
{
	if (dpcost(v, dpcost.cols - 1) < mincost)
	{
		mincost = dpcost(v, dpcost.cols - 1);
		minv = v;
	}
}

for (int u = dppath.cols - 1; u >= 0; u--)
{
	path[u] = minv;
	minv = dppath(minv, u);
}

上記コードで境界の空間的な連続性を与えるのが、以下の部分です

int jump = abs(disp(v, u) - disp(vv, u - 1));
int penalty = std::min(P1 * jump, P1 * P2);

隣接ノードとの視差の絶対差分に比例したペナルティを与えます
(ただし、ある程度の不連続性を許容するため差分がP2以上であればペナルティは一定とする)

ただこの計算、オーダーがO(w*h*h)となって結構重いです
これであってるんかな…

実装結果

前回紹介したデータセットを使って実験した結果が最初に見せた図です
入力の視差画像はデータセットから与えられます
路面上の視差d_R(v)はB-Splineではなく、直線で近似してます
まあ何となくあってそうに見える…

おわりに

Free Space Computationについて理解したこと、実装したことを紹介しました
今後はB-splineによる路面モデルの推定を実装したいです

EISATS

画像処理 車載

EISATS | Reinhard Klette: CCVというコンピュータビジョン向けのデータセットがありまして

EISATS
The .enpeda.. Image Sequence Analysis Test Site (EISATS) offers sets of image sequences for the purpose of comparative performance evaluation of stereo vision, optic flow, motion analysis, or further techniques in computer vision.

そのうちSET2は車載シーンが合成画像で作成されており
オプティカルフロー、視差、エゴモーション等の正解データが取得できます

ときどきOpenCV上で利用するのですが、毎度データの読み込み方法を忘れるので
自分へのメモとして各データを読み込んで表示するサンプルプログラムを作成しました
github.com

使用する際は前述のSET2から
SEQUENCE2のデータを適当なディレクトリにダウンロード・解凍して

/home/hoge/EISATS/SET2/S2/
├── colour-left-S2
├── disparityGT-S2
├── egoMotion
├── flowX-S2
├── flowY-S2

そのディレクトリをプログラムに渡します

./eisats-synthesized-sequences /home/hoge/EISATS/SET2/S2/

ソースファイルは1個だけなので全部載せておきます

#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <opencv2/opencv.hpp>

namespace {

	std::string numberString(int n)
	{
		std::stringstream ss;
		ss << std::setw(3) << std::setfill('0') << n;
		return ss.str();
	}

	cv::Mat readRaw(const std::string& filename)
	{
		std::ifstream ifs(filename, std::ios::binary);
		if (ifs.fail())
			return cv::Mat();

		std::string line;
		while (getline(ifs, line) && line[0] == '#') {}

		int width, height, depth;
		sscanf(line.c_str(), "%d %d %d", &width, &height, &depth);
		CV_Assert(depth == 32);

		cv::Mat img(height, width, CV_32F);
		ifs.read((char *)img.data, sizeof(float) * width * height);

		return img;
	}

	void readEgoMotion(const std::string& filename, cv::Matx33d& R, cv::Matx31d& t)
	{
		std::ifstream ifs(filename);
		if (ifs.fail())
			return;

		std::string line;
		while (getline(ifs, line) && line[0] == '#') {}

		sscanf(line.c_str(), "%lf %lf %lf %lf", &R(0, 0), &R(0, 1), &R(0, 2), &t(0)); getline(ifs, line);
		sscanf(line.c_str(), "%lf %lf %lf %lf", &R(1, 0), &R(1, 1), &R(1, 2), &t(1)); getline(ifs, line);
		sscanf(line.c_str(), "%lf %lf %lf %lf", &R(2, 0), &R(2, 1), &R(2, 2), &t(2));
	}

	void drawOpticalFlow(cv::Mat& img, const cv::Mat1f& flowx, const cv::Mat1f& flowy)
	{
		float maxnorm = 0;
		for (int i = 0; i < flowx.rows; ++i)
		{
			for (int j = 0; j < flowx.cols; ++j)
			{
				float x = flowx(i, j);
				float y = flowy(i, j);
				maxnorm = std::max(maxnorm, sqrtf(x * x + y * y));
			}
		}

		img.create(flowx.size(), CV_8UC3);
		for (int i = 0; i < flowx.rows; ++i)
		{
			for (int j = 0; j < flowx.cols; ++j)
			{
				float x = flowx(i, j);
				float y = flowy(i, j);

				// Convert flow angle to hue
				float angle = atan2f(y, x);
				if (angle < 0.f) angle += (float)(2 * CV_PI);
				angle /= (float)(2 * CV_PI);
				uchar hue = static_cast<uchar>(180 * angle);

				// Convert flow norm to saturation
				float norm = sqrtf(x * x + y * y) / maxnorm;
				uchar sat = static_cast<uchar>(255 * norm);

				img.at<cv::Vec3b>(i, j) = cv::Vec3b(hue, sat, 255);
			}
		}

		cv::cvtColor(img, img, cv::COLOR_HSV2BGR);
	}

} /* namespace */

int main(int argc, char *argv[])
{
	if (argc < 2)
	{
		std::cerr << "Usage: " << argv[0] << " [data path]" << std::endl;
		return -1;
	}

	std::string dataPath(argv[1]);

	for (int frameNo = 2; ; ++frameNo)
	{
		// Read left image
		cv::Mat leftImg = cv::imread(dataPath + "/colour-left-S2/img_c0_" + numberString(frameNo) + ".ppm");
		if (leftImg.empty())
			break;

		// Read optical flow map
		cv::Mat flowx = readRaw(dataPath + "/flowX-S2/flowU_from_" + numberString(frameNo - 1) + "_to_" + numberString(frameNo) + ".raw");
		cv::Mat flowy = readRaw(dataPath + "/flowY-S2/flowV_from_" + numberString(frameNo - 1) + "_to_" + numberString(frameNo) + ".raw");
		if (flowx.empty() || flowy.empty())
			break;

		// Read disparity map
		cv::Mat disp = readRaw(dataPath + "/disparityGT-S2/stereo_" + numberString(frameNo) + ".raw");
		if (disp.empty())
			break;

		// Read ego-motion
		cv::Matx33d R;
		cv::Matx31d t, r;
		readEgoMotion(dataPath + "/egoMotion/from_" + numberString(frameNo - 1) + "_to_" + numberString(frameNo) + ".txt", R, t);
		cv::Rodrigues(R, r);

		// Draw and display image
		cv::Mat flowImg, dispImg;
		drawOpticalFlow(flowImg, flowx, flowy);
		cv::normalize(disp, dispImg, 0.0, 1.0, cv::NORM_MINMAX);

		cv::putText(leftImg, "Rotation: " + std::to_string(r(0)) + " " + std::to_string(r(1)) + " " + std::to_string(r(2)),
			cv::Point(10, 20), 0, 0.5, cv::Scalar(0, 255, 255));
		cv::putText(leftImg, "Translation: " + std::to_string(t(0)) + " " + std::to_string(t(1)) + " " + std::to_string(t(2)),
			cv::Point(10, 40), 0, 0.5, cv::Scalar(0, 255, 255));

		cv::imshow("Left image", 16 * leftImg);
		cv::imshow("Flow image", flowImg);
		cv::imshow("Disp image", dispImg);

		char c = cv::waitKey(100);
		if (c == 27) // ESC
			break;
	}

	return 0;
}