仮眠プログラマーのつぶやき

自分がプログラムやっていて、思いついたことをつぶやいていきます。 2025年からzennに移行

自作ゲームやツール、ソースなどを公開しております。
①ポンコツ自動車シュライシュラー
DOWNLOAD
②流体力学ソース付き
汚いほうDOWNLOAD
綺麗なほうDOWNLOAD
③ミニスタヲズ
DOWNLOAD
④地下鉄でGO
DOWNLOAD
⑤ババドン
DOWNLOAD
⑥圧縮拳(ツール)
DOWNLOAD
⑦複写拳
DOWNLOAD
⑧布シミュレーション
DOWNLOAD
⑨minecraft巨大電卓地形データ
DOWNLOAD
⑩フリュードランダー
デジゲー博頒布α版
DOWNLOAD
⑪パズドラルート解析GPGPU版
DOWNLOAD
⑫ゲーム「流体de月面着陸」
DOWNLOAD

レイトレーシング

HLSLでリアルタイムレイトレーシング

のんびりと更新を続けて、今回でやっとライブドアブログ移転2回目の記事


まずHSPコンテストのことだが、今年は作成時間がとても少なく、締め切り直前にあたふたした感が否めない。
この前の記事で構想を温めていると自慢したはずの流体力学は、10/31日時点で全く自分の思い描いていた形に仕上げられないままであった。
今年の投稿は諦めるかと悩んだが、ブログの方でHSPで高速な演算ができることに期待をしているコメントを頂いていたことか後押しして、投稿することにした。
実際、自分にとってこの大きな発見を、一年も後に発表を引き延ばすことは耐えられるか危うかった。

投稿直前の10分でとりあえず動く形にして、付属として入れるはずだったリアルタイムレイトレーシングの方が完成度が高かったので、仕方なくそっちをメインということにしてサムネも作成した。
本当は「HSPでGPGPU」という名前にしたかったのだが、パソコンの時計を見たらもう58分になっていて、焦ってシフトキーを押す手が震えてしまい全部小文字になってしまった・・・
それとHSPでGPGPUは実はHLSLでGPGPUだった件

結局UPしたのは0:00を30秒くらいオーバーしていた・・・がなんとか受け付けてもらえたようだ。


そのような経緯があって、内容説明文や何から何まで完成度が1%となった。
本当はリアルタイムで数値流体力学シミュして、ゲームみたいに遊べるようにしたかったのだが。


そして11/18正午・・・、一次審査通過が決定!落ちた人には申し訳ないが、嬉しいけどこれじゃ不甲斐ない・・まぁ言い訳してもみっともない

来年こそは、時間たっぷりかけて作品を仕上げてやる!



また、参加者の作品を何十個か遊んでみたのだが、こういうゲーム自分もいつか作りたいと思ってたんだよねーというものが多く、各々自分の創造力を大切にしながらゲームを作っているなと言った感じだった。

やっぱりHSPコンテストの意義と言うか、最終目標は『HSPでこんなことが出来るのか!すげぇ!』という感動を初見に、あわよくば既存のユーザーにも与えるという事だと思う。

今年はそんな「すげぇ」と言わせてくれるような作品が多かったように思えるし、「将来すげぇものに化けそうだな」というような作品を作ってくる人も多かった。

私も少しながら、評価コメントの方で「すごい」と言わせられたので当面の目標は達成したかなと思っている。





さて本題だ

HLSLのピクセルシェーダを使って、512×512ピクセルのスクリーンにリアルタイムレイトレーシングをする方法だが、その手順を書いていこうと思う。
 
 前提としてHSPがインストールされていることと、easy3Dのver5.2.3.3が入っていることが条件。(なぜeasy3Dの過去バージョンを使っているかは後述)

今回もいきなりLEVEL100のモンスターを倒すようなマネはしないで、しっかりLEVEL1から、一番敷居の低いアルゴリズムで基礎の基礎からはじめていこう~!


【導入・・・的な】

レイトレのアルゴリズムといえば、


(1)視線は、スクリーン上のある画素を通って物体方向に向かう。この視線と交差する別の物体があるかどうか調べる。
(2)交差する物体が存在するなら、視線と物体との交点を求める。交差する物体が複数ある場合には、すべての物体について交点を求める。交点がない場合には背景の色とする。
(3)視線と交点との距離を求め、最も視点に近い物体を抽出する(隠面消去したことになる)。
(4)可視物体の輝度の計算をする。このとき、光源と交点を結ぶ線分と交差する物体を調べ、影の有無を判定する。
(5)反射・屈折がある物体なら、反射方向、屈折方向を求め、これらの方向を視線とみなして②③の処理を行い、屈折・反射して見える物体を抽出する。
以上の処理をスクリーン上のすべての画素について行う。
(5)の処理では、視線(レイ)と交差する物体が不透明物体の場合には、陰影計算によって物体の色を計算し、画素に塗る。一方、鏡面のような反射面に視線が交差した場合には、さらに反射方向を追跡する。また透明物体の場合には、屈折方向も追跡する。反射と屈折方向の両方を求める必要がある場合が多い。光線の反射や屈折方向を求めるには、その面の向き(すなわち法線)を求める必要がある。また、屈折方向を求めるには、法線方向のみだけではなく、その物体の屈折率が必要である。視線が反射と屈折の2本の線に分岐することから、追跡経路は2分木表現される。


http://www.rsch.tuis.ac.jp/~naka/naka/scola/member/chapter5_html/chapter5.html
に解説してあるとおりで間違いない。

だがこれはいわば
LEVEL50くらい手強いモンスターである。

ここからLEVEL1に下げて考えるには徹底的に簡略化する必要がある。


まず(1)は削る所がない。
(2)は、物体の形で交点を求めるプログラムが変わってくるが、ここは一番簡単な「球」の1種類でいいだろう。
背景色は最後の視線の方向で色と明るさ決定。(単色はダメ)
(3)ここも削る所がない。
(4)影の有無は面倒なので消去。(全て影なしとして計算)
(5)反射はあってもよいが、屈折は屈折率とかで面倒。光の吸収率も面倒なので100%光が反射するという設定でいこう。これで不気味な2分木とかいう経路演算は必要なくなる。これで相当簡単になるはず。

さらにレイトレではなぜかチェック模様の地面があるのが普通とされているのだが、地面と視線との当たり判定がコレまた面倒なので地面も取っ払い、球だけの世界で良しとすればどうだろう・・LEVEL1近くなったんじゃなかろうか





この方法をフローチャートにまとめるとこうなる。

ritrhrtyt

うん。かなり簡素化したな

ここで視線ベクトルと視点(位置)ベクトルは、それぞれ以下のベクトルtV、Mのこととして考えて欲しい。
4b8b9221.png

【テクスチャの準備】
GPGPUをやるということは、テクスチャを変数に見立てるということに他ならない。
ここで、レイトレで必要な変数と要素数を全てあげると、

・球の位置座標(x、y、z)と球の半径r = 4×球の数
・視点ベクトルM(x、y、z成分) = 3×全画素数
・視線ベクトルtV(tx、ty、tz成分) = 3×全画素数
tは上の図で言うカメラの位置から交点までの距離。
・視線と交差した球の識別IDを格納する変数 = 1×全画素数
・レンダリング用のバッファのr、g、b成分 = 3×全画素数


だ。

E3DCreateRenderTargetTexture命令で、テクスチャフォーマットはA32B32G32R32Fにしてバッファを作ると、
1ピクセルにR,G,B,Aの4つの情報を格納できる。さらに一つ一つの精度は単精度floatでレイトレには十分といえる精度だ。

さて上の変数をテクスチャ数がなるべく少なくなるように格納すると


・球の数だけ画素があるテクスチャ  ×1枚 : 球の位置x、y、z、半径rを格納
・512×512ピクセルのテクスチャ    ×3枚 : tV(tx,ty,tz)とIDで1枚(①)。M(x,y,z)で1枚(②)。レンダリングバッファのr、g、bで1枚(③)

ということで4枚のGPUバッファを用意すればいいことになる。


百聞は一見にしかずだ


32個の球の情報を格納したテクスチャは
kyu
こんなこんじだ。(大きさは64×64ピクセル)(節約しようと思えば8×8のバッファで事足りるが・・)


【フローチャート】
また、下の図は2枚+レンダリング用バッファのテクスチャを、レイトレーシングが完結するまでの各タイムステップごとにキャプチャしたものである。


a0a1
a2
a3

縦を①~③で区切って、横を(A)~(I)で区切ったとして

(A)は初期視線ベクトルを代入した時点での画像で、(I)は完成時の画像

①の列の画像は ・・まぁ書いてあるが・・ 視線ベクトル(x、y、z)成分がこの画像のr、g、bに対応してある。4つめの成分であるaの情報は割愛してある(表現できないし)




(A)の①
a0
フローチャートで言うと一番上
「視点、視線ベクトルxyz成分決定」
最初に飛ばす視線ベクトルtVを決定する。初期視点ベクトルMは全画素同値。
最初のtVは無限遠の空に飛んでいく視線なので
ここではtをできるだけ大きい数にとるのが理想(実際は65536を格納)
青地に、左が赤、上が緑っぽいグラデーションがかかっているのが分かるであろうか

赤はrつまりxの値が高いということを表している。
緑はgつまりyの値が高いということを表している。

赤い部分は、視線ベクトルが左方向に傾いていることを表している。
緑色の部分は、視線ベクトルが上方向に傾いていることを表している。
青はz方向なので手前から奥に向かう視線であることを表している。

マイナスの値は0にクランプされて表示されているが、クランプされているのは表示の方だけである。


(A)の②
a1
ここではまだ、球との当たり判定をしていないので、位置ベクトル=カメラの位置=(0,0,-7)で真っ暗な画像になっている。



(B)の①
b0
フローチャートで言うと2番目
「球と視線の交点を求める」
全ての球について当たり判定を行なう。

まず、一つ目の球でt1V(t1x,t1y,t1z)を求める。
これをr,g,b成分に代入するのだが、t1Vの解が求まらない場合(つまり交点がない場合)は代入を無視。
次に二つ目の球でt2V(t2x,t2y,t2z)を求める。
t1>t2なら代入。t1の値はr,g,b成分から求めればいい。(  t1=sqrt(r*r+g*g+b*b)  )

これを全ての球で行ない、最終的に一番近い交点への視線ベクトルが代入されたことになる。

そしてこの画像は全ての球との当たり判定後の状態。
球のところが黒く抜けているが、これは黒い部分に代入されているベクトル情報が例えば12.1×(x,y,z)みたくなっているので、無限遠に飛ぶ65536×(x,y,z)と比べて値がとても小さい → 相対的に黒く見える、という状態だ。(ここで(x,y,z)は単位ベクトルである)


(C)の②
c1
フローチャートで言うと3番目
「一番近い交点を視点ベクトルに代入」
交点がない場合、それまでの値は保持されたまま。
という訳で、視点ベクトルに代入された結果が(C)


(C)の①
c0
フローチャート4番目の「反射ベクトルを視点ベクトルに代入」
反射ベクトルの求め方はココで解説したとおり。





以降フローチャートの2~4が3回ループして、反射回数が増えていく・・・



(I)の③
i3
フローチャートの5番
「最後の視線の方向で背景色決定」
3回反射処理した結果、全ての画素の視線の先に球はなく、無限遠の背景が残るのみとなったら、視線ベクトルから背景色を決定する。

ここで、ライト(光源)ベクトルがあると便利である。
光源ベクトルと視線ベクトルの内積で1.0~-1.0の値が得られ、1.0が明るい水色の空、-1.0が暗い灰水色の空としてrgbを決定すれば意外とリアルな空になるはずである。





最後に少し応用で、新たに変数を作り反射した回数を記憶させておけば、反射による光の半減も再現できる。
あいにく視点(位置) ベクトルを格納しているテクスチャにはa成分だけ空きがある。
③の列のレンダリング画像は、実はその光の半減を取り入れている。





ひと通りフローチャート説明がこれで終了
カスタムシェーダは自分でくんでもらうとして 、とりあえずHSPでeasy3Dを介したHLSLによるリアルタイムレイトレーシングはこんな感じで動いている。

カスタムシェーダのソースとHSPのソースは
http://loda.jp/babadom/?id=114 
 からダウンロードできる。





だがこのソースははっきり言ってまだまだ改良の余地がありまくる。
今回の方法では、全てのドットで同じ処理をしないといけない制約があるので
例えば1ループ目ですでに球との交点がないドットでも、3ループ目まで当たり判定の処理をし続けないといけない。

だから全体的に無駄な処理が多い。
逆にまだまだ高速化が期待できるというわけだ。

どこまで挑戦できるか分からないが、飽くなき高速化への道はまだまだ続きそうだ。




ところでもしグラフィックボードが32bit浮動小数点の計算に対応していないとどうなるのか・・・
一応起動して、16bit浮動小数点のテクスチャで計算されるらしい。
だが少し計算にズレが生じるらしく、完成画像が少し荒くなる↓
msafr





あと前半の方でeasy3Dのバージョンが 5.2.3.3 を使用しないといけないと書いてあるが、
始めに「hspでgpgpu」を知り合いのパソコンで起動してもらおうとしたらなぜか起動しなく、その時easy3D.dllの最新バージョンを使っていた。
逆にインストール不要形式のバージョン 5.2.3.3で試してみたら動いたので、こちらを推奨しているというわけ

ちなみにその人のパソコンはwindows7でかなり最新のノートパソコンなのでDirectxは10以上だった気がするから、起動しなかった原因はよくわからない・・おちゃっこさん究明お願いします。



次回:流体力学のプログラムを作りたい!その2


PS
なんかコンテストTVに「hspでgpgpu」が紹介されてました!!審査員の方々ありがとうございます!!
あと2日で最終選考の結果が発表ですね。楽しみです

レイトレーシングアルゴリズム実験③

久しぶりに更新


夏休みに新しいノートパソコン買って、今まさにcore i5 450Mのスピードに感動しているところだ

前のノーパの8倍は早い!

最近暇ができてきたのでまたfortranでレイトレしてみた



パソコンも新しくしてスペックも上がり計算に2時間しかかからなかった


ce54e1ba.png




これもピクセルごとのRGBをfortranで出力したのをHSPで読み込んでスクリーン出力したものだ



相変わらず進歩しない風景だが一応ソースを載っけてみる



まずはfortranのソース↓(ちなみにコンパイルはg95っての使ってソースのファイル名はh.f90でした)






program IO
integer yoko,tate,kaiso,mate,jimenmate,i,k,m,mcnt1,htmate,mhtmate,hhpi,ii,kk,ll,h,jkl,jjk
real cb0,cb1,cb2,ca0,ca1,ca2,hsa0,hsa1,hsa2,hsb0,hsb1,hsb2,q0,q1,q2,qt2
real t,w,b,c,d4,vmr,vmg,vmb,vmi,kai,x,z,drtga,naiseki,light0,light1,light2
integer mtset(3,640,480)
real g(3,3)
real,allocatable,dimension(:,:) :: kyu
real,allocatable,dimension(:,:) :: ja
real,allocatable,dimension(:,:) :: jb1
real,allocatable,dimension(:,:) :: jb2
real,allocatable,dimension(:,:) :: j12
integer,allocatable,dimension(:,:,:) :: pset

print *,"grnd"
read *, jkl
allocate(ja(3,jkl))
allocate(jb1(3,jkl))
allocate(jb2(3,jkl))
allocate(j12(2,jkl))
jimenmate=jkl

print *,"circle?"
read *, jkl
allocate(kyu(5,jkl))
mate=jkl

print *,"アンチエイリアスx?"
read *, jkl
allocate(pset(3,jkl*640,480*jkl))

print *, "kaiso"
read *, kaiso

yoko=640*jkl
tate=480*jkl
t=-3.0
k=1
w=0.0
jjk=0
do i=1,mate
kyu(4,i)=0.5
kyu(1,i)=w
kyu(2,i)=4.9*mod(i,2)+0.27
kyu(3,i)=t
kyu(5,i)=kyu(4,i)**2
if (mod(i,2).eq.1) then
w=w+0.6
jjk=jjk+1
if (jjk.eq.k) then
jjk=0
w=-0.30*k
k=k+1
t=t+0.99
endif
endif

enddo
w=0.0
t=0.0
jjk=0
k=0


do i=1,jimenmate
ja(1,i)=0.01*rand()*2001-10.0
ja(2,i)=0.01*rand()*601+3.5
ja(3,i)=0.01*rand()*1500+9.25
jb1(1,i)=0.01*rand()*1000-4.95
jb1(2,i)=0.01*rand()*1000-4.95
jb1(3,i)=0.01*rand()*1000-4.95
jb2(1,i)=0.01*rand()*1000-4.95
jb2(2,i)=0.01*rand()*1000-4.95
jb2(3,i)=0.01*rand()*1000-4.95
j12(1,i)=0.7+0.01*rand()*130
j12(2,i)=0.7+0.01*rand()*130
w=sqrt(jb1(1,i)**2.0+jb1(2,i)**2.0+jb1(3,i)**2.0)
jb1(1,i)=jb1(1,i)/w
jb1(2,i)=jb1(2,i)/w
jb1(3,i)=jb1(3,i)/w
w=sqrt(jb2(1,i)**2.0+jb2(2,i)**2.0+jb2(3,i)**2.0)
jb2(1,i)=jb2(1,i)/w
jb2(2,i)=jb2(2,i)/w
jb2(3,i)=jb2(3,i)/w
enddo

light0=1.0*rand()*100-30.0
light1=1.0+rand()*50+rand()*(3.0+rand()*180)
light2=1.0*rand()*100-50.0
w=sqrt(light0**2+light1**2+light2**2)
light0=light0/w
light1=light1/w
light2=light2/w


do k=1,yoko
mcnt1=k
do ii=1,tate
ca0=0.0
ca1=3.0
ca2=-7.0
cb0=-0.004*yoko/jkl+0.008*mcnt1/jkl
cb1=0.004*tate/jkl-0.008*ii/jkl
cb2=9.3
w=sqrt(cb0*cb0+cb1*cb1+cb2*cb2)
cb0=cb0/w
cb1=cb1/w
cb2=cb2/w


hsa0=0.0
hsa1=0.0
hsa2=0.0
hsb0=0.0
hsb1=0.0
hsb2=0.0
mhtmate=-1
vmi=1.0
m=0
do hhpi=1,kaiso
if (m.eq.0) then
htmate=-1
t=99999990.0
do i=1,mate
if (mhtmate.ne.i) then


q0=ca0-kyu(1,i)
q1=ca1-kyu(2,i)
q2=ca2-kyu(3,i)
qt2=q0*q0+q1*q1+q2*q2
b=q0*cb0+q1*cb1+q2*cb2
c=qt2-kyu(5,i)
d4=b*b-c
if (d4.ge.0.0) then
kai=-sqrt(d4)-b
endif
if ((d4.gt.0.0) .and. (kai.gt.0.0) .and. (t.gt.kai)) then
t=kai
htmate=i
endif


endif
enddo
do i=1,jimenmate
if (mhtmate.ne.(i+10000)) then


q0=ca0-ja(1,i)
q1=ca1-ja(2,i)
q2=ca2-ja(3,i)
g(1,1)=cb0
g(2,1)=jb1(1,i)
g(3,1)=jb2(1,i)
g(1,2)=cb1
g(2,2)=jb1(2,i)
g(3,2)=jb2(2,i)
g(1,3)=cb2
g(2,3)=jb1(3,i)
g(3,3)=jb2(3,i)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga.eq.0.0) then
else
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga

if ((kai.lt.j12(1,i)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,i)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
kai=-kai/drtga
if ((kai.gt.0.0).and.(t.gt.kai)) then
t=kai
htmate=i+10000
endif
endif
endif
endif


endif
enddo
if (cb1.lt.0.0) then


x=cb0*ca1/abs(cb1)
z=cb2*ca1/abs(cb1)
kai=sqrt(x*x+ca1*ca1+z*z)
x=x+ca0
z=z+ca2

if (kai.lt.t) then
if (mod(abs(int(x*2.0)+int(z*2.0)),2) .eq. 1) then
vmr=10.0
vmg=10.0
vmb=93.0
else
vmr=155.0
vmg=123.0
vmb=190.0
endif
w=abs(x)+abs(z)
if (w.lt.400.0) then
vmr=vmr*(400.0-w)/400.0+3.0
vmb=vmb*(400.0-w)/400.0+4.0
vmg=vmg*(400.0-w)/400.0+11.0
else
vmr=3.0
vmb=4.0
vmg=11.0
endif
ca0=x
ca1=0.0
ca2=z
x=1.0
do kk=1,mate
q0=ca0-kyu(1,kk)
q1=-kyu(2,kk)
q2=ca2-kyu(3,kk)
qt2=q0*q0+q1*q1+q2*q2
b=q0*light0+q1*light1+q2*light2
c=qt2-kyu(5,kk)
if ((b*b-c).gt.(0.0)) then
x=0.0
endif
enddo

do kk=1,jimenmate
q0=ca0-ja(1,kk)
q1=-ja(2,kk)
q2=ca2-ja(3,kk)
g(1,1)=light0
g(2,1)=jb1(1,kk)
g(3,1)=jb2(1,kk)
g(1,2)=light1
g(2,2)=jb1(2,kk)
g(3,2)=jb2(2,kk)
g(1,3)=light2
g(2,3)=jb1(3,kk)
g(3,3)=jb2(3,kk)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga .ne. 0.0) then
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga
if ((kai.lt.j12(1,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
if ((kai/drtga) .lt. 0.0) then
x=0.0
endif
endif
endif
endif
enddo


if (x.eq.1.0) then
naiseki=cb0*light0-cb1*light1+cb2*light2
vmr=vmr+(1.3+naiseki)*55.0
if (naiseki.gt.0.95) then
vmr=vmr+(naiseki-0.95)*1730.0
endif
vmg=vmg+(1.3+naiseki)*35.0
if (naiseki.gt.0.95) then
vmg=vmg+(naiseki-0.95)*900.0
endif
vmb=vmb+(1.3+naiseki)*17.16
if (naiseki.gt.0.95) then
vmb=vmb+(naiseki-0.95)*450.0
endif
else
vmr=vmr/1.7
vmg=vmr/1.7
vmb=vmr/1.7
endif
htmate=-1
endif


endif
mhtmate=htmate
if (htmate.eq.-1) then
m=1
else


if (htmate.lt.10000) then
hsb0=(ca0-kyu(1,htmate)+t*cb0)/kyu(4,htmate)
hsb1=(ca1-kyu(2,htmate)+t*cb1)/kyu(4,htmate)
hsb2=(ca2-kyu(3,htmate)+t*cb2)/kyu(4,htmate)
hsa0=kyu(1,htmate)+hsb0*kyu(4,htmate)
hsa1=kyu(2,htmate)+hsb1*kyu(4,htmate)
hsa2=kyu(3,htmate)+hsb2*kyu(4,htmate)
endif

if (htmate.ge.10000) then
htmate=htmate-10000
hsa0=ca0+t*cb0
hsa1=ca1+t*cb1
hsa2=ca2+t*cb2
hsb0=jb1(2,htmate)*jb2(3,htmate)-jb1(3,htmate)*jb2(2,htmate)
hsb1=jb1(3,htmate)*jb2(1,htmate)-jb1(1,htmate)*jb2(3,htmate)
hsb2=jb1(1,htmate)*jb2(2,htmate)-jb1(2,htmate)*jb2(1,htmate)
w=sqrt(hsb0*hsb0+hsb1*hsb1+hsb2*hsb2)
hsb0=hsb0/w
hsb1=hsb1/w
hsb2=hsb2/w
htmate=htmate+10000
endif


naiseki=-(hsb0*cb0+hsb1*cb1+hsb2*cb2)
ca0=hsa0
ca1=hsa1
ca2=hsa2
cb0=2.0*naiseki*hsb0+cb0
cb1=2.0*naiseki*hsb1+cb1
cb2=2.0*naiseki*hsb2+cb2
vmi=vmi+0.7

endif
endif

enddo


if (cb1.ge.0.0) then

naiseki=(cb0*light0+cb1*light1+cb2*light2)
vmr=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmr=vmr+(naiseki-0.96)*8200.0
endif
vmg=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmg=vmg+(naiseki-0.96)*8200.0
endif
vmb=(1.5+naiseki)*190.0
if (naiseki.gt.0.96) then
vmb=vmb+(naiseki-0.96)*8200.0
endif

endif


ll=int(vmr/vmi)
if (ll.lt.256) then
pset(3,mcnt1,ii)=ll
else
pset(3,mcnt1,ii)=255
endif

ll=int(vmg/vmi)
if (ll.lt.256) then
pset(2,mcnt1,ii)=ll
else
pset(2,mcnt1,ii)=255
endif

ll=int(vmb/vmi)
if (ll.lt.256) then
pset(1,mcnt1,ii)=ll
else
pset(1,mcnt1,ii)=255
endif


vmr=0.0
vmg=0.0
vmb=0.0
enddo
enddo



open(11,file = 'file1.txt')
open(12,file = 'file2.txt')
open(13,file = 'file3.txt')
open(14,file = 'file4.txt')
i=0
h=0
k=tate/jkl+1
m=0
ii=0
kk=1
ll=0
do i=1,yoko
do h=1,tate
do m=1,3
mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)=mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)+pset(m,i,h)
enddo
enddo
enddo
m=0
h=0
i=0


hhpi=yoko*tate/jkl/jkl*3/4
do i=1,yoko*tate/jkl/jkl*3
h=mod(i,3)+1
m=mod(i,4)+1
if (h.eq.1) then
ii=mod(ii,yoko/jkl)+1
if (ii.eq.1) then
k=k-1
endif
endif
ll=ll+(mtset(h,ii,k)/jkl/jkl)*kk
kk=kk*256
if (m.eq.4) then
write(11+(i-1)/hhpi,*) ll
ll=0
kk=1
endif
enddo
close(11)
close(12)
close(13)
close(14)
print *,"end111111111111111111111"
stop
end





私は実行するときコマンドプロンプトから実行しているが、多分これを実行すると最初に

grnd?

と入力を求められる

ここに「四角形の板」の枚数を入力する。鏡みたいなやつだ。当然反射する。

上の生成画像では確か10って設定だったかな?


次に

circle?

と入力を求められる

ここに「」の個数を入力する

上の画像では確か15000くらい。


アンチエイリアス?

にはアンチエイリアスなしで1を、2,3,4となるごとにきめ細かくなっていくが計算時間が2乗倍になる

上の画像では確か4



kaiso?

は、最高反射計算回数。まぁ4以上の数字ならたいてい変な画像はできない

上の画像は確か19とか



これでenterキーを押せば計算開始、計算が終了する頃に「end111111111111111111111」っていう計算終了のお知らせメッセージが表示される

また出力ファイルで4つのテキストファイルが出力される


そこに全ピクセルのRGB情報が詰まっている




今度はHSPでスクリーン出力

そのソースは↓










onexit*enf
repeat 4
ww=3-cnt
exist ""+ww+""
if strsize=-1:IIIIIIIIIIDDDDDDD=ww:bsave ""+ww+"",IIIIIIIIIIDDDDDDD:break
if cnt=3:end
loop

screen 0,640,480:boxf
dim i,640*480*3/4
mref i,66
sdim k,16

h=640*480*3/4/4*IIIIIIIIIIDDDDDDD
repeat 1,IIIIIIIIIIDDDDDDD
sdim a,1
notesel a
noteload "file"+(cnt+1)+".txt"
repeat 640*480*3/4/4
noteget k,cnt
i.h=int(k)
h+
if cnt\900=0:await 1:redraw
loop
loop

bmpsave ""+IIIIIIIIIIDDDDDDD+".bmp"
delete ""+IIIIIIIIIIDDDDDDD+""
if IIIIIIIIIIDDDDDDD=1{
repeat 999999
w=0
repeat 4
exist ""+cnt+""
if strsize!-1:w+
loop
if w=0:break
wait 200
loop
screen 0,640,480:boxf
repeat 4
buffer 1+cnt,640,480:picload ""+cnt+".bmp"
loop
gsel 0
repeat 4
pos 0,0
gmode 5,640,480,256
pos 0,0
gcopy 1+cnt,0,0,640,480
loop
bmpsave "5.bmp"
delete "0.bmp"
delete "1.bmp"
delete "2.bmp"
delete "3.bmp"
stop
}else{
end
}
*enf
delete ""+IIIIIIIIIIDDDDDDD+""
end



相変わらず超汚い

こんなの誰も理解できんだろうな・・そのうち自分も分からなくなるよこれ



実はこのプログラムDualコア4スレッド対応プログラム

でもHSPはシングルスレッドなのでやっぱ4回起動する必要があるw


1スレッドで普通にプログラム実行すると1分くらいかかるので折角のcore i5なので「並列プログラム」を組んでみたというわけ


注意点はなるべく素早く4回起動しないとバグる!(ぇ



しかし別に原理的にはただ単に4つ起動させているだけなので1コアのCPUでも全然動くw


で、4つ起動させて全部の計算が終わると残り3つのプログラムは自動終了して1つのスクリーンに完成画像が表示されて「5.bmp」ってファイルに出力される。


これでやっと上の画像が完成

おしまい






っと、今のところfortranとHSPの組み合わせ研究はこんな感じで進んでいる。


今度はfortranでムービーの作成でもやってみようかと考えている!


でRGBデータ→ムービーの変換のところでHSPを使うって感じで

レイトレーシングアルゴリズム実験②

レイトレ実験はその②ですが、レイトレの記事はもうこれで5つ目です・・なんかタイトルのつけ方ミスった??



実は新しいプログラム言語 「fortran」 に手を出してみました!


これは計算に特化した言語でよく科学技術用の計算で用いられる、いわば四則演算がすごく速く行える言語です!


HSPでレイトレーシングをやるとすごく時間がかかるのでfortranに移植してみました。

当然「できあがる画像」はまったく一緒です!

同じアルゴリズムですから・・・



fortranで作ったCG↓

b114ac31.png



 

こんなにたくさんのオブジェクト(球)があるとHSPでは何時間計算かかるか分かりません


がfortranでは20分くらいで作れました。


なんだかよくわからないグラフィックになってしまいましたが、球が規則正しく敷き詰められているだけです。

下と上に1層ずつ、って



見ててきづいたのですが、画面の球にはさまれた部分(空と地面)の形は、本来なら長方形になるはずが、中心の部分がふっくらして見えません??

きっと目の錯覚と言うやつですね・・・・


どうでもいいです・・・・




ところでfortranに手を出したわけは・・

これからHSPをメインに使わなくなると言うわけではなく、あくまでレイトレーシングの計算をHSPで完成させたあと、fortranで動かすというカンジでやっていくためです。


やはりHSPはデバックが簡単だしやスクリプト画面もみやすいので。

それに一番は、もうそれに慣れてしまったところが大きいですかね




で一番自分が気になる速度比較ですが、計算速度はなんとHSPの30倍くらい早かったです!!


こんど結果をグラフにしてみたいと思います。




しかしfortranで計算が早くなったのはいいが画面出力のやり方がわからないので、計算結果→画面出力をほとんどHSPでやる羽目に・・


まずfortranでレイトレ演算結果の数値(RGB)をtxtに出力させ、それをHSPで読み込み、読み込んだ文字列を数値に変換しpsetで画面に打っていくという、なんともややこしく回りくどい方法を取らざるを得なくなり、結局その変換にすごい計算時間がかかるというオチであった・・orz



そんなこんなで4、5分くらいの計算ならHSPで全部やったほうが結果的に早いことが判明


もしかしたらfortranだけで画面出力できる方法があるのかもしれないがあまり書籍もないし、自分の能力も限界に近いので、数値→画像変換はHSPに頼ることにした。


やっぱわたしゃHSPから離れられんなぁ~




ところで↑のCG見て思った
ドットが荒い!

640×480の計算結果を640×480の画面に出力させてるから、1/1でドットが目立ち細かいところは何がなんだか分からなくなっている。


なので次に作る画像はアンチエイリアシングをかけてみようと思った。



擬似的なアンチエイリアシングの技法として目的の2~3倍程度の大きい画像を最初に作り、それを縮小してドットを滑らかにする、というものがある。

2倍の例だと

1280×960の計算結果を640×480の画面に出力させる、など



これは倍率を上げれば上げるほどドットが滑らかになるがその2乗に比例して演算時間がかかってしまう。



多分だが、本当のアンチエイリアシングは浮動小数点かなんか使って↑とは全く別のアルゴリズムでドットを滑らかにしているんだろうが

ま、レイトレーシングでアンチエイリアスするアルゴリズムなんて少し考えただけでも恐ろしく難しいので擬似的な方のアンチエイリアスでやることにした。




さて前回の記事で、四角板もレイトレーシングで描画できるのではと予言しましたが、やってみたら正解!

板の全反射なので「鏡」みたいな役割をするんじゃないかと思ってたら、その通りだった




以下のCGは3倍の倍率でやった擬似アンチエイリアシングの、四角形板&球混合バージョン

(これは最初からHSPで計算)

c745ce0e.png


うぉぉ・・

これには感動した!!(泣)


四角形の板は、球より難しいアルゴリズムであったがちゃんとピカピカ光ってる!


アンチエイリアシングも思いのほか効きが良くてうれしいー!


うんしばらくこれをデスクトップの壁紙にしよう!





四角形の板のアルゴリズムは次で解説します。

あと途中から語尾が変わってたのはまぁ文才がないってことで・・m(_ _ )m

レイトレーシングアルゴリズム実験①

まずレイトレの四角形の当たり判定を行なう前によくある球の当たり判定をおさらいしてみます



絶対必要なのが視線方向ベクトルと視線位置ベクトル


方向ベクトルは視線がどこの方向を向いているか、x、y、z、成分の3つの値で設定します

位置ベクトルはx、y、z、の位置をあらわしているベクトルです



視線方向をVベクトル、視線位置をMベクトルとし、判定をする球の中心位置ベクトルをQベクトル、球の単位ベクトルをP、半径をrとします


※VとPは単位ベクトル

4b8b9221.png




あとは数Bでやったハズの交点を求める計算をします。


球と線分の交点はベクトルVに変数tを掛けてtVとし、それが位置Qから見てrの距離に存在できるかを考えます。



上の図から Q+rP=M+tV


という式が立てられます

注意したいのはPベクトルの向きが決定されてないのでx、y、zの成分を3つ比べて連立方程式を解く

という方法は使えません。


なので |P|=1.0 という条件をうまく使って解くことを考えます

P^2 (Pの2乗) は1なので


r^2=(rP)^2=(M+tV-Q)^2  ・・・・・①


という式が立てられます。

MとQはすでにx、y、zの成分が決定しているので新しいベクトルHをつくり


H=M-Q


として①に代入&変形し


(H+tV)^2-r^2=0    ・・・・・・②


と t の2次方程式を作り、解の公式を使って t を求めます。

②を展開し


t^2 + 2tH・V + H^2 - r^2 = 0    ・・・・③


これで、ニーエーブンノ・・・の公式が使えますね


a=1

b=2H・V

c=H^2-r^2


bはベクトル同士の掛け算が入っているので内積を使います

H=(x1,y1,z1)

V=(x2,y2,z2)

のとき

H・V=x1*x2 + y1*y2 + z1*z2

です。


またbが2の倍数なので

b'=b/2  とおき


エーブンノマイナスビーダッシュ・・・の公式を使って


t=( -b'±√b'^2-ac ) / a


ここで t は0、1つまたは2つの解をもつことが判明します。


2つの場合は交点が2つ、すなわち球に視線が刺さっているような状態です。

1つの場合は接していると言えます。

0の場合は交点がありません。



また解が虚数でなくても t は正でなければいけないし

交点が2つの場合は

視線は球の奥じゃなく手前の交点で反射してますから解の公式の±の部分は-であると決定できます。

t は図の通り視線の長さです。


さてここまでが思考上の筋道です。

数学が絡むアルゴリズムを考える場合は、プログラムをいきなり書き始める前にこのように思考上で筋道を作ることが絶対必要です。





では、これらをプログラムでやってみましょう。



計算するに当たり初期状態で必要な変数は

ベクトルV、M、Qのx、y、z成分の値を記憶している変数です。

その変数はこのブログでは仮にVx,Vy,Vz,Mx・・・・・・とします。





Hx=Mx-Qx

Hy=My-Qy

Hz=Mz-Qz

b=Hx*Vx + Hy*Vy + Hz*Vz           ;内積の計算で解の公式のbに出力

c=Hx*Hx + Hy*Hy + Hz*Hz - r*r       ;ベクトルの距離の計算で解の公式のcに出力

D4=b*b-c                      ;解の公式の4分のDに出力

if D4<0.0:gosub*交点なし処理          ;ルートの中がマイナスのとき条件分岐

t=-b-sqrt(d4)

if t<0.0:gosub*交点なし処理           ;交点はあるが視点の後ろにある場合条件分岐

gosub*交点あり処理               ;視線と球が当たった時の処理、交点の位置はM+tV








以上が当たり判定部分の計算プログラムです。

紙の上で計算したときより、いらない計算や式は全部省いているので相当短くっています。


結果的にPベクトルは使わないので最初に設定しておく必要もなくプログラム上では変数を作る必要がありません。



交点がない場合はまた別の球と当たり判定をします

全ての球とぶつからなく、視線が地面or空まで伸びる場合はその色を、画面のドットの色として出力します。




交点がある場合は、入射してくる光線のベクトルを決定するためにまた別の計算を行わなければいけません。

その計算とは法線の計算です。

入射ベクトルを求めるには以下の図のように法線ベクトルを求める必要があります。

f8b0024f.png


法線ベクトルは球の中心から交点までのベクトルを単位ベクトルに変換すれば完成です。

-入射ベクトルを求める計算は図のように


-入射ベク = cos(θ) × 法線ベク × 2 + 視線ベク


ですぐ求まります。


cos(θ)は法線ベクと視線ベクの内積に等しいです。(ただし両方単位ベクトルでないといけない)



そして入射ベクトルが求まれば、視点位置=交点、視線ベクトル=-入射ベクトル、としてカメラの位置と方向を改めて決定し、再度当たり判定処理をすればいつかは地面or空に当たるときが来ますので、それまで何回もループさせます。

(処理が重い場合は反射の回数制限をする場合もある)

レイトレーシング試作CG実験

あれからいろいろプログラムを高速化するべく改善してみようとしましたがまだまだ自分のレイトレーシングの技術の足りなくて四苦八苦…


だけど、HSPでこんな綺麗なCGを作れるとは…と改めてレイトレーシングの凄さに痛感しました!


自分のイメージ通りに作れると感動しますよ!



自分で絵が描けなくても、または絵が超絶に下手でもコンピューターに描いてもらえばいいのです!
これぞプログラマーの特権!


なんとプログラム技術を伸ばせば絵の技術も同時にあがると!(ドット絵は含まれないが)



ところで何回かプログラムいじってて気がついたのですが、なぜレイトレーシングと言うと球と地面だけなのか…


グーグルの画像検索でレイトレーシングと検索するとほとんどが球の反射



でも、球しかない世界はつまらなすぎる…


レイトレの描画アルゴリズムを簡単に言うと、視点に入る光線を逆探知して物体との当り判定をし当たってたら反射する前の光線を逆探知…また当り判定に戻る…を繰り返します


物体との当たり判定にはベクトルを使って計算している訳だが、これはつまりベクトルで表現できる図形ならなんでも大丈夫ってことなんじゃ?
私の憶測ですけど


ベクトルで表現できる図形と言えばもうなんでもあり、ただの平面や歪んだ面、当然球も、面を組み合わせれば立体も、そして円柱を曲げて作ったような立体のリングなど…
ただよくあるレイトレでわざわざ球にする意味は確かにあり、有限のポリゴンで表現できない形だしレイトレの特徴を一番良く表現できてる図形だからなのでしょう


が、理論上平面やリングとかもレンダリング可能なはずだと私は思います!


そうと決まれば、早速実験!

平面のベクトル方程式

P=a+sB+tC (大文字がベクトル)

で平面のレンダリングをやってみようと思います!

プロフィール

toropippi

記事検索
アクセスカウンター

    QRコード
    QRコード
    • ライブドアブログ