Chango的数学Shader世界（十七）RayTrace三维分形（二）—

目的：

移植Inigo Quilez的可视化三维Julia的算法到ue4。

本篇仅考虑造型，忽略表面法线计算和光照模型等。

忽略分形相关细节数学原理。

参考：

最近墙越来越高，我也没法放链接了。

1.Inigo Quilez在ShaderToy上的Julia，网址后缀MsfGRr

2.维度：数学漫步第一季的复数（下）

3.Ray Tracing Quaternion Julia Sets on the GPU——Keenan Crane，University of Illinois at Urbana-Champaign

pdf,非论文。

观察：

我们意识到难点有二：

1.完整的JuliaSet没有解析式，只有递推式。

若要判断某一点是否可由初始集合经过k次迭代能被包含，消耗也是巨大而不现实的。

因此我们似乎既没办法直接RayTrace，也没办法RayMarch出一个Julia云。

（下面是a=0时的2D情况，演示了2次迭代，这个特殊情况我们能判断点是否在Julia内，因为Julia此时为unit disk）

2.标准Julia集极其小，看不清。

由于是z^2+a，“可见部分”就跟单位圆差不多大小。我们希望搞一个actor，能随意缩放，位移Julia。但又不改变Julia的那些公式（因为理论公式现成，要改麻烦）。

分析：

1.问题一：根据距离场RayMarch到Julia表面

简而言之，虽无Julia表面解析式，但有其距离场估计的表达式。

距离场是一个数值场S=F(x,y,z)，空间中每点数值为距离Julia表面最近的距离。

伪代码

bool intersectJulia(...)
{float h = 1.0;float t = 0.0;for( int i=0; i<最大March次数; i++ ){if( h<0.0001（碰到表面） ||t>最远距离 （没有碰到表面） ){break;}h = MarchStep( ro+rd*t, ...);t += h;}if( t<最远距离) { return true; }return false;
}

2.问题二：移动ro位置

Julia所有表达式不变，对传入的视线射线的端点位置调整。

ro本来是相机位置camPos。

缩放：对camPos缩放，越小，说明我们越接近单位圆去ray这个Julia，Julia越大。
移动:camPos-JuliaPos

步骤

1.计算ro

注意到RayMarchFunc也被连了进去。它以*1的形式参与ro的计算，而主函数以ro为参，保证了其usf内的函数先于主函数编译。

//RayMarchFunc.usfreturn 1;
}float4 qsqr( float4 a ) // square a quaterion
{return float4( a.x*a.x - a.y*a.y - a.z*a.z - a.w*a.w,2.0*a.x*a.y,2.0*a.x*a.z,2.0*a.x*a.w );
}float map(float3 p,float4 c )
{float4 z = float4(p,0.0);float md2 = 1.0;float mz2 = dot(z,z);float4 trap = float4(abs(z.xyz),dot(z,z));float n = 1.0;for( int i=0; i<11; i++ ){// dz -> 2·z·dz, meaning |dz| -> 2·|z|·|dz|// Now we take thr 2.0 out of the loop and do it at the end with an exp2md2 *= mz2;// z  -> z^2 + cz = qsqr(z) + c;  trap = min( trap, float4(abs(z.xyz),dot(z,z)) );mz2 = dot(z,z);if(mz2>4.0) break;n += 1.0;}return 0.25*sqrt(mz2/md2)*exp2(-n)*log(mz2);  // d = 0.5·|z|·log|z| / |dz|
}bool intersectJulia(float3 ro, float3 rd ,float4 c)
{float4 res;float resT = -1.0;float maxd = 10.0;float h = 1.0;float t = 0.0;for( int i=0; i<300; i++ ){if( h<0.0001||t>maxd ) break;h = map( ro+rd*t, c );t += h;}if( t<maxd ) { resT=t;return true; }return false;

2.计算rd

和上篇一样

3.主函数

其中c就是julia递推式z^2+a的a。

4.结合蓝图Actor

控制材质参数，包括Julia位置等。

5.结果

在ue4默认相机下，正常工作。

我注意到，当我使用ue4影视镜头时，由于其FOV并非90度，所以移动起来Julia的位置和大小不正常。这和上面的推导一致。

我还把Inigo Quilez计算法线的部分移植进来。这部分代码对我来说水平太高且不具普遍性，放弃理解了。先只写了一个diffuse，看起来不是很好看。对比最近的分形宣传视频，果然分形好看还是在光照和材质上。

当然，IQ大神自己是能信手拈来好看的光照模型的，并且作了AA。这是shaderToy原图：

但他的代码太过硬核，不适合学习。所以我决定找找其他代码清晰，而各方面都考虑得比较全面得Shader学习一下。

这是他的光照部分的代码：

// sky
{
float co = clamp( dot(-rd,nor), 0.0, 1.0 );
vec3 ref = reflect( rd, nor );
//float sha = softshadow( pos+0.0005*nor, ref, 0.001, 4.0, c );
float sha = occ;
sha *= smoothstep( -0.1, 0.1, ref.y );
float fre = 0.1 + 0.9*pow(1.0-co,5.0);col  = mate*0.3*vec3(0.8,0.9,1.0)*(0.6+0.4*nor.y)*occ;
col +=  2.0*0.3*vec3(0.8,0.9,1.0)*(0.6+0.4*nor.y)*sha*fre;
}// sun
{
const vec3 lig = sun;
float dif = clamp( dot( lig, nor ), 0.0, 1.0 );
float sha = softshadow( pos, lig, 0.001, 64.0, c );
vec3 hal = normalize( -rd+lig );
float co = clamp( dot(hal,lig), 0.0, 1.0 );
float fre = 0.04 + 0.96*pow(1.0-co,5.0);
float spe = pow(clamp(dot(hal,nor), 0.0, 1.0 ), 32.0 );
col += mate*3.5*vec3(1.00,0.90,0.70)*dif*sha;
col +=  7.0*3.5*vec3(1.00,0.90,0.70)*spe*dif*sha*fre;
}// extra fill
{
const vec3 lig = vec3( -0.707, 0.000, -0.707 );
float dif = clamp(0.5+0.5*dot(lig,nor), 0.0, 1.0 );
col += mate* 1.5*vec3(0.14,0.14,0.14)*dif*occ;
}// fake SSS
{
float fre = clamp( 1.+dot(rd,nor), 0.0, 1.0 );
col += mate* mate*0.6*fre*fre*(0.2+0.8*occ);
}

结语

之后参考比较规整的光照模型，将Julia Set表面整好看点...

1.技术差，就会丑。

2.在非学术的Shader实现里，我经常看到Inigo Quilez。我最佩服的就是他独立思考的能力。以参考中的那篇ShaderToy所附带的技术博客为例，他抛开常见的z=a+bi+cj+dk形式的复数，自创z=a+bi+cj形式复数，并规定基向量的运算规则，十分大胆。针对实际问题，更改思考框架的操作可以说是我的梦想。

3.在Ray Tracing Quaternion Julia Sets on the GPU中，作者十分风趣，各个细节都十分详细（就像把你当傻瓜一样），甚至充满半页纸的图，结尾直接上代码，你不理解都不行。与满篇废话又不着重点，吹上天又没东西的gp论文比，我太喜欢了。

作者的伊利诺伊大学厄巴纳-香槟分校似乎挺牛的。