碎片着色机上复杂Mobius变换的优化

Question

我正在开发自己的图形引擎，以呈现各种分形(例如，我的视频这里 )，我目前正在优化朱莉娅集匹配的代码(更多细节请参见这个问题和我的项目 )。在片段着色器中，我使用以下函数：

vec3 JuliaMatingLoop(dvec2 z)
{
    ...

    for (int k = some_n; k >= 0; --k)
    {
        // z = z^2
        z = dcproj(c_2(z));
        
        // Mobius Transformation: (ma[k] * z + mb[k]) / (mc[k] * z + md[k])
        z = dcproj(dc_div(cd_mult(ma[k], z) + mb[k], dc_mult(mc[k], z) + md[k]));
    }
    
    ...
}

在阅读了这之后，我意识到我正在进行代码中的Mobius转换，而且(从数学上讲)我可以使用矩阵来完成相同的操作。但是，a、b、c和d常数都是复数(在我的代码中表示为mak、mbk、mck和mdk )，而GLSL矩阵中的元素只包含实数(而不是vec2)。那么，对于我的问题:是否有一种方法来优化这些Mobius变换，使用GLSL中的矩阵？或者优化这部分代码的其他方法？

助手函数(我需要为这个部分使用双倍函数，所以我不能通过切换到使用浮点数来优化)：

// Squaring
dvec2 c_2(dvec2 c)
{
    return dvec2(c.x*c.x - c.y*c.y, 2*c.x*c.y);
}
// Multiplying
dvec2 dc_mult(dvec2 a, dvec2 b)
{
    return dvec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}
// Dividing
dvec2 dc_div(dvec2 a, dvec2 b)
{
    double x = b.x * b.x + b.y * b.y;
    return vec2((a.x * b.x + a.y * b.y) / x, (b.x * a.y - a.x * b.y) / x);
    
}
// Riemann Projecting
dvec2 dcproj(dvec2 c)
{
    if (!isinf(c.x) && !isinf(c.y) && !isnan(c.x) && !isnan(c.y))
        return c;
        
    return dvec2(infinity, 0);
}

Answer 1

我不确定这是否有帮助，但你可以用矩阵来做复杂的算术。

如果你把一个复数z看成是一个实的二向量，它的成分是Re(z)，那么Im(z)。

A*z + B ~ (Re(A)  -Im(A) ) * (Re(z)) + (Re(B))
          (Im(A)   Re(A) )   (Im(z))   (Im(B))

你当然想

(A*z + B) / (C*z + D)

如果你计算

A*z+b as (x)
         (y)

C*z+d as (x')
         (y')

那么你要找的答案是

inv( x' -y') * ( x)
   ( y' x' )   ( y)
i.e
    (1/(x'*x'+y'*y')) * (x'  y') * (x)
                        (-y' x')   (y)

但是，需要注意的一点是，在这些公式中，就像在代码中一样，除法的实现并不像以前那样可靠。问题在于评估b.x * b.x + b.y * b.y。这可能会溢出到无穷大，或者下流到0，即使除法的结果可能是相当合理的。一种常用的方法是史密斯的方法，如这里，如果你搜索“健壮的复杂除法”，你会发现最近的工作。这类事情通常无关紧要，但如果你要迭代到无穷大，它可能会产生影响。