プチコン3号 32日目 球を描く(2)

pc32-1.jpg

今回は、前回の続きで、上の画面のような前回よりもリアルな球を描いてみます。

前回は球を正面から光で照らした場合の陰影を描画しました。
しかし、球を正面から光で照らすシーンは、日常ではあまり出会いません。
もっとも一般的なのは、満月でしょうか・・・

今回は日常的に目にする、斜め上からの光で照らされている球を描画します。
また、光が物体の表面で鏡面反射されて生じる、ハイライトも描き込みます。

前回説明した「光の当たる角度」と「環境光」、さらに今回の「ハイライト」という3つの成分の合成でシェーディングを表現する手法は、「フォンの反射モデル」と呼ばれる古典的な手法です。
(ちなみに「フォン」は発明者の名前です。)

今回のプログラムでは、光のあたる角度は、球の正面から照らす場合を回転角0、球の真上から照らす場合を回転角π/2(90度)とし、X軸の周りの回転角だけで表します。左右へは傾けません。
これによって、画像が左右対称になり、全ピクセルの半分の計算で描画を終えることができます。
(左右に傾けても、対称軸が傾くだけで対称にはなりますが、ピクセルの位置を計算しづらくなります。)

光を斜め上方向に移動するには、25日目で解説した回転の座標変換を使います。

下図では、図左は前回のプログラムのように上(Z軸方向)から光が当たっています。
図右は、光が角度θだけ傾いています。
このとき、球面で一番明るくなるのは、真上からθだけ傾いた点Bになります。
では点Aの明るさはどうなるかというと、図左で真上から-θだけ傾いた点Cの明るさと同じになります。

pc32-3.jpg

以上は、球の中心の周りに、球も光も含めた全体が角度θだけ回転したと考えれば、直観的にも理解できると思います。

また、ハイライトはその強さのみを、0から1の範囲で与えることにします。

ハイライトの出現位置は、球の最も明るい場所ではなく、光の角度の半分の角度の位置に出ます。
たとえば

・光の方向が0度のときはハイライトの位置も0度(円の中心)
・光の方向が90度(上方)のときは、ハイライトの位置は45度(下図参照)

となります。

pc32-4.jpg

ハイライトの光の分布にもいろいろなモデルがありますが、要は鋭いピークがありつつ滑らかに変化していればもっともらしく見えます。
今回は「半分の角度で光を当てたときの明るさ」を500乗して作っています(Phong分布相当)。
また、反射光の色は物体の色ではなく、白色としました。

次の図は、今回のプログラムで角度とハイライトのパラメータをいろいろ変えてみた画像です。
上段は光の角度を左から0度(正面)、30度、60度、90度(真上)の4通りに変えています。
下段は、角度は60度のまま、ハイライトの係数を0.1、0.3、0.5、0.7と変えています。

pc32-2.jpg

プログラムは以下の通りです。
傾いた光に関する部分を紫、ハイライトに関する部分をオレンジにしています。

CLS:GCLS
T=MAINCNT
DRAWSPHERE 100,20,170,192,192,128,0.2,PI()/4,0.4
PRINT (MAINCNT-T)/60;" SECONDS ELAPSED"

DEF DRAWSPHERE X0,Y0,SIZE,R0,G0,B0,GL,T,SP
VAR R,RR,X,XX,Y,YY,Z,ZZ,YMAX,NZ,NY,C,R1,G1,B1
R=SIZE/2+0.5
RR=R*R
X0=X0+R
Y0=Y0+R
X=-0.1
WHILE R>=X
Y=-0.1
XX=X*X
YY=Y*Y
YMAX=RR-XX
ZZ=YMAX-YY
WHILE ZZ > 0
Z=SQR(ZZ)
ROTATE Y/R,Z/R,T OUT NY,NZ
NZ=MAX(0,NZ)
C=SHADING(NZ,1-GL,GL)
RGB_ADD 0,0,0,C,R0,G0,B0 OUT R1,G1,B1

ROTATE Y/R,Z/R,T/2 OUT NY,NZ
NZ=MAX(0,NZ)
C=SPECULAR(NZ,SP,0.002)
RGB_ADD R1,G1,B1,C,255,255,255 OUT R1,G1,B1
C=RGB2(255,R1,G1,B1)
GPSET X0-X,Y0-Y,C
GPSET X0+X,Y0-Y,C

ROTATE -Y/R,Z/R,T OUT NY,NZ
NZ=MAX(0,NZ)
C=SHADING(NZ,1-GL,GL)
RGB_ADD 0,0,0,C,R0,G0,B0 OUT R1,G1,B1

ROTATE -Y/R,Z/R,T/2 OUT NY,NZ
NZ=MAX(0,NZ)
C=SPECULAR(NZ,0.4,0.002)
RGB_ADD R1,G1,B1,C,255,255,255 OUT R1,G1,B1
C=RGB2(255,R1,G1,B1)
GPSET X0+X,Y0+Y,C
GPSET X0-X,Y0+Y,C
Y=Y+1
YY=Y*Y
ZZ=YMAX-YY
WEND
X=X+1
WEND
END

DEF SHADING(VAL,DIF,AMB)
RETURN AMB+DIF*VAL
END

DEF SPECULAR(VAL,SP_I,SP_SIZE)
RETURN SP_I*POW(VAL,1/SP_SIZE)
END

DEF RGB_ADD R0,G0,B0,K,R1,G1,B1 OUT R2,G2,B2
R2=MIN(255,R0+K*R1)
G2=MIN(255,G0+K*G1)
B2=MIN(255,B0+K*B1)
END

DEF RGB2(A,R,G,B)
IF RND(8) < (R AND 7) THEN R=R+8
IF RND(8) < (G AND 7) THEN G=G+8
IF RND(8) < (B AND 7) THEN B=B+8
RETURN RGB(A,R,G,B)
END

DEF ROTATE X,Y,T OUT X1,Y1
X1=COS(T)*X - SIN(T)*Y
Y1=SIN(T)*X + COS(T)*Y
END

回転角は変数Tに入っています。
回転の処理(手続きROTATE)は25日目に使ったものと同じです。

ハイライトの計算については、まず光がT/2の角度で入ってきた場合のピクセルの反射の強さを求めます。
そして、反射の強さを関数SPECULARでハイライトに変換します。
ハイライトは強さのパラメータが変数SPに入っています。
SPECULARは、反射の強さをVAL、ハイライトの強さをSP_I、ハイライトの広がり具合をSP_SIZEで与えます。
最後のSP_SIZEは小さくするほどハイライトも小さくなります。
ただし、実際にハイライト部分のピクセル数などを指定しているわけではなく、VALを1/SP_SIZE乗するだけです。
値としては、0.01以下が良いでしょう。
今回は0.002にしています。

画像の明暗は左右対称になりますが、形状自体は上下にも対象です。
そこで、ループの1回ごとに、上半分のピクセルの明るさと下半分のピクセルの明るさの計算を行っています。
ほぼ同じコードがループの前半と後半にありますが、前半が上半分、後半が下半分に関する計算です。

このプログラムで前回と同じ大きさの球(直径170ピクセル)を描くのに3.5秒かかります。
前回は全体の1/8のピクセルの明るさだけを計算すれば済みましたが、今回は全体の半分の計算が必要で、かつ計算自体も増えているためです。

プチコン3号 31日目 球を描く(1)

pc31-1.jpg

前回紹介したレイトレーシングの手法は、時間がかかりすぎるのでゲームの画像表示に使うには向いていません。
今回は、球を一つだけ、レイトレーシングよりも高速に描画するプログラムを作ってみます。
上の画面が実行結果で、直径170ピクセルの球を0.6秒弱で描画しています。
アニメーションに使うのは難しいですが、スプライトやBGの作成には使うことができるでしょう。

球を球らしく見せるには、陰影が重要です。
光の方向や強さなどに応じてピクセルの明るさを決めることを「シェーディング」と言います。

下図は球面が上方から来る光に照らされている様子を示しています。
(以下、球面の光の当たらない下半分は無視します。)
地球が太陽に照らされている様子を思い浮かべれば分かりやすいでしょう。
このときシェーディングは基本的には以下のようにして行えます。

まず、球面のごく小さな一部を考えます。
これは、ある傾きを持った平面にほぼ等しいと考えられます。
この小さな平面の明るさは、光が来る方向(図では垂直方向)に対して平面が直角(図では水平)になっているとき、最大になります。
平面が傾くにつれて暗くなっていき、面が垂直になったときに明るさはゼロになります。
この明るさは、ある長さの棒を平面に対して垂直に立てるとき、その棒の両端の座標のZ値の差分(図の赤い矢印)に比例します。

pc31-2.jpg

プログラムにするために、上記を数式に直してみます。
図では上から光が来ていますが、この上方向をZ軸とします。
そして、球面の中心を原点とします。
球面上の点は、半径をrとすると

 x2+y2+z2=r2

を満たします。

このとき、シェーディングに必要な値である赤い矢印の長さは、実は球の半径とzの値だけで決まります。
直感的にはちょっとわかりにくいかもしれませんが、zが最大となるのは図の球面の一番高い部分(中央部分)で、このときx=y=0でz=rとなり、図の球面の両端部分ではz=0となるのはすぐわかると思います。

球面に垂直に棒を立てると、その棒は必ず球の中心を指します。
従って

・棒を斜辺とし、水平および垂直な辺を持つ直角三角形
・球面の中心から棒を立てた点までの線分を斜辺とし、水平および垂直な辺を持つ直角三角形

は相似の関係になります。
すると、

・棒の長さと、赤い矢印の長さの比
・球面の半径と、棒を立てた点のZ座標値の比

が等しいことが言えるので、棒の長さが1のとき、赤い矢印の長さは(棒を立てた点のZ座標値/球面の半径)となります。

以上を用いて、球を光源側から見た場合(図では上方向から見た場合)の画像を描いてみます。
先に書いたとおり、球面の方程式は

 x2+y2+z2=r2

です。
従って半径rの球面を描くには、x,yを変化させて、それに対応するzの値を

 z=SQR(r2-x2-y2)

で求め、zの値に応じて(x,y)のピクセルの明るさを変化させます。

プログラムは以下のようになります。

CLS:GCLS
T=MAINCNT
DRAWSPHERE 100,20,170,255,255,176,0.2
PRINT (MAINCNT-T)/60;" SECONDS ELAPSED"

DEF DRAWSPHERE X0,Y0,SIZE,R0,G0,B0,GL
VAR R,RR,X,XX,Y,YY,Z,ZZ,YMAX,NZ,NY,C,R1,G1,B1
R=SIZE/2+0.5
RR=R*R
X0=X0+R
Y0=Y0+R
X=-0.08
WHILE R>=X
Y=X
XX=X*X
YY=XX
YMAX=RR-XX
ZZ=YMAX-YY
WHILE ZZ > 0
Z=SQR(ZZ)
NZ=Z/R
C=SHADING(NZ,1-GL,GL)
RGB_ADD 0,0,0,C,R0,G0,B0 OUT R1,G1,B1
C=RGB2(255,R1,G1,B1)
GPSET X0+X,Y0+Y,C
  GPSET X0-X,Y0+Y,C
GPSET X0+X,Y0-Y,C
GPSET X0-X,Y0-Y,C
  GPSET X0+Y,Y0+X,C
  GPSET X0-Y,Y0+X,C
GPSET X0+Y,Y0-X,C
GPSET X0-Y,Y0-X,C
  Y=Y+1
YY=Y*Y
ZZ=YMAX-YY
WEND
X=X+1
WEND
END

DEF SHADING(VAL,DIF,AMB)
RETURN AMB+DIF*VAL
END

DEF RGB_ADD R0,G0,B0,K,R1,G1,B1 OUT R2,G2,B2
R2=MIN(255,R0+K*R1)
G2=MIN(255,G0+K*G1)
B2=MIN(255,B0+K*B1)
END

DEF RGB2(A,R,G,B)
IF RND(8) < (R AND 7) THEN R=R+8
IF RND(8) < (G AND 7) THEN G=G+8
IF RND(8) < (B AND 7) THEN B=B+8
RETURN RGB(A,R,G,B)
END

手続きDRAWSPHEREが描画処理の本体です。
パラメータは

・左端のX座標
・上端のY座標
・直径
・色(R,G,B)
・環境光係数

となっています。

球のサイズは半径ではなく直径で指定しています。
球は一辺が直径に等しい正方形の領域に描画されます。

最後のパラメータ「環境光係数」には、通常0.1~0.2程度を指定します。
このパラメータは球全体の明るさを底上げします。
1.0を指定すると、陰影のない塗りつぶした円が描画されます。

DRAWSPHEREの中身は2重のループになっています。
外側のループがX座標を1ずつ増やし、内側のループがY座標を1ずつ増やします。
(X0,Y0)が中心の座標です。

内側のループの中で、GPSETを8回連続して行っています。
今回描画する球では、明るさの値は原点を中心とした同心円上で同じ値になります。
従って下図のように、座標を入れ替えたり符号を反転することによって、一度の計算で8か所のピクセルの値が決まります。
全体の1/8を計算すれば十分なので、ループは図の網掛けした部分だけをカバーしています。

pc31-3.jpg

RGB_ADD手続きは、その名の通り色の加算(と定数倍)を行う手続きです。
色(R0,G0,B0)に色(R1,G1,B1)をK倍したものを加算し、結果を(R2,G2,B2)に代入します。

RGB2関数は前回も説明しましたが、疑似的に中間色を生成して色の変化を滑らかに見せます。

なお、ピクセルが球の内部か外部かを判定する際、境界値付近での判定を意図的に制御するために、変数の初期値等を小数点以下のオーダーで幾何学的な値よりも大きくしたり小さくしたりしています。
たとえば

・半径Rが直径の1/2よりも0.5大きい
・Xの初期値が0ではなく-0.08

などです。
いずれも幾何学的には、0.5や-0.1は不要ですが、特に径が小さい円を描くとき、形状に微妙な違いが出てきます。
掲載したプログラムは、直径が16ピクセルや14ピクセルの形状が私の好みのもの(下図)になるように調整しています。
(その代わり、直径15ピクセルを指定しても16ピクセルになってしまいますが・・・)

PC31-4.png

今回のプログラムは、球を正面から光で照らした場合の陰影を計算しています。
ちょうど、満月のようなものです。

次回はこのプログラムを少し拡張して、前回のレイトレーシングの例で示したような、斜め方向から照らした場合の陰影を計算してみます。

プチコン3号 30日目 レイトレーシングでCGを描いてみる

pc30-1.jpg

今回は、コンピュータグラフィックスの古典的な手法であるレイトレーシング法のプログラムを紹介します。
これは、下記のリンク先にある、”Tiny Raytracer”というJavaScriptで書かれたプログラムをプチコン3号用に書き直したものです。

Gabriel Gambetta – Tiny Raytracer

レイトレーシングの原理は、物体から視点に届く光を計算するために、視点への光の通り道を逆算していくというものです。
そして、通り道の先に物体があったら、その物体に光が反射していると考えて、反射の先の通り道をさらに辿って行きます。
下の図はWikipediaの引用ですが、イメージがわきやすいと思います。

Ray trace diagram.svg
Ray trace diagram” by HenrikOwn work. Licensed under GFDL via Wikimedia Commons.

プログラムは140行あまりと短いので、まずは全体を掲載します。
なお今回は移植ですので、Miiverseへのアップロードは行いません。

OPTION STRICT

DEF READ_ARRAY A,SZ,LBL
VAR I
RESTORE LBL
FOR I=0 TO SZ
READ A[I]
NEXT
END

DEF DOT(A,B)
RETURN A[0]*B[0]+A[1]*B[1]+A[2]*B[2]
END

DEF A_MINUS_BK A,B,K,C
C[0]=A[0]-B[0]*K
C[1]=A[1]-B[1]*K
C[2]=A[2]-B[2]*K
END

VAR W=240 'IMAGE SIZE
VAR T,C

DIM SPHERES[36]
VAR SPECULAR=6
VAR REFLECT=7
@SPHERES
'IF YOU EDIT W VALUE, CHANGE "240" TO NEW W
DATA 240, 0,-240,0, 9,9,0, 240,2
DATA   1, 0,   0,3, 9,0,0, 240,3
DATA   1,-2,   1,4, 0,9,0,   9,4
DATA   1, 2,   1,4, 0,0,9, 240,5
READ_ARRAY SPHERES,35,"@SPHERES"

VAR AMBIENT_LIGHT=2

DIM LIGHTS[4]
@LIGHTS
DATA 8,2,2,0
READ_ARRAY LIGHTS,3,"@LIGHTS"

DEF CLOSEST_INTERSECTION(B,D,TMIN,TMAX)
 VAR _A,_V,_Q,_B,_D,_R,_F,_9Q
DIM J[3], SC[3]'SPHERE CENTER
T=W
_A=2*DOT(D,D)
_V=0
FOR _Q=0 TO 3 'NUM OF SPHERES
_9Q=_Q*9
_R   =SPHERES[_9Q]
SC[0]=SPHERES[_9Q+1]
SC[1]=SPHERES[_9Q+2]
SC[2]=SPHERES[_9Q+3]
A_MINUS_BK B,SC,1,J
_B=-2*DOT(J,D)
_D=_B*_B - 2*_A*(DOT(J,J)-_R*_R)
IF _D > 0 THEN
_D=SQR(_D)

_F=(_B-_D)/_A
IF TMIN<_F && _F<TMAX && _F<T THEN
_V=_9Q+1
T=_F
ENDIF

_F=(_B+_D)/_A
IF TMIN<_F && _F<TMAX && _F<T THEN
_V=_9Q+1
T=_F
ENDIF
ENDIF
NEXT
RETURN _V
END

DEF TRACE_RAY(B,D,TMIN,TMAX,DEPTH)
 VAR I,_S,_N,_I,_L,_U,_K,_4I
DIM N[3],X[3],SC[3],L[3],LC[3],M[3]
_S=CLOSEST_INTERSECTION(B,D,TMIN,TMAX)
IF _S == 0 THEN RETURN 0
SC[0]=SPHERES[_S]
SC[1]=SPHERES[_S+1]
SC[2]=SPHERES[_S+2]

A_MINUS_BK B,D,-T,X
A_MINUS_BK X,SC,1,N
_N=DOT(N,N)
_I=AMBIENT_LIGHT
FOR I=0 TO 0 'NUM OF LIGHTS-1
_4I=I*4
_U   =LIGHTS[_4I]
LC[0]=LIGHTS[_4I+1]
LC[1]=LIGHTS[_4I+2]
LC[2]=LIGHTS[_4I+3]
A_MINUS_BK LC,X,1,L
_K=DOT(N,L)
A_MINUS_BK L,N,2*_K/_N,M
IF CLOSEST_INTERSECTION(X,L,1/W,1)== 0 THEN
_I = _I + _U * ( MAX(0, _K/SQR(DOT(L,L)*_N)) + MAX(0,
POW(DOT(M,D)/SQR(DOT(M,M)*DOT(D,D)),SPHERES[_S+SPECULAR])))
ENDIF
NEXT

VAR LOCAL_COLOR=SPHERES[_S+3+C]*_I*2.8
VAR REF=SPHERES[_S+REFLECT]/9
VAR P[3]
DEPTH=DEPTH-1
IF DEPTH < 0 THEN
RETURN LOCAL_COLOR
ELSE
A_MINUS_BK D,N,2*DOT(N,D)/_N,P
RETURN TRACE_RAY(X,P,1/W,W,DEPTH)*REF + LOCAL_COLOR*(1-REF)
ENDIF
END

GCLS:CLS
VAR X,Y,H
VAR TSTART=MAINCNT
DIM CAM[3],POS[3]
DIM PX[3]
CAM[0]=0:CAM[1]=1:CAM[2]=0
H=W/2
FOR Y=H TO -H STEP -1
FOR X=-H TO H
FOR C=0 TO 2
POS[0]=X/W
POS[1]=Y/W
POS[2]=1
PX[C]=TRACE_RAY(CAM,POS,1,W,2)
NEXT
GPSET 80+X+H,H-Y,RGB2(255,PX[0],PX[1],PX[2])
LOCATE 0,0:?X+H,H-Y
NEXT
NEXT
?(MAINCNT-TSTART)/60;" SECONDS ELAPSED"

DEF RGB2(A,R,G,B)
 IF RND(8) < (R AND 7) THEN R=R+8
IF RND(8) < (G AND 7) THEN G=G+8
IF RND(8) < (B AND 7) THEN B=B+8
RETURN RGB(A,R,G,B)
END

移植はなるべくオリジナルの変数名等を変えずに行いました。
ただ、JavaScriptは変数名等で大文字と小文字を区別しますので、移植にあたっては小文字の変数名は先頭に”_”を付けて区別しています。

このプログラムの中身については、以下に概略を書いておきますが、きちんと理解するには高校程度の幾何(ベクトル)の知識が必要です。
オリジナルのほうの解説も参考にしてください(英語ですが)。

変数Wは、生成する画像の一辺のピクセル数です。
小さくすれば、画像も小さくなりますが、計算時間も短くなります。
この数値を変更する場合は、後述の球のデータ中の”240″という部分も、Wの値に一致させてください。

関数DOTは、2つのベクトルの内積(dot product)を計算します。
2つのベクトルa, bの内積は、|a||b| cosθ(|v|はvの長さ、θは2つのベクトルがなす角)になります。
しかし、実際には三角関数を使わずに計算することができます。(プログラムリストを見てください。)

内積は、方向の一致度合いを計算していると考えてください。
長さ1のベクトル同士では、完全に一致したとき1、直角の関係にあるとき0、完全に逆向きのとき-1になります。

配列SPHERESは、球に関するデータです。
1つの球について9つのパラメータがあります。
その内訳は、

・半径
・中心座標(X,Y,Z)
・色(R,G,B)
・輝き
・反射率

となっています。
色および反射率は、それぞれ0…9の範囲で指定してください。

光源の情報は配列LIGHTSに入っています。
1つの光源のパラメータは4つで、強度と光源の座標X,Y,Zです。
他に、環境光を表すAMBIENT_LIGHTというパラメータもあります。
光源の強度とAMBIENT_LIGHTは、合計して10になるように設定してください。

関数CLOSEST_INTERSECTIONは、点Bから方向Dへ向かう直線が最初にぶつかる球と、その交点を求めます。
どの球にもぶつからない場合は0、ぶつかる場合は正の数Sを返します(SPHERES[S]が、ぶつかった球の中心のX座標になります)。
交点を表すパラメータはグローバル変数Tに入ります。
交点は、B+D*Tになります。
TMINとTMAXは、ぶつかる点を探す範囲(Tの最小値と最大値)を指定します。

レイトレーシングの本体である関数TRACE_RAYは、各変数の意味を下図に示します。

pc30-4.jpg

一番最後にある関数RGB2は、プチコン3号の標準のRGB関数の代わりに使います。
プチコン3号では、RGBにそれぞれ8ビットの値を指定できますが、実際には上位5ビットしか使われません。
そのため、球面のようになめらかに階調が変化する画像では、色数が少ないことによる縞模様(バンディング)が目立ってしまいます。

RGB2は、乱数によって擬似的に中間色を生成します。
通常は、下位3ビットが0~7のどの値でも同じ色になってしまいますが、RGB2では(下位3ビットの値/8)の確率で1階調明るい色を出力します。
少ない階調で擬似的に中間色を表現する方法はディザリングといい、今回のものは超簡易版ですが、誤差拡散法などより高画質なアルゴリズムが各種提案されています。

次の図は、上が標準のRGB関数を使った場合、下がRGB2関数を使った場合です。

pc30-2.jpg
pc30-3.jpg