在片段着色器上优化复杂的莫比乌斯变换
Optimizing Complex Mobius Transformations on a Fragment Shader
我正在开发自己的图形引擎来渲染各种分形(如我的视频 here for example), and I'm currently working on optimizing my code for Julia Set matings (see this question and my project 了解更多详细信息)。在片段着色器中,我使用这个函数:
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]));
}
...
}
并且在阅读 this 之后,我意识到我在这段代码中进行莫比乌斯变换,并且(从数学上讲)我可以使用矩阵来完成相同的操作。但是,a、b、c 和 d 常数都是复数(表示为ma[k], mb[k], mc[k],和 md[k] 在我的代码中),而 GLSL 矩阵中的元素只包含实数(而不是 vec2)。所以我的问题是:有没有办法在 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);
}
我不确定这是否有帮助,但是是的,您可以通过矩阵进行复杂的算术运算。
如果将复数 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,即使除法的结果可能非常合理。一个常用的解决方法是 Smith 的方法,例如 here,如果您搜索 'robust complex division',您会找到更多最近的工作。通常这种事情无关紧要,但如果你迭代到无穷大,它可能会有所不同。
我正在开发自己的图形引擎来渲染各种分形(如我的视频 here for example), and I'm currently working on optimizing my code for Julia Set matings (see this question and my project 了解更多详细信息)。在片段着色器中,我使用这个函数:
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]));
}
...
}
并且在阅读 this 之后,我意识到我在这段代码中进行莫比乌斯变换,并且(从数学上讲)我可以使用矩阵来完成相同的操作。但是,a、b、c 和 d 常数都是复数(表示为ma[k], mb[k], mc[k],和 md[k] 在我的代码中),而 GLSL 矩阵中的元素只包含实数(而不是 vec2)。所以我的问题是:有没有办法在 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);
}
我不确定这是否有帮助,但是是的,您可以通过矩阵进行复杂的算术运算。
如果将复数 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,即使除法的结果可能非常合理。一个常用的解决方法是 Smith 的方法,例如 here,如果您搜索 'robust complex division',您会找到更多最近的工作。通常这种事情无关紧要,但如果你迭代到无穷大,它可能会有所不同。