GLSL MinGW 域 macOS OpenGL C 中级 Linux

GPU 上的域着色方法

Khaldoon Ghanem

4.86/5 (19投票s)

2012年11月11日

MIT

8分钟阅读

55556

1618

本文介绍了如何使用GPU上的域着色法来可视化单个复变量的复值函数。

Sample Image

引言

为了可视化单个复变量的复函数，我们需要一个四维空间！然而，通过将复平面上的值编码为颜色，我们能够在二维空间中可视化这些函数。这种方法被称为域着色法。

在本文中，我们简要介绍了域着色法。然后，我们解释了如何在GPU上实现它，以及解析函数公式和在复平面中导航等细节。

除了这个程序在可视化复函数方面的实用性之外，它还提供了一个简单而富有启发性的使用着色器的示例。实际上，这个程序是我第一次接触着色器，但结果对我来说非常有趣，所以我认为值得分享。

在开始之前，我们先列出程序的功能，这样我们就知道我们的目标是什么。程序应该允许以下功能：

使用常见的数学表达式指定函数的公式。
指定多个函数并它们之间进行切换。
在复平面中进行缩放和平移。
在不同的颜色映射之间切换

域着色法

为了可视化函数 $f: \mathbb{C} \rightarrow \mathbb{C}: w = f(z)$

用某种颜色映射或方案覆盖复平面。即为每个点 $w$ 赋予颜色。
对于定义域复平面上的每一点 $z$ ，计算 $w = f(z)$ ，然后用 $w$ 的颜色来着色 $z$ 。

举例说明

作为一个简单的例子，假设我们要可视化 $z^2$ 函数在四个点 $1, -1, i, -i$ 上的值。我们给值域平面上的这些点中的每一个赋予一个独特的颜色。然后我们计算这些点的函数值，并根据结果着色： $f(1)=f(-1)=1$ （绿色）和 $f(i)=f(-i)=-1$ （黄色）。右边的图像称为颜色映射，左边的图像是 $z^2$ 的绘图。
有关该方法和解释所得图表的详细说明，请参阅这篇非常好的文章：使用域着色法可视化复解析函数

在GPU上实现

复平面域的着色是“尴尬并行”的，因为平面上的每个点都是独立处理的。为了在此问题中利用GPU，我们将复数点映射到屏幕像素，并使用片段着色器对其进行着色。流程大纲如下：

使用图形API创建一个覆盖整个屏幕的矩形。矩形的片段将由图形管线的固定功能生成。这些片段就是屏幕的最终像素，因为只有一个图元，并且它没有被裁剪。
片段着色器为每个片段执行。在片段着色器中：
- 使用从主程序传递过来的变换矩阵，将片段的屏幕坐标转换为复数点 $z$ 。
- 在该点计算函数 $f(z) = w$ 。
- 根据结果点 $w$ 处的颜色映射确定片段的颜色。

遵循这些步骤，片段着色器的main函数在GLSL中看起来会是这样：

out vec4 gl_FragColor;

uniform mat3 windowToComplex;
vec2 f(vec2 z);
vec4 complexToColor(vec2 w);            

void main()
{

    vec2 z = vec2(windowToComplex * vec3(gl_FragCoord.x, gl_FragCoord.y, 1.0));

    vec2 w = f(z);

    gl_FragColor = complexToColor(w);
}

关于代码的一些说明

vec2 数据类型用于表示复数。
windowToComplex 是一个3x3矩阵，用于将片段（像素）的坐标转换为复平面。它作为uniform从应用程序传递给着色器，并在用户在复平面中导航时更改。
vec2 f(vec2 z) 是我们要可视化的复函数。该函数的实现不是固定的，它将根据用户输入动态构建。
vec4 complexToColor(vec2 w) 根据颜色映射返回与复数 $w$ 相关的颜色。

以下三个小节分别讨论windowToComplex、f 和 complexToColor。

复平面导航

假设窗口的宽度为w，高度为h，我们希望看到以x0+iy0为中心、最小维度为side的复平面区域（另一个维度将由窗口的纵横比决定），那么窗口和感兴趣区域之间的变换矩阵windowToComplex是什么？

从上图可以看出，变换可以分三步完成：

将窗口中心移到0,0，即平移-w/2, -h/2。
缩放窗口，使其最小维度等于side，即如果纵横比小于1则缩放side/w，如果纵横比大于1则缩放side/h。
将窗口中心与区域中心对齐，即平移x0, y0。

因此，最终的变换矩阵就是对应于每个步骤的矩阵的乘积（当然，会用齐次坐标表示以允许平移）。构建好矩阵后，我们将其作为一个uniform变量发送到着色器。

//Caution: This code is inside application not the shader!
void updateView()
{
    float ratio = w / h;

    if (ratio < 1)
        scale = side / w;
    else
        scale = side / h;

    m[0] = scale;
    m[1] = 0;
    m[2] = scale * (-w / 2) + x0;
    m[3] = 0;
    m[4] = scale;
    m[5] = scale * (-h / 2) + y0;
    m[6] = 0;
    m[7] = 0;
    m[8] = 1;

    //Send matrix to shader
    glUniformMatrix3fv(loc, GL_TRUE, 1, m);
}

显然，平移只需操作x0和y0即可完成，而缩放则是通过操作side来完成的。

计算复函数及其公式解析

首先，由于GLSL不支持复数，我们实现复数的常见运算：

float cAbs(vec2 z)
{
    return sqrt(z.x * z.x + z.y * z.y);
}

float cArg(vec2 z)
{
    return atan(z.y, z.x);
}

vec2 cMul(vec2 z1, vec2 z2)
{
    return vec2(z1.x * z2.x - z1.y * z2.y, z1.x * z2.y + z1.y * z2.x);
}

vec2 cDiv(vec2 z1, vec2 z2)
{
    if (z2.y == 0)
    {
        return vec2(z1.x / z2.x, z1.y / z2.x);
    }
    if (z2.x == 0)
    {
        return vec2(z1.y / z2.y, -(z1.x / z2.y));
    }

    float r = z2.x * z2.x + z2.y * z2.y;
    return vec2((z1.x * z2.x + z1.y * z2.y) / r,
            (z1.y * z2.x - z1.x * z2.y) / r);
}

vec2 cPow(vec2 z, int k)
{
    vec2 res = vec2(1.0, 0);
    if (k >= 0)
    {
        for (; k > 0; k--)
            res = cMul(res, z);
    }
    else
    {
        for (; k < 0; k++)
            res = cDiv(res, z);
    }
    return res;
}

vec2 cSin(vec2 z)
{
    return vec2(sin(z.x) * cosh(z.y), cos(z.x) * sinh(z.y));
}

vec2 cCos(vec2 z)
{
    return vec2(cos(z.x) * cosh(z.y), sin(z.x) * sinh(z.y));
}

vec2 cTan(vec2 z)
{
    float cx = cos(z.x);
    float shy = sinh(z.y);
    float temp = cx * cx + shy * shy;
    return vec2((sin(z.x) * cx) / temp, (shy * cosh(z.y)) / temp);
}

vec2 cExp(vec2 z)
{
    return exp(z.x) * vec2(cos(z.y), sin(z.y));
}

vec2 cLog(vec2 z)
{
    return vec2(log(cAbs(z)), cArg(z));
}

如引言所述，我们希望使用常见的数学表达式来指定复函数。类似于“`(z-2)^2*(z+1-2i)*(z+2+2i)/z^3+exp(sin(z))`”这样的表达式，应该在着色器中被解析为“`cDiv(cMul(cMul(cPow((z - vec2(2.0, 0.0)), 2), (z + vec2(1.0, 0.0) - vec2(0.0,2.0))), (z + vec2(2.0, 0.0) + vec2(0.0,2.0))), cPow(z, 3))+cExp(cSin(z))`”。

现在，由于着色器是在运行时编译的，我们不需要函数的中间表示，我们的任务就简化为将输入字符串转换为相应的GLSL代码，并将其编译为函数vec2 f(vec2 z)的实现。

为了执行转换，我们使用了lex和yacc工具。用于转换的语法规则和动作在文件paraser.l和parser.y中。该语法允许用户以如下方式指定复函数：

数字仅使用十进制表示法。
字母“z”和“Z”用于表示复变量。
字母“i”可以加在任何数字前面，使其成为虚数。例如：5i，-3.5i，i，0.04i。
运算符“+”、“-”、“/”、“*”的含义、结合性和优先级与c/c++语言相同。
运算符“^”是幂运算符，左结合，优先级最高。只允许整数（正数或负数）幂。
数字可以加在复变量前面表示乘法，例如：-4.5z。提示：这比-4.5*z快，因为前者被翻译为-4.5*z，而后者被翻译为cMul(vec2(5.0, 0.0), z)。
可以使用函数“sin”、“cos”、“log”、“exp”。

颜色映射

目标是为每个复数赋予一种颜色。原则上，您可以使用任何图像来覆盖复平面。然而，由此产生的函数图通常不容易解释。此外，图像无法覆盖整个无限平面，因此我们需要一种更系统的方法。

对于复数 $w$ ，我们将其辐角 $\arg(w)$ 编码为从平滑颜色序列（颜色渐变）中选择的色调。同时，我们将其绝对值 $|w|$ 编码为亮度。

您可以选择不同的颜色渐变来编码辐角。上图显示了两种常见的渐变。第一种在三种颜色之间插值：黑色、红色和黄色，而第二种则跨越了整个光谱。第一种渐变有助于确定复函数的零点和极点的阶。

对于亮度，一种选择是将零绝对值赋为纯黑色，将大于阈值的绝对值赋为纯白色，并在两者之间进行插值。另一种选择是使用绝对值对数 $\log_2|w|$ 的小数部分。取对数的好处是可以在极点附近平衡爆炸效应，同时取小数部分可以使无限的绝对值与介于零和一之间的亮度值之间产生周期性映射。

我们使用上述两种辐角编码颜色渐变，结合对数亮度方法，实现了两种颜色映射。

// First color map       
float PI = atan(1.0)*4.0;
float cAbs(vec2 z);
float cArg(vec2 z);
vec4 complexToColor(vec2 w)
{
    //Compute color hue
    float phi = cArg(w);
    if (phi < 0.0)
        phi += 2.0 * PI; // Restrict to interval [0,2PI]
    phi = degrees(phi);
    vec4 hue = vec4(0.0);
    vec4 c1 = vec4(0.0, 0.0, 0.0, 1.0); //Black
    vec4 c2 = vec4(1.0, 0.0, 0.0, 1.0); //Red
    vec4 c3 = vec4(1.0, 1.0, 0.0, 1.0); //Yellow
    //In the upper half of the plane, interploate between black and red
    if (phi >= 0.0 && phi < 180.0)
    {
        float s = (phi) / 180.0;
        hue = c2 * s + (1.0 - s) * c1;
    }
    //In the lower half of the plane, interploate between red and yellow
    else if (phi >= 180.0 && phi < 360.0)
    {
        float s = (phi - 180.0) / 180.0;
        hue = c3 * s + (1.0 - s) * c2;
    }
    
    //Compute brightness
    float r = cAbs(w);
    float brightness = fract(log2(r));
        
    return brightness * hue;
}

// Second color map         
vec4 HSVtoRGB(float h, float s, float v)
{
    //Convert between the HSV and RGB color model.
    //Taken from http://www.cs.rit.edu/~ncs/color/t_convert.html and rewritten for GLSL    
    int i;
    float f, p, q, t;
    vec4 RGB;

    if (s == 0.0)
    {
        // achromatic (grey)
        RGB.x = RGB.y = RGB.z = v;
        RGB.w = 1.0;
        return RGB;
    }

    h /= 60.0; // sector 0 to 5
    i = int(floor(h));
    f = h - float(i); // factorial part of h
    p = v * (1.0 - s);
    q = v * (1.0 - s * f);
    t = v * (1.0 - s * (1.0 - f));

   if(i==0)
    {
        RGB.x = v;
        RGB.y = t;
        RGB.z = p;
    }
    else if(i==1)
    {
        RGB.x = q;
        RGB.y = v;
        RGB.z = p;
    }
    else if(i==2)
    {
        RGB.x = p;
        RGB.y = v;
        RGB.z = t;
    }
    else if(i==3)
    {
        RGB.x = p;
        RGB.y = q;
        RGB.z = v;
    }
    else if(i==4)
    {
        RGB.x = t;
        RGB.y = p;
        RGB.z = v;
    }
    else if(i==5)
    {
        RGB.x = v;
        RGB.y = p;
        RGB.z = q;
    }
    else
    {
        RGB.x = 1.0;
        RGB.y = 1.0;
        RGB.z = 1.0;
    }
    RGB.w = 1.0;
    return RGB;
}

float PI = atan(1.0)*4.0;
float cAbs(vec2 z);
float cArg(vec2 z);

vec4 complexToColor(vec2 w)
{
    
    float r = cAbs(w);
    float brightness = fract(log2(r));
    
    //Compute color hue on color wheel
    float phi = cArg(w);
    if (phi < 0.0)
        phi += 2.0 * PI;
    float hue = degrees(phi);
        
    return HSVtoRGB(hue, 1.0, brightness); //Full saturation (1.0) is used
}

在函数和颜色映射之间切换

在读取用户输入的函数列表后，我们将每个函数翻译成GLSL着色器。我们还加载了用于颜色映射的着色器。然后编译所有着色器。现在，要可视化特定的函数和特定的颜色映射，我们只需要创建一个新的程序（着色程序），附加适当的着色器，链接该程序并使用它。当然，每次更改程序时，我们需要重新获取uniform变量windowToComplex的位置。

编译

该程序使用c99编写。它应该可以在Linux、Mac OS X和Windows（使用MinGW）上编译并运行。
在Linux和Mac上，需要安装以下软件：Glut库、Glew库、lex工具、yacc工具。
在Windows上（使用MinGW），不需要安装任何东西，因为它们已经包含在“Win32”目录中。
要运行该程序，计算机应具备至少支持GLSL 1.2版本的GPU和驱动程序。
已知问题：在Windows上，鼠标滚轮不起作用。这是Glut库在Windows上的一个问题。

使用

程序从“stdin”逐行读取输入。
每一行指定一个要渲染的复函数。
如果您使用键盘输入，请记住按“Ctrl+D”结束输入。
可以使用命令行重定向运算符“|”和“<”将输入重定向到其他源。
例如，您可以使用以下命令从文件“./res/misc.data”读取函数：./coloring.exe < ./res/misc.data
一个有趣的用法是使用该程序来观察泰勒级数如何随着项数的增加而收敛到其对应的函数。文件exp_tylor.data列出了指数函数的前20项泰勒级数。

最后说明

这个项目是旨在利用GPU可视化不同类复函数的一个更大项目的一部分。

分形，如Mandelbrot集和Newton分形。
单个复变量的复函数（本文）。
定义在单位球面上的复函数，如球谐函数。
在3D空间中定义的复函数，使用Marching Tetrahedra方法。

欢迎访问项目主页，查看报告和代码。