请教这份代码的数学！

Question

#define pi 3.141592653589793

// Wythoff symbol
vec4 wythoff_symbol = vec4(3,2,5,0);

// Construct spherical triangle from pqr
mat3 construct_tri(vec3 pqr){
    vec3 a = pi/pqr;
    vec3 c = cos(a); float sp=sin(a.x);
    vec3 l1=vec3(1,0,0);
    vec3 l2=vec3(-c.x, sp, 0);
    float x3=-c.z;
    float y3=-(c.y+c.x*c.z)/sp;
    float z3=sqrt(1.0-x3*x3-y3*y3);
    vec3 l3=vec3(x3,y3,z3);
    return mat3(l1,l2,l3);
}

// Cross product intersection
vec3 intersect(vec3 a, vec3 b){return normalize(cross(a,b));}

// Bisector of two great circles
vec3 bisect(vec3 l1, vec3 l2){return intersect(cross(l1,l2),0.5*(l1+l2));}

// Flip point into triangle (magic number 5 iterations)
void flip_into_tri(mat3 tri, inout vec3 x, out mat3 M){
    vec3 d=x*tri;
    M=mat3(1.0);
    for(int k=0;k<5;k++){
        if(min(d.x,min(d.y,d.z))>=0.0) break;
        for(int j=0;j<3;j++){
            if(d[j]<0.0){
                vec3 tj=tri[j];
                M=mat3(reflect(M[0],tj),reflect(M[1],tj),reflect(M[2],tj));
                x=reflect(x,tj);
                d=x*tri;
            }
        }
    }
}

// Minimal check_domain
void check_domain(mat3 tri, vec3 x, float type,
                  out vec3 tri_vert, out vec3 tri_region,
                  out vec3 face_normal, out vec3 edge){
    vec3 p0=intersect(tri[1],tri[2]);
    vec3 p1=intersect(tri[2],tri[0]);
    vec3 p2=intersect(tri[0],tri[1]);
    if(type==0.0) tri_vert=p2;
    else {
          vec3 l_b1=bisect(tri[2],tri[0]);
          if(type==1.0) tri_vert=intersect(l_b1,tri[1]);
          else{
              vec3 l_b2=bisect(tri[0],tri[1]);
              tri_vert=intersect(l_b1,l_b2);
          }
    }
    vec3 l_a0=intersect(tri_vert,tri[0]);
    vec3 l_a1=intersect(tri_vert,tri[1]);
    vec3 l_a2=intersect(tri_vert,tri[2]);
    
    float d0=dot(x,l_a0), d1=dot(x,l_a1), d2=dot(x,l_a2);
    
    if(d1<0.0 && d2>=0.0){
        tri_region=vec3(1,0,0); 
        edge=abs(d1)<abs(d2)?l_a1:l_a2;
    }
    else if(d2<0.0 && d0>=0.0){
        tri_region=vec3(0,1,0); 
        edge=abs(d2)<abs(d0)?l_a2:l_a0;
    }
    else{
        tri_region=vec3(0,0,1); 
        edge=abs(d0)<abs(d1)?l_a0:l_a1;
    }
    face_normal=mat3(p0,p1,p2)*tri_region;
}

// Wythoff wrapper
void wythoff(vec4 ws, vec3 x, out mat3 tri, out mat3 M,
             out vec3 tri_vert, out vec3 tri_region,
             out vec3 face_normal, out vec3 edge){
    vec3 pqr=ws.xyz; float type=ws.w;
    tri=construct_tri(pqr);
    flip_into_tri(tri,x,M);
    check_domain(tri,x,type,tri_vert,tri_region,face_normal,edge);
    face_normal=face_normal*M; tri_vert=tri_vert*M; edge=edge*M;
}

// Minimal distance function
float map(vec3 pos){
    vec3 x=normalize(pos);
    mat3 tri,M;
    vec3 tri_vert,tri_region,face_normal,edge;
    wythoff(wythoff_symbol,x,tri,M,tri_vert,tri_region,face_normal,edge);
    return dot(pos,face_normal)-dot(face_normal,tri_vert);
}

// Normal calc
vec3 calc_normal(vec3 p){vec3 eps=vec3(0.001,0,0);
    return normalize(vec3(
        map(p+eps.xyy)-map(p-eps.xyy),
        map(p+eps.yxy)-map(p-eps.yxy),
        map(p+eps.yyx)-map(p-eps.yyx)
    ));
}

// Raymarch
vec2 cast_ray(vec3 ro, vec3 rd){
    float t=0.0,h=8.0;
    for(int i=0;i<25;i++){
        h=map(ro+rd*t); t+=h; if(abs(h)<0.01) break;
    }
    return vec2(t,1.0);
}

// Shade
vec3 shade(vec3 ro, vec3 rd){
    vec2 tm=cast_ray(ro,rd);
    if(tm.x>5.) return vec3(0.3);
    vec3 p=ro+rd*tm.x; vec3 x=normalize(p);
    mat3 tri,M; vec3 tri_vert,tri_region,face_normal,edge;
    wythoff(wythoff_symbol,x,tri,M,tri_vert,tri_region,face_normal,edge);
    vec3 n=calc_normal(p);
    
    return (face_normal*0.5+0.5);
}

// Main
void mainImage(out vec4 fragColor, in vec2 fragCoord){
    vec2 uv=(fragCoord-iResolution.xy*0.5)/iResolution.y;
    float t=iTime;
    float r=3.;
    vec3 tgt=vec3(0), cpos=vec3(r*cos(t),1.,r*sin(t));
    vec3 up=vec3(0,1,0);
    vec3 rz=normalize(tgt-cpos), rx=normalize(cross(rz,up)), ry=cross(rx,rz);
    mat3 Rview=mat3(rx,ry,rz);
    vec3 rd=Rview*normalize(vec3(uv,1.));
    vec3 color=shade(cpos,rd);
    fragColor=vec4(color,1.0);
}

Answer 1

晚上好🌙🌙🌙
本答案参考通义千问

这段代码是一个基于 Wythoff 构造 的三维图形生成算法，用于创建 非欧几何结构（如双曲空间中的多面体）。它通过一系列数学运算和几何变换来实现图形的构建、渲染与光照计算。

---

一、代码中涉及的数学原理

1. Wythoff 符号

Wythoff 符号是用于描述 正多面体或更一般的对称结构 的一种符号，通常表示为 p q r，其中 p、q、r 是整数，代表三角形的边角对应的反射次数。

• 在本代码中，vec4 wythoff_symbol = vec4(3,2,5,0); 表示一个 Wythoff 符号为 (3,2,5,0) 的结构。
• 这种结构通常用于构造 球面三角形 或 双曲三角形，并利用反射对称性进行扩展。

2. 构造球面三角形

函数 construct_tri(vec3 pqr) 使用 三角形角度 来构造一个单位球面上的三角形。

数学公式：

• 将角度转换为弧度：a = pi / pqr
• 计算各边的角度余弦值：c = cos(a)
• 构造三个向量 l1, l2, l3，表示三角形的三个顶点方向：

vec3 l1 = vec3(1, 0, 0);
vec3 l2 = vec3(-c.x, sp, 0); // sp = sin(a.x)
vec3 l3 = vec3(x3, y3, z3);

其中 x3, y3, z3 是通过三角形的几何约束计算得到的坐标。

重点： 这个三角形是球面上的一个 正三角形，用于后续的反射对称操作。

---

3. 交点计算（cross product intersection）

函数 intersect(vec3 a, vec3 b) 返回两个向量的 叉积归一化 向量，这表示两向量所定义的平面的法线方向。

重点： 叉积可以用来找到两个向量之间的垂直方向，常用于几何中的法线计算。

---

4. 角平分线（bisect）

函数 bisect(vec3 l1, vec3 l2) 返回两个向量的角平分线方向，其方法是：

• 计算两个向量的 和 和 叉积。
• 然后求出它们的 交点。

重点： 角平分线用于确定三角形内部的对称轴，是 Wythoff 构造的重要步骤。

---

5. 翻转到三角形内（flip_into_tri）

函数 flip_into_tri(mat3 tri, inout vec3 x, out mat3 M) 用于将一个点 x 投影到三角形内部，通过多次反射操作。

数学原理：

• 如果点 x 不在三角形内部（即 d.x < 0 或其他坐标小于 0），则用三角形的边进行反射。
• 每次反射都会更新变换矩阵 M，以记录当前的反射状态。

重点： 这是一个典型的 反射群 操作，用于生成对称结构。

---

6. 检查区域（check_domain）

函数 check_domain(...) 用于判断点 x 所在的三角形区域，并返回：

• tri_vert：三角形的顶点。
• tri_region：点所在的区域（例如顶点、边、面）。
• face_normal：面的法线方向。
• edge：边的方向。

重点： 通过点与三角形边的夹角判断点的位置，这是图形学中常用的点-面关系判定方法。

---

7. Wythoff 构造封装函数

函数 wythoff(...) 封装了所有上述操作，用于生成最终的 几何结构，并通过 map(...) 函数计算距离场。

---

8. 距离场（map）

函数 map(vec3 pos) 用于计算点 pos 到目标结构的距离。

公式：

return dot(pos, face_normal) - dot(face_normal, tri_vert);

重点： 这是一个标准的 距离场计算方式，用于 raymarching 渲染。

---

9. 法线计算（calc_normal）

函数 calc_normal(...) 使用 有限差分法 计算点 p 处的法线向量。

重点： 法线是渲染光照效果的基础，通过微小位移计算梯度得到。

---

10. 光线追踪（raymarch）

函数 cast_ray(...) 使用 raymarching 算法，沿着光线方向逐步推进，直到遇到结构表面。

重点： Raymarching 是一种适用于隐式曲面（如距离场）的高效渲染技术。

---

11. 着色（shade）

函数 shade(...) 负责根据光线与结构的交互，计算颜色。

重点： 通过 face_normal 和 tri_region 来决定颜色，实现了基本的光照效果。

---

二、关键数学流程总结（有序列表）

1. 定义 Wythoff 符号

   • vec4 wythoff_symbol = vec4(3,2,5,0);
   • 用于指定对称结构的类型。

2. 构造球面三角形

   • 使用 construct_tri(...) 生成三角形的三个方向向量。

3. 计算交点和角平分线

   • 使用 intersect(...) 和 bisect(...) 来辅助几何分析。

4. 反射翻转到三角形内部

   • 通过 flip_into_tri(...) 将点 x 映射到三角形内部。

5. 区域检测与法线计算

   • 使用 check_domain(...) 判断点所在区域，并计算法线。

6. 生成距离场

   • 通过 map(...) 计算点到结构的距离。

7. 法线计算

   • 使用 calc_normal(...) 计算点法线。

8. 光线追踪（Raymarching）

   • 通过 cast_ray(...) 实现光线与结构的交互。

9. 着色

   • 根据法线和区域信息计算颜色。

---

三、修改建议（如有需要）

如果你希望该代码能够运行在 WebGL 2.0 中，请注意以下几点：

1. 使用 #version 300 es 替代 #define pi ...

   • 在 WebGL 中不能使用宏定义，应直接使用 const float pi = 3.141592653589793;

2. 避免使用 for(int k=0;k<5;k++)

   • 在 GLSL 中，for 循环必须使用常量表达式作为条件，否则会报错。
   • 修改为：

for(int k=0; k < 5; k++) {
    ...
}

1. 确保所有变量都初始化

   • 在 flip_into_tri(...) 中，M 必须被正确初始化，否则可能引发错误。

2. 完整代码片段修复建议（部分）

const float pi = 3.141592653589793;

// Wythoff symbol
vec4 wythoff_symbol = vec4(3.0, 2.0, 5.0, 0.0);

// Construct spherical triangle from pqr
mat3 construct_tri(vec3 pqr){
    vec3 a = pi/pqr;
    vec3 c = cos(a); 
    float sp = sin(a.x);
    vec3 l1 = vec3(1.0, 0.0, 0.0);
    vec3 l2 = vec3(-c.x, sp, 0.0);
    float x3 = -c.z;
    float y3 = -(c.y + c.x * c.z) / sp;
    float z3 = sqrt(1.0 - x3*x3 - y3*y3);
    vec3 l3 = vec3(x3, y3, z3);
    return mat3(l1, l2, l3);
}

// Cross product intersection
vec3 intersect(vec3 a, vec3 b){ return normalize(cross(a,b)); }

// Bisector of two great circles
vec3 bisect(vec3 l1, vec3 l2){ return intersect(cross(l1,l2), 0.5*(l1+l2)); }

// Flip point into triangle (magic number 5 iterations)
void flip_into_tri(mat3 tri, inout vec3 x, out mat3 M){
    vec3 d = x * tri;
    M = mat3(1.0);
    for(int k=0; k < 5; k++){
        if(min(d.x, min(d.y, d.z)) >= 0.0) break;
        for(int j=0; j < 3; j++){
            if(d[j] < 0.0){
                vec3 tj = tri[j];
                M = mat3(reflect(M[0], tj), reflect(M[1], tj), reflect(M[2], tj));
                x = reflect(x, tj);
                d = x * tri;
            }
        }
    }
}

---

四、总结

这段代码的核心思想是通过 Wythoff 构造 生成复杂的几何结构，并通过 raymarching 和 距离场 技术进行渲染。它结合了 反射对称性、球面几何、法线计算、光照模型 等多种图形学技术。

重点： 它是基于 非欧几何 的图形生成算法，适合用于生成双曲空间中的复杂结构。