はじめに
Halideに関するリンクは以下の記事にまとめてあります．
Halideによる画像処理まとめ
Halideとは，画像処理専用のプログラミング言語，つまりドメイン固有言語です．
ドメイン固有言語は，なんでもかけるC，C++，Java，Pythonといったプログラミング言語とは違って，特定の作業に特化したプログラミング言語です．
例えば，SQLなどもドメイン固有言語と言っても良いでしょう．
特定の作業に特化すると，簡単な記述で複雑な処理を実現できたり通常のコンパイラでは出来ない最適化が可能になったります．
一方で，チューリング完全である必要も無いため，その言語では書けないプログラムが出てくることすらありえます．
Halideは関数型言語ですが，チューリング完全ではありません．つまり，プログラムのすべては記述することが出来ません．
そのため，C++にラップする形で使われています．
プログラミング言語Halideは，まるでライブラリをリンクして関数やクラスを使うかのようにして動かします．
つまり，ヘッダをincludeしてとライブラリをリンクして使います．
コンパイラを使ってgcc,iccなどといったコンパイラのようなhalideccのようなものはありません．（※ 似たようなことは出来ますが）
詳細は以下のリンク．静的コンパイルvsJITコンパイルにあります．
Halideの何がすごいのかというと，簡単な記述で，画像処理の高速化に卓越した人よりも高速なコードが生成可能になる点です．
画像処理の高速化のためには，並列化，ベクトル化以外にも様々なめんどくさくて泥臭い作業があります．
それを，Halideの画像処理コンパイラは，自動的にいろいろやってくれます．
もっとも典型的な画像処理の高速化方法は，「並列化」してマルチコアで動かし，「ベクトル化」してSIMD命令を発行し，そして，「ブロック化」して処理を局所化することでキャッシュ効率を最大化するといった，３つの処理が上げられます．
このような最適化を試みるコードを書くと，画像をX×Yのタイルに区切って，それをSIMDイントリンシックでベクトル化して，タイルを並列化して実行して，順序が最適になるようにスケジューリングを書いて．．．といったコードを書く必要があり，単純な３ｘ３の平均値フィルタですら膨大な量の長さのコードになってしまいます．
例えば，上記のようにコードを限界までハンドチューンして展開した３ｘ３のボックスフィルタは以下のようになります．
void blur_3x3(const Image& in, Image& blury)
{
__m128i one_third = _mm_set1_epi16(21846);
#pragma omp parallel for
for (int yTile = 0; yTile < in.height(); yTile += 32)
{
__m128i a, b, c, sum, avg;
__m128i blurx[(256 / 8)*(32 + 2)]; // allocate tile blurx array
for (int xTile = 0; xTile < in.width(); xTile += 256)
{
__m128i *blurxPtr = blurx;
for (int y = -1; y < 32 + 1; y++)
{
const uint16_t *inPtr = &(in[yTile + y][xTile]);
for (int x = 0; x < 256; x += 8)
{
a = _mm_loadu_si128((__m128i*)(inPtr - 1));
b = _mm_loadu_si128((__m128i*)(inPtr + 1));
c = _mm_load_si128((__m128i*)(inPtr));
sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
avg = _mm_mulhi_epi16(sum, one_third);
_mm_store_si128(blurxPtr++, avg);
inPtr += 8;
}
}
blurxPtr = blurx;
for (int y = 0; y < 32; y++)
{
__m128i *outPtr = (__m128i *)(&(blury[yTile + y][xTile]));
for (int x = 0; x < 256; x += 8)
{
a = _mm_load_si128(blurxPtr + (2 * 256) / 8);
b = _mm_load_si128(blurxPtr + 256 / 8);
c = _mm_load_si128(blurxPtr++);
sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
avg = _mm_mulhi_epi16(sum, one_third);
_mm_store_si128(outPtr++, avg);
}
}
}
}
}
この最適化するときは，
- どのループを並列化するか？
- どの単位をベクトル化するか？
- タイリングする単位を何×何にするのか？
- タイリングしたブロックをどの順番で実行するのか？
- あまった画素はどのようにするのか？
- あまらないようにどのように画素をパディングするのか？
- AVXを使うのかSSEまでにするのかNeonがいるのか？
など考慮するべき条件がたくさんあります．
しかも，この条件は，プログラム中のifやswitchで切り替えることでは実現できず，無数の条件を持つコードとして書き起こす必要があります．
どのパターンが最も速いかを試すために，すべてのコードを書き起こしていては時間が全く足りません．
それがHalideを使うと以下の通りになります．
まず，３×３の平均値フィルタを，縦・横の変数分離型の形のフィルタに分解した定義をします．
次に，そのフィルタをタイル分割し，ベクトル化し，並列に実行するというスケジュールを定義して実行します．
Func blur_3x3(Func input)
{
Func blur_x, blur_y;
Var x, y, xi, yi;
// The algorithm - no storage or order
blur_x(x, y) = (input(x-1, y) + input(x, y) + input(x+1, y))/3;
blur_y(x, y) = (blur_x(x, y-1) + blur_x(x, y) + blur_x(x, y+1))/3;
// The schedule - defines order, locality; implies storage
blur_y.tile(x, y, xi, yi, 256, 32)
.vectorize(xi, 8).parallel(y);
blur_x.compute_at(blur_y, x).vectorize(x, 8);
return blur_y;
}
まず，アルゴリズムを定義して何が計算されるのか定義して，そしてスケジュールで何がいつ計算されるのか設計するこのようにアルゴリズムと計算法ができるだけ分離するように分離した言語がHalideです．
上記プログラムの例だと，
Func blur_x, blur_y;
がアルゴリズムの定義
.tale, .vectorize, parallel
などのメソッドがスケジュールの定義です．
このメソッドに任意に数字を入れれば様々な条件のコードが生成できます．
このように，アルゴリズム設計と計算順序の最適化が分離されることで，最適なコードを様々なパターンで生成可能になり，最適なコード生成が総当りで試すことが可能になります．
一方で，最適化をがんばるとコンパイルに丸一日かかることすらありえます．
出た当時の論文は，それくらいかかると書いてあります．
ただし現在はだいぶ改善しているようです．
ライブラリとの違い
ユーザとしては，画像処理ライブラリとして提供されているライブラリがあれば，こんなコンパイラはあまり出番がないのじゃないか？という考え方もあるでしょう．
しかし，プログラミング言語として実装されているという点は，処理速度に強力なアドバンテージがあります．
例えば，画像全体を二乗したかったとしましょう．
OpenCVであれば以下のようにかけば良いでしょう．
Mat a(1024,1024,CV_32F);
Mat dest;
multiply(a, a, dest);
このコードは，最適に書いた関数さえ使えば画像処理コンパイラをあえて使う必要はどこにもありません．
ここで，ほんの少しだけ処理を複雑にして見ましょう．
例えば，画像全体を二乗して何かの値を足したかったかったとしましょう．
OpenCVであれば以下のようにかけば良いでしょう．
Mat a(1024,1024,CV_32F);
Mat b(1024,1024,CV_32F);
Mat dest;
multiply(a, a, a);
add(a, b, dest);
これは，たとえmultiplyやaddのそれぞれの関数が限界まで最適化されていたとしても，最適な演算にはなりえません．
なぜなら，multiplyとaddの関数それぞれで，画像全体を2度も走査しているからです．
本来なら一回の走査で二乗して足せばいいためです．
更にはsimdのFMA命令を使ったらもっと早くなります．
ライブラリは，単純な加算，乗算などの基本演算は用意してあるものの，その組み合わせは，数は爆発に増加するため用意されていることが少なく，基本演算を並べて発行するだけにとどまります．
このあたりは，MatlabやPythonがどうにも速く計算が出来ないという理由とも共通しています．
しかし，Halideを使えばこのような状況でも，数式を定義するだけで簡単に，そして高速な演算を書くことが出来ます．
例えば入力画像を二乗して20を足すコードを書くとします．
計算の結果，
- OpenCVは73 ms
- Halideは9.8ms
と圧倒的な違いがありました．
計算機環境は，Xeon X5690 3.47GHz Dual CPU (12core 24thread)です．
使ったコードは以下のようになります．
#include "Halide.h"
#include <stdio.h>
using namespace Halide;
using namespace cv;
// Support code for loading pngs.
#include "halide_image_io.h"
using namespace Halide::Tools;
int main(int argc, char **argv)
{
//Halide part//////////////////////////////////////////////////////
Image<uint8_t> input = load_image("lenna-gray-8kx8k.png");
Var x("x"), y("y"), xi("xi"), yi("yi");
//Algorithm
Func out("out");
out(x, y) = input(x, y)*input(x, y) + 20;
Image<uint8_t> result(input.width() - 2, input.height() - 2);
//Schedule
result.set_min(1, 1);
out.compute_root();
out.tile(x, y, xi, yi, 256, 32).vectorize(xi, 4).parallel(y);
//main loop
for (int i = 0; i < 10; i++)
{
std::cout << i << ": ";
int64 start = cv::getTickCount();
out.realize(result);
int64 end = cv::getTickCount();
std::cout << (end - start) * 1000 / cv::getTickFrequency() << " [ms]" << std::endl;
}
//OpenCV part //////////////////////////////////////////////////////
Mat a = imread("lenna-gray-8kx8k.png", 0);
for (int i = 0; i < 10; i++)
{
std::cout << i << ":";
Mat b = Mat::zeros(a.size(), CV_8UC3);
int64 start = cv::getTickCount();
multiply(a, a, b);
add(b, 20, b);
int64 end = cv::getTickCount();
std::cout << (end - start) * 1000 / cv::getTickFrequency() << " [ms]" << std::endl;
}
return 0;
}
vs OpenCV
Halideでコンパイルした３ｘ３の平均値フィルタの性能を比べてみましょう．
入力画像は，8192×8192の巨大なグレイスケール画像です．
比較対象となるOpenCVの関数はblur関数を使いました．これは，インテグラルイメージやSIMD演算の最適化など十分なチューンがされています．
結果は，
- Halide: 72ms
- OpenCV: 123ms
で，Halideのほうが高速でした．
計算機環境は，Xeon X5690 3.47GHz Dual CPU (12core 24thread)です．
テストしたコードは以下のようになっています．
#include "Halide.h"
#include <stdio.h>
using namespace Halide;
using namespace cv;
// Support code for loading pngs.
#include "halide_image_io.h"
using namespace Halide::Tools;
int main(int argc, char **argv)
{
//Halide part//////////////////////////////////////////////////////
Image<uint8_t> input = load_image("lenna-gray-8kx8k.png");
Var x("x"), y("y"), xi("xi"), yi("yi");
//Algorithm
Func input_16("input_16");
input_16(x, y) = cast<uint16_t>(input(x, y));
Func blur_x("blur_x");
blur_x(x, y) = (input_16(x - 1, y) + input_16(x, y) + input_16(x + 1, y));
Func blur_y("blur_y");
blur_y(x, y) = (blur_x(x, y - 1) + blur_x(x, y) + blur_x(x, y + 1)) / 9;
Func output("output");
output(x, y) = cast<uint8_t>(blur_y(x, y));
Image<uint8_t> result(input.width() - 2, input.height() - 2);
//Schedule
result.set_min(1, 1);
blur_y.compute_root();
blur_y.tile(x, y, xi, yi, 256, 32).vectorize(xi, 4).parallel(y);
blur_x.compute_at(blur_y, x).vectorize(x, 4);
//main loop
for (int i = 0; i < 10; i++)
{
std::cout << i << ": ";
int64 start = cv::getTickCount();
output.realize(result);
int64 end = cv::getTickCount();
std::cout << (end - start) * 1000 / cv::getTickFrequency() << "[ms]" << std::endl;
}
//OpenCV part //////////////////////////////////////////////////////
Mat a = imread("lenna-gray-8kx8k.png", 0);
for (int i = 0; i < 10; i++)
{
std::cout << i << ":";
Mat b = Mat::zeros(a.size(), CV_8UC3);
int64 start = cv::getTickCount();
blur(a, b, Size(3, 3));
int64 end = cv::getTickCount();
std::cout << (end - start) * 1000 / cv::getTickFrequency() << "[ms]" << std::endl;
}
printf("Success!\n");
return 0;
}
Halideチュートリアルへのショートコメント
Halideを始めるに当たり，あまり多くの情報はありませんが，公式サイトにチュートリアルがあります．
このチュートリアルの一部を翻訳要約しました．
Halideチュートリアルのリンクはこちら．
00 イントロダクション
プログラミングを覚えるためには，自分の手を汚して，動かして覚えましょう！
01 Func, Vars, Exprsからはじめよう
1回目の最初のチュートリアルだけ全訳＋コメントです．
このチュートリアルでは，まずは，コードをコピペしてコンパイルできることを確認しましょう．
windowsユーザ用へのコメント：Windows用のバイナリを落としたら，Debug, Releaseにパスを通して，include,libもVisual studio等のパスを通してください．Halide.lib,Halide.dllとincludeのファイルの中身があればこれで動きます．
使用するコードは以下のコードです．
以下では，それぞれについて解説していきます．
#include "Halide.h"
#include <stdio.h>
int main(int argc, char **argv)
{
Halide::Func gradient;
Halide::Var x, y;
Halide::Expr e = x + y;
gradient(x, y) = e;
Halide::Image<int32_t> output = gradient.realize(800, 600);
for (int j = 0; j < output.height(); j++) {
for (int i = 0; i < output.width(); i++) {
// We can access a pixel of an Image object using similar
// syntax to defining and using functions.
if (output(i, j) != i + j) {
printf("Something went wrong!\n"
"Pixel %d, %d was supposed to be %d, but instead it's %d\n",
i, j, i + j, output(i, j));
return -1;
}
}
}
printf("Success!\n");
}
Halideでは，Halide.hをインクルードすることで，基本的には関連する全てのincludeファイルを読み込むことが出来ます．
#include "Halide.h"
Funcオブジェクトは，パイプラインステージというものを表しています．
これは純粋関数で，各画素がどの値を持つべきかを定めてます．
このオブジェクトを，計算結果の画像と思って差し支えないでしょう．
このサンプルプログラムは，グレイスケール画像に対して，左上から順番にの対角方向にグラディエーションとなる値を書き込む，シングルステージの画像処理パイプラインです．
Halide::Func gradient;
VarオブジェクトはFuncで使用される，もしくは定義される変数群の名前です．
これら単体だけでは特に意味はありません．
典型的には，x や yといったx-y座標を表す変数名などを使います．
Halide::Var x, y;
Funcは，変数やその他の関数からなるExprとして定義されます．
これは，任意の整数座標上で成り立つ関数として定義されています．
ここでは，x+yの値を持つExprが定義されており，またVarsのオブジェクトには演算子のオーバーロードが適切にされています．
x+yのような式はExprオブジェクトになります．
Halide::Expr e = x + y;
次に，画素x,y上で，画像がExpr eの値を持つように，Funcオブジェクトに定義を追加します．
左辺はFuncと変数Varを持つ定義です．
右辺はExprオブジェクトで同様の変数Varを持ちます．
gradient(x, y) = e;
実はこれらの定義は，以下のようにもかけます．
gradient(x, y) = x + y;
実際は，こちらのほうが一般的な書き方ですが，要素を網羅するために，中間記述であるExprをあえて示しました．
このように定義されたFuncのコードは，実際はまだ出力画像の計算自体をしていません．
この段階では，Func，Expr，Varは画像処理パイプラインの構造を表現するようにメモリに入っているだけです．
これは，メタプログラミングです．
このC++のプログラムは Halideプログラムをメモリ中に構築します．実際の画素データの計算は次で行います．
Funcを実現するためにrealizeメソッドを呼びましょう．
これは，定義したパイプラインを実装するコードをJITコンパイルします．
（この場合はFuncは1段しかないパイプラインです．）
そのとき，HalideにFuncを評価する定義域（上記で定義されたxとyの範囲）を教える必要があります．
それは出力画像の解像度にもなります．
また，Halide.hにより，基本的な画像構造のためのテンプレートを用いた型が提供されます．
さて，800×600の画像を作るとしたら以下のようになります．
Halide::Image<int32_t> output = gradient.realize(800, 600);
Halideは型推論機能を持ちます．
Varオブジェクトは32bit整数であるため，Exprオブジェクト'x+y'もまた32ビット整数です．
また，gradient関数は32bit画像です．
Halideの型とキャストのルールはCと同様です．
なお，このrealizeは，関数を順次実行してくれるスケジューリングメソットです．
より詳細なスケジューラは次以降のチュートリアルで紹介されます．
最後のForループの２重ループは実際の出力を確認するようのループです．
printfデバッグ結果で確認しましょう．
for (int j = 0; j < output.height(); j++) {
for (int i = 0; i < output.width(); i++) {
// We can access a pixel of an Image object using similar
// syntax to defining and using functions.
if (output(i, j) != i + j) {
printf("Something went wrong!\n"
"Pixel %d, %d was supposed to be %d, but instead it's %d\n",
i, j, i+j, output(i, j));
return -1;
}
}
}
これですべてが完了しました．
Funcを定義して，realizeをコールする．これで，画像を生成するためのマシンコードが生成され，そして実行されました．
printf("Success!\n");
02 画像を処理してみよう
画像を読み込んで，輝度値を1.5倍にするサンプルです．
Imageは画像入出力のためのオブジェクトです．
画像IOのために，libpngが必要です．
libpngとzlib（libpngのコンパイルに必要）をダウンロードしてコンパイルしてください．
両者ともにinclude, libをパスに通した後，プロジェクトファイルにlibpngとzlibのライブラリをリンクするように記述してください．
また，読み込む画像を"images/rgb.png"でハードコーディングしてあるため必要に応じて変えてください．
03 Halideで生成されたコードを確認しよう
実際にHalideでコンパイルされて生成したコードを確認する方法を知りましょう．
04 トレース，print，print_whenによるデバッグ
Halide専用printfデバッグを覚えましょう．
05 Vectorize, parallelize, unroll, tileによるスケジューリング
高速なバイナリを動かす仕組みです．
グラディエント関数をいくつかの方法で定義・スケジューリングします．
どの順番で画素が計算されていくかに注目してください．
このメソッドで，実際のスケジュールを表示します
this.print_loop_nest();
06 任意の定義域でFuncを実現しよう
(0,0)から始まるわけではない，任意の矩形で定義域を定義する方法を学びましょう
07 多段パイプラインを作ってみよう
生産者と消費者モデルによる分散強調型の計算の実行順序を学びます．
これが分かって初めて，冒頭の平均値フィルタが作れます．
08 多段パイプラインのスケジューリングを覚えよう
上のスケジュールです．
09 Funcが複数回実行される場合や定義がアップデートされる場合やリダクション処理を覚えよう
何度も実行されたり，定義が変わったり，リダクション処理がされる場合です．
10 ahead-of-time (AOT)を覚えよう(パート１)
ahead-of-time (AOT)コンパイラ→事前コンパイラ→事前にコンパイルするコンパイラ→つまりgcc,vc,iccのようないつものコンパイラ．
対義語は，javaなどのJITコンパイラ．
HalideはJITコンパイラとして動作しているが，静的ライブラリをコンパイルして生成するようにしても使えますよということを学ぶところ．
10 ahead-of-time (AOT)を覚えよう(パート２)
上に引き続き，コンパイルしたものをリンクして使ってみようという話し．
11 クロスプラットフォーム
Halideはクロスコンパイラとしてもちゃんとうごくよって話し．x86はもちろんのこと，ARMやGPUまで．
FPGAは下にあるdarkroomという言語のほうがいい．同じDSLだけど設計指針が違う．消費電力の最小化．
12 GPU
上に引き続きGPUが使えますよという話し．
GPUのCUDAカーネルのプログラミングも若干関数型言語っぽい仕様が入ってますが，もっと使いやすくなってますよ．
13 タプル
タプルって便利ですよね．なんでC++にデフォルトで無いんでしょ．
関数型言語は，Listとかタプルがないと書きにくくてしょうががありません．
14 Halideの型システム
Halideの型，ここで公開するから覚えてね．
Expr u8 = cast<uint8_t>(x);
Expr u16 = cast<uint16_t>(x);
Expr u32 = cast<uint32_t>(x);
Expr u64 = cast<uint64_t>(x);
Expr s8 = cast<int8_t>(x);
Expr s16 = cast<int16_t>(x);
Expr s32 = cast<int32_t>(x);
Expr s64 = cast<int64_t>(x);
Expr f32 = cast<float>(x);
Expr f64 = cast<double>(x);
デザインパターン：ジェネレータ（パート１）
デザインパターンのジェネレータとしてどうやってHalideの言語で記述するかが書いてあります．
再利用を促します．
デザインパターン：ジェネレータ（パート２）
上の続き
RGBの画像とそのメモリレイアウト（パート１）
カラー画像のデータはRGBRGBRGB...と並んでいるのが当たり前だと思っていませんか？
状況に応じて事前に並び替えをしたほうが高速に動作することが多々あります．
それを覚えましょう．
RGBの画像とそのメモリレイアウト（パート２）
その２です．
まとめ
画像処理プログラミング言語Halideの紹介をしました．
（vs OpenCVアドベントカレンダーの一部なのに，vs OpenCVのおまけ感．．．）
あと，画像処理専用のプログラミング言語の研究を始めました．
並列化やSIMD化といったチューニングをチューニングに終わらせずに研究にしたい人や専用言語の設計をしたい人，一緒にやりましょう！twitter @fukushima1981
日本語情報ページ
- https://tech.d-itlab.co.jp/programming/567/
- http://derivecv.tumblr.com/post/70538301269
- http://nebuta.hatenablog.com/entry/2013/08/04/135932
その他の画像処理専用プログラミング言語
Darkroom FPGAやASIC用の画像処理言語．