メモブログ

行列計算の並列処理でCUDAを使ったものはありますが、OpenCLを使ったものはどうかなと探してみると、clMathというものがありました。
clMathは、現在clBLASとclFFTがあり、行列計算や高速フーリエ変換をOpenCLを使ってできるようです。

今回は、clBLASの使い方について簡単にメモっておきます。

OpenCLはインストール済みと仮定します。
githubにclMathのデータがあります。

[1] https://github.com/clMathLibraries

ユーザ登録が必要かどうかはわかりませんが、自分の場合はgithubにユーザ登録して、clBLASのDownload ZIPをダウンロードしました。

clblas-develop.zipを適当なフォルダに解凍して、CMake(cmake-gui)を立ち上げます。

clBLAS-develop/srcフォルダをソースフォルダに設定し、適当にclBLAS-develop/buildフォルダを作成してそこをビルドフォルダに設定します。
Configureを押すと、ACMLが見つからないというエラーと、Boostが見つからないというエラーがでていました。

ACMLとBoostは、サンプルプログラム？等に必要な雰囲気でなくてもかまわない気がしましたが、一応入れることにしました。

ACMLは、下記のページからacml5.3.1-win64.exe PGI fortran用をダウンロードしてインストールしました。

[2] http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/acml-downloads-resources/

ACML_ROOTにF:/AMD/acml5.3.1/win64といったインストールフォルダを指定すればACMLのエラーは消えます。
Boostは、省略します。

Generateを押せば、buildフォルダにclBLAS.slnができていますので、ビルドすればエラーなく成功しました。

ライブラリを使用するときは、clBLAS-develop/src及びclBLAS-develop/src/includeにあるヘッダファイルと、clBLAS-develop/build/library/ReleaseフォルダにあるclBLAS.libとclBLAS.dllがあれば良さそうです。

一番基本となる行列同士の掛け算は、[3]の下のほうにサンプルソースコードがあります。

[3] https://github.com/clMathLibraries/clBLAS

APIの説明は、[4]にあります。

[4] http://clmathlibraries.github.io/clBLAS/

[3]では、clblasSgemmという関数を使って行列の掛け算を計算しています。

Sはsingle precision、gemmはGeneral matrix-matrix multiplicationの略で、double型で計算をしたい場合はclblasDgemmという関数になります。

で、引数がかなりわかりにくかったです。APIドキュメントには書かれていますが、ひとつの関数で下記のような計算を行うことができるようになっています。


[4]より引用

α、βは複素数のバージョンもあるのですが、上記の2関数は実数用のものになっています。
C=αAB+βCで、AやBは引数によって転置させて計算することも可能となっています。
つまり、普通掛け算というとC=ABだと思いますが、これを行うためにはα=1.0、β=0.0に設定しないとおかしな結果になってしまいます。

それ以外は、引数の説明に気をつけながら値を渡せば特に問題はないかと思います。

実行結果

自作行列掛け算と、clBLASを使用した場合(clblasDgemm)の計算時間を比較してみました。
1024x1024の行列A,B（値は-1.0～1.0の一様乱数）を掛けて、行列Cに代入します。

自作行列掛け算の時間：9629[ms]
clBLAS(転送時間含む)：55[ms]

自作行列クラスでデータ生成して、それを一時配列にコピーして、GPU上にメモリを作成して、CPU->GPUにデータ転送して、掛け算計算して、GPU->CPUにデータ転送して、自作クラスにデータコピーして結果を返すまでの時間が上記の結果でした。
一応計算結果をファイルに出力して比較してみましたが、どちらも同じ結果になっていたので合っていると思います。
無駄なデータ転送時間を入れてもこの差はさすが並列処理ですね。まぁ比較対象がひどすぎるというのもありますが。

とりあえず後はGEMV - General matrix-Vector multiplicationとかを試してみようとは思います。
○○solveとかあるので逆行列とかも計算できそうな気がしますが、まだやり方はわかっていません。

前回の続きです。

シャドウマッピングの手法で、影領域であれば加算しない、という方法で試してみました。

実行結果



一応前回の結果と比べると、遮られているところが少し暗くはなっているのですが、これであっているんでしょうかね。

前回の続きです。
slideshareというサイトでPPTを共有できるようですね。そのうち試してみたいです。



実行結果

物体が手前にあるTk&ltT1場合を赤色、物体がスポットライト内T1&ltTk&ltT2を緑色、物体がスポットライトの外側にある場合はそのままレンダリングした結果です（物体の色は0にしています）。



物体に当たったらそこで積分をやめますので、物体があるところは暗くなります。



遮られて暗くなるのはシャドウボリュームとか使えばいいんでしょうかね。シャドウマッピングでもできそうですがレイを進める度にデプス値の比較などは面倒そうです。



暗めのレンダリング結果＋スポットライトです。それっぽいような何か薄いような。

今はdt=0.005ぐらいにして計算しているので、かなり時間がかかっています。参考サイトの方では、点光源の場合ですが数値計算ではなく積分の式を解くことで一気に計算できています。
以前参考にしたときは式をスルーしていたのですが、今見ると理解できるので少しはましになっているんでしょうかね。

ちなみに参考サイトで2次元FEMを使った変形プログラムが公開されていました。
背景をRGBA=(0,0,0,0)にしたRLE圧縮なしのTGAファイルを読み込ませたらちゃんと変形して面白かったです。
以前のBox2Dの三角形分割もこれぐらいきれいにできれば良さそうですね。

http://blog.mmacklin.com/2012/06/27/2d-fem/

スポットライトのネタを二つ程．文字いっぱいの図を貼り付けるぐらいならPPT上げた方が良いですかね･･･．

1．円錐とレイの衝突判定

円柱とレイの衝突判定は[2]に載っていたのですが，円錐は載っていなかったので計算してみました．





結構計算量多そうなので，調べればもっといい方法がありそうですね．

実行結果

奥の面までの距離t2と手前の面までの距離t1を使って，t2-t1を表示しています．
衝突しない場合は-1にしています．



厚いところは白っぽくなっているので一応あっていそうです。
次はレイマーチングでボリュームレンダリング的なことでしょうかね。

2．スポットライトの指数及び減衰項をもう少しわかりやすく設定する方法

スポットライトの広がりに対する指数pや，点光源もそうですが減衰項k0, k1, k2をどれくらいに設定すれば良いのか直感的にはわからず，試行錯誤的に良さそうな数値を設定していたので，もう少しわかりやすい設定方法を考えてみました．

入力は，スポットライトの円錐の高さhと，広がりθ，始点の光の強さC0と，底面の中心の光の強さC1と，円錐の横の端に来たときの減衰率C2とします．







謎なパラメータは極力無くしたかったのですが，仕方がないですかね．
4つの値を決めるために入力が5+1個になってしまっていますが，まぁ適当に試行錯誤するよりかはわかりやすいのではないかと思います。

実行結果

h=2.0, θ=15°，C0=2.0, C1=1.0, C2=0.01, R=0.9で計算した結果，p=132.84, k0=0.5, k1=0.025, k2=0.1125になりました．実際どういう明るさになるかは次回テストですかね．





[1] http://www.densu.jp/pinpo/kukan4.pdf
[2] Eric Lengyel著，狩野智英訳：「ゲームプログラミングのための3Dグラフィックス数学」，株式会社ボーンデジタル

参考文献[1]などに、複素数の屈折率が出てきます。
複素数ってだけで何か難しそうに感じて困ったものです。

というわけで、演算子のオーバーロードの勉強も兼ねて複素数のクラスを作ることにしました。STLか何かにあった気がするので自分で実装する意味はありませんが。

[2]を参考にして実装した結果、以下のような感じになりました。

MyComplex.h

上記ソースコードで何か問題があっても責任はとりません。
たぶん問題ないと思うのですが、オーバーロードとか慣れていないので変なところがあるかもしれません。

本題のフレネルの式の計算ですが、まず入射角から屈折角を計算する必要があります。
スネルの法則[3]を使いますが、n2/n1って複素数の時どうなるのかよくわかりませんでしたが、とりあえず割り算してその複素数の大きさ|n2/n1|でいいようでした。

次は反射率の計算ですが、p波とかs波とかよくわかってないですがとりあえず複素数のまま計算します。

double getReflectRate(CMyComplex n1, CMyComplex n2, const double LH)
{
    CMyComplex value;
    // 入射角（ラジアン）
    double alpha = acos(LH);
    double sin1, sin2, beta, rp, rs, r;
    // 屈折角（ラジアン）
    sin1 = sin(alpha);
    sin2 = sin1*(n1/n2).Length;
    beta = asin(sin2);

    // rp
    // 分子: n1cosB-n2cosA
    // 分母: n1cosB+n2cosA
    // 計算
    value = ((n1*cos(beta))-(n2*cos(alpha)))/((n1*cos(beta))+(n2*cos(alpha)));
    rp = value.Length*value.Length;

    // rs
    // 分子: n1cosA-n2cosB 
    // 分母: n1cosA+n2cosB 
    // 計算
    value = (n1*cos(alpha)-n2*cos(beta))/(n1*cos(alpha)+n2*cos(beta));
    rs = value.Length*value.Length;

    // 偏光を考えない場合、平均
    r = 0.5*(rs+rp);
    return r;
}

演算子のオーバーロード便利ですね。ほぼ式そのままで計算できます。

実行結果

[1]に載っていた0.6-3.6iのグラフが下記のようになりました。



複素数が混じるとR+T=1にならない場合が多いので何か間違っているかもしれません。
上のグラフはとりあえずT=1.0-Rにしてごまかしています。

後はおまけでHDRビューアをガンマ補正できるように修正しました。
以下は異方性クック・トランスのテスト結果です。



あぴみくさんは元のデータが素晴らしいので例題としては逆にわかりにくいのかもしれないですね。
適当な計算でもきれいに見えてしまうという…。

[1] http://research.tri-ace.com/Data/BasicOfReflectance.ppt
[2] http://www.kairo-nyumon.com/electric_circuits.html
[3] http://en.wikipedia.org/wiki/Fresnel_equations

昨日の続きです。

やはり負の値も扱いたかったので、-1～1の場合、保存時には0～2にシフトして保存、読込時には-1～1の範囲に戻すことにしました。
シフトする量はヘッダのコメントなどがあればそこに入れたかったのですが、無駄が無さそうでしたので別ファイルで読書きすることにしました。

実行結果



オレオレ仕様になってしまいましたが、-1～1の範囲のように一応負の値でも取り扱えるようになりました。
その他、大きい画像は適当に縮小したり、マウスの中ボタンを押すとそこの画素値を表示するようにしたり、キーボード入力で0～1の範囲に設定したりと、とりあえず必要な機能は実装できたかなと思います。

Luminance HDRを使ったのでいいのですが、とりあえず線形変換だけでいいのですぐ開いて確認して閉じるための自分用ビューアを作ることにしました。
対応ファイルフォーマットはRGBEだけで、読込み／書込みも[1]のプログラムを使用させていただいております。

このプログラムを作成中に気付いたのですが、RGBEって負の値は扱えないんですね。浮動小数点＝float型だから負の値でもなんでも扱えるだろうと思って負の値を出力してました…。球面調和関数とかもHDRで負の値出力していたような気がするので確認してみないと…。
（rgbe.txtにはちゃんと0からって書いてました。ファイルフォーマット理解せずに適当に使ってると駄目ですね…。）

実行結果

[2]のサイトの画像でテストさせていただいております。



通常0～1の範囲ですが、0～21.75の間のデータになっているので、最初は暗いですが、マウスドラッグで範囲を調整すれば線形変換で一応暗いところ、明るいところの確認だけはできるようなものです。

[1] http://www.graphics.cornell.edu/online/formats/rgbe/
[2] http://www.anyhere.com/gward/hdrenc/pages/originals.html

そろそろ本来の目的であるボリュームレンダリングに入りたいところですが、今回も横道のアンサンブル平均です。
海面ネタを何回かやってたのですが、コースティクスの計算をしていなかったので、試してみることにしました。
海面にはFFTで作ってポリゴン出力したものを使います。あまり海面っぽくないですね。



ある視点からレンダリングして、深度値を見やすく出力すると下図のようになります。



光源からそれぞれのピクセルへのレイと法線ベクトルから、スネルの法則を使って屈折ベクトルを計算します。
このとき、RGBそれぞれに屈折率を変えて３本のレイに分割します。
そして、適当に決めた床平面（凹凸なし）との交点を計算して、RGBそれぞれに1/3の色を足しこみます。
全てのピクセルに対してこの処理を行った結果（1回のアンサンブル平均）が下図のようになります。



屈折率ですが、差が小さいと分光現象が現れなかったので、かなり差をつけてますので物理的には正しくありません。（R:1.33, G:1.38, B:1.44）

800回のアンサンブル平均で下図のようになりました。



これは床平面に集光する光の計算しかしていないので、海底や水面のレンダリングなどと合わせるともう少しそれっぽくなるのではないかなーと思います。

接続部をもう少し何とかできないかなということで短時間フーリエ変換のようなものを試してみます。

画像等に離散フーリエ変換を適用する場合に問題になるのが、境界で不連続になる点です。
本来、左端から右端までが1周期の信号で、これが繰り返される仮定になっていますので、左端の画素値と右端の画素値が異なっている時点で不連続ですし、同じ画素値でも勾配的に不連続になってしまったりしてあまりきれいにつながらなくなります。
前回は、境界部分にガウスフィルタをかけて平滑化しましたが、短時間フーリエ変換を用いて滑らかに繋がらないか試してみます。

まずは1次元でテストしてみます。

原信号256点を、2つ並べて512点にします。真ん中が不連続点になります。



0～255、128～383、256～511の3つの区間の256点の信号に分けて窓関数（今回はハミング窓）を掛けます。

それぞれをFFTして、ローパスフィルタやハイパスフィルタを掛けます。
短時間フーリエ変換の長所？ですが、信号ごとに異なるフィルタを掛けることが出来ます。今回は、真ん中の信号（不連続点を含む信号）だけ、強めにローパスフィルタを掛けてなるべく滑らかな信号にしてしまいます。
そして、それぞれをIFFTします。

3つの信号を全て足し合わして、真ん中の256点だけ取り出します[1]。
（参考サイトでは、逆ハミング窓を掛けてもう一度窓関数を掛けたりしています[2][3]）
比較として、512点に対してFFT、フィルタリング(全体に強いローパスフィルタ)、IFFTを行った場合を並べてみます。

高周波成分200個削除した場合（左：全体、右：短時間フーリエ変換）
　　　

高周波成分244個削除した場合（左：全体、右：短時間フーリエ変換）
　　　

全体にローパスフィルタを掛けた場合、境界部(中心と端)が多少振動してしまっています。
短時間フーリエ変換の場合は、一応多少はまともになっているように見えます。

次は2次元（画像）で試してみます。
シームレステクスチャにしたい画像を2x2の4枚並べます。



今回は、これを128x128の画像ごとに、64ピクセルずつずらしながら、計7x7枚の部分領域に分割しました。
そして、それぞれに2次元のハミング窓[4]を掛けて、FFT→フィルタリング(ハイパス＋ローパス)→IFFTして結果を足し合わせます。
フィルタリングは、不連続となる十字の部分により強いローパスフィルタを掛けるようにしています。

ハミング窓を掛けた部分画像の例



フィルタリング→IFFTした画像



結果をすべて足し合わせた画像



この中で、シームレステクスチャとして使用するのは真ん中の領域だけですので、それを切り出します。



ついでにこの画像に対してもFFT→フィルタリング→IFFTして、明るさを調整して、2x2に並べてみます。



残念ながら前回の結果とあまり変わりませんでした。

[1] http://www.cqpub.co.jp/hanbai/books/30/30911/30911_6syo.pdf
[2] http://junzo.sakura.ne.jp/c_sharp22/c_sharp22.htm
[3] http://www.me.cs.scitec.kobe-u.ac.jp/~takigu/jugyou/shingou/0210/ensyuu/Fourier.pdf
[4] http://www.kkaneko.com/seminar/20031128/abe_seminar_031128.ppt