Elimination of points in the set

[nessalc]

In an attempt to optimize my own simple program, I need to eliminate points inside the set. The main cardioid and period-2 bulb tests are in there (for initial points), but I still wind up maxing the iteration count for just as many points (actually, slightly more) as I skip with that method. I saw something about "orbit traps" and thought that would be my ticket, but that seems to indicate a different method of rendering, rather than a way to eliminate points that lie within the set, and I haven't found any tutorial on implementing orbit traps (for that purpose) if they are what I need.

I've seen references here to "Pauldelbrot et al"'s method(s) of eliminating points in the set, but search as I might, I can't find the method being referred to. Any ideas?

[claude]

Here's what I do:

Code:

// quadratic Mandelbrot set distance estimator (exterior + interior)
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef unsigned int N;
typedef long double R;
typedef long double complex C;

R cabs2(C z) {
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

int wucleus(C *z0, C c, N period) {
  R eps = 1e-10;
  R er2 = 16;
  C z = *z0;
  for (N j = 0; j < 64; ++j) {
    C dz = 1.0;
    for (N k = 0; k < period; ++k) {
      dz = 2.0 * dz * z;
      z = z * z + c;
    }
    R z2 = cabs(z);
    if (! (z2 < er2)) {
      break;
    }
    z = *z0 - (z - *z0) / (dz - 1.0);
    R e = cabs(z - *z0);
    *z0 = z;
    if (e < eps) {
      return 1;
    }
  }
  return 0;
}

R interior_distance(C z, C c, N period, R pixel_size) {
  if (wucleus(&z, c, period)) {
    C dz = 1.0;
    C dzdz = 0.0;
    C dc = 0.0;
    C dcdz = 0.0;
    for (N j = 0; j < period; ++j) {
      dcdz = 2.0 * (z * dcdz + dz * dc);
      dc = 2.0 * z * dc + 1.0;
      dzdz = 2.0 * (dz * dz + z * dzdz);
      dz = 2.0 * z * dz;
      z = z * z + c;
    }
    return (1.0 - cabs2(dz)) / (cabs(dcdz + dzdz * dc / (1.0 - dz)) * pixel_size);
  }
  return -1.0;
}

int main(int argc, char **argv) {
  N width = atoi(argv[1]);
  N height = atoi(argv[2]);
  C c0 = atof(argv[3]) + atof(argv[4]) * I;
  R r = atof(argv[5]);
  N maxiters = atoi(argv[6]);
  R er2 = 65536.0 * 65536.0;
  R pixel_size = r * 2.0 / height;
  printf("P5\n%u %u\n255\n", width, height);
  unsigned char *image = calloc(1, width * height);
  #pragma omp parallel for
  for (N y = 0; y < height; ++y) {
    fprintf(stderr, "%ud\r", y);
    for (N x = 0; x < width; ++x) {
      R py = height/2.0 - (y + 0.5);
      R px = (x + 0.5) - width/2.0;
      C c = c0 + (px + py * I) * pixel_size;
      C z = 0.0;
      C d = 0.0;
      R minz2 = er2;
      R de = -1.0;
      for (N i = 1; i < maxiters; ++i) {
        d = 2.0 * z * d + 1.0;
        z = z * z + c;
        R z2 = cabs2(z);
        if (! (z2 < er2)) {
          de = sqrt(z2) * log(z2) / (cabs(d) * pixel_size);
          break;
        }
        if (z2 < minz2) {
          minz2 = z2;
          de = interior_distance(z, c, i, pixel_size);
          if (de >= 0.0) {
            break;
          }
        }
      }
      de = fmin(fmax(256.0 * tanh(de), 0.0), 255.0);
      image[y * width + x] = de;
    }
  }
  fwrite(image, width * height, 1, stdout);
  return 0;
}

For each likely candidate period (where |z_n| is smaller than any previous |z_n|), find a point ("wucleus", stupid name I chose...) in the periodic attractor with Newton's method - then you can use it to compute the interior distance estimate (if Newton's method succeeded) - if Newton's method doesn't succeed then it won't be interior to a component of that period.

EDIT: usage instructions:

Code:

gcc -std=c99 -Wall -pedantic -Wextra -O3 -march=native -fopenmp -o mandelbrot mandelbrot.c -lm
./mandelbrot 1920 1080 -0.75 0 1.25 1000 > output.pgm