root/trunk/scene/fractal.d

module scene.fractal;

import scene.common, base.math, base.math: bit_isnan;

version(GNU) {
  import gcc.builtins;
  alias __builtin_logf logf;
} else {
  alias log logf;
}

float ipow(int I)(float f) {
  version(GNU) {
    return __builtin_powif(f, I);
  } else {
    static if (I == 0) return 1;
    else static if (I == 1) return f;
    else static if ((I%2) == 1) return sqr(ipow!(I/2)(f))*f;
    else return sqr(ipow!(I/2)(f));
  }
}

string ctPowExpand(string text) {
  string buf;
  while (true) {
    int pos = text.ctFind("^");
    if (pos == -1) return buf ~ text;
    char var = text[pos-1];
    buf ~= text[0 .. pos-1];
    int num = pos+1;
    while (text[num] >= '0' && text[num] <= '9')
      num ++;
    buf ~= "z"~text[pos+1 .. num]~"."~var;
    text = text[num .. $];
  }
}

float fn(vec c, int iters = 20) {
  auto z = c;
  // vec dt = vec(1, 0, 0);
  FP r = z.length, dl = 1;
  // auto ÎŽ = Spherical(0, 0, 1);
  int i;
  for (; i < iters; ++i) {
    auto r7 = ipow!(7)(r);
    dl = dl*r7*8 + 1; // Will this work, I wonder?
    /**
      z.x = |z|^8 * cos(asin(z.z / |z| )*8) * cos(atan2(z.y, z.x)*8) + c.x
      z.y = |z|^8 * cos(asin(z.z / |z| )*8) * sin(atan2(z.y, z.x)*8) + c.y
      z.z = |z|^8 * sin(asin(z.z / |z| )*8) + c.y
      Maxima codegen GO
    **/
    // Funs
    {
      vec z2 = void, z4 = void, z6 = void, z8 = void, z10 = void, z12 = void, z14 = void, z16 = void;
      vec z3 = void, z5 = void, z7 = void, z9 = void, z11 = void, z13 = void, z15 = void;
      // OH NOES HER NUMERIC COPROCESSORS
      for (int k = 0; k < vec.real_size; ++k) {
        z2.real_field[k] = z.real_field[k] * z.real_field[k];
        z3.real_field[k] = z2.real_field[k] * z.real_field[k];
        z4.real_field[k] = z2.real_field[k] * z2.real_field[k];
        z5.real_field[k] = z4.real_field[k] * z.real_field[k];
        z6.real_field[k] = z2.real_field[k] * z4.real_field[k];
        z7.real_field[k] = z6.real_field[k] * z.real_field[k];
        z8.real_field[k] = z4.real_field[k] * z4.real_field[k];
      }
      for (int k = 0; k < 2; ++k) { // z only needed up to ^8
        z9.real_field[k] = z8.real_field[k] * z.real_field[k];
        z10.real_field[k] = z8.real_field[k] * z2.real_field[k];
        z11.real_field[k] = z10.real_field[k] * z.real_field[k];
        z12.real_field[k] = z8.real_field[k] * z4.real_field[k];
        z13.real_field[k] = z12.real_field[k] * z.real_field[k];
        z14.real_field[k] = z8.real_field[k] * z6.real_field[k];
        z15.real_field[k] = z14.real_field[k] * z.real_field[k];
        z16.real_field[k] = z8.real_field[k] * z8.real_field[k];
      }
      vec nz = void;
      // I wish I had smart goo.
      // trigsimp(trigexpand(sqrt(x*x+y*y+z*z)^8*cos(asin(z/sqrt(x*x+y*y+z*z))*8)*cos(atan2(y, x)*8)));
      nz.x = mixin(ctPowExpand("((y^8-28*x^2*y^6+70*x^4*y^4-28*x^6*y^2+x^8)*z^8+(-28*y^10+756*x^2*y^8-1176*x^4*y^6-1176*x^6*y^4+756*x^8*y^2-28*x^10)*z^6+
(70*y^12-1820*x^2*y^10+1050*x^4*y^8+5880*x^6*y^6+1050*x^8*y^4-1820*x^10*y^2+70*x^12)*z^4+
(-28*y^14+700*x^2*y^12+308*x^4*y^10-2772*x^6*y^8-2772*x^8*y^6+308*x^10*y^4+700*x^12*y^2-28*x^14)*z^2+y^16-24*x^2*y^14-36*x^4*y^12+88*x^6*y^10+198*x^8*y^8+88*x^10*y^6-36*x^12*y^4
-24*x^14*y^2+x^16)/(y^8+4*x^2*y^6+6*x^4*y^4+4*x^6*y^2+x^8)")) + c.x;
      // trigsimp(trigexpand(sqrt(x*x+y*y+z*z)^8*sin(asin(z/sqrt(x*x+y*y+z*z))*8)*cos(atan2(y, x)*8)));
      nz.y = mixin(ctPowExpand("-((8*z.x*y^7-56*x^3*y^5+56*x^5*y^3-8*x^7*z.y)*z^8+(-224*z.x*y^9+1344*x^3*y^7-1344*x^7*y^3+224*x^9*z.y)*z^6+
(560*z.x*y^11-2800*x^3*y^9-3360*x^5*y^7+3360*x^7*y^5+2800*x^9*y^3-560*x^11*z.y)*z^4+(-224*z.x*y^13+896*x^3*y^11+2464*x^5*y^9-2464*x^9*y^5-896*x^11*y^3+224*x^13*z.y)*z^2+8*z.x*y^15-
24*x^3*y^13-120*x^5*y^11-88*x^7*y^9+88*x^9*y^7+120*x^11*y^5+24*x^13*y^3-8*x^15*z.y)/(y^8+4*x^2*y^6+6*x^4*y^4+4*x^6*y^2+x^8)")) + c.y;
      // trigsimp(trigexpand(sqrt(x*x+y*y+z*z)^8*sin(asin(z/sqrt(x*x+y*y+z*z))*8)));
      nz.z = mixin(ctPowExpand("sqrt(y^2+x^2)*(-8*z^7+(56*y^2+56*x^2)*z^5
      +(-56*y^4-112*x^2*y^2-56*x^4)*z^3+(8*y^6+24*x^2*y^4+24*x^4*y^2+8*x^6)*z.z)")) + c.z;
      z = nz;
    }
    r = z.length;
    if (r > 16) {
      break;
    }
  }
  // logln(c, ": estimated distance ", 0.5 * r * logf(r) / dl, " at ", i);
  return 0.5 * r * logf(r) / dl;
}

final class Fractal : SceneObject {
  float ε; int iters; bool adapt;
  static SceneObject fractal(ref string s) {
    float ε = 0.00001; int iters = 20; bool adapt;
    mixin(MatchExpr("s: [([[auto$#adapt=true$] $ε$|$#adapt=false$][, $iters$])]"));
    return new Fractal(ε, adapt, iters);
  }
  static this() {
    parse_fns["fractal"] = &fractal;
    builtin["fractal"] = true;
    help["fractal"] = ": stuff";
  }
  this(float ε, bool adapt, int iters) { this.ε = ε; this.adapt = adapt; this.iters = iters; }
  Fractal dup() { return new Fractal(ε, adapt, iters); }
  struct FractBuffer {
    /*void next(ref Hit hit) {
    }*/
  }
  mixin ThreadSlabSet!(FractBuffer);
  vec getGrad(vec v, float localε) {
    auto p = fn(v);
    return vec(
      fn(v + vec(localε,0,0), iters) - p,
      fn(v + vec(0,localε,0), iters) - p,
      fn(v + vec(0,0,localε), iters) - p
    );
  }
  override {
    string toString() { return "fractal"; }
    void iterate(void delegate(ref SceneObject) dg) { }
    void init(ref void* threadp) {
      auto c = cache_init(threadp);
    }
    void setup(ref Hit hit, ref Ray ray, ref void* threadp) {
      auto c = cache(threadp);
      // hit a sphere centered around origin radius 1
      auto psq = ray.pos.lensq, dsd = ray.dir.lensq,
        b = (-ray.pos).dot(ray.dir) / dsd, disc = b*b + (2 - psq) / dsd;
      if (disc < 0) return;
      
      auto d = sqrt(disc), t1 = b - d, t2 = b + d;
      FP[2] bracket = void;
      bracket[0] = t1;
      bracket[1] = t2;
      if (bracket[1] <= HitTolerance) return;
      auto cur = bracket[0];
      if (cur < 0) cur = 0;
      int x;
      int max_steps = iters * 2; // arbitrary choice
      while (true) {
        ++x;
        auto dist = fn(ray.at(cur), iters);
        cur += dist;
        float localε = ε;
        /*logln("Local ε at ", cur, " => ", ray.at(cur), ": ", ray.radius(cur), " * ", localε,
          "; ray angle ", ray.angle,
          ", ray origin ", ray.cam_origin,
          ", distance ", (ray.cam_origin - ray.at(cur)).length
        );*/
        if (adapt)
          localε *= ray.radius(cur);
        
        if (localε !> float.epsilon)
          localε = float.epsilon; // just to be sure
        
        if (dist < localε) {
          auto norm = getGrad(ray.at(cur), localε).normalized;
          
          if (norm.dot(ray.dir) > 0 && max_steps--) {
            cur += localε - dist; // definitely move away at least this much
            continue; // we're pointed AWAY. Not stopping.
          }
          // don't produce a hit if we're clearly moving away
          if (max_steps) hit.set(cur, true, norm);
          return;
        }
        if (cur !< bracket[1]) return; // left again
      }
    }
    Bound genBound() {
      return Bound(vec(-2), vec(2)); // I hope
    }
  }
}