Compare commits

...

1 commit

Author SHA1 Message Date
Bruno BELANYI fdb0d3d55d WIP: old changes rebased on top 2020-12-18 16:35:24 +01:00
4 changed files with 120 additions and 43 deletions

View file

@ -17,8 +17,14 @@ sources = [
'src/ppm.c',
]
cc = meson.get_compiler('c')
deps = [
cc.find_library('m', required : false),
]
executable(
'buddhabrot',
sources : sources,
c_args : '-D_XOPEN_SOURCE',
c_args : '-D_XOPEN_SOURCE=500',
dependencies : deps,
)

View file

@ -12,12 +12,15 @@
#include "image.h"
struct data {
float min_x;
float max_x;
float min_y;
float max_y;
struct window {
float min_x;
float max_x;
float min_y;
float max_y;
} window;
size_t max_iter;
size_t repeats;
size_t h;
size_t w;
@ -37,32 +40,104 @@ static size_t data_index(size_t h, size_t w, const struct data *d) {
return h * d->w + w;
}
static float scale(float min, float max, size_t x, size_t width) {
return min + x * (max - min) / (width - 1);
}
static size_t inv_scale(float min, float max, float x, size_t width) {
return (x - min) / (max - min) * (width - 1);
}
static bool compute_orbit(struct my_complex z,
struct orbit *o, size_t max_iter) {
float r = 0.0;
float i = 0.0;
for (o->len = 0; o->len < max_iter; ++o->len) {
o->values[o->len].r = r;
o->values[o->len].i = i;
float xtemp = r * r - i * i + z.r;
i = 2 * r * i + z.i;
r = xtemp;
if (r * r + i * i >= 4)
return true; // Diverges
o->values[o->len].r = r;
o->values[o->len].i = i;
return false; // Diverges
}
return false;
return true;
}
static float orbit_probability(const struct window *w, const struct orbit *o) {
for (size_t i = 0; i < o->len; ++i) {
if (w->min_x <= o->values[i].r && o->values[i].r <= w->max_x
&& w->min_y <= o->values[i].i && o->values[i].i <= w->max_y)
return 1;
}
return FLT_MIN; // Cannot be 0 because of sampler
}
static float transition_probability(const struct orbit *src,
const struct orbit *dst, size_t max_iter) {
return (1 - (max_iter - src->len) / max_iter)
/ (1 - (max_iter - dst->len) / max_iter);
}
static bool mutate_near(const struct orbit *current,
struct orbit *new, size_t max_iter) {
struct my_complex next = current->values[1]; // Always the first one after 0
const double magnifictation = 1; // FIXME: proportional to magnification
const double r1 = (1.f / magnifictation) * 0.0001;
const double r2 = (1.f / magnifictation) * 0.1;
const double phi = M_PI * 2 * random() / LONG_MAX;
const double r = r2 * exp(-log(r2 / r1) * random() / LONG_MAX);
next.r += r * cos(phi);
next.i += r * sin(phi);
return compute_orbit(next, new, max_iter);
}
static bool mutate_any(struct orbit *new, size_t max_iter) {
struct my_complex next = {
2 - 4. * random() / LONG_MAX,
2 - 4. * random() / LONG_MAX,
};
while (next.r * next.r + next.i * next.i >= 4) {
next.r = 2 - 4. * random() / LONG_MAX;
next.i = 2 - 4. * random() / LONG_MAX;
}
return compute_orbit(next, new, max_iter);
}
static void mutate(const struct data *d, struct orbit *current) {
struct orbit *new =
calloc(1, d->max_iter * sizeof(*new->values) + sizeof(*new));
if (!new)
err(EXIT_FAILURE, "could not allocate orbit data");
// Only consider values that diverge
if (random() < LONG_MAX / 5)
while (!mutate_any(new, d->max_iter))
continue;
else
while (!mutate_near(current, new, d->max_iter))
continue;
const float prob_cur = orbit_probability(&d->window, current);
const float prob_new = orbit_probability(&d->window, new);
const float t_new_cur = transition_probability(new, current, d->max_iter);
const float t_cur_new = transition_probability(current, new, d->max_iter);
const float alpha = (prob_new * t_cur_new) / (prob_cur * t_new_cur);
if (alpha > (float)random() / LONG_MAX)
memcpy(current, new, d->max_iter * sizeof(*new->values) + sizeof(*new));
free(new);
}
static size_t inv_scale(float min, float max, float x, size_t width) {
return (x - min) / (max - min) * width;
}
static void compute_buddahbrot(struct data *d) {
@ -71,25 +146,19 @@ static void compute_buddahbrot(struct data *d) {
if (!o)
err(EXIT_FAILURE, "could not allocate orbit data");
for (size_t i = 0; i < d->h; ++i) {
for (size_t j = 0; j < d->w; ++j) {
struct my_complex z = {
scale(d->min_x, d->max_x, j, d->w),
scale(d->min_y, d->max_y, i, d->h),
};
// Only keep if it diverges
if (!compute_orbit(z, o, d->max_iter))
continue;
while (true)
if (mutate_any(o, d->max_iter))
break; // Initialize with correct value
for (size_t i = 0; i < o->len; ++i) {
size_t w = inv_scale(d->min_x, d->max_x,
o->values[i].r, d->w);
size_t h = inv_scale(d->min_y, d->max_y,
o->values[i].i, d->h);
if (w >= d->w || h >= d->h)
continue; // We've diverged, keep going just in case
d->buf[data_index(h, w, d)] += 1;
}
for (size_t i = 0; i < d->repeats; ++i) {
mutate(d, o); // Computes its orbit
for (size_t i = 0; i < o->len; ++i) {
size_t w = inv_scale(d->window.min_x, d->window.max_x,
o->values[i].r, d->w);
size_t h = inv_scale(d->window.min_y, d->window.max_y,
o->values[i].i, d->h);
d->buf[data_index(h, w, d)] += 1;
}
}
@ -108,7 +177,7 @@ static void fillout_image(const struct data *d, struct image *i) {
}
}
void buddhabrot(struct image *i, size_t max_iter) {
void buddhabrot(struct image *i, size_t max_iter, size_t repeats) {
struct data *d = calloc(1, sizeof(*d->buf) * i->w * i->h + sizeof(*d));
if (!d)
err(EXIT_FAILURE, "could not allocate buddhabrot data");
@ -116,12 +185,13 @@ void buddhabrot(struct image *i, size_t max_iter) {
d->h = i->h;
d->w = i->w;
d->max_iter = max_iter;
d->repeats = repeats;
// FIXME" make it user-definable
d->min_x = -2.5;
d->max_x = 1.0;
d->min_y = -1.0;
d->max_y = 1.0;
d->window.min_x = -2.5;
d->window.max_x = 1.0;
d->window.min_y = -1.0;
d->window.max_y = 1.0;
compute_buddahbrot(d);
fillout_image(d, i);

View file

@ -5,6 +5,6 @@
struct image; // Forward declaration
void buddhabrot(struct image *i, size_t max_iter);
void buddhabrot(struct image *i, size_t max_iter, size_t repeats);
#endif /* !BUDDHABROT_H */

View file

@ -4,6 +4,7 @@
#include "buddhabrot.h"
#include "image.h"
#include "buddhabrot.h"
#include "mandelbrot.h"
#include "options.h"
#include "ppm.h"
@ -17,7 +18,7 @@ int main(int argc, char *argv[]) {
switch (opt.render) {
case BUDDHABROT:
buddhabrot(image, opt.max_iter);
buddhabrot(image, opt.max_iter, 10000); // FIXME: allow user option
break;
case MANDELBROT:
mandelbrot(image, opt.max_iter);