mandelbrot

Commit: 4fd96a0677500f10d13679550fc5dc437851179d
Parent: 82c38a88350ed138682a761145f5df6adbae5c9b
Author: Randy Palamar
Date:   Sun, 26 May 2024 18:37:06 -0600

first pass attempt at approximating high depth zooms

This is implemented based on Perturbation Theory. The code has
some issues but will be expanded upon soon.

Diffstat:
M
frag.glsl
 | 
40
+++++++++++++++++++++++++++++++---------
M
main.c
 | 
124
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------

2 files changed, 139 insertions(+), 25 deletions(-)
diff --git a/frag.glsl b/frag.glsl
@@ -7,6 +7,9 @@ uniform uvec2 u_screen_dim;
 uniform vec2  u_top_left;
 uniform vec2  u_bottom_right;
 
+uniform vec2 u_z_n[300];
+uniform bool u_use_approx = false;
+
 const int   iterations    = 300;
 const float escape_radius = 4.0;
 
@@ -84,26 +87,45 @@ vec2 map_mandelbrot(vec2 v)
 	return vec2(u_top_left.x, u_bottom_right.y) + v * scale;
 }
 
+vec2 complex_mult(vec2 a, vec2 b)
+{
+	vec2 result;
+	result.x = a.x * a.x - b.y * b.y;
+	result.y = 2 * a.x * b.y;
+	return result;
+}
+
 void main()
 {
 	vec2 vf = gl_FragCoord.xy / u_screen_dim.xy;
-	vec2 xy0 = map_mandelbrot(vf);
+	vec2 d  = map_mandelbrot(vf);
 
 	int i;
+	vec2 z = u_use_approx ? u_z_n[0] : d;
+	vec2 e = vec2(0);
 	float xx = 0, yy = 0;
-	vec2 xy = xy0;
-	for (i = 0; i < iterations && xx + yy < escape_radius; i++) {
-		xx = xy.x * xy.x;
-		yy = xy.y * xy.y;
-		xy = vec2(xx - yy + xy0.x, 2 * xy.x * xy.y + xy0.y);
+
+	if (u_use_approx) {
+		for (i = 0; i < iterations && xx + yy < escape_radius; i++) {
+			xx = z.x * z.x;
+			yy = z.y * z.y;
+			e = complex_mult(2 * u_z_n[i], e) + complex_mult(e, e) + d;
+			z = u_z_n[i] + e;
+		}
+	} else {
+		for (i = 0; i < iterations && xx + yy < escape_radius; i++) {
+			xx = z.x * z.x;
+			yy = z.y * z.y;
+			z  = vec2(xx - yy + d.x, 2 * z.x * z.y + d.y);
+		}
 	}
 
 	/* extra iterations to reduce error in escape calculation */
 	/* fun value j = 5 */
 	for (int j = 0; j < 2; j++) {
-		xx = xy.x * xy.x;
-		yy = xy.y * xy.y;
-		xy = vec2(xx - yy + xy0.x, 2 * xy.x * xy.y + xy0.y);
+		xx = z.x * z.x;
+		yy = z.y * z.y;
+		z  = vec2(xx - yy + d.x, 2 * z.x * z.y + d.y);
 	}
 
 	float zmag = sqrt(xx + yy);
diff --git a/main.c b/main.c
@@ -29,9 +29,12 @@ typedef struct { u8 *beg, *end; } Arena;
 
 #define ABS(x)           ((x) <= 0 ? -(x) : (x))
 #define BETWEEN(x, l, h) ((x) >= (l) && (x) <= (h))
+#define ARRAY_COUNT(a)   (sizeof(a) / sizeof(*a))
 
 #define EPS 0.00001f
 
+#define MAX_ITERATIONS 300
+
 #include "util.c"
 
 #ifdef __unix__
@@ -50,13 +53,28 @@ static struct {
 	u32 vao, vbo;
 	i32 pid;
 	i32 height, width;
-	i32 u_screen_dim, u_top_left, u_bottom_right;
+	union {
+		struct {
+			i32 screen_dim, z_n;
+			i32 top_left, bottom_right;
+			i32 use_approx;
+		};
+		i32 E[5];
+	} uniforms;
 	v2  dP;
 	Rect boundary;
 	f32 zoom;
+	v2 *z_n;
 	Colour clear_colour;
 } g_glctx;
 
+static const char *uniform_names[] = {
+	"u_screen_dim",
+	"u_z_n",
+	"u_top_left",
+	"u_bottom_right",
+	"u_use_approx",
+};
 
 static Rect default_boundary = {
 	.top_left     = (v2){ .x = -2.5, .y =  1.5 },
@@ -74,6 +92,21 @@ move_rect(Rect r, v2 delta)
 }
 
 static v2
+rect_center(Rect r)
+{
+	return (v2) {
+		.x = r.bottom_right.x - r.top_left.x,
+		.y = r.top_left.y - r.bottom_right.y,
+	};
+}
+
+static f32
+magnitude_v2(v2 v)
+{
+	return v.x * v.x + v.y * v.y;
+}
+
+static v2
 sub_v2(v2 a, v2 b)
 {
 	return (v2){ .x = a.x - b.x, .y = a.y - b.y };
@@ -170,6 +203,19 @@ scroll_callback(GLFWwindow *win, f64 xdelta, f64 ydelta)
 }
 
 static void
+recalculate_z_n(void)
+{
+	v2 *z_n = g_glctx.z_n;
+	z_n[0] = rect_center(g_glctx.boundary);
+	for (u32 i = 1; i < MAX_ITERATIONS; i++) {
+		f32 xx = z_n[i - 1].x * z_n[i - 1].x;
+		f32 yy = z_n[i - 1].y * z_n[i - 1].y;
+		z_n[i].x = xx - yy + z_n[0].x;
+		z_n[i].y = 2 * z_n[i - 1].x * z_n[i - 1].y + z_n[0].y;
+	}
+}
+
+static void
 mouse_button_callback(GLFWwindow *win, i32 btn, i32 act, i32 mod)
 {
 #if 0
@@ -187,6 +233,10 @@ mouse_button_callback(GLFWwindow *win, i32 btn, i32 act, i32 mod)
 	if (btn == GLFW_MOUSE_BUTTON_RIGHT && act == GLFW_PRESS) {
 		g_glctx.boundary = default_boundary;
 		g_glctx.dP = (v2){0};
+		g_glctx.zoom = 1.0;
+		recalculate_z_n();
+		glUniform2fv(g_glctx.uniforms.z_n, MAX_ITERATIONS,
+		             (f32 *)g_glctx.z_n);
 	}
 }
 
@@ -234,7 +284,6 @@ spawn_window(void)
 	glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 6);
 	glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
 
-	g_glctx.u_screen_dim = -1;
 	g_glctx.window = glfwCreateWindow(g_glctx.width, g_glctx.height,
 	                                  "Mandelbrot Viewer", NULL, NULL);
 	if (g_glctx.window == NULL)
@@ -321,15 +370,45 @@ get_arena(void)
 static void
 validate_uniforms(void)
 {
-	g_glctx.u_screen_dim = glGetUniformLocation(g_glctx.pid, "u_screen_dim");
-	g_glctx.u_top_left   = glGetUniformLocation(g_glctx.pid, "u_top_left");
-	g_glctx.u_bottom_right   = glGetUniformLocation(g_glctx.pid, "u_bottom_right");
+	for (u32 i = 0; i < ARRAY_COUNT(g_glctx.uniforms.E); i++) {
+		i32 uid = glGetUniformLocation(g_glctx.pid, uniform_names[i]);
+		g_glctx.uniforms.E[i] = uid;
+	}
+}
+
+#if 0
+static v2
+map_screen_to_rect(Rect bounds, v2 screen_pos)
+{
+	ASSERT(BETWEEN(screen_pos.x, 0, 1.0) && BETWEEN(screen_pos.y, 0, 1.0));
+
+	v2 size   = {
+		.x = ABS(bounds.bottom_right.x - bounds.top_left.x),
+		.y = ABS(bounds.top_left.y - bounds.bottom_right.y),
+	};
+
+	return (v2){
+		.x = bounds.top_left.x + size.x * screen_pos.x,
+		.y = bounds.bottom_right.y + size.y * screen_pos.y,
+	};
+}
+
+static void
+print_motion(void)
+{
+	v2 dP  = g_glctx.dP;
+	v2 ddP = g_glctx.ddP;
+	printf("dP  = { .x = %0.02f, .y = %0.02f }\n", dP.x, dP.y);
+	printf("ddP = { .x = %0.02f, .y = %0.02f }\n", ddP.x, ddP.y);
 }
+#endif
 
 int
 main(void)
 {
 	Arena memory = get_arena();
+	g_glctx.z_n = alloc(&memory, v2, MAX_ITERATIONS);
+	g_glctx.zoom = 1.0;
 
 	g_glctx.boundary = default_boundary;
 
@@ -364,14 +443,15 @@ main(void)
 		f32 current_time = glfwGetTime();
 		f32 dt = current_time - last_time;
 		last_time = current_time;
-		if (++fcount > 1000) {
-			printf("FPS: %0.03f | dt = %0.03f [ms]\n"
-			       "bounds = { { %0.04f, %0.04f }, { %0.04f, %0.04f } }\n",
-			       1 / dt, dt * 1e3,
-			       g_glctx.boundary.top_left.x,
-			       g_glctx.boundary.top_left.y,
-			       g_glctx.boundary.bottom_right.x,
-			       g_glctx.boundary.bottom_right.y);
+
+		i32 ua = 1.0e-8 > magnitude_v2(sub_v2(g_glctx.boundary.top_left,
+		                                      g_glctx.boundary.bottom_right));
+
+		if (++fcount > 300) {
+			v2 bound_cent = rect_center(g_glctx.boundary);
+			printf("FPS: %0.03f | dt = %0.03f [ms] | approx = %d\n"
+			       "Center: <%f, %f>\n",
+			       1 / dt, dt * 1e3, ua, bound_cent.x, bound_cent.y);
 			fcount = 0;
 		}
 
@@ -398,9 +478,21 @@ main(void)
 		delta.y = v.y * dt;
 		g_glctx.boundary = move_rect(g_glctx.boundary, delta);
 
-		glUniform2fv(g_glctx.u_top_left, 1, (f32 *)&g_glctx.boundary.top_left);
-		glUniform2fv(g_glctx.u_bottom_right, 1, (f32 *)&g_glctx.boundary.bottom_right);
-		glUniform2ui(g_glctx.u_screen_dim, g_glctx.width, g_glctx.height);
+		if (ua) {
+			recalculate_z_n();
+			glUniform2fv(g_glctx.uniforms.z_n, MAX_ITERATIONS,
+			             (f32 *)g_glctx.z_n);
+		}
+
+		glUniform2fv(g_glctx.uniforms.top_left, 1,
+		             (f32 *)&g_glctx.boundary.top_left);
+		glUniform2fv(g_glctx.uniforms.bottom_right, 1,
+		             (f32 *)&g_glctx.boundary.bottom_right);
+		glUniform2ui(g_glctx.uniforms.screen_dim,
+		             g_glctx.width, g_glctx.height);
+
+		glUniform1i(g_glctx.uniforms.use_approx, ua);
+
 		clear_colour(g_glctx.clear_colour);
 		glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 		glfwSwapBuffers(g_glctx.window);