同样场景，Taichi版比p5.js版大概快了20倍左右（！）
现在还只支持线段与射线相交，而且衰减函数也是拍脑袋想的，而且可能是因为我写的有误差再加上一些微妙的差别，两个版本的结果看起来有些不同，不过大体上我已经知道了我想知道的这两者之间的数量级的差距 XD
JavaScript版：
https://editor.p5js.org/nitroglycerine/sketches/yDq3_-hY2
Taichi版：
# -*- coding: utf-8 -*-
import taichi as ti
# Javascript: 0.1 FPS
# i7-6820HQ CPU: 2.70 FPS
# Quadro M1000M CUDA: 2.03 FPS
ti.init(arch=ti.cpu)
# 把场景写进kernel导致编译时间太长，故关掉高级优化
ti.core.toggle_advanced_optimization(False)
W = 720; H = 480;
MAGNIFY_DISP = 4
color_buffer = ti.Vector(3, dt=ti.f32, shape=(W, H));
gui = ti.GUI('Haha', (W, H))
# Taichi's coordinate system is different from Processing's
#
# Taichi Processing
#
# ^ (0, 480) +------------> (720,0
# | |
# | |
# +----------> (720, 0) V (0, 480)
@ti.func
def Perp(v):
return ti.Vector([-v[1], v[0], 0.0])
# 返回：[t, tact]
@ti.func
def IntersectLineWithRay(l, r):
p1 = l[0]; p2 = l[1]; o = r[0]; d = ti.normalized(r[1]);
p1p2 = p2 - p1
n = ti.normalized(Perp(p1p2))
op1 = p1 - o;
dist = ti.dot(n, op1);
ddn = ti.dot(n, d);
t = dist / ddn;
tact = ti.Vector([0.0, 0.0, 0.0]);
if t > 0:
tact = o + t * d
t0 = ti.dot((tact - p1), p1p2)
t1 = ti.dot((p2 - tact), p1p2)
completion = t0 / (t0+t1)
if completion >= 0 and completion <= 1:
pass
else:
t = -1
return t, tact
@ti.kernel
def Clear():
for i, j in color_buffer:
color_buffer[i, j] = ti.Vector([0,0,0]);
@ti.func
def Diffuse(l, d):
n = Perp(l[1]-l[0]);
if (ti.dot(n, d) > 0): n = n*(-1)
theta = ti.random(ti.f32) * 3.14159*2
ret = ti.Vector([ti.cos(theta), ti.sin(theta), 0.0]);
for x in range(0, 4):
theta = ti.random(ti.f32) * 3.14159*2
ret = ti.Vector([ti.cos(theta), ti.sin(theta), 0.0]);
if ti.dot(n, ret) < 0:
break
return ret
NSAMP = 4;
x0 = 0.2*W; x1 = 0.25*W; x2 = 0.3*W; y1 = 0.7*H; y2 = 0.55*H; y0 = 0.85*H;
g1 = [ti.Vector([x1,y1,0]), ti.Vector([x2,y1,0])]
g2 = [ti.Vector([x2,y1,0]), ti.Vector([x2,y2,0])]
g3 = [ti.Vector([x0,y2,0]), ti.Vector([x2,y2,0])]
g4 = [ti.Vector([x0,y2,0]), ti.Vector([x0,y0,0])]
g5 = [ti.Vector([x0,y0,0]), ti.Vector([x2,y0,0])]
#####
#
# ###
# #
#####
@ti.func
def IntersectG(o, d, t_, this_intensity_, done_):
d_out = d;
t = t_;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(g1, [o, d]);
if t1 > 0 and t>t1:
t = t1;
this_intensity = ti.Vector([0.5, 0.0, 0.0]);
done = True
t1, _ = IntersectLineWithRay(g2, [o, d]);
if t1 > 0 and t>t1: t = t1; d_out = Diffuse(g2, d); done = False;
t1, _ = IntersectLineWithRay(g3, [o, d]);
if t1 > 0 and t>t1: t = t1; this_intensity = ti.Vector([0.5, 0.0, 0.0]); done = True
t1, _ = IntersectLineWithRay(g4, [o, d]);
if t1 > 0 and t>t1: t = t1; d_out = Diffuse(g4, d); done = False
t1, _ = IntersectLineWithRay(g5, [o, d]);
if t1 > 0 and t>t1: t = t1; this_intensity = ti.Vector([0.5, 0.0, 0.0]); done = True
return t, d_out, this_intensity, done
#
# #
# #
#####
# #
x0 = 0.325*W; x1=0.375*W; x2=0.425*W; y1=0.7*H; y2=0.55*H; y0=0.85*H; y3=0.65*H;
a1 = [ti.Vector([x1,y0,0]), ti.Vector([x2,y1,0])]
a2 = [ti.Vector([x1,y0,0]), ti.Vector([x0,y1,0])]
a3 = [ti.Vector([x0,y3,0]), ti.Vector([x2,y3,0])]
a4 = [ti.Vector([x0,y1,0]), ti.Vector([x0,y2,0])]
a5 = [ti.Vector([x2,y1,0]), ti.Vector([x2,y2,0])]
@ti.func
def IntersectA(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(a1, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(a1, d); done = False;
t1, _ = IntersectLineWithRay(a2, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(a2, d); done = False;
t1, _ = IntersectLineWithRay(a3, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.0, 0.8, 0.0]); done = True
t1, _ = IntersectLineWithRay(a4, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(a4, d); done = False;
t1, _ = IntersectLineWithRay(a5, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(a5, d); done = False;
return t, d_out, this_intensity, done;
# #
## ##
# # #
# #
# #
x0=0.45*W; x1=0.5*W; x2=0.55*W; y1=0.7*H; y2=0.55*H; y0=0.85*H;
m1 = [ti.Vector([x0,y0,0]), ti.Vector([x0,y2,0])]
m2 = [ti.Vector([x0,y0,0]), ti.Vector([x1,y2,0])]
m3 = [ti.Vector([x2,y0,0]), ti.Vector([x1,y2,0])]
m4 = [ti.Vector([x2,y0,0]), ti.Vector([x2,y2,0])]
@ti.func
def IntersectM(o, d, t_, this_intensity_, done_):
d_out = d;
t = t_;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(m1, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(m1, d); done = False;
t1, _ = IntersectLineWithRay(m2, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.4, 0.4, 0.0]); done = True
t1, _ = IntersectLineWithRay(m3, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.4, 0.4, 0.0]); done = True
t1, _ = IntersectLineWithRay(m4, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(m4, d); done = False;
return t, d_out, this_intensity, done;
#####
#
#####
#
#####
x0=0.575*W; x1=0.625*W; x2=0.675*W; y1=0.7*H; y2=0.55*H; y0=0.85*H;
e1 = [ti.Vector([x0,y0,0]), ti.Vector([x0,y2,0])]
e2 = [ti.Vector([x0,y0,0]), ti.Vector([x2,y0,0])]
e3 = [ti.Vector([x0,y2,0]), ti.Vector([x2,y2,0])]
e4 = [ti.Vector([x0,y1,0]), ti.Vector([x1,y1,0])]
@ti.func
def IntersectE(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(e1, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(e1, d); done = False;
t1, _ = IntersectLineWithRay(e2, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0, 0.5, 0.5]); done = True
t1, _ = IntersectLineWithRay(e3, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0, 0.5, 0.5]); done = True
t1, _ = IntersectLineWithRay(e4, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0, 0.5, 0.5]); done = True
return t, d_out, this_intensity, done;
#####
#
#####
#
#####
x0=0.70*W; x1=0.75*W; x2=0.80*W; y1=0.7*H; y2=0.55*H; y0=0.85*H;
s1 = [ti.Vector([x0,y0,0]), ti.Vector([x2,y0,0])]
s2 = [ti.Vector([x0,y0,0]), ti.Vector([x0,y1,0])]
s3 = [ti.Vector([x0,y1,0]), ti.Vector([x2,y1,0])]
s4 = [ti.Vector([x2,y1,0]), ti.Vector([x2,y2,0])]
s5 = [ti.Vector([x0,y2,0]), ti.Vector([x2,y2,0])]
@ti.func
def IntersectS(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_
done = done_
t1, _ = IntersectLineWithRay(s1, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.5, 0.2, 0.2]); done = True
t1, _ = IntersectLineWithRay(s2, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(s2, d); done = False;
t1, _ = IntersectLineWithRay(s3, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.5, 0.2, 0.2]); done = True
t1, _ = IntersectLineWithRay(s4, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(s4, d); done = False;
t1, _ = IntersectLineWithRay(s5, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.5, 0.2, 0.2]); done = True
return t, d_out, this_intensity, done;
#####
#
#####
#
#####
x0=0.25*W; x1=0.3*W; x2=0.35*W; y1=0.3*H; y2=0.15*H; y0=0.45*H;
v21 = [ti.Vector([x0,y0,0]), ti.Vector([x2,y0,0])]
v22 = [ti.Vector([x2,y0,0]), ti.Vector([x2,y1,0])]
v23 = [ti.Vector([x0,y1,0]), ti.Vector([x2,y1,0])]
v24 = [ti.Vector([x0,y1,0]), ti.Vector([x0,y2,0])]
v25 = [ti.Vector([x0,y2,0]), ti.Vector([x2,y2,0])]
@ti.func
def Intersect2(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(v21, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v21, d); done = False;
t1, _ = IntersectLineWithRay(v22, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v22, d); done = False;
t1, _ = IntersectLineWithRay(v23, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v23, d); done = False;
t1, _ = IntersectLineWithRay(v24, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v24, d); done = False;
t1, _ = IntersectLineWithRay(v25, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v25, d); done = False;
return t, d_out, this_intensity, done;
#####
# #
# #
# #
#####
x0=0.45*W; x1=0.5*W; x2=0.55*W; y1=0.3*H; y2=0.15*H; y0=0.45*H;
v01 = [ti.Vector([x0,y0,0]), ti.Vector([x2,y0,0])]
v02 = [ti.Vector([x2,y0,0]), ti.Vector([x2,y2,0])]
v03 = [ti.Vector([x2,y2,0]), ti.Vector([x0,y2,0])]
v04 = [ti.Vector([x0,y0,0]), ti.Vector([x0,y2,0])]
@ti.func
def Intersect0(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(v01, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.2, 0.2, 0.4]); done = True
t1, _ = IntersectLineWithRay(v02, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v02, d); done = False;
t1, _ = IntersectLineWithRay(v03, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.2, 0.2, 0.4]); done = True
t1, _ = IntersectLineWithRay(v04, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v04, d); done = False;
return t, d_out, this_intensity, done;
#
##
# #
#
#####
x0=0.65*W; x1=0.7*W; x2=0.75*W; y1=0.3*H; y2=0.15*H; y0=0.45*H;
v11 = [ti.Vector([x1,y0,0]), ti.Vector([x1-0.05*W,y0-0.075*H,0])]
v12 = [ti.Vector([x1,y0,0]), ti.Vector([x1,y2,0])]
v13 = [ti.Vector([x2,y2,0]), ti.Vector([x0,y2,0])]
@ti.func
def Intersect1(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(v11, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.4, 0.4, 0.8]); done = True
t1, _ = IntersectLineWithRay(v12, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v12, d); done = False;
t1, _ = IntersectLineWithRay(v13, [o, d]);
if t1 > 0 and t > t1: t = t1; d_out = Diffuse(v13, d); done = False;
return t, d_out, this_intensity, done;
##############
# #
# #
# #
# #
##############
vw1 = [ti.Vector([0.15*W,0.5*H,0]), ti.Vector([0.15*W,0.1*H,0])]
vw2 = [ti.Vector([0.85*W,0.5*H,0]), ti.Vector([0.85*W,0.1*H,0])]
vw3 = [ti.Vector([0.15*W,0.05*H,0]), ti.Vector([0.85*W,0.05*H,0])]
vw4 = [ti.Vector([0.15*W,0.95*H,0]), ti.Vector([0.85*W,0.95*H,0])]
vw5 = [ti.Vector([0.15*W,0.9*H,0]), ti.Vector([0.15*W,0.5*H,0])]
vw6 = [ti.Vector([0.85*W,0.9*H,0]), ti.Vector([0.85*W,0.5*H,0])]
@ti.func
def IntersectW(o, d, t_, this_intensity_, done_):
t = t_;
d_out = d;
this_intensity = this_intensity_;
done = done_
t1, _ = IntersectLineWithRay(vw1, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.6, 0.6, 1.4]); done = True
t1, _ = IntersectLineWithRay(vw2, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.6, 0.6, 1.4]); done = True
t1, _ = IntersectLineWithRay(vw3, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.6, 0.6, 1.4]); done = True
t1, _ = IntersectLineWithRay(vw4, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.4, 0.8, 0.8]); done = True
t1, _ = IntersectLineWithRay(vw5, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.1, 0.2, 0.2]); done = True
t1, _ = IntersectLineWithRay(vw6, [o, d]);
if t1 > 0 and t > t1: t = t1; this_intensity = ti.Vector([0.1, 0.2, 0.2]); done = True
return t, d_out, this_intensity, done;
@ti.kernel
def Fill():
for i, j in color_buffer:
intensity = ti.Vector([0.0, 0.0, 0.0])
samp = 0
while samp < NSAMP: # 每个像素采这么多次
o = ti.Vector([i*1.0, j*1.0, 0.0]); # 必须是三维的不能是二维的？
theta = ti.random(ti.f32) * 3.14159*2
d = ti.Vector([ti.cos(theta), ti.sin(theta), 0.0])
# obj = None 不能设为none？
t_total = 0.0 # 光走过的距离
level = 0
while level < 3:
d_orig = d;
t = 1000000000.0
# 求交点
# this_intensity = ti.Vector([0,0,0]) # 错误：整型的不能赋浮点值
this_intensity = ti.Vector([0.0, 0.0, 0.0])
done = False
t, d, this_intensity, done = IntersectG(o, d_orig, t, this_intensity, done); # G
t, d, this_intensity, done = IntersectA(o, d_orig, t, this_intensity, done); # A
t, d, this_intensity, done = IntersectM(o, d_orig, t, this_intensity, done); # M
t, d, this_intensity, done = IntersectE(o, d_orig, t, this_intensity, done); # E
t, d, this_intensity, done = IntersectS(o, d_orig, t, this_intensity, done); # S
t, d, this_intensity, done = Intersect2(o, d_orig, t, this_intensity, done); # 2
t, d, this_intensity, done = Intersect0(o, d_orig, t, this_intensity, done); # 0
t, d, this_intensity, done = Intersect1(o, d_orig, t, this_intensity, done); # 1
t, d, this_intensity, done = IntersectW(o, d_orig, t, this_intensity, done); # Walls
# 总结一下求交点结果决定是再发射间接光还是退出
if t != 1000000000.0:
t_total += t
damp = 1 / (1 + ti.exp(t_total/400))
EPS = 0.001;
o = o + d_orig * (t - EPS);
if done: intensity += this_intensity * damp; break
level += 1 # 一次递归结束了
samp += 1; # 一次采样结束了
color_buffer[i, j] += intensity * 1.0 / NSAMP
g_frame_count = 0;
while True:
while gui.get_event(ti.GUI.PRESS):
if gui.event.key == ti.GUI.ESCAPE:
exit()
#Clear();
Fill();
g_frame_count += 1;
img = (color_buffer).to_numpy() * (1.0 / (g_frame_count)) * MAGNIFY_DISP;
gui.set_image(img);
gui.show()