root/trunk/dmpp/main.d

module main;
version(freetype) import SDL_ttf;
import std.math, std.string, std.stdio;
import tools.mersenne, tools.base;
import qd, SDL_ttf;

alias Mersenne Random;

int SDL_SaveBMP(SDL_Surface *surface, char []file) {
  return SDL_SaveBMP_RW(surface, SDL_RWFromFile(toStringz(file), "wb\0".ptr), 1);
}

struct rect(T, U) {
  T x, y;
  U w, h;
  string toString() { return format("rect (", x, ", ", y, ")[", w, ", ", h, "]"); }
}

bool maybe() { return ((rand&1)==0); }
T max(T)(T a, T b) { return (a>b?a:b); }

struct arguments {
  string[string] args;
  bool has(string s) { return (s in args)!=null; }
  string opIndex(string id) { return args[id]; }
  string need(string what) {
    if (!(what in args)) throw new Exception("Required parameter "~what~" not provided");
    return args[what];
  }
  void opCall(string[] args) {
    string last;
    void push(string what) { if (!last.length) return; this.args[last]=what; last=""; }
    foreach (arg; args) {
      if (!arg.length) continue;
      if (arg[0]=='-') { push(""); last=arg[1..$]; }
      else push(arg);
    }
    push("");
  }
}

struct point(T) { T x, y; }
struct pair(T) { T a, b; }

struct converter(T, U) {
  private rect!(T, T) origin;
  private rect!(U, U) target;
  private point!(real) ratio;
  void opCall(rect!(T,T) start, rect!(U,U) end) {
    origin=start; target=end;
    ratio.x=(target.w*1.0)/(origin.w*1.0);
    ratio.y=(target.h*1.0)/(origin.h*1.0);
  }
  complex!(U) toFn(point!(T) p) {
    complex!(U) np = void;
    np.re=cast(U)((p.x-origin.x)*ratio.x+target.x);
    np.im=cast(U)((p.y-origin.y)*ratio.y+target.y);
    return np;
  }
  point!(T) undo(complex!(U) p) {
    point!(T) up = void; /// undone point
    up.x=cast(T)((p.re-target.x)/ratio.x+origin.x);
    up.y=cast(T)((p.im-target.y)/ratio.y+origin.y);
    return up;
  }
  bool isIn(complex!(U) c) {
    static bool between(U v, U b, U t) { return ((v<b)&&(v>t)) || ((v<t)&&(v>b)); }
    return (between(c.re, target.x, target.x+target.w) && between(c.im, target.y, target.y+target.h));
  }
}

struct complex(T) {
  T re, im;
  static complex!(T) opCall(T re, T im) { complex!(T) result = void; result.re = re; result.im = im; return result; }
  complex!(T) opMul(complex!(T) b) {
    complex result = void;
    result.re=cast(T)(re*b.re-im*b.im);
    result.im=cast(T)(re*b.im+im*b.re);
    return result;
  }
  complex!(T) opAdd(complex!(T) b) {
    complex!(T) result = void;
    result.re=cast(T)(re+b.re); result.im=cast(T)(im+b.im);
    return result;
  }
  complex!(T) opDiv(complex!(T) b) {
    complex!(T) result = void;
    T sqs=cast(T)(b.re*b.re+b.im*b.im);
    result.re=cast(T)((re*b.re+im*b.im)/sqs); result.im=cast(T)((im*b.re-re*b.im)/sqs);
    return result;
  }
  void opAddAssign(complex!(T) b) {
    re=cast(T)(re+b.re); im=cast(T)(im+b.im);
  }
  void opMulAssign(complex!(T) b) {
    T tre=cast(T)(re*b.re-im*b.im);
    im=cast(T)(im*b.re+re*b.im);
    re=tre;
  }
}

alias float myval;
alias complex!(myval) mycomplex;

void swap(T)(inout T a, inout T b) { T c=a; a=b; b=c; }

struct FuncWrapper {
  myval a, b, c, d, e, f;
  void function(inout myval re, inout myval im) fn;
  myval[3] col;
  public myval p;
  static FuncWrapper opCall(myval _a, myval _b, myval _c, myval _d, myval _e, myval _f, void function(inout myval re, inout myval im) _fn, myval[] _col, myval _p) {
    FuncWrapper fw; with(fw) {
      fw.a=_a; fw.b=_b; fw.c=_c; fw.d=_d; fw.e=_e; fw.f=_f;
      fw.fn=_fn;
      col[]=_col[];
      p=_p;
    }
    return fw;
  }
  static FuncWrapper opCall(void function(inout myval re, inout myval im) _fn, uint seed, int factor=1) {
    auto rnd=new Random; rnd.seed(seed);
    myval getval() {
      auto v=rnd();
      myval res=(((v*1.0)/typeof(v).max));
      /*res=res>0.5?res*res:-1.0-res*res;
      res*=1.5;
      return res;*/
      return (res-0.5)*2.0;
    }
    myval[3] col; foreach (inout c; col) c=(rnd()*1.0)/(1.0*typeof(rnd()).max);
    static void normalize(inout myval r, inout myval g, inout myval b, myval to) {
      myval factor=to/(r*0.299+g*0.587+b*0.114);
      r*=factor; g*=factor; b*=factor;
    }
    static void boost(inout myval r, inout myval g, inout myval b, myval by) {
      static void _boost(inout myval lowest, inout myval highest, myval factor) {
        lowest=lowest*factor;
        highest=1-(1-highest)*factor;
      }
      /// w as in _walue_
      static ubyte max(myval[3] w) { ubyte m=0; if (w[1]>w[0]) { m=1; if (w[2]>w[1]) m=2; } else if (w[2]>w[0]) m=2; return m; }
      static ubyte min(myval[3] w) { ubyte m=0; if (w[1]<w[0]) { m=1; if (w[2]<w[1]) m=2; } else if (w[2]<w[0]) m=2; return m; }
      auto high=max([r,g,b]); auto low=min([r,g,b]);
      _boost(*[&r,&g,&b][min([r,g,b])], *[&r,&g,&b][max([r,g,b])], by);
    }
    normalize(col[0], col[1], col[2], 1.0);
    boost(col[0], col[1], col[2], 0.5);
    if (factor==-1) factor=(rnd()%2)+10;
    auto res=opCall(getval, getval, getval, getval, getval, getval, _fn, col, factor);
    return res;
  }
  void opCall(ref mycomplex p) {
    myval nre=a*p.re+b*p.im+c;
    myval nim=d*p.re+e*p.im+f;
    fn(nre, nim);
    p.re=nre; p.im=nim;
  }
}

void function(inout myval re, inout myval im) fns[];
char[][] fnames;
myval rsq(myval re, myval im) { return re*re+im*im; }
/// This is actually faster than doing a switch. I benchmarked it.
void setupFns() {
  fnames~="00-Linear";
  fns~=function void(inout myval re, inout myval im) {};
  
  fnames~="01-Sinusoidal";
  fns~=function void(inout myval re, inout myval im) { re=sin(re); im=sin(im); };
  
  fnames~="02-Spherical";
  fns~=function void(inout myval re, inout myval im) { myval rs=re*re+im*im; re/=rs; im/=rs; };
  
  fnames~="03-Swirl";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re, im)); myval a=atan2(im, re); a+=r; re=r*cos(a); im=r*sin(a); };
  
  fnames~="04-Horseshoe";
  fns~=function void(inout myval re, inout myval im)
  {
    myval a=atan2(re, im);
    myval c1=sin(a); myval c2=cos(a);
    myval nre=c1*re - c2*im;
    im=c2*re + c1*im;
    re=nre;
  };
  
  fnames~="05-Polar";
  fns~=function void(inout myval re, inout myval im)
  { myval r=1.0/sqrt(rsq(re, im)); re=atan2(im, re)/PI; im=r-1.0; };
    
  fnames~="06-Handkerchief";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); re=r*sin(a+r); im=r*cos(a-r); };
  
  fnames~="07-Heart";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); a*=r; re=r*sin(a); im=-r*cos(a); };
  
  fnames~="08-Disc";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re, im)); myval a=atan2(im, re); a/=PI; r*=PI; re=a*sin(r); im=a*cos(r); };
  
  fnames~="09-Spiral";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); re=(cos(a)+sin(r))/r; im=(sin(a)-cos(r))/r; };
  
  fnames~="10-Hyperbolic";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); re=sin(a)/r; im=cos(a)*r; };
  
  fnames~="11-Diamond";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); re=sin(a)*cos(r); im=cos(a)*sin(r); };
  
  fnames~="12-Ex";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); re=r*pow(sin(a+r), 3); im=r*pow(cos(a-r), 3); };
  
  fnames~="13-Julia";
  fns~=function void(inout myval re, inout myval im)
  { myval r=sqrt(rsq(re,im)); myval a=atan2(im, re); r=sqrt(r); a/=2.0; re=r*cos(a); im=r*sin(a); };
  
  writefln("Set up ", fns.length, " functions.");
}

int lookup(FuncWrapper f) {
  foreach (i, fn; fns) if (fn is f.fn) return i;
  writefln("Lookup failed. Throwing SIGFPE.");
  int e=0; e=1/e; throw new Exception("Lookup failed");
}

void logfn (FuncWrapper f) {
  writefln("Function ", fnames[lookup(f)], ":\n  x = p.x * ", f.a, " + p.y * ", f.b, " + ", f.c,
    ";\n  y = p.x * ", f.d, " + p.y * ", f.e, " + ", f.f, ";"c);
  writefln("Color: R", f.col[0], " G", f.col[1], " B", f.col[2]);
}

void logfn(FuncWrapper[] fs) {
  foreach (fw; fs) logfn(fw);
}

import std.c.stdlib:malloc, free;
final class fractal {
  alias rect!(myval, myval) myrect;
  static struct pixel { myval[3] col; }
  pixel[] buffer; int w, h;
  FuncWrapper[] fws;
  myrect range;
  converter!(ushort, myval) conv;
  ulong itercount;
  void cleanup() {
    foreach (inout col; buffer) with(col) { foreach (inout c; col) c=0; }
    conv(rect!(ushort, ushort)(cast(ushort)0, cast(ushort)0, cast(ushort)w, cast(ushort)h), range); // Init
    itercount=0;
  }
  this(int _w, int _h, FuncWrapper[] _fws, myrect _range) {
    w=_w; h=_h;
    buffer=(cast(pixel*)malloc(w*h*pixel.sizeof))[0..w*h];
    fws=_fws; range=_range;
    logfn(fws);
    cleanup;
  }
  ~this() { free(buffer.ptr); buffer=null; }
  bool addpixel(mycomplex cp, myval[3] pixcol) {
    if (!conv.isIn(cp)) return false;
    auto pp=conv.undo(cp);
    auto index=pp.x*h+pp.y;
    // if ((index<0)||(index>=buffer.length)) *(cast(int*)null)=0;
    for (int i = 0; i < 3; ++i) buffer[index].col[i] += pixcol[i];
    return true;
  }
  void iterate(ulong count) {
    /// Prescale the probability values to 0..maxint range
    if (!fws.length) return;
    ushort[] probability; ushort running_max=0;
    float psum=0.0; foreach (fw; fws) psum+=fw.p;
    // any rationale for this? float rndscale=(ushort.max*0.5)/psum;
    float rndscale=ushort.max/psum;
    ushort probsum=0;
    foreach (fw; fws) {
      running_max+=cast(ushort)(fw.p*rndscale);
      probability~=running_max; probsum+=running_max;
    }
    assert(probability[$-1] > ushort.max * 0.9); // guard against overflow
    mycomplex n;
    void resetpoint() {
      //n.re=0; n.im=0;
      n.re=((rand()*1.0)/typeof(rand()).max)*2.0-1.0;
      n.im=((rand()*1.0)/typeof(rand()).max)*2.0-1.0;
      //for (int i=0; i<20; ++i) fws[rnd.next%fws.length](n);
    }
    resetpoint;
    myval[3] col; foreach (inout c; col) c=0;
    /// Initialize the point
    int e=0;
    while (count) {
      /// 500000 is close to 524.288, which is 2^19
      if (!(count&524287)) writefln(count, "...");
      auto r = rand(), prerandom = cast(ushort) r; r >>= 16;
      // things are done twice to make full use of the uint
      if (prerandom >= probability[$-1]) continue; // highly unlikely
      int choice=0; while (probability[choice]<prerandom) choice++;
      auto func = fws[choice]; func(n);
      for (int i = 0; i < 3; ++i) col[i] = (col[i] + func.col[i]) / 2f;
      if (addpixel(n, col)) itercount++;
      prerandom = cast(ushort) r;
      if (prerandom >= probability[$-1]) continue; // still highly unlikely
      choice=0; while (probability[choice]<prerandom) choice++;
      func = fws[choice]; func(n);
      for (int i = 0; i < 3; ++i) col[i] = (col[i] + func.col[i]) / 2f;
      if (addpixel(n, col)) itercount++;
      count -= 2;
    }
  }
  Area render(config c, Area target=null) {
    auto dg={
      myval ipp=(itercount*1.0)/((buffer.length)*1.0); // iterations per pixel
      size_t i=0;
      for (ushort x=0; x<w; x+=c.dispscale) { if (!(x&255)) writefln("Column ", x);
        for (ushort y=0; y<h; y+=c.dispscale) {
          //myval tf(myval v) { return pow(log(v*c.brightness+1.0), c.gamma); } // Transform function
          myval tf(myval v) { return pow(log(v+1) / log(sqrt(2f))*c.brightness, c.gamma); } // Transform function
          myval factor=1.0/ipp; // Scaling factor
          static ubyte cap(uint v, ubyte c) { if (v>c) return c; return cast(ubyte)v; } // Self-explanatory
          ubyte[3] fcol; // "f" colors
          foreach (ix, cv; buffer[i].col) fcol[ix]=cap(cast(uint)(tf(cv*factor)*/*bness**/256.0), 255);
          auto col=rgb(fcol);
          pset(x/c.dispscale, y/c.dispscale, (c.invert?col.invert:col));
          i+=c.dispscale;
        }
        i+=h*(c.dispscale-1);
      }
    };
    if (target) screen.With(target) = dg();
    else screen.With(c.res[0]/c.dispscale, c.res[1]/c.dispscale) = { dg(); target = display; };
    return target;
  }
}

FuncWrapper[] seedgen(uint seed, char[] argstr) {
  FuncWrapper[] fws;
  auto rand=new Random; rand.seed(seed);
  auto amount=rand()%6;
  if (argstr.length) {
    foreach (fract; split(argstr, ",")) {
      while (fract[0]==' ') fract=fract[1..$]; while (fract[$-1]==' ') fract=fract[0..$-1];
      auto params=split(fract, " "); assert(params.length==7, "Must be seven params to fractal");
      float[6] coeffs; foreach (i, str; params[1..$]) coeffs[i]=atof(str);
      auto fn=FuncWrapper(fns[cast(uint)atoi(params[0])], rand());
      with (fn) { a=coeffs[0]; b=coeffs[1]; c=coeffs[2]; d=coeffs[3]; e=coeffs[4]; f=coeffs[5]; }