2013年12月2日月曜日

『プログラミングの魔導書 Vol.3』の予約を開始しました!


Vol.2の発売から2年ほど経過してしまいましたが、本日、『プログラミングの魔導書 Vol.3』の予約を開始しました!


今回のテーマは「並行、並列、分散」です。
一般的なプログラミングでとくに難しいと言われるこの分野について、各分野の有識者の方々に執筆していただきました!

本書のコンテンツ、および執筆者は以下のとおりです:
  1. 序文(熊崎 宏樹)
  2. Lock-free入門(熊崎 宏樹)
  3. OpenACC(藤田 典久)
  4. ErlangとScalaにおけるアクターモデルの紹介(幾田 雅仁)
  5. C#の非同期処理(岩永 信之)
  6. Real World STM ~作って学ぶSTM~(石井 大海)
  7. データ並列への招待(shelarcy)
  8. 合成可能なメッセージパッシング ~Concurrent ML の紹介~(小笠原 啓)
  9. コルーチンスタイルプログラミング(高橋 晶)
  10. 画像検索入門(miyabiarts)

ご予約は以下のページから承っております。

書籍版は受注生産のみとなっておりますので、お早めにご予約ください!

2012年12月21日金曜日

Curl Noise追記

GPGPU Advent Calender@21日目ですが、完全にネタ切れのため前回書き忘れていた部分を補足してお茶を濁します!
by jirohcl

まず、前回説明し忘れた重要なnoise関数。
float noise(float2 p,__global float2* rt,__global unsigned char* rti)
{
 float2 floor_p = floor(p);//float2(floor(p.x),floor(p.y))
 int i=(int)floor_p.x, j=(int)floor_p.y;
 float2 n00=rt[hash_idx(i,j,rti)];
 float2 n10=rt[hash_idx(i+1,j,rti)];
 float2 n01=rt[hash_idx(i,j+1,rti)];
 float2 n11=rt[hash_idx(i+1,j+1,rti)];

 float2 pf = p - floor_p;
 float sx= pow(pf.x,3)*(10-pf.x*(15-pf.x*6));
      float sy= pow(pf.y,3)*(10-pf.y*(15-pf.y*6));

 return bilerp(   dot(pf,n00)                   , dot(pf-(float2)(1.f,0.f),n10),
                         dot(pf-(float2)(0.f,1.f),n01) , dot(pf-(float2)(1.f,1.f),n11),
                         sx, sy);
}

となっており、パーティクルの座標と4点少しずらした座標でbilerpして求めます。
またランダムテーブル自体も

//host
float spin_t = 0.5f*0.03f/0.1f*time_;
for(unsigned int i=0; i<RANDOM_SIZE; ++i){
 float theta=rnd_table_spin_[i]*spin_t;
 float c=std::cos(theta),s=std::sin(theta);
 rnd_table_[i][0]= c*rnd_table_orig_[i][0]+s*rnd_table_orig_[i][1];
 rnd_table_[i][1]=-s*rnd_table_orig_[i][0]+c*rnd_table_orig_[i][1];
}


このように円を書くように更新していく事で、うずまくような動きになります。(つまりCurl)
float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
               ^^^^^^^

ノイズゲインを上げることによってこの動きは顕著になります。(前回の↑部分がゲイン)

といった所で前回の説明忘れ部分は終了です。
次にパフォーマンス比較を忘れていたので行いました。
100k particles
GPU          :   10ms
CPU single  : 183ms
CPU omp    :   93ms

CPU singleはシングルスレッド。CPU ompはOpenMPを用いた並列化の結果になります。
大体想定どおりの結果です。
CPU ompが若干遅いですが2coreのアレなCPUを使ってるので現行の4core辺りなら40-60ms位に落ち着くのではないでしょうか。

Curl Noise自体が並列化し易く、かつ入力データが多い例なのでOpenCLでの高速化が上手くはまるアルゴリズムでしたのでこのような(理想的な)結果に落ち着きましたが、何でもかんでもGPGPUで高速化できるものではありません。
例えば入力が多く、並列化が容易であっても分岐が多くなるとGPUでの高速化は上手くはまらないのでアルゴリズム審美眼を磨いて、「これはGPUでうまくいく!」といった勘を養うことがGPGPUでは重要ではないかなぁといった所で僕のノープランgdgd GPGPU Advent Calenderを終了しようと思います。
atnd立ち上げた山田さんに多謝!

2012年12月13日木曜日

Curl Noise with OpenCL

どうも再び@jirohclです。そしてGPGPU Advent Calender@12日目です!
人数が少なく、下手をすればもう一度まわってくるらしいのですが完全にネタも尽きてはたしてどうなる事やら。
まぁ愚痴はアンドロメダ星雲にでもすっ飛ばしておくとして、今回はこんなものを実装してみましたネタです。

Curl Noiseと聞いてピンと来る方も多いと思いますが、今をときめくスクウェアエニックス様が作ったAgni's Philosophyでも使われていたパーティクルの技法です。
主な特徴としては、パーティクル同士の判定を行わないために超高速、剛体との衝突がなんかそれっぽく見える、といった所でしょうか。

これをどうOpenCLに持っていくかというと、パーティクルの更新が都合よく並列化に向いていたためupdate関数をカーネル化してやりました。
__kernel void update(__global float2*        in,
__global float2* rx,
__global float * t,
__global float2* rt,
__global unsigned char* rti){
int const id = get_global_id(0);
float2 mid = in[id];
mid = mid+0.5f*curl_vel(in[id],rx,t[0],rt,rti)*0.04f;
in[id] += curl_vel(mid,rx,t[0],rt,rti)*0.04f;
}

ここで特筆する点は乱数テーブル(rt)をホストから渡している所でしょうか。
この乱数テーブルは円状に配置された点をswapしてテーブル化してあります。

次に
float2 curl_vel(float2 x,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
float const delta_x = (1e-4);
float2 v;
v.x = -(potential((float2)(x.x,x.y+delta_x),rx,t,rt,rti)
- potential((float2)(x.x,x.y-delta_x),rx,t,rt,rti))/(2*delta_x);
v.y = (potential((float2)(x.x+delta_x,x.y),rx,t,rt,rti)
- potential((float2)(x.x-delta_x, x.y),rx,t,rt,rti))/(2*delta_x);
return v;
}

速度の計算です、パーティクルの現在位置から微小距離ずらし、ポテンシャル関数を用いて速度を求めます。
float potential(float2 px,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
float2 wind_velocity = (float2)(0.1f,0.5f);
float const inv_env_sqr = 1/sqrt(0.17f);
float numer= inv_env_sqr*cross2(px,wind_velocity);
float denom= inv_env_sqr;
//---------------hard coded rigid X( ---
float2 dx=px-rx[0];//centrer
float psi=cross2(px,rx[1]);
float phi=((sqrt(dx.x*dx.x)+sqrt(dx.y*dx.y))-0.06f);//radius;
float m=1/(sqrt(phi)+1e-18);
numer+=m*psi;
denom+=m;
//--------------------------------------
float base_psi=numer/denom;

float d=phi/0.1f;
float g=smooth_step(1.1-px.x);
float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
return base_psi+turb;
}

そしてこのポテンシャル関数の中で、回転するノイズ関数、剛体や風等の周りからの力の付与を計算することで流体を表現するらしいです。

実行結果


実装に疲れたので説明はおざなりですが中々面白い技術なので皆様試してみてはいかがでしょうか?

Windows用の実行ファイルになり、かつ他の環境で動くかどうか怪しいですが一応デモも置いておきます。
※PCが壊れた等の責任は負いかねますので悪しからず

2012年12月4日火曜日

Hello Oldschool GPGPU World

お久しぶりです。@jirohclです。
特に書くことも思い浮かばず日が過ぎ、気がつけばまたAdvent Calenderの季節になりました。

今回、私が参加するのはこちら、GPGPU Advent Calender@3日目となります。
この記事は温故知新と聞こえの良い事を言いつつ、ネタが思い浮かばなかった苦肉の策、
タイトルの通り古臭いGPGPUをやってみようと思います。

古臭いGPGPUとは何かと言うと、プログラマブルシェーダ、特にピクセルシェーダを用いて汎用計算をやってみましょうというアレです。
これを知ることでCUDA/OpenCL等は裏側で一体何をしているのか押し知ることが出来、また効率的な実装の方法も見えてくるはず。(建前)

まず、GPGPUをするために何が必要かというとご存知の通り
1.デバイスへの入力
2.処理関数
3.ホストへの出力
これら3つ。

まず、1の入力はテクスチャを使うのがポピュラーな方法です。
次に2の処理関数に関してはフラグメントシェーダーで頑張ります。
3はglReadPixels等のレンダリング後の画像を読み込む関数を使うことで結果を取り出します。

ではまずデバイス側のコード(行列積)を見てみましょう。
//vs
uniform mat4 mvp;

in vec3 in_pos;
in vec2 in_tex;

out vec2 out_tex;

void main(void)
{
gl_Position = mvp * vec4(gl_Vertex.xyz,0.f);
out_tex = in_tex;
}

//fs
uniform int mat_size;

uniform sampler2D mat_1;
uniform sampler2D mat_2;

in vec2 texcoord;

void main (void)
{
float result = 0.f;
for(int i = 0;i<mat_size;++i)
{
result += texture2D( mat_1,vec2(i/(float)mat_size,texcoord.y) )
* texture2D( mat_2,vec2(texcoord.y,i/(float)mat_size) );
}
gl_FragColor = vec4(result,1.f);
}

バーテクスシェーダ自体は何も特別なことは行わず、フラグメントシェーダで積を演算しています。
ここでフラグメントシェーダでの演算ですが実は一つ大きな問題があり、
フラグメントシェーダー内での処理には時間制限があるために行列のサイズが大きくなると計算が正常に終了しなくなります。OpenCLを使ったことがある方はきっと身に覚えがあるかと思います。アレです。

では次にホストコードを、と思いましたが全部羅列すると酷いことになるので一部
//データ作成
GLuint mat_lhs;
::glGenTextures( 1, &mat_lhs );
::glBindTexture( GL_TEXTURE_2D, mat_lhs );
::glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,mat_size,mat_size,0,GL_LUMINANCE,GL_FLOAT,data);
::glBindTexture( GL_TEXTURE_2D, 0 );

//実行
::glEnableClientState(GL_VERTEX_ARRAY);
::glEnableClientState(GL_TEXTURE_COORD_ARRAY);
::glBindBuffer(GL_ARRAY_BUFFER, texed_quad);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, index_buffer);
::glDrawElements(GL_TRIANGLES,buffer_size/sizeof(GLuint),GL_UNSIGNED_INT,0);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
::glBindBuffer(GL_ARRAY_BUFFER, 0);
::glDisableClientState(GL_VERTEX_ARRAY);
::glDisableClientState(GL_TEXTURE_COORD_ARRAY);
::glFlush();
::glReadPixels(0,0,mat_size,mat_size,GL_RGBA,GL_FLOAT,result);

デバイスへのデータ転送、カーネルの実行する”だけ”でもこれだけの決まりごとの命令が必要になります。
また、この他にもシェーダの初期化、ユニフォーム、イン/アウトへの関連付け、等々面倒くさい作業が目白押しとなっております。

CUDA/OpenCLが難しそうと思った方、古くさいGPGPUを動かすためのコードを見ると、現在の環境は大分めぐまれてると思いませんか?
思った方は是非GPGPU Advent Calenderに登録して記事を書いてみましょう!

2012年3月26日月曜日

Sparse Voxel Octree(SVO)

大体月1更新の@jirohclです。
今回はここ最近ステマをしているSparse Voxel Octree(SVO)についてです。

SVOとはVoxelOctreeをメモリ効率を上げてマッピングする手法で、
基本的なボクセルデータvec_t rgbaの上位bitをアドレスとして設定します。
これによって空白ボクセル部にメモリを取られず、高効率でボクセルデータを格納することが出来ます。

svo_image

(64bit type):
X: [ID:2][node x:10][node y:4] //ID : leafの所持フラグ
Y:[node y:6][node z:10]
Z:[color r:8][color g:8]
W:[color b:8][color a:8]
※数値はbit

SVOはRaytraceとの親和性が高く、また法線マップなどを使わずに詳細なディティールを表現できるため、
tech6等の次世代のゲームエンジンで積極的に取り入れられようとしています。

しかし、この実装通りに行うと、OpenCLでレイトレースするためには再帰関数がサポートされない関係で
bool hit_geom_level0(ray r,float *t,__global uint4 * geom)
bool hit_geom_level1(ray r,float *t,__global uint4 * geom)
と自前で再帰定義、呼び出しを行わなくてはいけません。
mipmap levelがそこまで高くなければ問題ありませんが、純粋にOctreeを作っていては納得できる分割数に至るまでにかなり階層が深くなってしまいます。
そこでmipmap level N(最大レベル)の分割数をある程度上げ、階層を固定化することでこの問題はある程度解決します。

また、SVOを使用することで、レイとジオメトリとの交差判定で得た距離をLevel of Detailに用いることで比較的簡単にCone Tracingの効果を得られる事も一応記しておきます(まだしっかりと実験していないので詳細は無しで)。

SVO自体の可視化
http://www.youtube.com/watch?v=CJaN38CqHpY&context=C43a407aADvjVQa1PpcFM1-KZzH_fwDhkQS7JnlNktMbYClnWVAl0=
上データを使用した単純なパストレース
http://www.youtube.com/watch?v=M7n8Cy7hGVg&context=C46616bcADvjVQa1PpcFM1-KZzH_fwDnaz0vsZ61jyfCMC7vbypAM=

若干取り留めのない文章になってしまいましたが、この文ですこしでも「よし!SVOやる!」という方が増えてくれれば幸いです。
次回未定。

2012年2月24日金曜日

Raycasting

3D CG Programmingnの時間です。

さて、前回ボクセルシーンを生成した訳ですが。
今回はそのデータを用いてRayCastを行ってみました。

RayCastとは:
RayMarchingとほぼ同じ意味で使われますが、簡単に言えばボリュームデータに対し、線分を少しづつ進めて交差判定を行う方法の事です。
用途としては、1[煙等、透過するオブジェクトのレンダリング]
2[X線画像等の体積を持ったデータの可視化]
ですが、最近リアルタイムレンダリングの界隈でライティング計算に用いるようになってきました。

流行に乗るために簡単な実装を行いました。
video

で、問題のRaycasting部分
//in raycast function
//*r is ray
int const size = res * center.w;
float4 const mini = center - size*0.5f;

for(int i = 1;1;++i)
{
float4 const p = r.orig_ + r.dir_*i*center.w;
float4 const r_vox_pos = ( p - mini )/center.w;
int4 const r_vox_pos_i = convert_int(r_vox_pos);

//exit without hit
if(r_vox_pos_i.x < 0 || r_vox_pos_i.x >= res-1 ||
r_vox_pos_i.y < 0 || r_vox_pos_i.y >= res-1 ||
r_vox_pos_i.z < 0 || r_vox_pos_i.z >= res)
return false;

sphere sp;
sp.pos_ = convert_float(r_vox_pos_i)*center.w - size/2.f + center.w*0.5f;
sp.r_ = center.w*0.5f;
float t;

if( sphere_intersect(r,sp,&t) )
{
// hit actions
return true;
}

ここで交差した部分のボリュームデータが半透明な物であれば線分の開始点を書き換え(必要であればライティングも)て再びRaycastを行えばボリュームレンダリングになります。
また、今回は交差判定にbounding sphereを用いてますが、一般的にはAABBを使った方がクオリティの問題で良いと思います。

ではまた何かネタが出来たときに。

2012年2月1日水曜日

プログラミングの魔導書 Vol.3 始動

高橋 晶です。

プログラミングの魔導書 Vol.3のプロジェクトが本格始動しました。Vol.3のテーマは「並行/並列/分散」です。
2012年6月以降に発売する予定で進めています。

今回も、前回に負けないくらいのメンバーと記事を取り揃える予定ですので、発売まで楽しみにお待ちください。