GLSLでレイトレーシング

『メンガーのスポンジ』

今日は、昔作ったプログラムを修正したりした。

立方体をくりぬいて出来る、自己相似なフラクタル図形。

これを普通にポリゴンで描画しようと思ったら、
n段階のスポンジは20のn乗の立方体で出来ているから・・・段階数増やすとちょい厳しい。

描画方法

3D描画には、キャンバスに図形(ポリゴンとか)を描いて行く方法と、
キャンバスの各ピクセルから視線を逆向きに辿って色を決める方法がある。
今回は、後者のレイトレーシングを使う。

プログラム概略は下図の通り。

フラクタル図形の隙間(空間立方体)の大きさ計算して、効率的に視線を逆に辿っていこうという作戦。
図では赤→橙→黄→緑　と4段階で立体まで視線を辿っている。

GLSLにするとこんなかんじ
vertex shader
ここでは特に何かやってるわけじゃない

varying vec3 pos;
varying vec3 view;
void main(){
	vec3 trans=gl_ModelViewMatrix[3].xyz;
	view=gl_Vertex.xyz+vec3(dot(gl_NormalMatrix[0],trans),dot(gl_NormalMatrix[1],trans),dot(gl_NormalMatrix[2],trans));
	pos=gl_Vertex.xyz;
	gl_Position=ftransform();
}

fragment shader
レイトレーシングを行う部分

varying vec3 pos;
varying vec3 view;
uniform int N;
uniform float SIZE;
const int MAXLOOP=100;
int inblock(float x,float y,float z){
	if(x<0.||y<0.||z<0.||x>SIZE||y>SIZE||z>SIZE)return 0;
	float S=SIZE/3.;
	int retval=int(S);
	for(int n=0;n<N;n++){
		int cnt=(x<S||x>2.*S?1:0)+(y<S||y>2.*S?1:0)+(z<S||z>2.*S?1:0);
		if(cnt==0||cnt==1)return retval;
		x-=floor(x/S)*S;x*=3.;
		y-=floor(y/S)*S;y*=3.;
		z-=floor(z/S)*S;z*=3.;
		retval/=3;
	}
	return -1;
}

vec3 move(inout vec3 p,vec3 b,vec3 v,float s){
	vec3 t1=(b-p)/v;
	vec3 t2=t1+s/v;
	float M=SIZE*65536.;
	if(t1.x<=0.)t1.x=M;if(t1.y<=0.)t1.y=M;if(t1.z<=0.)t1.z=M;
	if(t2.x<=0.)t2.x=M;if(t2.y<=0.)t2.y=M;if(t2.z<=0.)t2.z=M;
	
	float min1=min(min(t1.x,t1.y),t1.z);
	float min2=min(min(t2.x,t2.y),t2.z);
	if(min1<min2){
		p+=min1*v;
		if(min1==t1.x){p.x=b.x;return vec3(-0.5,0,0);}
		else if(min1==t1.y){p.y=b.y;return vec3(0,-0.5,0);}
		else{p.z=b.z;return vec3(0,0,-0.5);}
	}else{
		p+=min2*v;
		if(min2==t2.x){p.x=b.x+s;return vec3(0.5,0,0);}
		else if(min2==t2.y){p.y=b.y+s;return vec3(0,0.5,0);}
		else{p.z=b.z+s;return vec3(0,0,0.5);}
	}
}
vec3 trace(vec3 p,vec3 v){
	vec3 bv=vec3(0,0,0);
	if(v.x>0.&&p.x<=0.){p-=p.x*v/v.x;p.x=0.;bv=vec3(0.5,0,0);}
	if(v.y>0.&&p.y<=0.){p-=p.y*v/v.y;p.y=0.;bv=vec3(0,0.5,0);}
	if(v.z>0.&&p.z<=0.){p-=p.z*v/v.z;p.z=0.;bv=vec3(0,0,0.5);}
	if(v.x<0.&&p.x>=1.){p+=(1.-p.x)*v/v.x;p.x=1.;bv=vec3(-0.5,0,0);}
	if(v.y<0.&&p.y>=1.){p+=(1.-p.y)*v/v.y;p.y=1.;bv=vec3(0,-0.5,0);}
	if(v.z<0.&&p.z>=1.){p+=(1.-p.z)*v/v.z;p.z=1.;bv=vec3(0,0,-0.5);}
	p*=SIZE;
	
	for(int i=0;i<MAXLOOP;i++){
		vec3 pp=p+bv;
		int size=inblock(pp.x,pp.y,pp.z);
		if(size==0)return vec3(0,0,0);
		if(size<0)return p/SIZE*(0.9+0.05*dot(vec3(1,2,3),bv));
		float s=float(size);
		vec3 b=floor((p+bv)/s)*s;
		bv=move(p,b,v,s);
	}
	return p/SIZE/2.;
}

void main(void){
	gl_FragColor.a=1.;
	gl_FragColor.rgb=trace((pos+vec3(1,1,1))/2.,normalize(view));
	if(gl_FragColor.a==0.)gl_FragColor.rgb=vec3(1,1,0);
}

inblock 再帰関数(for文に展開済)　フラクタルの形を決定付けている
move 隙間の境界まで視線を辿る関数
trace inblockとmoveを繰り返して視線を逆に辿って行く
場合分けが多くて嫌な感じ。

MacOSX付属のOpenGLShaderBuilderに読み込ませて、
Nに段階数(例 4)、SIZEに3のN乗(例 81)を設定してやれば描画できる。

他のボリュームレンダリングとかにも応用できたり。

P.S.
目指せ3日坊主回避