Chris_M_Thomasson

Crude Vector Field Port to C++...

Jan 20th, 2021
2,182
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /*
  2.     Spiral Fractal Experiment for Alf
  3.     By: Chris M. Thomasson
  4.     Version: Pre-Alpha (0.0.0)
  5.  
  6.  
  7.     This program is free software: you can redistribute it and/or modify
  8.     it under the terms of the GNU General Public License as published by
  9.     the Free Software Foundation, either version 3 of the License, or
  10.     (at your option) any later version.
  11.  
  12.     This program is distributed in the hope that it will be useful,
  13.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  14.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15.     GNU General Public License for more details.
  16.  
  17.     You should have received a copy of the GNU General Public License
  18.     along with this program.  If not, see <https://www.gnu.org/licenses/>.
  19.  
  20. *_____________________________________________________________*/
  21.  
  22.  
  23. #include <cstdlib>
  24. #include <cstdio>
  25. #include <cassert>
  26. #include <complex>
  27. #include <cmath>
  28.  
  29.  
  30. #define CT_WIDTH (1920 * 1)
  31. #define CT_HEIGHT (1080 * 1)
  32.  
  33. #define CT_PI 3.1415926535897932384626433832795
  34. #define CT_PI2 (CT_PI * 2)
  35.  
  36.  
  37. typedef std::complex<double> ct_complex;
  38.  
  39.  
  40. struct ct_rgb_meta
  41. {
  42.     double red;
  43.     double blue;
  44.     double green;
  45.     double alpha;
  46. };
  47.  
  48.  
  49. #define CT_RGB_META(mp_red, ct_blue, ct_green, ct_alpha) \
  50.     { (mp_red), (ct_blue), (ct_green), (ct_alpha) }
  51.  
  52.  
  53.  
  54. struct ct_rgb
  55. {
  56.     unsigned char r;
  57.     unsigned char g;
  58.     unsigned char b;
  59.     struct ct_rgb_meta meta; // meta data for this color
  60. };
  61.  
  62.  
  63. #define CT_RGB(mp_red, mp_green, mp_blue) \
  64.     { (mp_red), (mp_green), (mp_blue), CT_RGB_META(0., 0., 0., 0.) }
  65.  
  66.  
  67. struct ct_canvas
  68. {
  69.     unsigned long width;
  70.     unsigned long height;
  71.     struct ct_rgb* buf;
  72. };
  73.  
  74. bool
  75. ct_canvas_create(
  76.     struct ct_canvas* const self,
  77.     unsigned long width,
  78.     unsigned long height
  79. ) {
  80.     std::size_t size = width * height * sizeof(*self->buf);
  81.  
  82.     self->buf = (ct_rgb*)std::calloc(1, size);
  83.  
  84.     if (self->buf)
  85.     {
  86.         self->width = width;
  87.         self->height = height;
  88.  
  89.         return true;
  90.     }
  91.  
  92.     return false;
  93. }
  94.  
  95. void
  96. ct_canvas_destroy(
  97.     struct ct_canvas const* const self
  98. ) {
  99.     std::free(self->buf);
  100. }
  101.  
  102.  
  103.  
  104. double
  105. ct_canvas_log_density_post_processing_get_largest(
  106.     struct ct_canvas const* const self
  107. ) {
  108.     double largest = 0;
  109.  
  110.     std::size_t size = self->width * self->height;
  111.  
  112.     for (std::size_t i = 0; i < size; ++i)
  113.     {
  114.         struct ct_rgb const* c = self->buf + i;
  115.  
  116.         if (largest < c->meta.alpha)
  117.         {
  118.             largest = c->meta.alpha;
  119.         }
  120.     }
  121.  
  122.     return largest;
  123. }
  124.  
  125.  
  126. void
  127. ct_canvas_log_density_post_processing(
  128.     struct ct_canvas* const self
  129. ) {
  130.     double largest = ct_canvas_log_density_post_processing_get_largest(self);
  131.     std::size_t size = self->width * self->height;
  132.  
  133.     printf("largest = %f\n", largest);
  134.  
  135.     for (std::size_t i = 0; i < size; ++i)
  136.     {
  137.         struct ct_rgb* color = self->buf + i;
  138.  
  139.         if (color->meta.alpha > 0)
  140.         {
  141.             struct ct_rgb_meta c = color->meta;
  142.  
  143.             double dense = log10(c.alpha) / c.alpha;
  144.  
  145.             c.red *= dense;
  146.             if (c.red > 1.0) c.red = 1.0;
  147.  
  148.             c.green *= dense;
  149.             if (c.green > 1.0) c.green = 1.0;
  150.  
  151.             c.blue *= dense;
  152.             if (c.blue > 1.0) c.blue = 1.0;
  153.  
  154.             color->r = (unsigned char)(c.red * 255);
  155.             color->g = (unsigned char)(c.green * 255);
  156.             color->b = (unsigned char)(c.blue * 255);
  157.         }
  158.     }
  159. }
  160.  
  161.  
  162.  
  163.  
  164. bool
  165. ct_canvas_save_ppm(
  166.     struct ct_canvas const* const self,
  167.     char const* fname
  168. ) {
  169.     std::FILE* fout = std::fopen(fname, "wb");
  170.  
  171.     if (fout)
  172.     {
  173.         char const ppm_head[] =
  174.             "P6\n"
  175.             "# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)";
  176.  
  177.         std::fprintf(fout, "%s\n%lu %lu\n%u\n",
  178.             ppm_head,
  179.             self->width, self->height,
  180.             255U);
  181.  
  182.         std::size_t size = self->width * self->height;
  183.  
  184.         for (std::size_t i = 0; i < size; ++i)
  185.         {
  186.             struct ct_rgb* c = self->buf + i;
  187.             std::fprintf(fout, "%c%c%c", c->r, c->g, c->b);
  188.         }
  189.  
  190.         if (!std::fclose(fout))
  191.         {
  192.             return true;
  193.         }
  194.     }
  195.  
  196.     return false;
  197. }
  198.  
  199.  
  200.  
  201.  
  202. struct ct_axes
  203. {
  204.     double xmin;
  205.     double xmax;
  206.     double ymin;
  207.     double ymax;
  208. };
  209.  
  210. struct ct_axes
  211.     ct_axes_from_point(
  212.         ct_complex z,
  213.         double radius
  214.     ) {
  215.     struct ct_axes axes = {
  216.         z.real() - radius, z.real() + radius,
  217.         z.imag() - radius, z.imag() + radius
  218.     };
  219.  
  220.     return axes;
  221. }
  222.  
  223.  
  224. struct ct_plane
  225. {
  226.     struct ct_axes axes;
  227.     double xstep;
  228.     double ystep;
  229. };
  230.  
  231.  
  232. void
  233. ct_plane_init(
  234.     struct ct_plane* const self,
  235.     struct ct_axes const* axes,
  236.     unsigned long width,
  237.     unsigned long height
  238. ) {
  239.     self->axes = *axes;
  240.  
  241.     double awidth = self->axes.xmax - self->axes.xmin;
  242.     double aheight = self->axes.ymax - self->axes.ymin;
  243.  
  244.     assert(width > 0 && height > 0 && awidth > 0.0);
  245.  
  246.     double daspect = fabs((double)height / width);
  247.     double waspect = fabs(aheight / awidth);
  248.  
  249.     if (daspect > waspect)
  250.     {
  251.         double excess = aheight * (daspect / waspect - 1.0);
  252.         self->axes.ymax += excess / 2.0;
  253.         self->axes.ymin -= excess / 2.0;
  254.     }
  255.  
  256.     else if (daspect < waspect)
  257.     {
  258.         double excess = awidth * (waspect / daspect - 1.0);
  259.         self->axes.xmax += excess / 2.0;
  260.         self->axes.xmin -= excess / 2.0;
  261.     }
  262.  
  263.     self->xstep = (self->axes.xmax - self->axes.xmin) / width;
  264.     self->ystep = (self->axes.ymax - self->axes.ymin) / height;
  265. }
  266.  
  267.  
  268.  
  269. struct ct_plot
  270. {
  271.     struct ct_plane plane;
  272.     struct ct_canvas* canvas;
  273. };
  274.  
  275.  
  276. void
  277. ct_plot_init(
  278.     struct ct_plot* const self,
  279.     struct ct_axes const* axes,
  280.     struct ct_canvas* canvas
  281. ) {
  282.     ct_plane_init(&self->plane, axes, canvas->width - 1, canvas->height - 1);
  283.     self->canvas = canvas;
  284. }
  285.  
  286.  
  287. bool
  288. ct_plot_point(
  289.     struct ct_plot* const self,
  290.     ct_complex z,
  291.     struct ct_rgb const* color
  292. ) {
  293.     long x = (long)((z.real() - self->plane.axes.xmin) / self->plane.xstep);
  294.     long y = (long)((self->plane.axes.ymax - z.imag()) / self->plane.ystep);
  295.  
  296.     if (x > -1 && x < (long)self->canvas->width &&
  297.         y > -1 && y < (long)self->canvas->height)
  298.     {
  299.         // Now, we can convert to index.
  300.         std::size_t i = x + y * self->canvas->width;
  301.  
  302.         assert(i < self->canvas->height * self->canvas->width);
  303.  
  304.         self->canvas->buf[i] = *color;
  305.         return true;
  306.     }
  307.  
  308.     return false;
  309. }
  310.  
  311.  
  312.  
  313. bool
  314. ct_plot_add(
  315.     struct ct_plot* const self,
  316.     ct_complex z,
  317.     struct ct_rgb const* color
  318. ) {
  319.     long x = (long)((z.real() - self->plane.axes.xmin) / self->plane.xstep);
  320.     long y = (long)((self->plane.axes.ymax - z.imag()) / self->plane.ystep);
  321.  
  322.     if (x > -1 && x < (long)self->canvas->width &&
  323.         y > -1 && y < (long)self->canvas->height)
  324.     {
  325.         // Now, we can convert to index.
  326.         std::size_t i = x + y * self->canvas->width;
  327.  
  328.         assert(i < self->canvas->height * self->canvas->width);
  329.  
  330.         struct ct_rgb* exist = &self->canvas->buf[i];
  331.  
  332.         exist->meta.red += color->meta.red;
  333.         exist->meta.green += color->meta.green;
  334.         exist->meta.blue += color->meta.blue;
  335.         exist->meta.alpha += 1.0;
  336.  
  337.         return true;
  338.     }
  339.  
  340.     return true;
  341. }
  342.  
  343.  
  344. // slow, so what for now... ;^)
  345. void ct_circle(
  346.     struct ct_plot* const plot,
  347.     ct_complex c,
  348.     double radius,
  349.     unsigned int n
  350. ) {
  351.     double abase = CT_PI2 / n;
  352.  
  353.     for (unsigned int i = 0; i < n; ++i)
  354.     {
  355.         double angle = abase * i;
  356.  
  357.         ct_complex z = {
  358.             c.real() + cos(angle) * radius,
  359.             c.imag() + sin(angle) * radius
  360.         };
  361.  
  362.         struct ct_rgb rgb = { 255, 255, 255 };
  363.  
  364.         ct_plot_point(plot, z, &rgb);
  365.     }
  366. }
  367.  
  368.  
  369. void ct_line(
  370.     struct ct_plot* const plot,
  371.     ct_complex p0,
  372.     ct_complex p1,
  373.     unsigned int n
  374. ) {
  375.     ct_complex dif = p1 - p0;
  376.  
  377.     double normal_base = 1. / n;
  378.  
  379.     for (unsigned int i = 0; i < n; ++i)
  380.     {
  381.         double normal = normal_base * i;
  382.  
  383.         ct_complex z = p0 + dif * normal;
  384.  
  385.         struct ct_rgb rgb = { 255, 255, 0 };
  386.  
  387.         ct_plot_point(plot, z, &rgb);
  388.     }
  389. }
  390.  
  391.  
  392.  
  393.  
  394.  
  395. /* Vector Field
  396.    By: Chris M. Thomasson
  397. ______________________________________________________*/
  398. #include <vector>
  399.  
  400. struct ct_vfield_point
  401. {
  402.     ct_complex p;
  403.     double m;
  404.     ct_rgb color;
  405. };
  406.  
  407. typedef std::vector<ct_vfield_point> ct_vfield_vector;
  408.  
  409.  
  410. ct_complex
  411. ct_normalize(
  412.     ct_complex p
  413. ) {
  414.     double dis = std::abs(p);
  415.  
  416.     if (!std::isnan(dis) && dis > 0.0000001)
  417.     {
  418.         p /= dis;
  419.     }
  420.  
  421.     return p;
  422. }
  423.  
  424. ct_complex
  425. ct_vfield_normal(
  426.     ct_vfield_vector const& self,
  427.     ct_complex p,
  428.     double npow
  429. ) {
  430.     ct_complex g = { 0.0, 0.0 };
  431.  
  432.     const int imax = self.size();
  433.  
  434.     for (int i = 0; i < imax; ++i)
  435.     {
  436.         ct_complex dif = self[i].p - p;
  437.  
  438.         double sum = dif.real() * dif.real() + dif.imag() * dif.imag();
  439.         double mass = std::pow(sum, npow);
  440.  
  441.         if (!std::isnan(mass) && mass > 0.000001)
  442.         {
  443.             g.real(g.real() + self[i].m * dif.real() / mass);
  444.             g.imag(g.imag() + self[i].m * dif.imag() / mass);
  445.         }
  446.     }
  447.  
  448.     return ct_normalize(g);
  449. }
  450.  
  451.  
  452. double ct_pi_normal(
  453.     ct_complex p,
  454.     double start_angle = 0.0
  455. ) {
  456.     double angle = std::arg(p) + start_angle;
  457.     if (angle < 0.0) angle = angle + CT_PI2;
  458.     angle = angle / CT_PI2;
  459.     return angle;
  460. }
  461.  
  462.  
  463.  
  464. void ct_vfield_plot_line(
  465.     ct_vfield_vector const& self,
  466.     struct ct_plot* const plot,
  467.     ct_complex p0,
  468.     double power,
  469.     double astart,
  470.     unsigned int n,
  471.     ct_rgb const& color
  472. ) {
  473.     double radius = .028; // 1 / (n - 0.0) * CT_PI;// .01;
  474.     ct_complex z = p0;
  475.  
  476.     for (unsigned int i = 0; i < n; ++i)
  477.     {
  478.         //double ir = i / (n - 1.0);
  479.  
  480.         ct_complex vn = ct_vfield_normal(self, z, power);
  481.  
  482.  
  483.         double angle = astart + arg(vn);
  484.  
  485.         //double xxxx = ct_pi_normal(angle, 0);
  486.  
  487.         ct_complex zn = {
  488.             z.real() + std::cos(angle) * radius,
  489.             z.imag() + std::sin(angle) * radius
  490.         };
  491.  
  492.         //ct_rgb color_angle = color;// CT_RGBF(1, 1, 1);
  493.  
  494.        // plot.line(z, zn, 1.9, color_angle);
  495.  
  496.         ct_line(plot, z, zn, 128);
  497.  
  498.         z = zn;
  499.     }
  500. }
  501.  
  502.  
  503. void ct_vfield_plot_line_para(
  504.     ct_vfield_vector const& self,
  505.     struct ct_plot* const plot,
  506.     ct_complex p0,
  507.     double power,
  508.     double astart,
  509.     double astart_radius,
  510.     unsigned int n,
  511.     ct_rgb lcolor
  512. ) {
  513.     double abase = CT_PI2 / n;
  514.     double radius = astart_radius;
  515.  
  516.     //ct_rgb bcolor = { 1, 0, 0, { 0, 0, 0 } };
  517.  
  518.     ct_complex xxxSTART = ct_vfield_normal(self, p0, power);
  519.     double angleXXX = std::arg((p0 + p0 * xxxSTART) - p0);
  520.  
  521.     for (unsigned int i = 0; i < n; ++i)
  522.     {
  523.         //double ir = i / (n - 1.);
  524.  
  525.         double angle = abase * i + .0001 + angleXXX;
  526.  
  527.         ct_complex p = {
  528.             std::cos(angle) * radius + p0.real(),
  529.             std::sin(angle) * radius + p0.imag()
  530.         };
  531.  
  532.         ct_rgb xcolor = lcolor;
  533.  
  534.         ct_line(plot, p0, p, 256);
  535.  
  536.         ct_vfield_plot_line(self, plot, p, power, astart, 1024, xcolor);
  537.     }
  538. }
  539.  
  540.  
  541. // Plotter Lines
  542. void ct_vfield_plot_lines(
  543.     ct_vfield_vector const& self,
  544.     struct ct_plot* const plot,
  545.     double power,
  546.     double astart,
  547.     double astart_radius,
  548.     unsigned int n
  549. ) {
  550.     std::printf("\n\nct_vfield_plot_lines plotting...\n\n");
  551.  
  552.     for (unsigned int i = 0; i < self.size(); ++i)
  553.     {
  554.         //double ir = i / (self.size() - 1.);
  555.         if (self[i].m > 0.0) continue;
  556.         ct_rgb lcolor = self[i].color;
  557.         ct_vfield_plot_line_para(self, plot, self[i].p, power, astart, astart_radius, n, lcolor);
  558.         std::printf("line:(%u of %llu)\r", i, (unsigned long long)self.size() - 1);
  559.     }
  560. }
  561.  
  562.  
  563.  
  564.  
  565.  
  566.  
  567.  
  568.  
  569. /* Spiral Fractal
  570.    By: Chris M. Thomasson
  571. ______________________________________________________*/
  572.  
  573. void ct_spiral_fractal(
  574.     unsigned int ri,
  575.     unsigned int rn,
  576.     struct ct_plot* const plot,
  577.     ct_vfield_vector& vfield,
  578.     ct_complex p0,
  579.     ct_complex p1,
  580.     unsigned int n
  581. ) {
  582.     if (ri >= rn) return;
  583.  
  584.     ct_complex dif = p1 - p0;
  585.  
  586.     double sangle = std::arg(dif);
  587.     double sradius = std::abs(dif);
  588.  
  589.     double rbase = sradius / n;
  590.     double abase = CT_PI / 2 / n;
  591.  
  592.     for (unsigned long i = 0; i < n; ++i)
  593.     {
  594.         double radius = rbase * i;
  595.  
  596.         double angle = 0 + std::sin(abase * i) * CT_PI2 + sangle;
  597.         double anglenx = 0 + std::sin(abase * (i + 1)) * CT_PI2 + sangle;
  598.  
  599.         ct_complex s0_raw = { std::cos(angle), std::sin(angle) };
  600.         ct_complex s1_raw = { std::cos(anglenx), std::sin(anglenx) };
  601.  
  602.         ct_complex s0 = p0 + s0_raw * radius;
  603.         ct_complex s1 = p0 + s1_raw * (radius + rbase);
  604.  
  605.         ct_rgb color = { 0, 0, 0, { .7, .6, .5 } };
  606.  
  607.         ct_plot_add(plot, s0, &color);
  608.         ct_plot_add(plot, s1, &color);
  609.  
  610.         //ct_plot_add(plot, -s0, &color);
  611.        // ct_plot_add(plot, -s1, &color);
  612.  
  613.  
  614.         if (ri == rn - 1)
  615.         {
  616.             //ct_line(plot, s0, s1, 128);
  617.  
  618.             vfield.push_back({ s0, -1, { 1, 0, 0, { 0, 0, 0 } } });
  619.             //vfield.push_back({ s1, -1, { 1, 0, 0, { 0, 0, 0 } } });
  620.         }
  621.  
  622.  
  623.         if (!((i + 1) % 2))
  624.         {
  625.             ct_spiral_fractal(ri + 1, rn, plot, vfield, s0, s1, n * 2);
  626.         }
  627.  
  628.         else
  629.         {
  630.             ct_spiral_fractal(ri + 1, rn, plot, vfield, s1, s0, n * 2);
  631.         }
  632.     }
  633. }
  634.  
  635.  
  636.  
  637.  
  638. /* Entry
  639. ______________________________________________________*/
  640. int main(void)
  641. {
  642.     std::printf("\n\nChris M. Thomasson's Vector Field based on Spiral Plot\n\n");
  643.  
  644.  
  645.     struct ct_canvas canvas;
  646.  
  647.     if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT))
  648.     {
  649.         ct_complex plane_origin = { 0, 0 };
  650.         double plane_radius = 2.0;
  651.  
  652.         struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius);
  653.  
  654.         struct ct_plot plot;
  655.         ct_plot_init(&plot, &axes, &canvas);
  656.  
  657.         // Our unit circle
  658.         ct_circle(&plot, { 0, 0 }, 1.0, 2048);
  659.         ct_circle(&plot, { 2, 0 }, 2.0, 2048);
  660.         ct_circle(&plot, { -2, 0 }, 2.0, 2048);
  661.         ct_circle(&plot, { 0, 2 }, 2.0, 2048);
  662.         ct_circle(&plot, { 0, -2 }, 2.0, 2048);
  663.  
  664.         {
  665.             ct_vfield_vector vfield;
  666.  
  667.             std::printf("\n\nct_spiral_fractal plotting...\n\n");
  668.             ct_spiral_fractal(0, 3, &plot, vfield, { -1, 0 }, { 1, 0 }, 3);
  669.  
  670.             ct_vfield_plot_lines(vfield, &plot, 2., 1.2, 0.04, 6);
  671.         }
  672.  
  673.         ct_canvas_save_ppm(&canvas, "ct_plane.ppm");
  674.  
  675.         ct_canvas_destroy(&canvas);
  676.  
  677.         std::printf("\n\nPlot Complete! :^)\n\n");
  678.  
  679.         return EXIT_SUCCESS;
  680.     }
  681.  
  682.     return EXIT_FAILURE;
  683. }
RAW Paste Data