mandulia-0.5: src/rjulia.c

/*
Mandulia -- Mandelbrot/Julia explorer
Copyright (C) 2010  Claude Heiland-Allen <claudiusmaximus@goto10.org>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <math.h>
#include <stdlib.h>

#include "rjulia.h"

#ifndef PRECISION
#define PRECISION 1
#endif

#if PRECISION == 1
typedef float R;
#define cos cosf
#define sqrt sqrtf
#define log logf
#define log2 log2f
#define fmin fminf
#define fmax fmaxf
#define SCAN "%f"
#else
typedef double R;
#define SCAN "%lf"
#endif

static const int channels = 4;
static const R escapeRadius = 65536.0;
static const R escapeRadius2 = 65536.0 * 65536.0;
static const int escapeLimit = 64;

static const R sqrt5 = 2.23606797749979;
static const R sqrt6 = 2.449489742783178;
static const R sqrt7 = 2.6457513110645907;

struct point {
  unsigned char *pixel;
  R x;
  R y;
};

struct context {
  struct point *points;
  int width;
  int height;
};

struct context *julia_new(int width, int height) {
  if (width > 0 && height > 0) {
    struct context *ctx = calloc(1, sizeof(struct context));
    if (ctx) {
      ctx->points = calloc(width * height, sizeof(struct point));
      if (ctx->points) {
        ctx->width = width;
        ctx->height = height;
        return ctx;
      }
      free(ctx);
    }
  }
  return 0;
}

/*
void julia_delete(struct context *ctx) {
  if (ctx) {
    if (ctx->points) {
      free(ctx->points);
    }
    free(ctx);
  }
}
*/

static inline int min(int x, int y) {
  return x < y ? x : y;
}

static inline int max(int x, int y) {
  return x > y ? x : y;
}

static inline void colour(R v, unsigned char *r, unsigned char *g, unsigned char *b) {
  int rr = 128 * cos(sqrt5 * v) + 128;
  int gg = 128 * cos(sqrt6 * v) + 128;
  int bb = 128 * cos(sqrt7 * v) + 128;
  *r = min(max(rr, 0), 255);
  *g = min(max(gg, 0), 255);
  *b = min(max(bb, 0), 255);
}

static inline double mandelbrot1(double x, double y) {
  R px = 0, py = 0, cx = x, cy = y, px1, py1, px2, py2, pxy, d2;
  for (int i = 0; i < 16384; ++i) {
    // p_n+1 := p_n * p_n + c
    // d_n+1 := 2 * p_n * d_n + 1
    px2 = px * px;
    py2 = py * py;
    d2 = px2 + py2;
    if (d2 > escapeRadius2) {
      return i - log2(log2(d2)/log2(escapeRadius2));
    }
    pxy = px * py;
    px1 = px2 - py2 + cx;
    py1 = 2 * pxy + cy;
    px = px1;
    py = py1;
  }
  return -1;
}

void julia(struct context *ctx, unsigned char *image, double dcx, double dcy) {
  if (! ctx || ! image) {
    return;
  }
  struct point *points = ctx->points;
  const int width = ctx->width;
  const int height = ctx->height;
  const int stride = width * channels;
  const R cx = dcx;
  const R cy = dcy;
  const R k = 4;
  const R hw = height * k / width;
  const R wh = width  * k / height;
  const R sx = width < height ? k : wh;
  const R sy = width < height ? hw : k;
  const int ntotal = width * height;
  int npoints;
  { // initialize
    npoints = ntotal;
    struct point *p = points;
    for (int j = 0; j < height; ++j) {
      R y = j * sy / height - sy/2;
      for (int i = 0; i < width; ++i) {
        R x = i * sx / width - sx/2;
        p->pixel = image + j * stride + i * channels;
        p->x = x;
        p->y = y;
        p->pixel[0] = 0;
        p->pixel[1] = 0;
        p->pixel[2] = 0;
        p->pixel[3] = 255;
        ++p;
      }
    }
  }
  int ndone = 0;
  int n = 0;
  int escapes;
  R u = mandelbrot1(dcx, dcy);
  u = u > 0 ? u : 0;
  do {
    escapes = 0;
    { // iterate
      struct point *p = points;
      for (int k = 0; k < npoints; ++k) {
        R x = p->x;
        R y = p->y;
        int ok = 1;
        for (int e = 0; ok && (e < escapeLimit); ++e) {
          R zx2 = x * x - y * y;
          R zy2 = 2 * x * y;
          R zx = zx2 + cx;
          R zy = zy2 + cy;
          R zd2 = zx * zx + zy * zy;
          if (zd2 >= escapeRadius2) {
            R d = sqrt(zd2);
            // colourify pixel
            unsigned char *pixel = p->pixel;
            R v = (n+e) - log2(log2(d)/log2(escapeRadius));
            v /= 32;
            colour(u + v, &pixel[0], &pixel[1], &pixel[2]);
            R a = fmin(fmax(8 * (v - 0.1) * v, 0), 1);
            pixel[3] = a * 255;
            p->pixel = 0;
            ++escapes;
            ok = 0;
          } else {
            x = zx;
            y = zy;
          }
        }
        p->x = x;
        p->y = y;
        ++p;
      }
      n += escapeLimit;
    }
    { // compact memory
      struct point *src = points;
      struct point *dst = points;
      int ncopied = 0;
      for (int k = 0; k < npoints; ++k) {
        if (! (src->pixel)) {
          ++src;
          ++ndone;
        } else {
          unsigned char *pixel = src->pixel;
          R x = src->x;
          R y = src->y;
          dst->pixel = pixel;
          dst->x = x;
          dst->y = y;
          ++src;
          ++dst;
          ++ncopied;
        }
      }
      npoints = ncopied;
    }
  } while(escapes);
}