基于WebGL2的曼德博查看器实现级数近似/微扰理论时图像失真

我基于WebGL2和JavaScript开发曼德博集合查看器，正尝试从基础层面实现级数近似功能。目前我根据迭代次数、所用参考点（当前取复平面视口中心）的不同，会得到异常、失真、不完整的渲染图像。
我的数学基础不算深厚，所以完全有可能我在某处的计算出了问题。我知道参考点选择、系数支持的无失真最大迭代次数等是独立的问题，但我目前只想先让基础概念跑通。现在我仅在原点附近(-0.5, 0.0)的参考点跳过1次迭代，我认为结果不该失真到这个程度。
如果有人能指出并解释代码中导致该问题的逻辑错误，我将不胜感激。我会附上所有相关代码，包括fragment着色器和JavaScript代码。

fragment.glsl代码 #

#version 300 es
precision highp float;

in vec2 uv;

uniform sampler2D palette;
uniform sampler2D orbit;

layout(std140) uniform blockData {
  vec2 center;
  float scale;
  float aspect;
  float radius;
  int iterations;
  int width;
  mat4x2 coefficients;
  vec2 zStart;
};

out vec4 fragColor;

ivec2 uvFromIndex(int i) {
  int col = i % width;
  int row = i / width;
  return ivec2(col, row);
}

vec2 cxMul(vec2 a, vec2 b) {
  return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

vec2 cxPow(vec2 a, float n) {
  float angle = atan(a.y, a.x);
  float r = length(a);
  float real = pow(r, n) * cos(n * angle);
  float im = pow(r, n) * sin(n * angle);
  return vec2(real, im);
}

float iterate() {
  vec2 dc = scale * vec2(aspect * (uv.x - 0.5), (uv.y - 0.5)) + center;
  vec2 z = zStart;
  vec2 dz = cxMul(coefficients[0], dc) + cxMul(coefficients[1], cxPow(dc, 2.)) + cxMul(coefficients[2], cxPow(dc, 3.)) + cxMul(coefficients[3], cxPow(dc, 4.));

  for(int i = 0; i < iterations; i++) {
    dz = cxMul(2. * z + dz, dz) + dc;
    if(length(dz) > radius)
      return float(i) + 1. - log(log(length(dz))) / log(2.);
    z = texelFetch(orbit, uvFromIndex(i), 0).rg;
  }
  return 0.;
}

void main() {
  float n = iterate();
  fragColor = texture(palette, vec2(n / float(iterations), 0.5));
}

app.ts代码 #

import PicoGL, { DrawCall } from "picogl";

import { Complex } from "./complex";
import palettes from "./palettes";
import shaders from "./shaders";
import UniformBufferWrapper from "./UniformBufferWrapper";
import { getReferenceOrbit, getSeriesCoefficients } from "./utils";

const digitCodeRegex = /^Digit(\d)$/;

const defaultValues = {
  center: [-0.5, 0] as [number, number],
  scale: 3,
  radius: 2,
  iterations: 256,
  width: 16,
};

export default (function main() {
  const canvas = document.querySelector("canvas") as HTMLCanvasElement;

  const app = PicoGL.createApp(canvas).resize(
    canvas.clientWidth,
    canvas.clientHeight
  );

  const quadPositions = app.createVertexBuffer(
    PicoGL.FLOAT,
    2,
    new Float32Array([
      -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0,
    ])
  );

  const quadArray = app
    .createVertexArray()
    .vertexAttributeBuffer(0, quadPositions);

  const uniformBuffer = new UniformBufferWrapper(
    app,
    [
      PicoGL.FLOAT_VEC2,
      PicoGL.FLOAT,
      PicoGL.FLOAT,
      PicoGL.FLOAT,
      PicoGL.INT,
      PicoGL.INT,
      PicoGL.FLOAT_MAT4x2,
      PicoGL.FLOAT_VEC2,
    ],
    [
      "center",
      "scale",
      "aspect",
      "radius",
      "iterations",
      "width",
      "coefficients",
      "zStart",
    ]
  );

  let drawCall: DrawCall;
  let image: HTMLImageElement;

  let z = new Complex(0);
  let coefficients = getSeriesCoefficients(
    z,
    new Complex(...defaultValues.center),
    // defaultValues.iterations
    1
  );
  console.log("coefficients:", coefficients);
  let zStart = z.toNumberArray();

  uniformBuffer
    .set("center", new Float32Array(defaultValues.center))
    .set("scale", new Float32Array([defaultValues.scale]))
    .set("aspect", new Float32Array([canvas.width / canvas.height]))
    .set("radius", new Float32Array([defaultValues.radius]))
    .set("iterations", new Int32Array([defaultValues.iterations]))
    .set("width", new Int32Array([defaultValues.width]))
    .set("coefficients", new Float32Array(coefficients))
    .set("zStart", new Float32Array(zStart))
    .update();

  let orbit = getReferenceOrbit(
    z,
    new Complex(...defaultValues.center),
    defaultValues.iterations
  );

  const orbitTexture = app.createTexture2D(
    orbit,
    defaultValues.width,
    defaultValues.width,
    {
      internalFormat: PicoGL.RG32F,
      wrapS: PicoGL.CLAMP_TO_EDGE,
      wrapT: PicoGL.CLAMP_TO_EDGE,
      magFilter: PicoGL.NEAREST,
      minFilter: PicoGL.NEAREST,
    }
  );

  app
    .createPrograms([shaders.vertex, shaders.fragment])
    .then(function ([program]) {
      image = new Image();

      image.onload = function () {
        const paletteTexture = app.createTexture2D(image, {
          magFilter: PicoGL.NEAREST,
          minFilter: PicoGL.NEAREST,
        });

        drawCall = app
          .createDrawCall(program, quadArray)
          .texture("palette", paletteTexture)
          .texture("orbit", orbitTexture)
          .uniformBlock("blockData", uniformBuffer.getBuffer());

        drawCall.draw();
      };

      image.src = palettes[0];
    });

  window.onresize = function () {
    app.resize(canvas.clientWidth, canvas.clientHeight);

    uniformBuffer
      .set("aspect", new Float32Array([canvas.width / canvas.height]))
      .update();

    drawCall.draw();
  };

  canvas.onclick = function (e) {
    updateScreenBounds(e, 0.8);
  };

  canvas.oncontextmenu = function (e) {
    e.preventDefault();
    updateScreenBounds(e, 1.25);
  };

  let result: RegExpExecArray | null;

  window.onkeydown = function (e) {
    if ((result = digitCodeRegex.exec(e.code)) !== null) {
      const index = +result[1];
      image.src = palettes[index];
    }
  };

  function updateScreenBounds(position: MousePosition, scaleFactor: number) {
    const [center, scale, aspect] = uniformBuffer.get(
      "center",
      "scale",
      "aspect"
    );

    const mouse = new Float32Array([
      aspect[0] * scale[0] * (position.x / canvas.width - 0.5),
      scale[0] * ((canvas.height - position.y) / canvas.height - 0.5),
    ]);

    const newCenter = [
      center[0] + mouse[0] * 0.5,
      center[1] + mouse[1] * 0.5,
    ] as [number, number];

    z = new Complex(0);
    coefficients = getSeriesCoefficients(
      z,
      new Complex(...newCenter),
      // defaultValues.iterations
      1
    );
    console.log("coefficients:", coefficients);
    zStart = z.toNumberArray();

    uniformBuffer
      .set("center", new Float32Array(newCenter))
      .set("scale", new Float32Array([scale[0] * scaleFactor]))
      .set("coefficients", new Float32Array(coefficients))
      .set("zStart", new Float32Array(zStart))
      .update();

    orbit = getReferenceOrbit(
      z,
      new Complex(...newCenter),
      defaultValues.iterations
    );

    orbitTexture.data(orbit);

    drawCall.draw();
  }
})();

type MousePosition = { x: number; y: number };

utils.ts代码 #

import { Complex } from "./complex";

export function getReferenceOrbit(
  z: Complex,
  c: Complex,
  iterations: number,
  radius = 2
) {
  const orbit = new Float32Array(iterations * 2);
  for (let i = 0; i < iterations * 2; i += 2) {
    iterate(z, c);
    if (z.abs().toNumber() > radius)
      console.warn(`${c} is not a valid reference point!`);
    orbit[i] = z.re.toNumber();
    orbit[i + 1] = z.im.toNumber();
  }
  return orbit;
}

export function getSeriesCoefficients(
  z: Complex,
  c: Complex,
  iterations: number,
  radius = 2
) {
  const coefficients = [
    new Complex(1),
    new Complex(0),
    new Complex(0),
    new Complex(0),
  ];
  for (let i = 0; i < iterations; i++) {
    iterate(z, c);
    if (z.abs().toNumber() > radius)
      console.warn(`${c} is not a valid reference point!`);
    const [A, B, C, D] = [...coefficients];
    coefficients[0] = A.mul(z.mul(2)).add(1);
    coefficients[1] = B.mul(z.mul(2)).add(A.mul(A));
    coefficients[2] = C.mul(z.mul(2)).add(B.mul(A.mul(2)));
    coefficients[3] = D.mul(z.mul(2))
      .add(C.mul(A.mul(2)))
      .add(B.mul(B));
  }
  return coefficients.flatMap((x) => x.toNumberArray());
}

function iterate(z: Complex, c: Complex) {
  const x = z.mul(z).add(c);
  z.re = x.re;
  z.im = x.im;
}

complex.ts代码 #

import Decimal from "decimal.js";

export namespace Complex {
  export type Constructor = typeof Complex;
  export type Value = [Decimal.Value, Decimal.Value?] | [Complex];
}

export class Complex {
  re: Decimal;
  im: Decimal;

  constructor(...n: Complex.Value) {
    if (n[0] instanceof Complex) {
      this.re = n[0].re;
      this.im = n[0].im;
    } else {
      this.re = new Decimal(n[0]);
      this.im = new Decimal(n[1] ?? 0);
    }
  }

  add(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(this.re.add(c.re), this.im.add(c.im));
  }

  sub(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(this.re.sub(c.re), this.im.sub(c.im));
  }

  mul(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(
      this.re.mul(c.re).sub(this.im.mul(c.im)),
      this.re.mul(c.im).add(this.im.mul(c.re))
    );
  }

  div(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(
      this.re
        .mul(c.re)
        .add(this.im.mul(c.im))
        .div(c.re.mul(c.re).add(c.im.mul(c.im))),
      this.im
        .mul(c.re)
        .sub(this.re.mul(c.im))
        .div(c.re.mul(c.re).add(c.im.mul(c.im)))
    );
  }

  abs() {
    return this.re.mul(this.re).add(this.im.mul(this.im)).sqrt();
  }

  toNumberArray() {
    return [this.re.toNumber(), this.im.toNumber()];
  }
}

核心错误1：dc变量定义错误 #

级数近似的输入是当前像素点和参考点的复坐标差值，你当前的dc计算多了+ center，变成了像素点的绝对复坐标，完全不符合级数近似的输入要求。
修复：把fragment着色器里的dc定义改为

vec2 dc = scale * vec2(aspect * (uv.x - 0.5), (uv.y - 0.5));

核心错误2：迭代逻辑错误 #

你当前的循环是先更新dz，再读取参考轨道的z值，导致每一步计算dz用的是上一步的z值，和系数计算的逻辑不匹配。另外逃逸判断的是z+dz的总长度，不是dz的长度，你现在直接判断dz的长度完全错误。
修复：调整循环逻辑：

for(int i = 0; i < iterations; i++) {
  // 先取当前步的参考z值
  z = texelFetch(orbit, uvFromIndex(i), 0).rg;
  // 计算当前z总长度，判断逃逸
  vec2 zTotal = z + dz;
  if(length(zTotal) > radius)
    return float(i) + 1. - log(log(length(zTotal))) / log(2.);
  // 更新dz
  dz = cxMul(2.0 * z + dz, dz) + dc;
}
// 迭代完还没逃逸的属于集合内，返回迭代次数而不是0
return float(iterations);

其他适配修改 #

1. 参考轨道生成时要重新初始化z为0，你当前传入的z是已经经过系数迭代的变量，会导致轨道和系数计算的迭代步数错位。
2. 你用的是4阶泰勒展开，仅适合视口范围很小的高倍缩放场景，默认scale=3的视口太大，高阶项误差会非常明显，测试时建议先把scale缩小到0.1以内验证逻辑正确性。

内容的提问来源于stack exchange，提问作者badross92