为基于 WebGL2 的 Mandelbrot 查看器应用程序实现系列 approximation/perturbation 理论

Question

我已经开始使用 WebGL2 和 JavaScript 构建一个 mandelbrot 查看器应用程序，并且正在尝试在基本级别上实现级数逼近。我目前正在获取 weird/distored/incomplete 图像，具体取决于迭代次数和我使用的参考点（目前只是复平面中视口的中心）。

我的数学背景不深，所以我完全有可能（很明显）在某处搞砸了我的计算。我知道在选择参考点、在不扭曲图像的情况下可以用系数完成的迭代次数等方面存在不同的问题，但我只是想看看基本概念是否有效。目前我只跳过 1 次迭代，参考点靠近原点 (-0.5, 0.0)，我不认为结果应该如此扭曲。

如果有人能够看到并向我解释代码中的逻辑哪里出了问题导致了这个结果，我将不胜感激。我将尝试包含所有相关代码，包括片段着色器和 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()];
  }
}

Answer 1

我解决了我自己的问题。系数和当前“xn”(z) 值在 1 次迭代中不同步。重新阅读 this paper 帮助我意识到我正在跳过 x 的迭代，这会导致 delta n 计算失败。我现在对如何对屏幕中心以外的其他参考点应用相同的逻辑有了更好的了解（在计算 delta 0 时应该只需要添加 x0 到中心的距离）并且可以开始探索优化，例如选择好的参考点并检测故障，计算用系数迭代多少次等

在广泛的缩放级别插入 x sub 0 系数后，分形看起来很正常，并且当应用超出该系数时，在更深的缩放级别上变得越来越准确，就像它应该的那样。相关固定码：

fragment.glsl

// ...
float iterate() {
  vec2 d0 = scale * vec2(aspect * (uv.x - 0.5), (uv.y - 0.5)); // + center - x0;
  vec2 dn = cxMul(coefficients[0], d0) + cxMul(coefficients[1], cxPow(d0, 2.)) + cxMul(coefficients[2], cxPow(d0, 3.)) + cxMul(coefficients[3], cxPow(d0, 4.));
  vec2 xn;

  for(int i = 0; i < iterations; i++) {
    xn = texelFetch(orbit, uvFromIndex(i), 0).xy;
    dn = cxMul(2. * xn + dn, dn) + d0;
    if(length(dn) > radius)
      return float(i) + 1. - log(log(length(dn))) / log(2.);
  }
  return 0.;
}
// ...

app.ts

// ...
const defaultValues = {
  center: [-0.5, 0] as [number, number],
  scale: 2.5,
  radius: 2,
  iterations: 100,
  width: 10,
  coefficients: [1, 0, 0, 0, 0, 0, 0, 0],
};

(function main() {
// ...

为基于 WebGL2 的 Mandelbrot 查看器应用程序实现系列 approximation/perturbation 理论

Implementing series approximation/perturbation theory for WebGL2 based Mandelbrot viewer application

fractals

mandelbrot

webgl

typescript

webgl2