SHARE
TWEET

Power Stack Mandelbulb

Chris_M_Thomasson Nov 19th, 2019 138 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /*
  2.     Experimental Power Stack Mandelbulb for Owen
  3.     By: Chris M. Thomasson
  4. *_____________________________________________________________*/
  5.  
  6.  
  7. #include <stdlib.h>
  8. #include <stdio.h>
  9. #include <assert.h>
  10. #include <complex.h>
  11. #include <tgmath.h>
  12. #include <stdbool.h>
  13.  
  14.  
  15. #define CT_WIDTH 256    // width of plane
  16. #define CT_HEIGHT 256   // height of plane
  17. #define CT_FRAMES 256   // slices of the mbulb
  18. #define CT_ITERS 64     // iterations per pixel
  19.  
  20.  
  21. /* The Canvas
  22. ___________________________________*/
  23. struct ct_rgb
  24. {
  25.     unsigned char r;
  26.     unsigned char g;
  27.     unsigned char b;
  28. };
  29.  
  30.  
  31. struct ct_canvas
  32. {
  33.     unsigned long width;
  34.     unsigned long height;
  35.     struct ct_rgb* buf;
  36. };
  37.  
  38. bool
  39. ct_canvas_create(
  40.     struct ct_canvas* const self,
  41.     unsigned long width,
  42.     unsigned long height
  43. ){
  44.     size_t size = width * height * sizeof(*self->buf);
  45.  
  46.     self->buf = calloc(1, size);
  47.  
  48.     if (self->buf)
  49.     {
  50.         self->width = width;
  51.         self->height = height;
  52.  
  53.         return true;
  54.     }
  55.  
  56.     return false;
  57. }
  58.  
  59. void
  60. ct_canvas_destroy(
  61.     struct ct_canvas const* const self
  62. ){
  63.     free(self->buf);
  64. }
  65.  
  66. bool
  67. ct_canvas_save_ppm(
  68.     struct ct_canvas const* const self,
  69.     char const* fname
  70. ){
  71.     FILE* fout = fopen(fname, "wb");
  72.  
  73.     if (fout)
  74.     {
  75.         char const ppm_head[] =
  76.             "P6\n"
  77.             "# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)";
  78.  
  79.         fprintf(fout, "%s\n%lu %lu\n%u\n",
  80.                 ppm_head,
  81.                 self->width, self->height,
  82.                 255U);
  83.  
  84.         size_t size = self->width * self->height;
  85.  
  86.         for (size_t i = 0; i < size; ++i)
  87.         {
  88.             //unsigned int c = self->buf[i];
  89.             struct ct_rgb* c = self->buf + i;
  90.  
  91.             fprintf(fout, "%c%c%c", c->r, c->g, c->b);
  92.         }
  93.  
  94.         if (! fclose(fout))
  95.         {
  96.             return true;
  97.         }
  98.     }
  99.  
  100.     return false;
  101. }
  102.  
  103.  
  104.  
  105. /* The Axes
  106. ___________________________________*/
  107. struct ct_axes
  108. {
  109.     double xmin;
  110.     double xmax;
  111.     double ymin;
  112.     double ymax;
  113. };
  114.  
  115. struct ct_axes
  116. ct_axes_from_point(
  117.     double complex z,
  118.     double radius
  119. ){
  120.     struct ct_axes axes = {
  121.         creal(z) - radius, creal(z) + radius,
  122.         cimag(z) - radius, cimag(z) + radius
  123.     };
  124.  
  125.     return axes;
  126. }
  127.  
  128.  
  129.  
  130. /* The Plane
  131. ___________________________________*/
  132. struct ct_plane
  133. {
  134.     struct ct_axes axes;
  135.     double xstep;
  136.     double ystep;
  137. };
  138.  
  139.  
  140. void
  141. ct_plane_init(
  142.     struct ct_plane* const self,
  143.     struct ct_axes const* axes,
  144.     unsigned long width,
  145.     unsigned long height
  146. ){
  147.     self->axes = *axes;
  148.  
  149.     double awidth = self->axes.xmax - self->axes.xmin;
  150.     double aheight = self->axes.ymax - self->axes.ymin;
  151.  
  152.     assert(width > 0 && height > 0 && awidth > 0.0);
  153.  
  154.     double daspect = fabs((double)height / width);
  155.     double waspect = fabs(aheight / awidth);
  156.  
  157.     if (daspect > waspect)
  158.     {
  159.         double excess = aheight * (daspect / waspect - 1.0);
  160.         self->axes.ymax += excess / 2.0;
  161.         self->axes.ymin -= excess / 2.0;
  162.     }
  163.  
  164.     else if (daspect < waspect)
  165.     {
  166.         double excess = awidth * (waspect / daspect - 1.0);
  167.         self->axes.xmax += excess / 2.0;
  168.         self->axes.xmin -= excess / 2.0;
  169.     }
  170.  
  171.     self->xstep = (self->axes.xmax - self->axes.xmin) / width;
  172.     self->ystep = (self->axes.ymax - self->axes.ymin) / height;
  173. }
  174.  
  175.  
  176.  
  177. /* The Plot
  178. ___________________________________*/
  179. struct ct_plot
  180. {
  181.     struct ct_plane plane;
  182.     struct ct_canvas* canvas;
  183. };
  184.  
  185.  
  186. void
  187. ct_plot_init(
  188.     struct ct_plot* const self,
  189.     struct ct_axes const* axes,
  190.     struct ct_canvas* canvas
  191. ){
  192.     ct_plane_init(&self->plane, axes, canvas->width - 1, canvas->height - 1);
  193.     self->canvas = canvas;
  194. }
  195.  
  196.  
  197. bool
  198. ct_plot_add(
  199.     struct ct_plot* const self,
  200.     double complex z,
  201.     struct ct_rgb const* color
  202. ){
  203.     long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep;
  204.     long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep;
  205.  
  206.     if (x > -1 && x < (long)self->canvas->width &&
  207.         y > -1 && y < (long)self->canvas->height)
  208.     {
  209.         // Now, we can convert to index.
  210.         size_t i = x + y * self->canvas->width;
  211.  
  212.         assert(i < self->canvas->height * self->canvas->width);
  213.  
  214.         struct ct_rgb exist = self->canvas->buf[i];
  215.  
  216.         exist.r += 1;
  217.  
  218.         self->canvas->buf[i] = exist;
  219.         return true;
  220.     }
  221.  
  222.     return true;
  223. }
  224.  
  225.  
  226. bool
  227. ct_plot_pixel(
  228.     struct ct_plot* const self,
  229.     long x,
  230.     long y,
  231.     struct ct_rgb const* color
  232. ){
  233.     if (x > -1 && x < (long)self->canvas->width &&
  234.         y > -1 && y < (long)self->canvas->height)
  235.     {
  236.         // Now, we can convert to index.
  237.         size_t i = x + y * self->canvas->width;
  238.  
  239.         assert(i < self->canvas->height * self->canvas->width);
  240.  
  241.         self->canvas->buf[i] = *color;
  242.         return true;
  243.     }
  244.  
  245.     return false;
  246. }
  247.  
  248.  
  249. void
  250. ct_plot_clear(
  251.     struct ct_plot* const self,
  252.     struct ct_rgb const* color
  253. ){
  254.     size_t size = self->canvas->height * self->canvas->width;
  255.  
  256.     for (size_t i = 0; i < size; ++i)
  257.     {
  258.         self->canvas->buf[i] = *color;
  259.     }
  260. }
  261.  
  262.  
  263. double complex
  264. ct_project_to_xy(
  265.     struct ct_plot* const self,
  266.     long x,
  267.     long y
  268. ){
  269.     double complex p =
  270.         (self->plane.axes.xmin + x * self->plane.xstep) +
  271.         (self->plane.axes.ymax - y * self->plane.ystep) * I;
  272.  
  273.     return p;
  274. }
  275.  
  276.  
  277. bool
  278. ct_plot_point(
  279.     struct ct_plot* const self,
  280.     double complex z,
  281.     struct ct_rgb const* color
  282. ){
  283.     long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep;
  284.     long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep;
  285.  
  286.     if (x > -1 && x < (long)self->canvas->width &&
  287.         y > -1 && y < (long)self->canvas->height)
  288.     {
  289.         // Now, we can convert to index.
  290.         size_t i = x + y * self->canvas->width;
  291.  
  292.         assert(i < self->canvas->height * self->canvas->width);
  293.  
  294.         self->canvas->buf[i] = *color;
  295.         return true;
  296.     }
  297.  
  298.     return false;
  299. }
  300.  
  301.  
  302. // slow, so what for now... ;^)
  303. void ct_circle(
  304.     struct ct_plot* const plot,
  305.     double complex c,
  306.     double radius,
  307.     unsigned int n
  308. ){
  309.     double abase = 6.2831853071 / n;
  310.  
  311.     for (unsigned int i = 0; i < n; ++i)
  312.     {
  313.         double angle = abase * i;
  314.  
  315.         double complex z =
  316.             (creal(c) + cos(angle) * radius) +
  317.             (cimag(c) + sin(angle) * radius) * I;
  318.  
  319.         struct ct_rgb rgb = { 255, 255, 255 };
  320.  
  321.         ct_plot_point(plot, z, &rgb);
  322.     }
  323. }
  324.  
  325.  
  326.  
  327. /* Simple vec3 for just what I need
  328. ___________________________________*/
  329. struct ct_vec3
  330. {
  331.     double x;
  332.     double y;
  333.     double z;
  334. };
  335.  
  336.  
  337. double ct_vev3_length(
  338.     struct ct_vec3 const p
  339. ){
  340.     double sum = p.x * p.x + p.y * p.y + p.z * p.z;
  341.     return sqrt(sum);
  342. }
  343.  
  344.  
  345. struct ct_vec3
  346. ct_vec3_add(
  347.     struct ct_vec3 const p,
  348.     struct ct_vec3 const addend
  349. ){
  350.     struct ct_vec3 sum = {
  351.         p.x + addend.x,
  352.         p.y + addend.y,
  353.         p.z + addend.z
  354.     };
  355.  
  356.     return sum;
  357. }
  358.  
  359.  
  360.  
  361. /* The Power Stack Mandelbulb
  362. ___________________________________*/
  363. void
  364. ct_iterate_mbulb_pixel(
  365.     struct ct_plot* const plot,
  366.     struct ct_vec3 const z0,
  367.     struct ct_vec3 const c0,
  368.     double power,
  369.     long x,
  370.     long y,
  371.     unsigned int n
  372. ){
  373.     struct ct_vec3 zs = z0;
  374.     struct ct_vec3 cs = c0;
  375.  
  376.     for (unsigned long i = 0; i < n; ++i)
  377.     {
  378.         double r = ct_vev3_length(zs);
  379.         double rpower = pow(r, power);
  380.  
  381.         double angle0 = atan2(zs.z, sqrt(zs.x * zs.x + zs.y * zs.y));
  382.         double angle1 = atan2(zs.y, zs.x);
  383.  
  384.         struct ct_vec3 znew;
  385.         znew.x = rpower * cos(power * angle1) * cos(power * angle0);
  386.         znew.y = rpower * sin(power * angle1) * cos(power * angle0);
  387.         znew.z = rpower * sin(power * angle0);
  388.  
  389.         zs = ct_vec3_add(znew, cs);
  390.  
  391.         if (r > 2.0)
  392.         {
  393.             return;
  394.         }
  395.     }
  396.  
  397.     struct ct_rgb rgb = { 0, 255, 0 };
  398.     ct_plot_pixel(plot, x, y, &rgb);
  399. }
  400.  
  401.  
  402. // Gain the field...
  403. void
  404. ct_iterate_mbulb_plane(
  405.     struct ct_plot* const plot,
  406.     double zaxis,
  407.     double power,
  408.     unsigned int n
  409. ){
  410.     for (long y = 0; y < plot->canvas->height; ++y)
  411.     {
  412.         for (long x = 0; x < plot->canvas->width; ++x)
  413.         {
  414.             double complex cz = ct_project_to_xy(plot, x, y);
  415.  
  416.             // Well, our per slice coords...
  417.             struct ct_vec3 z;
  418.  
  419.             z.x = creal(cz);
  420.             z.y = cimag(cz);
  421.             z.z = zaxis;
  422.  
  423.             // Iterate the slice
  424.             ct_iterate_mbulb_pixel(plot, z, z, power, x, y, n);
  425.         }
  426.  
  427.         if (! (y % (plot->canvas->height / 10)))
  428.         {
  429.             printf("ct_iterate_mbulb_plane: %lu of %lu\r",
  430.                    y + 1, plot->canvas->height);
  431.  
  432.             fflush(stdout);
  433.         }
  434.     }
  435.  
  436.     printf("ct_iterate_mbulb_plane: %lu of %lu\n\n",
  437.            plot->canvas->height, plot->canvas->height);
  438.  
  439.     fflush(stdout);
  440. }
  441.  
  442.  
  443. /* The Animation Driver
  444. ___________________________________*/
  445. void ct_anime(
  446.     struct ct_plot* const plot,
  447.     unsigned long frame_n,
  448.     unsigned long n
  449. ){
  450.     assert(frame_n > 0);
  451.  
  452.     double zmin = 0.;
  453.     double zmax = 1.;
  454.     double zdif = zmax - zmin;
  455.     double zbase = zdif / (frame_n - 1.);
  456.  
  457.     double pmin = 2.;
  458.     double pmax = 10.;
  459.     double pdif = pmax - pmin;
  460.     double pbase = pdif / (frame_n - 1.);
  461.  
  462.     char fname[256] = { '\0' };
  463.  
  464.     for (unsigned long frame_i = 0; frame_i < frame_n; ++frame_i)
  465.     {
  466.         double zaxis = zmin + zbase * frame_i;
  467.         double power = pmin + pbase * frame_i;
  468.  
  469.         sprintf(fname, "ct_anime_%lu.ppm", frame_i);
  470.  
  471.         printf("ct_anime: %lu of %lu, zaxis(%lf), power(%lf), (%s)\n",
  472.            frame_i, frame_n, zaxis, power, fname);
  473.  
  474.         fflush(stdout);
  475.  
  476.         struct ct_rgb black = { 0, 0, 0 };
  477.         ct_plot_clear(plot, &black);
  478.  
  479.         ct_iterate_mbulb_plane(plot, zaxis, power, n);
  480.  
  481.         ct_canvas_save_ppm(plot->canvas, fname);
  482.     }
  483.  
  484.     printf("\n\nct_anime complete! %lu frames: \n", frame_n);
  485.     fflush(stdout);
  486. }
  487.  
  488.  
  489.  
  490.  
  491. void ct_test_plane(
  492.     struct ct_plot* const plot
  493. ){
  494.     printf("Generating...\n\n");
  495.     fflush(stdout);
  496.  
  497.     ct_anime(plot, CT_FRAMES, CT_ITERS);
  498. }
  499.  
  500.  
  501. int main(void)
  502. {
  503.     struct ct_canvas canvas;
  504.  
  505.     if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT))
  506.     {
  507.         double complex plane_origin = 0+0*I;
  508.         double plane_radius = 1.5;
  509.  
  510.         struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius);
  511.  
  512.         struct ct_plot plot;
  513.         ct_plot_init(&plot, &axes, &canvas);
  514.  
  515.  
  516.         ct_test_plane(&plot);
  517.  
  518.  
  519.         // Our unit circle
  520.         ct_circle(&plot, 0+0*I, 1.0, 2048);
  521.         ct_circle(&plot, 2+0*I, 2.0, 2048);
  522.         ct_circle(&plot, -2+0*I, 2.0, 2048);
  523.         ct_circle(&plot, 0+2*I, 2.0, 2048);
  524.         ct_circle(&plot, 0-2*I, 2.0, 2048);
  525.  
  526.         ct_canvas_save_ppm(&canvas, "ct_plane.ppm");
  527.  
  528.         ct_canvas_destroy(&canvas);
  529.  
  530.         return EXIT_SUCCESS;
  531.     }
  532.  
  533.     return EXIT_FAILURE;
  534. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top