読者です 読者をやめる 読者になる 読者になる

おぺんcv

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

Visual Studio 2015でHalideを始める

Halideという画像処理向けのプログラミング言語を始めるべく
Visual Studio 2015でHalideの環境構築をしてみました
今回はその手順について紹介します

はじめに

HalideとはMITが開発している画像処理向けの言語です
公式ページやブログなどの紹介によると以下の特徴があるようです

  • C++内で使用するDSL(domain-specific language)である
  • 画像処理を簡単に記述できる
  • 色んなターゲット(x86、ARM、CUDA等)向けに高速化できる

簡単に処理を記述できて、尚且つ一度書いたコードを色んなターゲットで
高速に動かせるとしたら、とても魅力的ですね
リポジトリのREADMEによるとWindows(VS2015)にも対応しているようなので、
試しに導入してみることにしました

参考:

やること

  • Visual Studio 2015でHalideをビルドする
  • Halidを使った簡単なプログラムをビルド・実行する

事前に必要なもの

  • Visual Studio 2015 Update 3 以降
  • CMake
  • libjpegとlibpng (画像ファイルを読み書きするtutorialを動かす場合)

また、コマンドラインでの作業が主体になるので
お好きなコンソールエミュレータを立ち上げておきましょう

ソースの入手

Halid

GitHubからHalidのリポジトリをクローンします
ここではD:\halideというディレクトリにソースを置くことにします

mkdir D:\halide
cd D:\halide
git clone https://github.com/halide/Halide.git
LLVMとClang

HalidのビルドにはLLVMとClangが必要になります
LLVM Download PageからLLVM source codeとClang source codeとダウンロードします
2017年4月22日の時点で最新バージョンは4.0.0でした
解凍してHalidのソースと同じディレクトリに置きます

$pwd
D:\halide

$ls
Halide/  cfe-4.0.0.src/  llvm-4.0.0.src/

ビルド

LLVMのビルド

ソースと同じ場所にllvm-buildというディレクトリを作成しビルドを行います
まずはcmake

$mkdir D:\halide\llvm-build
$cd D:\halide\llvm-build
$cmake -DCMAKE_INSTALL_PREFIX=../llvm-install -DLLVM_ENABLE_TERMINFO=OFF -DLLVM_TARGETS_TO_BUILD=X86;ARM;NVPTX;AArch64;Mips;Hexagon -DLLVM_ENABLE_ASSERTIONS=ON -DLLVM_BUILD_32_BITS=OFF -DCMAKE_BUILD_TYPE=Release ../llvm-4.0.0.src -G "Visual Studio 14 Win64"

次にビルド

MSBuild.exe /m /t:Build /p:Configuration=Release .\INSTALL.vcxproj

私の場合、MSBuild.exeへのパスが通ってなかったので環境変数の設定から
以下のディレクトリをPathに追加しました

C:Program Files (x86)\MSBuild\14.0\Bin

ビルドが終わるまでしばらく待ちます
CPU使用率が100%になってイイ感じ(?)です

Clangのビルド

LLVM同様、cmakeとビルドを行います

$mkdir D:\halide\clang-build
$cd D:\halide\clang-build
$cmake -DCMAKE_INSTALL_PREFIX=../llvm-install -DLLVM_ENABLE_TERMINFO=OFF -DLLVM_TARGETS_TO_BUILD=X86;ARM;NVPTX;AArch64;Mips;Hexagon -DLLVM_ENABLE_ASSERTIONS=ON -DLLVM_BUILD_32_BITS=OFF -DCMAKE_BUILD_TYPE=Release ../cfe-4.0.0.src -G "Visual Studio 14 Win64"
MSBuild.exe /m /t:Build /p:Configuration=Release .\INSTALL.vcxproj
Halideのビルド

最後にHalideをビルドします

$mkdir D:\halide\halide-build
$cd D:\halide\halide-build
$cmake -DLLVM_DIR=../llvm-install/lib/cmake/llvm -DLLVM_VERSION=37 -DCMAKE_BUILD_TYPE=Release -G "Visual Studio 14 Win64" ../halide
MSBuild.exe /m /t:Build /p:Configuration=Release .\ALL_BUILD.vcxproj

ビルドが完了するとhalide-build以下にHalideがインストールされています
私の場合、無事ビルド完了…と思いきやtutorialコードのビルドで失敗してしまいました
それについては余談で触れたいと思います

Halideのプロジェクトを作成してみる

ビルドできたらtutorialを実行するも良し、ですが
ここでは自分でHalideのプロジェクトを作成してみましょう
プロジェクト用ディレクトリを作成し、CMakeLists.txtとソースコードを置きます

halide_sample
   |- CMakeLists.txt
   |- main.cpp
CMake

CMakeLists.txtの中身はこんな感じです
Halide_DIRにHalideのインストール先を設定し
そこからヘッダとライブラリのパスを設定します
Halide.libもリンクしておきます

cmake_minimum_required(VERSION 2.8)

set(Halide_DIR D:/halide/halide-build)

include_directories(${Halide_DIR}/include)
link_directories(${Halide_DIR}/lib/Release)

file(GLOB srcs ./*.cpp ./*.h*)
add_executable(sample ${srcs})

target_link_libraries(sample Halide)
ソース

main.cppの中身はこんな感じです
内容はtutorialのlesson_01_basics.cppから引っ張ってきました
Halide::FuncとHalide::Varで関数を記述し、realizeで呼び出すって感じでしょうか

#include <Halide.h>
#include <iostream>

int main()
{
	Halide::Func gradient;
	Halide::Var x, y;
	gradient(x, y) = x + y;
	Halide::Buffer<int32_t> output = gradient.realize(3, 3);
	
	for (int i = 0; i < output.height(); i++)
	{
		for (int j = 0; j < output.width(); j++)
		{
			std::cout << output(i, j) << " ";
		}
		std::cout << std::endl;
	}

	return 0;
}
プロジェクトの作成とビルド

cmakeを使ってプロジェクトの作成とビルドを行います

$mkdir build
$cd build
$cmake .. -G "Visual Studio 14 Win64"
cmake --build . --config Release
パスの設定

実行にはHalide.dllが必要になりますので
環境変数のPathにdllの場所を指定しておきましょう

D:\halide\halide-build\bin\Release
実行結果

ちゃんと (x, y) = x + y という結果になってました

0 1 2
1 2 3
2 3 4

おわりに

今回はVisual Studio 2015でHalideの環境構築するところまでを紹介しました
今後はより複雑なアルゴリズムの記述や、高速化の方法について調べてみようと思います

【余談】tutorialのビルドにハマった話

私の場合、tutorialのビルド時に以下のようなビルドエラーに遭遇してしまいました

  • jpeg.libで未定義のシンボル__iob_funcが参照されている
  • jpeg.libで未定義のシンボルsprintfが参照されている
  • INT32が再定義されている

ビルドエラーは画像ファイルを読み書きするtutorial(lesson_02,07,09,12)で起きてました
Halide/tutorialのCMakeLists.txtを見ると、これらのサンプルはlibjpegとlibpngが入ってないと
ビルドされないはずなのですが、libjpegとlibpngなんて入れたっけなあ…
と思ってたら、どうやらAnaconda2(python環境)を入れた際に一緒に入ってたみたいです

__iob_funcの未定義

これについては、ググると回避方法が見つかりました
Halide/tools/halide_image_io.hで__iob_funcを定義してやります

FILE _iob[] = { *stdin, *stdout, *stderr };
extern "C" FILE * __cdecl __iob_func(void)
{
    return _iob;
}

参考:

sprintfの未定義

なんでsprintfが未定義なの?ってのが謎ですが
どっかに定義されてるだろってことで、extern宣言したら回避できました(適当)
これもhalide_image_io.hに書いてます

extern "C" int sprintf ( char * str, const char * format, ... );
NT32の再定義

どうやらINT32がlibjpgのjmorecfg.hと、Windows Kits\8.1\Include\shared\basetsd.hで
両方定義されちゃってるようです
jmorecfg.hをみると、INT32はQGLOBAL_Hが定義されていない場合に定義されるようです
そこで、halide_image_io.hでjpeglib.hをインクルードする前にQGLOBAL_Hを定義してやります

#ifndef HALIDE_NO_JPEG
#define QGLOBAL_H
#include "jpeglib.h"
#undef QGLOBAL_H
#endif

以上でビルドエラーを回避することができましたが、かなり適当なのでご参考までに

Stixel Computationを実装してみた multi-layer編

はじめに

以前Stixel Computationを実装してみたで紹介したStixelですが
今回はその発展形である、Multi-Layered Stixelを実装してみました

概要

Stixelとは、Daimlerが研究しているステレオカメラによる障害物検出の手法です
画像を幅5~10ピクセル程度の細長い棒に区切って、そこからさらに障害物の領域を検出します
以前紹介したStixel計算手法では、区切られた1つの列からは1つのStixel(最も手前側にあるもの)
しか検出できませんでしたが、Multi-Layered Stixelでは複数のStixelを検出可能になりました

デモ

物体上に出てる棒状のものがStixelです
以前のStixelとは違い、ガードレールの奥にいる車両のStixelも検出できています

www.youtube.com

アルゴリズム解説

ざっくり説明するとMulti-Layered Stixelの計算は
画像の各列を路面(ground)、物体(object)、空(sky)のいずれかのセグメントに分類する問題です
列を3つのセグメントに分割する組み合わせは色々と考えられますが、各セグメントに対し

  • Data項:セグメントがあるクラスに属しているとき期待される視差値と実際の視差値の誤差
  • Prior項:セグメントの割り当て方に関する先験情報

からなるコスト関数を定義し、これらの和が最小となる分割を求めます
最小化問題を効率的に解くために動的計画法が用いられます

さらなる詳細は[1]の3 The Multi-Layered Stixel Worldを参照下さい
数式の展開や実装方法など、かなり細かく解説してあります

おわりに

参考文献は著者の博士論文のようで、非常に読み応えがありました
なかなか内容が理解できず、3回くらいの挫折を経て何とか実装に至りました
また実装ができても当初は処理がめちゃめちゃ重く
inline化やテーブル化などの高速化を頑張りました

Stixel検出後の処理として、Stixelの追跡やグループ化などの処理があるようです
これらもいつか実装できればと思っています

3次元のcv::Matを作る

はじめに

小ネタです
cv::Matが多次元に対応しているのを知ったのでメモ

作り方

多次元Matのコンストラクタがこちら(mat.hppより抜粋)
ndimsで次元数を、sizesで各次元の長さを、typeで要素の型を指定します

    /** @overload
    @param ndims Array dimensionality.
    @param sizes Array of integers specifying an n-dimensional array shape.
    @param type Array type. Use CV_8UC1, ..., CV_64FC4 to create 1-4 channel matrices, or
    CV_8UC(n), ..., CV_64FC(n) to create multi-channel (up to CV_CN_MAX channels) matrices.
    */
    Mat(int ndims, const int* sizes, int type);

サンプル

2×3×4のint行列を作ってみました

#include <opencv2/opencv.hpp>
#include <iostream>

int main()
{
	const int ni = 2, nj = 3, nk = 4;
	const int ndims = 3;
	const int sizes[ndims] = { ni, nj, nk };

	cv::Mat m(ndims, sizes, CV_32S);
	
	for (int i = 0; i < ni; i++)
		for (int j = 0; j < nj; j++)
			for (int k = 0; k < nk; k++)
				m.at<int>(i, j, k) = i + j + k;

	for (int i = 0; i < ni; i++)
		for (int j = 0; j < nj; j++)
			for (int k = 0; k < nk; k++)
				std::cout << i << "+" << j << "+" << k << " = " << m.at<int>(i, j, k) << std::endl;

	// std::cout << m << std::endl; // OpenCV Error: Assertion failed (m.dims <= 2)

	return 0;
}

(i, j, k)で要素にアクセスできるのが素晴らしいです
ちなみに、coutで表示しようとしたらassertに引っかかりました
dims <= 2までしか表示できないようです

おわりに

今までcv::Matは2次元しか対応してないと勝手に思っていて
3次元のデータを扱う際はcv::Matではなく生ポインタつかってました
早く気付いていれば良かったと思う今日この頃…

kd-treeを実装してみた

はじめに

kd-treeを実装してみました
最近仕事でよく使うので勉強がてら

kd-treeとは

最近傍探索を効率的にできるデータ構造です
kd木 - Wikipedia

kd-treeが使えるライブラリとしてはFLANNPCLが有名どころでしょうか

ソースコード

以下に公開してあります
github.com

kd-treeの構築と、以下の探索機能を実装してます

  • 最近傍探索 (Nearest neighbor search)
  • k-最近傍探索 (K-nearest neighbor search)
  • 半径内に含まれる近傍の探索 (Radius search)

あと、ヘッダ1個includeするだけで使えるのでお手軽です

アルゴリズム

kd-treeの構築

wikipediaに載ってる疑似コードが分かりやすかったので引用させていただきます

function kdtree (list of points pointList, int depth)
{
    if pointList is empty
        return nil;
    else
    {
        // 深さに応じて軸を選択し、軸が順次選択されるようにする
        var int axis := depth mod k;

        // 点のリストをソートし、中央値の点を選択する
        select median from pointList;

        // ノードを作成し、部分木を構築する
        var tree_node node;
        node.location := median;
        node.leftChild := kdtree(points in pointList before median, depth+1);
        node.rightChild := kdtree(points in pointList after median, depth+1);
        return node;
    }
}

一方私の書いたコードはこんな↓感じです

Node* buildRecursive(int* indices, int npoints, int depth)
{
	if (npoints <= 0)
		return nullptr;

	const int axis = depth % PointT::DIM;
	const int mid = (npoints - 1) / 2;

	std::nth_element(indices, indices + mid, indices + npoints, [&](int lhs, int rhs)
	{
		return points_[lhs][axis] < points_[rhs][axis];
	});

	Node* node = new Node();
	node->idx = indices[mid];
	node->axis = axis;

	node->next[0] = buildRecursive(indices, mid, depth + 1);
	node->next[1] = buildRecursive(indices + mid + 1, npoints - mid - 1, depth + 1);

	return node;
}

疑似コードとほとんど同じです
私の場合は点そのものではなく、点へのインデックスをノードに保持しています
(中央値を取得するためにnth_elementを初めて使った…)

最近傍探索

k-d treeが構築できたらlet's最近傍探索
ググってヒットしたこちらの資料をもとに実装しました
(URLを見るとスタンフォード大の宿題らしい)
https://web.stanford.edu/class/cs106l/handouts/assignment-3-kdtree.pdf

原理は2分探索と同じで、「探している値は木の半分のこっち側にあるよね」
ってのを繰り返しながら左か右の子ノードに降りていきます
その過程で現在ノードとクエリとの距離を計算し最近傍の値を更新します
葉までたどり着いたら、自分の兄弟ノード(親から見て自分が左だったら、右のノード)の方に
最近傍の可能性がないかチェックして、ある場合は兄弟ノードの方も探索します

void nnSearchRecursive(const PointT& query, const Node* node, int *guess, double *minDist) const
{
	if (node == nullptr)
		return;

	const PointT& train = points_[node->idx];

	const double dist = distance(query, train);
	if (dist < *minDist)
	{
		*minDist = dist;
		*guess = node->idx;
	}

	const int axis = node->axis;
	const int dir = query[axis] < train[axis] ? 0 : 1;
	nnSearchRecursive(query, node->next[dir], guess, minDist);

	const double diff = fabs(query[axis] - train[axis]);
	if (diff < *minDist)
		nnSearchRecursive(query, node->next[!dir], guess, minDist);
}

K-近傍探索や半径探索も探索の流れは同じです

デモ

最近傍探索/k-最近傍探索/半径探索を実行した結果です
赤い点がクエリ(画像の中心座標)、青い点が近傍点になります

おわりに

kd-treeを実装してみました
実装に当たり、再帰関数をいっぱい書きました
こういう場合は再帰を使うとシンプルに書けるんだなーって思いました

ARM NEONの使い方 除算編

四則演算ラスト!
今回は除算編です

アーキテクチャによる違い

32bit ARM

32bit ARMには除算を直接行うSIMD命令がないため
a/bという除算を行う場合、まず逆数を計算するSIMD命令によってbの逆数を計算し
これをaと掛け算することによって除算を表現します

64bit ARM

64bit ARMには除算を行うSIMD命令が追加されています

逆数の計算

vrecpe[q]_f32(va)

32bit浮動小数ベクタvaの各要素の逆数の近似値を計算します
近似値なのでCPU(IEEE754)よりも精度が低いです

サンプル

32bit浮動小数ベクタva, vbの除算をva*(vbの逆数)という形でやってみます
vaとvbの各レーンを同じ値に設定しているので、除算後の値が1になるのが望ましい結果です

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

static void print_vector(float32x4_t value)
{
	float lane[4];
	vst1q_f32(lane, value);
	for (int i = 0; i < 4; i++)
		printf("lane[%d]: %f\n", i, lane[i]);
}

int main()
{
	float a[4] = { 1, 2, 3, 4 };
	float b[4] = { 1, 2, 3, 4 };

	float32x4_t va = vld1q_f32(a);
	float32x4_t vb = vld1q_f32(b);

	float32x4_t vrecipb = vrecpeq_f32(vb); // approximate reciprocal of vb
	float32x4_t vc = vmulq_f32(va, vrecipb);

	printf("va:\n");
	print_vector(va);

	printf("\nvb:\n");
	print_vector(vb);

	printf("\n1/vb:\n");
	print_vector(vrecipb);

	printf("\nva * (1/vb):\n");
	print_vector(vc);

	return 0;
}
実行結果
va:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

vb:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

1/vb:
lane[0]: 0.998047
lane[1]: 0.499023
lane[2]: 0.333008
lane[3]: 0.249512

va * (1/vb):
lane[0]: 0.998047
lane[1]: 0.998047
lane[2]: 0.999023
lane[3]: 0.998047

近似値なので結構精度低いなーという印象
まあ用途次第でしょうけど

ベクタ逆数のステップ

よくわからないタイトルですが
先ほどの逆数の精度を良くする方法です

vrecps[q]_f32(va)

ARM Information Centerの説明が以下です

VRECPS (ベクタ逆数のステップ)は、あるベクタの要素を対応する別のベクタの要素で乗算し、その結果を 2 から減算して、最終的な結果をデスティネーションベクタの要素に返します。

最初この説明を見たときはなんのこっちゃと思いましたが
ニュートン法で逆数の近似を求める際に使うようです

ある数 d とその逆数の現在の近似値 x_nがあるとき、ニュートン法によって
 x_{n+1} = x_n(2 - dx_n)
という反復を繰り返すことで近似値が 1/dに収束します
この (2 - dx_n)を計算してくれるのがvrecpsというわけ

サンプル

先ほどの逆数のサンプルにニュートン法の反復計算を追加してみました

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

static void print_vector(float32x4_t value)
{
	float lane[4];
	vst1q_f32(lane, value);
	for (int i = 0; i < 4; i++)
		printf("lane[%d]: %f\n", i, lane[i]);
}

int main()
{
	float a[4] = { 1, 2, 3, 4 };
	float b[4] = { 1, 2, 3, 4 };

	float32x4_t va = vld1q_f32(a);
	float32x4_t vb = vld1q_f32(b);

	float32x4_t vrecipb = vrecpeq_f32(vb); // approximate reciprocal of vb
	float32x4_t vc = vmulq_f32(va, vrecipb);

	printf("va:\n");
	print_vector(va);

	printf("\nvb:\n");
	print_vector(vb);

	printf("\n1/vb:\n");
	print_vector(vrecipb);

	printf("\nva * (1/vb):\n");
	print_vector(vc);

        //  Newton-Raphson iteration
	int steps = 3;
	for (int i = 0; i < steps; i++)
	{
		printf("\nNewton step: %d\n", i + 1);

		float32x4_t vtmp = vrecpsq_f32(vb, vrecipb);
		vrecipb = vmulq_f32(vrecipb, vtmp);
		vc = vmulq_f32(va, vrecipb);

		printf("\n1/vb:\n");
		print_vector(vrecipb);

		printf("\nva * (1/vb):\n");
		print_vector(vc);
	}

	return 0;
}
実行結果
va:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

vb:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

1/vb:
lane[0]: 0.998047
lane[1]: 0.499023
lane[2]: 0.333008
lane[3]: 0.249512

va * (1/vb):
lane[0]: 0.998047
lane[1]: 0.998047
lane[2]: 0.999023
lane[3]: 0.998047

Newton step: 1

1/vb:
lane[0]: 0.999996
lane[1]: 0.499998
lane[2]: 0.333333
lane[3]: 0.249999

va * (1/vb):
lane[0]: 0.999996
lane[1]: 0.999996
lane[2]: 0.999999
lane[3]: 0.999996

Newton step: 2

1/vb:
lane[0]: 1.000000
lane[1]: 0.500000
lane[2]: 0.333333
lane[3]: 0.250000

va * (1/vb):
lane[0]: 1.000000
lane[1]: 1.000000
lane[2]: 1.000000
lane[3]: 1.000000

Newton step: 3

1/vb:
lane[0]: 1.000000
lane[1]: 0.500000
lane[2]: 0.333333
lane[3]: 0.250000

va * (1/vb):
lane[0]: 1.000000
lane[1]: 1.000000
lane[2]: 1.000000
lane[3]: 1.000000

(定量的に評価したわけではありませんが)
ステップを進めると精度が良くなっていくのがわかります

除算 (64bit ARMのみ)

vdiv[q]_f32(va, vb)

32bit浮動小数ベクタvaとvbの除算を行います

サンプル

vdivq_f32を使用して32bit浮動小数ベクタva, vbの除算をやってみます
コンパイルする際は64bit用のコンパイラ(aarch64-linux-gnu-g++など)を使ってね

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

static void print_vector(float32x4_t value)
{
	float lane[4];
	vst1q_f32(lane, value);
	for (int i = 0; i < 4; i++)
		printf("lane[%d]: %f\n", i, lane[i]);
}

int main()
{
	float a[4] = { 1, 2, 3, 4 };
	float b[4] = { 1, 2, 3, 4 };

	float32x4_t va = vld1q_f32(a);
	float32x4_t vb = vld1q_f32(b);

	float32x4_t vc = vdivq_f32(va, vb);

	printf("va:\n");
	print_vector(va);

	printf("\nvb:\n");
	print_vector(vb);

	printf("\nva / vb:\n");
	print_vector(vc);

	return 0;
}
実行結果
va:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

vb:
lane[0]: 1.000000
lane[1]: 2.000000
lane[2]: 3.000000
lane[3]: 4.000000

va / vb:
lane[0]: 1.000000
lane[1]: 1.000000
lane[2]: 1.000000
lane[3]: 1.000000

そのまま除算ができるとありがたいですね
また64bit ARM(ARMv8)では浮動小数SIMD演算がIEEE754準拠になったため
結果もCPUと変わらないはずです

参考

NEONの32bitと64bitの違いが解説されてます
今回の記事を書くにあたりとても参考になりました
community.arm.com

次回

次回は初期化編です
ARM NEONの使い方 初期化編

Stixel Computationを実装してみた

6D-Visionについて調べてみたで紹介したStixelですが
ちまちまと実装を重ね、ようやくひと段落したので公開することにしました

デモ

物体上に出てる棒のようなものがStixelです
入力の視差はOpenCVのStereo SGBMで計算しています

youtu.be


youtu.be

使い方

README.mdに記載しましたので
興味があれば使ってやって下さい

アルゴリズム解説

執筆中…!

参考文献

  • [1] D. Pfeiffer, U. Franke: “Efficient Representation of Traffic Scenes by means of Dynamic Stixels”, IEEE Intelligent Vehicles Symposium IV 2010, San Diego, CA, Juni 2010
  • [2] H. Badino, U. Franke, and D. Pfeiffer, “The stixel world - a compact medium level representation of the 3d-world,” in DAGM Symposium, (Jena, Germany), September 2009.

ARM NEONの使い方 乗算編

今回は乗算編です
特に変わった内容はありません

乗算

vmul[q]_<type>(va, vb)

64bit(qが付く場合は128bit)のベクタvaとvbの乗算

サンプル

符号付き16bit整数ベクタの乗算をしてみます

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

int main()
{
	int16_t a[4] = { 1, 2, 3, 4 };
	int16_t b[4] = { 5, 6, 7, 8 };

	int16x4_t va = vld1_s16(a);
	int16x4_t vb = vld1_s16(b);

	int16x4_t vc = vmul_s16(va, vb);

	int16_t c[4];
	vst1_s16(c, vc);
	for (int i = 0; i < 4; i++)
		printf("c[%d]: %d\n", i, c[i]);

	return 0;
}
実行結果
c[0]: 5
c[1]: 12
c[2]: 21
c[3]: 32

整数を使っていて桁あふれを考慮する場合はvmull(Vector long multiply)を使いましょう

積和演算

みんな大好き(?)積和演算
内積とか計算するときに重宝しますね

vmla[q]_<type>(va, vb, vc)

64bit(qが付く場合は128bit)のベクタva、vb、vcについて

va + vb * vc

を計算します

サンプル

符号付き16bit整数のベクタの積和演算です

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

int main()
{
	int16_t a[4] = { 1, 1, 1, 1 };
	int16_t b[4] = { 2, 2, 2, 2 };
	int16_t c[4] = { 3, 3, 3, 3 };

	int16x4_t va = vld1_s16(a);
	int16x4_t vb = vld1_s16(b);
	int16x4_t vc = vld1_s16(c);

	int16x4_t vd = vmla_s16(va, vb, vc);

	int16_t d[4];
	vst1_s16(d, vd);
	for (int i = 0; i < 4; i++)
		printf("d[%d]: %d\n", i, d[i]);

	return 0;
}
実行結果
d[0]: 7
d[1]: 7
d[2]: 7
d[3]: 7

doubling multiplyってなんぞや

乗算に関しては先ほどのmulとmlaを使うことがほとんどだと思いますが
命令一覧を見ると doubling multiply というのがあったので使ってみました

サンプル

vqdmull(Vector saturating doubling long multiply)を使って乗算をしてみると…

#include <stdio.h>
#include <stdint.h>
#include <arm_neon.h>

int main()
{
	int16_t a[4] = { 1, 2, 3, 4 };
	int16_t b[4] = { 5, 6, 7, 8 };

	int16x4_t va = vld1_s16(a);
	int16x4_t vb = vld1_s16(b);

	int32x4_t vc = vqdmull_s16(va, vb);

	int32_t c[4];
	vst1q_s32(c, vc);
	for (int i = 0; i < 4; i++)
		printf("c[%d]: %d\n", i, c[i]);

	return 0;
}
実行結果
c[0]: 10
c[1]: 24
c[2]: 42
c[3]: 64

結果が2倍になった
どうやら乗算してさらに2倍する演算のようです
使い道はあるのかな…

次回

次回は除算編です
ARM NEONの使い方 除算編