shithub: fractals

ref: 24cb1dcb30000df8f556468ef6b27a899e6a5db7
dir: /gen.c/

#include <u.h>
#include <libc.h>
#include <draw.h>
#include <memdraw.h>
#include "complex.h"
#include "gen.h"

int is∈(Complex *c);
double len(Complex *c)
{
	Complex r;
	r.r = c->r * c->r;
	r.r = r.r < 0.00001 && r.r > 0.00001 ? 0.00001 : r.r;
	r.i = c->i * c->i;
	r.i = r.i < 0.00001 && r.i > 0.00001 ? 0.00001 : r.i;
	double ret = r.r + r.i;
	ret = ret < 3200000. ? ret : 3200000.;
	return sqrt(ret);
};

float lerp(float a, float b, float s)
{
	return s*b + (1.-s)*a;
}

Memimage*
generate(int x, int y, SetParams params, DrawParams dparams)
{
	if (x < 8 || y < 8) {
		werrstr("invalid image size");
		return nil;
	}
	
	Rectangle r;
	r.min.x = r.min.y = 0;
	r.max.x = x;
	r.max.y = y;
	Memimage *img = allocmemimage(r, XRGB32);
	if (!img) {
		return nil;
	}
	
	for (int nx = 0; nx < x; nx++) {
		for (int ny = 0; ny < y; ny++) {
			int valid, red, green, blue;
			Complex c;
			
			float uvx = (float)nx/(float)x;
			float uvy = (float)ny/(float)y;
			uvx -= .5;
			uvy -= .5;
			uvx *= dparams.scale;
			uvy *= dparams.scale;
			uvx += dparams.offsetx;
			uvy += dparams.offsety;
			
			c.r = uvx;
			c.i = uvy;
			
			valid = calc(&c, c, params);
			if (!valid)
				goto err;
			
			if (is∈(&c)) {
				red = 0;
				green = 0;
				blue = (int)(len(&c) * 255.);
			} else {
				float fmix = (float)valid/(float)params.iterations;
				red = lerp(0., 255., pow(fmix, 0.5));
				green = lerp(255., 0., pow(fmix, 0.2)) * 0.5;
				blue = 0;
			}
			
			if (red > 255)
				red = 255;
			if (green > 255)
				green = 255;
			if (blue > 255)
				blue = 255;
			
			int arri = (ny*x + nx) * 4;
			img->data->bdata[arri+0] = blue;
			img->data->bdata[arri+1] = green;
			img->data->bdata[arri+2] = red;
			img->data->bdata[arri+3] = 255;
		}
	}
	
	return img;
	
err:
	freememimage(img);
	return nil;
}

int
calc(Complex *out, Complex c, SetParams params)
{
	Complex t, f;
	
	if (params.iterations < 1) {
		werrstr("iterations < 1");
		return 0;
	}
	
	if (params.julia) {
		/* calculate julia */
		t = c;
		f = params.c;
	} else {
		/* calculate mandelbrot */
		t = params.c;
		f = c;
	}
	
	int iter = -1;
	
	for (int i = 0; i < params.iterations; i++) {
		if (t.r >= 2. || t.i >= 2.) {
			iter = i;
			break;
		}
		t = cpow2(&t);
		t = cadd(&t, &f);
	}
	
	*out = t;
	return iter > 0 ? iter : 1;
}

int
∈(Complex c, SetParams params)
{
	Complex r;
	int valid = calc(&r, c, params);
	if (!valid)
		return -1;
	return is∈(&r);
}

int
is∈(Complex *c)
{
	return (c->r <= 2. && c->i <= 2.);
}