/* "Experimental Power Stack Mandelbulb for Owen" * By: Chris M. Thomasson * Source: comp.lang.c digest, Dec 5th 2019. * Topic: A 3d volumentric renderer [sic] of my PowerBulb... */ #include #include #include #include #include #include #define CT_WIDTH 256 // width of plane #define CT_HEIGHT 256 // height of plane #define CT_FRAMES 256 // slices of the mbulb #define CT_ITERS 64 // iterations per pixel /* the canvas */ struct ct_rgb { unsigned char r, g, b; }; struct ct_canvas { unsigned long width; unsigned long height; struct ct_rgb* buf; }; static bool ct_canvas_create(struct ct_canvas* const self, unsigned long width, unsigned long height) { size_t size = width * height * sizeof(*self->buf); self->buf = calloc(1, size); if (self->buf) { self->width = width; self->height = height; return true; } return false; } static void ct_canvas_destroy( struct ct_canvas const* const self) { free(self->buf); } static bool ct_canvas_save_ppm(struct ct_canvas const* const self, char const* fname) { FILE* fout = fopen(fname, "wb"); if (fout) { char const ppm_head[] = "P6\n" "# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)"; fprintf(fout, "%s\n%lu %lu\n%u\n", ppm_head, self->width, self->height, 255U); size_t size = self->width * self->height; for (size_t i = 0; i < size; ++i) { // unsigned int c = self->buf[i]; struct ct_rgb* c = self->buf + i; fprintf(fout, "%c%c%c", c->r, c->g, c->b); } if (!fclose(fout)) return true; } return false; } /* the axes */ struct ct_axes { double xmin, xmax, ymin, ymax; }; static struct ct_axes ct_axes_from_point(double complex z, double radius) { struct ct_axes axes = { creal(z) - radius, creal(z) + radius, cimag(z) - radius, cimag(z) + radius }; return axes; } /* The Plane */ struct ct_plane { struct ct_axes axes; double xstep; double ystep; }; static void ct_plane_init(struct ct_plane* const self, struct ct_axes const* axes, unsigned long width, unsigned long height) { self->axes = *axes; double awidth = self->axes.xmax - self->axes.xmin; double aheight = self->axes.ymax - self->axes.ymin; assert(width > 0 && height > 0 && awidth > 0.0); double daspect = fabs((double)height / (double)width); double waspect = fabs(aheight / awidth); if (daspect > waspect) { double excess = aheight * (daspect / waspect - 1.0); self->axes.ymax += excess / 2.0; self->axes.ymin -= excess / 2.0; } else if (daspect < waspect) { double excess = awidth * (waspect / daspect - 1.0); self->axes.xmax += excess / 2.0; self->axes.xmin -= excess / 2.0; } self->xstep = (double)(self->axes.xmax - self->axes.xmin) / (double)width; self->ystep = (double)(self->axes.ymax - self->axes.ymin) / (double)height; return; } /* the plot */ struct ct_plot { struct ct_plane plane; struct ct_canvas* canvas; }; static void ct_plot_init(struct ct_plot* const self, struct ct_axes const* axes, struct ct_canvas* canvas) { ct_plane_init(&self->plane, axes, canvas->width - 1, canvas->height - 1); self->canvas = canvas; } #if (0) static bool ct_plot_add(struct ct_plot* const self, double complex z, struct ct_rgb const* color) { long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep; long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep; if (x > -1 && x < (long)self->canvas->width && y > -1 && y < (long)self->canvas->height) { // Now, we can convert to index. size_t i = x + y * self->canvas->width; assert(i < self->canvas->height * self->canvas->width); struct ct_rgb exist = self->canvas->buf[i]; exist.r += 1; self->canvas->buf[i] = exist; return true; } return true; } #endif static bool ct_plot_pixel(struct ct_plot* const self, long x, long y, struct ct_rgb const* color) { if (x > -1 && x < (long)self->canvas->width && y > -1 && y < (long)self->canvas->height) { // Now, we can convert to index. size_t i = x + y * self->canvas->width; assert(i < self->canvas->height * self->canvas->width); self->canvas->buf[i] = *color; return true; } return false; } static void ct_plot_clear(struct ct_plot* const self, struct ct_rgb const* color) { size_t size = self->canvas->height * self->canvas->width; for (size_t i = 0; i < size; ++i) self->canvas->buf[i] = *color; return; } static double complex ct_project_to_xy(struct ct_plot* const self, long x, long y) { double complex p = (self->plane.axes.xmin + x * self->plane.xstep) + (self->plane.axes.ymax - y * self->plane.ystep) * I; return p; } static bool ct_plot_point(struct ct_plot* const self, double complex z, struct ct_rgb const* color) { long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep; long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep; if (x > -1 && x < (long)self->canvas->width && y > -1 && y < (long)self->canvas->height) { // Now, we can convert to index. size_t i = x + y * self->canvas->width; assert(i < self->canvas->height * self->canvas->width); self->canvas->buf[i] = *color; return true; } return false; } /* slow, so what for now... ;^) */ static void ct_circle(struct ct_plot* const plot, double complex c, double radius, unsigned int n) { double abase = 6.2831853071 / n; for (unsigned int i = 0; i < n; ++i) { double angle = abase * i; double complex z = (creal(c) + cos(angle) * radius) + (cimag(c) + sin(angle) * radius) * I; struct ct_rgb rgb = { 255, 255, 255 }; ct_plot_point(plot, z, &rgb); } return; } /* simple vec3 for just what I need */ struct ct_vec3 { double x, y, z; }; static double ct_vev3_length(struct ct_vec3 const p) { double sum = p.x * p.x + p.y * p.y + p.z * p.z; return sqrt(sum); } static struct ct_vec3 ct_vec3_add(struct ct_vec3 const p, struct ct_vec3 const addend) { struct ct_vec3 sum = { p.x + addend.x, p.y + addend.y, p.z + addend.z }; return sum; } /* the Power Stack Mandelbulb */ static void ct_iterate_mbulb_pixel(struct ct_plot* const plot, struct ct_vec3 const z0, struct ct_vec3 const c0, double power, long x, long y, unsigned int n) { struct ct_vec3 zs = z0; struct ct_vec3 cs = c0; for (unsigned long i = 0; i < n; ++i) { double r = ct_vev3_length(zs); double rpower = pow(r, power); double angle0 = atan2(zs.z, sqrt(zs.x * zs.x + zs.y * zs.y)); double angle1 = atan2(zs.y, zs.x); struct ct_vec3 znew; znew.x = rpower * cos(power * angle1) * cos(power * angle0); znew.y = rpower * sin(power * angle1) * cos(power * angle0); znew.z = rpower * sin(power * angle0); zs = ct_vec3_add(znew, cs); if (r > 2.0) return; } struct ct_rgb rgb = { 0, 255, 0 }; ct_plot_pixel(plot, x, y, &rgb); return; } /* gain the field... */ static void ct_iterate_mbulb_plane(struct ct_plot* const plot, double zaxis, double power, unsigned int n) { for (long y = 0; y < plot->canvas->height; ++y) { for (long x = 0; x < plot->canvas->width; ++x) { double complex cz = ct_project_to_xy(plot, x, y); // Well, our per slice coords... struct ct_vec3 z; z.x = creal(cz); z.y = cimag(cz); z.z = zaxis; // Iterate the slice ct_iterate_mbulb_pixel(plot, z, z, power, x, y, n); } if (! (y % (plot->canvas->height / 10))) { printf("ct_iterate_mbulb_plane: %lu of %lu\r", y + 1, plot->canvas->height); fflush(stdout); } } printf("ct_iterate_mbulb_plane: %lu of %lu\n\n", plot->canvas->height, plot->canvas->height); fflush(stdout); return; } /* The Animation Driver */ static void ct_anime(struct ct_plot* const plot, unsigned long frame_n, unsigned long n) { assert(frame_n > 0); double zmin = 0.; double zmax = 1.; double zdif = zmax - zmin; double zbase = zdif / ((double)frame_n - 1.); double pmin = 2.; double pmax = 10.; double pdif = pmax - pmin; double pbase = pdif / ((double)frame_n - 1.); char fname[256] = { '\0' }; for (unsigned long frame_i = 0; frame_i < frame_n; ++frame_i) { double zaxis = zmin + zbase * (double)frame_i; double power = pmin + pbase * (double)frame_i; sprintf(fname, "ct_anime_%03lu.ppm", frame_i); printf("ct_anime: %lu of %lu, zaxis(%lf), power(%lf), (%s)\n", frame_i, frame_n, zaxis, power, fname); fflush(stdout); struct ct_rgb black = { 0, 0, 0 }; ct_plot_clear(plot, &black); ct_iterate_mbulb_plane(plot, zaxis, power, n); ct_canvas_save_ppm(plot->canvas, fname); } printf("\n\nct_anime complete! %lu frames: \n", frame_n); fflush(stdout); return; } static void ct_test_plane(struct ct_plot* const plot) { printf("Generating...\n\n"); fflush(stdout); ct_anime(plot, CT_FRAMES, CT_ITERS); return; } int main(void) { struct ct_canvas canvas; if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT)) { double complex plane_origin = 0+0*I; double plane_radius = 1.5; struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius); struct ct_plot plot; ct_plot_init(&plot, &axes, &canvas); ct_test_plane(&plot); // Our unit circle ct_circle(&plot, 0+0*I, 1.0, 2048); ct_circle(&plot, 2+0*I, 2.0, 2048); ct_circle(&plot, -2+0*I, 2.0, 2048); ct_circle(&plot, 0+2*I, 2.0, 2048); ct_circle(&plot, 0-2*I, 2.0, 2048); ct_canvas_save_ppm(&canvas, "ct_canvas.ppm"); ct_canvas_destroy(&canvas); return EXIT_SUCCESS; } return EXIT_FAILURE; }