code:
# Created by inigo quilez - iq/2020
# License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
#
# The Julia set of f(z) = z³ + c, as rendered for the Youtube
# video called "Geodes": https://www.shadertoy.com/view/3llyzl
#
# I simplified a few things, reduced the number of GI bounces
# and did some temporal reprojection to keep it more or less
# real-time while looking similar to the one in the video.
#
# Explanations:
# https://iquilezles.org/www/articles/distancefractals/distancefractals.htm
# https://iquilezles.org/www/articles/orbittraps3d/orbittraps3d.htm
#
# Related shaders:
#
# Julia - Quaternion 1 : https://www.shadertoy.com/view/MsfGRr
# Julia - Quaternion 2 : https://www.shadertoy.com/view/lsl3W2
# Julia - Quaternion 3 : https://www.shadertoy.com/view/3tsyzl
import taichi as ti
import taichi_glsl as ts
import time
ti.init(ti.cpu)
iResolution = ti.Vector([1024, 512])
pixels = ti.Vector.field(4, ti.f32, (iResolution.x, iResolution.y))
iTime = ti.field(ti.f32, ())
iFrame = ti.field(ti.i32, ())
iChannel0 = ti.Vector.field(4, ti.f32, (1024, 1024))
@ti.func
def xy(v):
return ti.Vector([v.x, v.y])
@ti.func
def xz(v):
return ti.Vector([v.x, v.z])
@ti.func
def xxx(v):
return ti.Vector([v.x, v.x, v.x])
@ti.func
def xyy(v):
return ti.Vector([v.x, v.y, v.y])
@ti.func
def xyz(v):
return ti.Vector([v.x, v.y, v.z])
@ti.func
def yxy(v):
return ti.Vector([v.y, v.x, v.y])
@ti.func
def yyx(v):
return ti.Vector([v.y, v.y, v.x])
@ti.func
def texture(sampler, P):
P = ts.clamp(P, 0.0, 1.0)
return sampler[ti.cast(sampler.shape * P, ti.i32)]
@ti.func
def textureLod(sampler, P, lod):
return texture(sampler, P)
@ti.func
def texelFetch(sampler, P, lod):
return sampler[P]
@ti.kernel
def makeIChannel0():
for texCoord_x, texCoord_y in iChannel0:
iChannel0[texCoord_x, texCoord_y] = ti.Vector([
ti.sin(0.5 * texCoord_x),
1.0 - ti.sin(0.5 * texCoord_x),
0.0,
1.0,
])
# disable TRAPs to see just the set
TRAPS = True
# disable CUT to see the whole set
CUT = True
kNumIte = 200
kPrecis = 0.00025
kC = ti.Vector([-2, 6, 15, -6]) / 22.0
kFocLen = 3.0
if TRAPS:
kBSRad = 2.0
else:
kBSRad = 1.2
kNumBounces = 3
# --------------------------------------
# oldschool rand() from Visual Studio
# --------------------------------------
seed = ti.field(ti.i32, ())
seed[None] = 1
@ti.func
def srand(s):
seed[None] = s
@ti.func
def rand():
seed[None] = seed[None] * 0x343fd + 0x269ec3
return (seed[None] >> 16) & 32767
@ti.func
def frand():
return ti.cast(rand(), ti.f32) / 32767.0
# --------------------------------------
# hash to initialize the random seed (copied from Hugo Elias)
# --------------------------------------
@ti.func
def hash(n):
n = (n << 13) ^ n
return n * (n * n * 15731 + 789221) + 1376312589
#--------------------------------------------------------------------------------
# http://amietia.com/lambertnotangent.html
#--------------------------------------------------------------------------------
@ti.func
def cosineDirection(nor):
u = frand() * 2.0 - 1.0
a = frand() * 6.28318531
return (nor + ti.Vector([
ti.sqrt(1.0 - u * u) * ti.cos(a),
ti.sqrt(1.0 - u * u) * ti.sin(a),
u,
])).normalized()
#--------------------------------------------------------------------------------
# quaternion manipulation
#--------------------------------------------------------------------------------
@ti.func
def qSquare(q):
return ti.Vector([
q.x * q.x - q.y * q.y - q.z * q.z - q.w * q.w,
2.0 * q.x * q.y,
2.0 * q.x * q.z,
2.0 * q.x * q.w,
])
@ti.func
def qCube(q):
q2 = q * q
return ti.Vector([
q.x * (q2.x - 3.0 * q2.y - 3.0 * q2.z - 3.0 * q2.w),
q.y * (3.0 * q2.x - q2.y - q2.z - q2.w),
q.z * (3.0 * q2.x - q2.y - q2.z - q2.w),
q.w * (3.0 * q2.x - q2.y - q2.z - q2.w),
])
@ti.func
def qLength2(q):
return q.dot(q)
#--------------------------------------------------------------------------------
# ray-sphere intersection
# http://iquilezles.org/www/articles/intersectors/intersectors.htm
#--------------------------------------------------------------------------------
@ti.func
def iSphere(ro, rd, rad):
result = ti.Vector([-1.0, 0.0])
b = ro.dot(rd)
c = ro.dot(ro) - rad * rad
h = b * b - c
if h >= 0.0:
h = ti.sqrt(h)
result = ti.Vector([-b - h, -b + h])
return result
#--------------------------------------------------------------------------------
# build camera rotation matrix
#--------------------------------------------------------------------------------
@ti.func
def setCamera(ro, ta, cr):
cw = (ta - ro).normalized()
cp = ti.Vector([ti.sin(cr), ti.cos(cr), 0.0])
cu = (cw.cross(cp)).normalized()
cv = (cu.cross(cw)).normalized()
return ti.Matrix.rows([cu, cv, cw])
#--------------------------------------------------------------------------------
# SDF of the Julia set z³+c
# https://iquilezles.org/www/articles/distancefractals/distancefractals.htm
#--------------------------------------------------------------------------------
@ti.func
def mapFunc(p):
z = ti.Vector([p.x, p.y, p.z, 0.0])
dz2 = 1.0
m2 = 0.0
n = 0.0
if ti.static(TRAPS):
o = 1e10
for i in ti.static(range(kNumIte)):
# z' = 3z² -> |z'|² = 9|z²|²
dz2 *= 9.0 * qLength2(qSquare(z))
# z = z³ + c
z = qCube(z) + kC
# stop under divergence
m2 = qLength2(z)
# orbit trapping : https://iquilezles.org/www/articles/orbittraps3d/orbittraps3d.htm
if ti.static(TRAPS):
o = ti.min(o, (xz(z) - ti.Vector([0.45, 0.55])).norm() - 0.1)
# exit condition
if m2 > 256.0:
break
n += 1.0
# sdf(z) = log|z|·|z|/|dz| : https://iquilezles.org/www/articles/distancefractals/distancefractals.htm
d = 0.25 * ti.log(m2) * ti.sqrt(m2 / dz2)
if ti.static(TRAPS):
d = ti.min(o, d)
if ti.static(CUT):
d = ti.max(d, p.y)
return ti.Vector([d, n])
#--------------------------------------------------------------------------------
# Compute Normal to SDF
#--------------------------------------------------------------------------------
#if 1
# https://iquilezles.org/www/articles/normalsSDF/normalsSDF.htm
@ti.func
def calcNormal1(pos):
e = ti.Vector([1.0, -1.0]) * 0.5773 * kPrecis
return (xyy(e) * mapFunc(pos + xyy(e)).x +
yyx(e) * mapFunc(pos + yyx(e)).x +
yxy(e) * mapFunc(pos + yxy(e)).x +
xxx(e) * mapFunc(pos + xxx(e)).x).normalized()
# else
# https://iquilezles.org/www/articles/juliasets3d/juliasets3d.htm
@ti.func
def calcNormal(p):
result = ti.Vector([0.0, 1.0, 0.0])
# if ti.static(TRAPS):
# the code below only works for the actual Julia set, not the traps
z = ti.Vector([p.x, p.y, p.z, 0.0])
# identity derivative
J = ti.Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
])
for i in ti.static(range(kNumIte)):
# f(q) = q³ + c =
# x = x²x - 3y²x - 3z²x - 3w²x + c.x
# y = 3x²y - y²y - z²y - w²y + c.y
# z = 3x²z - y²z - z²z - w²z + c.z
# w = 3x²w - y²w - z²w - w²w + c.w
#
# Jacobian, J(f(q)) =
# 3(x²-y²-z²-w²) 6xy 6xz 6xw
# -6xy 3x²-3y²-z²-w² -2yz -2yw
# -6xz -2yz 3x2-y²-3z²-w² -2zw
# -6xw -2yw -2zw 3x²-y²-z²-3w²
k1 = 6.0 * z.x * z.y
k2 = 6.0 * z.x * z.z
k3 = 6.0 * z.x * z.w
k4 = 2.0 * z.y * z.z
k5 = 2.0 * z.y * z.w
k6 = 2.0 * z.z * z.w
sx = z.x * z.x
sy = z.y * z.y
sz = z.z * z.z
sw = z.w * z.w
mx = 3.0 * sx - 3.0 * sy - 3.0 * sz - 3.0 * sw
my = 3.0 * sx - 3.0 * sy - sz - sw
mz = 3.0 * sx - sy - 3.0 * sz - sw
mw = 3.0 * sx - sy - sz - 3.0 * sw
# chain rule of jacobians
J = J * ti.Matrix([
[mx, -k1, -k2, -k3],
[k1, my, -k4, -k5],
[k2, -k4, mz, -k6],
[k3, -k5, -k6, mw],
])
# q = q³ + c
z = qCube(z) + kC
# exit condition
if z.dot(z) > 256.0:
break
if p.y <= 0.0:
result = xyz(J@z).normalized()
return result
#--------------------------------------------------------------------------------
# ray-scene intersection
#--------------------------------------------------------------------------------
@ti.func
def raycast(ro, rd):
result = ti.Vector([-2.0, 0.0])
tmax = 7.0
tmin = kPrecis
# intersect clipping plane
if ti.static(CUT):
kSplit = 0.01
tpS = (kSplit - ro.y) / rd.y
if tpS > 0.0:
if ro.y > kSplit:
tmin = ti.max(tmin, tpS)
else:
tmax = ti.min(tmax, tpS)
# intersect lower clipping plane
if ti.static(True):
tpF = (-0.8 - ro.y) / rd.y
if tpF > 0.0:
tmax = ti.min(tmax, tpF)
# intersect bounding sphere
hasResult = False
if ti.static(True):
bv = iSphere(ro, rd, kBSRad)
if bv.y < 0.0:
hasResult = True
else:
tmin = ti.max(tmin, bv.x)
tmax = ti.min(tmax, bv.y)
if not hasResult:
# raymarch
res = ti.Vector([-1.0, 0.0])
t = tmin
lt = 0.0
lh = 0.0
for i in range(1024):
res = mapFunc(ro + rd * t)
if res.x < kPrecis:
break
lt = t
lh = res.x
if ti.static(not TRAPS):
t += ti.min(res.x, 0.2)
else:
t += ti.min(res.x, 0.01) * (0.5 * 0.5 * frand())
if t > tmax:
break
# linear interpolation for better isosurface
if lt > 0.0001 and res.x < 0.0:
t = lt - lh * (t - lt) / (res.x - lh)
if t < tmax:
res.x = t
else:
res.x = -1.0
result = res
return result
#--------------------------------------------------------------------------------
# color of the surface
#--------------------------------------------------------------------------------
@ti.func
def colorSurface(pos, nor, tn):
col = 0.5 + 0.5 * ti.cos(ti.log(tn.y) / ti.log(2.0)
* 0.9 + 3.5 + ti.Vector([0.0, 0.6, 1.0]))
if pos.y > 0.0:
col = ts.mix(col, ti.Vector([1.0, 0.0, 0.0]), 0.2)
inside = ts.smoothstep(14.0, 15.0, tn.y)
col *= ti.Vector([0.45, 0.42, 0.40]) + \
ti.Vector([0.55, 0.58, 0.60]) * inside
col = ts.mix(col * col * (3.0 - 2.0 * col), col, inside)
col = ts.mix(ts.mix(col, ti.Vector(
col.dot(ti.Vector([0.3333, 0.0, 0.0]))), -0.4), col, inside)
return ts.clamp(col * 0.65, 0.0, 1.0)
#--------------------------------------------------------------------------------
# Render the scene through super simplified path-tracing
#--------------------------------------------------------------------------------
@ti.func
def render(fragCoord, ro, rd, resPos, resT):
result = ti.Vector([0.0, 0.0, 0.0])
colorMask = ti.Vector([1.0, 0.0, 0.0])
resT = 1e20
# path-tracing
for bounce in ti.static(range(kNumBounces)):
tn = raycast(ro, rd)
t = tn.x
if t < 0.0:
if bounce > 0:
result = colorMask * 1.65 * ts.step(0.0, rd.y)
break
else:
result = ti.Vector(
[ts.clamp(0.02 + 0.021 * rd.y, 0.0, 1.0), 0.0, 0.0])
break
else:
pos = ro + t * rd
nor = calcNormal(pos)
if bounce == 0:
resT = t
resPos = pos
colorMask *= colorSurface(pos, nor, tn)
rd = cosineDirection(nor)
ro = pos + nor * kPrecis
return result, resPos, resT
@ti.func
def mainImage(fragColor, fragCoord):
#-----------------------------------------------
# init random seed
#-----------------------------------------------
q = ti.cast(fragCoord, ti.i32)
srand(hash(q.x + hash(q.y + hash(1117 * iFrame[None]))))
#-----------------------------------------------
# camera
#-----------------------------------------------
an = 0.5 + iTime[None] * 0.03
ro = 2.0 * ti.Vector([ti.sin(an), 0.8, ti.cos(an)])
ta = ti.Vector([0.0, -0.1, 0.0])
if ti.static(CUT):
ta = ti.Vector([0.0, -0.3, 0.0])
cam = setCamera(ro, ta, 0.0)
#-----------------------------------------------
# ray direction
#-----------------------------------------------
p = (2.0 * fragCoord - xy(iResolution)) / iResolution.y
rd = (cam @ ti.Vector([p.x, p.y, kFocLen])).normalized()
#-----------------------------------------------
# render fractal
#-----------------------------------------------
pos = ti.Vector([0.0, 0.0, 0.0])
resT = 0.0
col, pos, resT = render(fragCoord, ro, rd, pos, resT)
#-----------------------------------------------
# reproject to previous frame and pull history
#-----------------------------------------------
oldCam = ti.Matrix.rows([
texelFetch(iChannel0, ti.Vector([0, 0]), 0),
texelFetch(iChannel0, ti.Vector([1, 0]), 0),
texelFetch(iChannel0, ti.Vector([2, 0]), 0),
])
# world space point
wpos = ti.Vector([pos.x, pos.y, pos.z, 1.0])
# convert to camera space (note inverse multiply)
cpos = oldCam @ wpos
# convert to NDC space (project)
npos = kFocLen * xy(cpos) / cpos.z
# convert to screen space
spos = 0.5 + 0.5 * npos * ti.Vector([iResolution.y / iResolution.x, 1.0])
# convert to raster space
rpos = spos * xy(iResolution)
# read color+depth from this point's previous screen location
ocolt = textureLod(iChannel0, spos, 0.0)
# if we consider the data contains the history for this point
if iFrame[None] > 0 and resT < 100.0 and (rpos.y > 1.5 or rpos.x > 3.5):
# blend with history (it's a IIR low pas filter really)
col = ts.mix(xyz(ocolt), col, 0.06)
# output
if q.y == 0 and q.x < 3:
# camera matrix in lower left three pixels, for next frame
if q.x == 0:
fragColor = ti.Vector([cam[0, 0], cam[0, 1], cam[0, 2], -
1.0 * ti.Vector([cam[0, 0], cam[0, 1], cam[0, 2]]).dot(ro)])
elif q.x == 1:
fragColor = ti.Vector([cam[1, 0], cam[1, 1], cam[1, 2], -
1.0 * ti.Vector([cam[1, 0], cam[1, 1], cam[1, 2]]).dot(ro)])
else:
fragColor = ti.Vector([cam[2, 0], cam[2, 1], cam[2, 2], -
1.0 * ti.Vector([cam[2, 0], cam[2, 1], cam[2, 2]]).dot(ro)])
else:
# color and depth
fragColor = ti.Vector([col.x, col.y, col.z, resT])
return fragColor
@ti.kernel
def fragmentShader():
for fragCoord_x, fragCoord_y in pixels:
pixels[fragCoord_x, fragCoord_y] = mainImage(
pixels[fragCoord_x, fragCoord_y],
ti.Vector([fragCoord_x, fragCoord_y])
)
if __name__ == "__main__":
gui = ti.GUI("Julia Quaternion 3", pixels.shape)
makeIChannel0()
iTime[None] = 0
iFrame[None] = 0
preTime = time.time()
while gui.running:
nowTime = time.time()
iTime[None] += (nowTime - preTime)
preTime = nowTime
iFrame[None] += 1
fragmentShader()
gui.set_image(pixels)
gui.show()
error:
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py", line 552, in <module>
fragmentShader()
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/kernel_impl.py", line 715, in wrapped
return primal(*args, **kwargs)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/kernel_impl.py", line 643, in __call__
key = self.ensure_compiled(*args)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/kernel_impl.py", line 629, in ensure_compiled
self.materialize(key=key, args=args, arg_features=arg_features)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/kernel_impl.py", line 475, in materialize
taichi_kernel = _ti_core.create_kernel(taichi_ast_generator,
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/kernel_impl.py", line 465, in taichi_ast_generator
transform_tree(tree, ctx)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/transform.py", line 8, in transform_tree
ASTTransformer()(ctx, tree)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 26, in __call__
raise e
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 23, in __call__
return method(ctx, node)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 484, in build_Module
build_stmt(ctx, stmt)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 26, in __call__
raise e
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 23, in __call__
return method(ctx, node)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 448, in build_FunctionDef
build_stmts(ctx, node.body)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 1083, in build_stmts
build_stmt(ctx, stmt)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 26, in __call__
raise e
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 23, in __call__
return method(ctx, node)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 911, in build_For
return ASTTransformer.build_struct_for(ctx,
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 829, in build_struct_for
build_stmts(ctx, node.body)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 1083, in build_stmts
build_stmt(ctx, stmt)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 26, in __call__
raise e
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 23, in __call__
return method(ctx, node)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer.py", line 69, in build_Assign
build_stmt(ctx, node.value)
File "/Users/vw/miniconda3/envs/taichi/lib/python3.8/site-packages/taichi/lang/ast/ast_transformer_utils.py", line 32, in __call__
raise TaichiCompilationError(msg)
taichi.lang.exception.TaichiCompilationError: On line 535 of file "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py":
pixels[fragCoord_x, fragCoord_y] = mainImage(
^^^^^^^^^^
pixels[fragCoord_x, fragCoord_y],
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ti.Vector([fragCoord_x, fragCoord_y])
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
On line 487 of file "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py":
col, pos, resT = render(fragCoord, ro, rd, pos, resT)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
On line 433 of file "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py":
tn = raycast(ro, rd)
^^^^^^^^^^^^^^^
On line 380 of file "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py":
res = mapFunc(ro + rd * t)
^^^^^^^^^^^^^^^^^^^^
On line 241 of file "/Users/vw/Projects/Taichi-Julia-Quaternion-3/main.py":
o = ti.min(o, (xz(z) - ti.Vector([0.45, 0.55])).norm() - 0.1)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: Invalid constant scalar expression: <class 'NoneType'>
这个报错没看懂,前面我好几个函数都这样写也没问题。然后这个函数改成不是ti.func,就没错了。我这里是错了什么?