crude vector field plotter experiment...

[Chris M. Thomasson]

Fwiw, here is a highly experimental and very crude test of my vector
field algorihtm that produces a file called ct_plane.ppm. The code was
quickly ported from pure C99. However, its going to make a C++
programmers eyes bleed! Using defines and C io, ect... My port was
basically just enough to get it to compile in C++11. I need to work on
it to make proper classes out of ct_canvas, ct_axes, ct_plane and ct_rgb.

I am just wondering if anybody else can get it to compile on their C++11
compiler of choice? Does it produce the ppm; can you view it? How many
warnings do you get?


Any advise?


Here is the code:
__________________________________
/*
Spiral Fractal Experiment for Alf
By: Chris M. Thomasson
Version: Pre-Alpha (0.0.0)


This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.

*_____________________________________________________________*/


#include <cstdlib>
#include <cstdio>
#include <cassert>
#include <complex>
#include <cmath>


#define CT_WIDTH (1920 * 1)
#define CT_HEIGHT (1080 * 1)

#define CT_PI 3.1415926535897932384626433832795
#define CT_PI2 (CT_PI * 2)


typedef std::complex<double> ct_complex;


struct ct_rgb_meta
{
double red;
double blue;
double green;
double alpha;
};


#define CT_RGB_META(mp_red, ct_blue, ct_green, ct_alpha) \
{ (mp_red), (ct_blue), (ct_green), (ct_alpha) }



struct ct_rgb
{
unsigned char r;
unsigned char g;
unsigned char b;
struct ct_rgb_meta meta; // meta data for this color
};


#define CT_RGB(mp_red, mp_green, mp_blue) \
{ (mp_red), (mp_green), (mp_blue), CT_RGB_META(0., 0., 0., 0.) }


struct ct_canvas
{
unsigned long width;
unsigned long height;
struct ct_rgb* buf;
};

bool
ct_canvas_create(
struct ct_canvas* const self,
unsigned long width,
unsigned long height
) {
std::size_t size = width * height * sizeof(*self->buf);

self->buf = (ct_rgb*)std::calloc(1, size);

if (self->buf)
{
self->width = width;
self->height = height;

return true;
}

return false;
}

void
ct_canvas_destroy(
struct ct_canvas const* const self
) {
std::free(self->buf);
}



double
ct_canvas_log_density_post_processing_get_largest(
struct ct_canvas const* const self
) {
double largest = 0;

std::size_t size = self->width * self->height;

for (std::size_t i = 0; i < size; ++i)
{
struct ct_rgb const* c = self->buf + i;

if (largest < c->meta.alpha)
{
largest = c->meta.alpha;
}
}

return largest;
}


void
ct_canvas_log_density_post_processing(
struct ct_canvas* const self
) {
double largest =
ct_canvas_log_density_post_processing_get_largest(self);
std::size_t size = self->width * self->height;

printf("largest = %f\n", largest);

for (std::size_t i = 0; i < size; ++i)
{
struct ct_rgb* color = self->buf + i;

if (color->meta.alpha > 0)
{
struct ct_rgb_meta c = color->meta;

double dense = log10(c.alpha) / c.alpha;

c.red *= dense;
if (c.red > 1.0) c.red = 1.0;

c.green *= dense;
if (c.green > 1.0) c.green = 1.0;

c.blue *= dense;
if (c.blue > 1.0) c.blue = 1.0;

color->r = (unsigned char)(c.red * 255);
color->g = (unsigned char)(c.green * 255);
color->b = (unsigned char)(c.blue * 255);
}
}
}




bool
ct_canvas_save_ppm(
struct ct_canvas const* const self,
char const* fname
) {
std::FILE* fout = std::fopen(fname, "wb");

if (fout)
{
char const ppm_head[] =
"P6\n"
"# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)";

std::fprintf(fout, "%s\n%lu %lu\n%u\n",
ppm_head,
self->width, self->height,
255U);

std::size_t size = self->width * self->height;

for (std::size_t i = 0; i < size; ++i)
{
struct ct_rgb* c = self->buf + i;
std::fprintf(fout, "%c%c%c", c->r, c->g, c->b);
}

if (! std::fclose(fout))
{
return true;
}
}

return false;
}




struct ct_axes
{
double xmin;
double xmax;
double ymin;
double ymax;
};

struct ct_axes
ct_axes_from_point(
ct_complex z,
double radius
) {
struct ct_axes axes = {
z.real() - radius, z.real() + radius,
z.imag() - radius, z.imag() + radius
};

return axes;
}


struct ct_plane
{
struct ct_axes axes;
double xstep;
double ystep;
};


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 / 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 = (self->axes.xmax - self->axes.xmin) / width;
self->ystep = (self->axes.ymax - self->axes.ymin) / height;
}



struct ct_plot
{
struct ct_plane plane;
struct ct_canvas* canvas;
};


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;
}


bool
ct_plot_point(
struct ct_plot* const self,
ct_complex z,
struct ct_rgb const* color
) {
long x = (long)((z.real() - self->plane.axes.xmin) /
self->plane.xstep);
long y = (long)((self->plane.axes.ymax - z.imag()) /
self->plane.ystep);

if (x > -1 && x < (long)self->canvas->width &&
y > -1 && y < (long)self->canvas->height)
{
// Now, we can convert to index.
std::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;
}



bool
ct_plot_add(
struct ct_plot* const self,
ct_complex z,
struct ct_rgb const* color
) {
long x = (long)((z.real() - self->plane.axes.xmin) /
self->plane.xstep);
long y = (long)((self->plane.axes.ymax - z.imag()) /
self->plane.ystep);

if (x > -1 && x < (long)self->canvas->width &&
y > -1 && y < (long)self->canvas->height)
{
// Now, we can convert to index.
std::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->meta.red += color->meta.red;
exist->meta.green += color->meta.green;
exist->meta.blue += color->meta.blue;
exist->meta.alpha += 1.0;

return true;
}

return true;
}


// slow, so what for now... ;^)
void ct_circle(
struct ct_plot* const plot,
ct_complex c,
double radius,
unsigned int n
) {
double abase = CT_PI2 / n;

for (unsigned int i = 0; i < n; ++i)
{
double angle = abase * i;

ct_complex z = {
c.real() + cos(angle) * radius,
c.imag() + sin(angle) * radius
};

struct ct_rgb rgb = { 255, 255, 255 };

ct_plot_point(plot, z, &rgb);
}
}


void ct_line(
struct ct_plot* const plot,
ct_complex p0,
ct_complex p1,
unsigned int n
) {
ct_complex dif = p1 - p0;

double normal_base = 1. / n;

for (unsigned int i = 0; i < n; ++i)
{
double normal = normal_base * i;

ct_complex z = p0 + dif * normal;

struct ct_rgb rgb = { 255, 255, 0 };

ct_plot_point(plot, z, &rgb);
}
}





/* Vector Field
By: Chris M. Thomasson
______________________________________________________*/
#include <vector>

struct ct_vfield_point
{
ct_complex p;
double m;
ct_rgb color;
};

typedef std::vector<ct_vfield_point> ct_vfield_vector;


ct_complex
ct_normalize(
ct_complex p
) {
double dis = std::abs(p);

if (! std::isnan(dis) && dis > 0.0000001)
{
p /= dis;
}

return p;
}

ct_complex
ct_vfield_normal(
ct_vfield_vector const& self,
ct_complex p,
double npow
) {
ct_complex g = { 0.0, 0.0 };

const int imax = self.size();

for (int i = 0; i < imax; ++i)
{
ct_complex dif = self[i].p - p;

double sum = dif.real() * dif.real() + dif.imag() * dif.imag();
double mass = std::pow(sum, npow);

if (! std::isnan(mass) && mass > 0.000001)
{
g.real(g.real() + self[i].m * dif.real() / mass);
g.imag(g.imag() + self[i].m * dif.imag() / mass);
}
}

return ct_normalize(g);
}


double ct_pi_normal(
ct_complex p,
double start_angle = 0.0
) {
double angle = std::arg(p) + start_angle;
if (angle < 0.0) angle = angle + CT_PI2;
angle = angle / CT_PI2;
return angle;
}



void ct_vfield_plot_line(
ct_vfield_vector const& self,
struct ct_plot* const plot,
ct_complex p0,
double power,
double astart,
unsigned int n,
ct_rgb const& color
) {
double radius = .028; // 1 / (n - 0.0) * CT_PI;// .01;
ct_complex z = p0;

for (unsigned int i = 0; i < n; ++i)
{
//double ir = i / (n - 1.0);

ct_complex vn = ct_vfield_normal(self, z, power);


double angle = astart + arg(vn);

//double xxxx = ct_pi_normal(angle, 0);

ct_complex zn = {
z.real() + std::cos(angle) * radius,
z.imag() + std::sin(angle) * radius
};

//ct_rgb color_angle = color;// CT_RGBF(1, 1, 1);

// plot.line(z, zn, 1.9, color_angle);

ct_line(plot, z, zn, 128);

z = zn;
}
}


void ct_vfield_plot_line_para(
ct_vfield_vector const& self,
struct ct_plot* const plot,
ct_complex p0,
double power,
double astart,
double astart_radius,
unsigned int n,
ct_rgb lcolor
) {
double abase = CT_PI2 / n;
double radius = astart_radius;

//ct_rgb bcolor = { 1, 0, 0, { 0, 0, 0 } };

ct_complex xxxSTART = ct_vfield_normal(self, p0, power);
double angleXXX = std::arg((p0 + p0 * xxxSTART) - p0);

for (unsigned int i = 0; i < n; ++i)
{
//double ir = i / (n - 1.);

double angle = abase * i + .0001 + angleXXX;

ct_complex p = {
std::cos(angle) * radius + p0.real(),
std::sin(angle) * radius + p0.imag()
};

ct_rgb xcolor = lcolor;

ct_line(plot, p0, p, 256);

ct_vfield_plot_line(self, plot, p, power, astart, 1024, xcolor);
}
}


// Plotter Lines
void ct_vfield_plot_lines(
ct_vfield_vector const& self,
struct ct_plot* const plot,
double power,
double astart,
double astart_radius,
unsigned int n
) {
std::printf("\n\nct_vfield_plot_lines plotting...\n\n");

for (unsigned int i = 0; i < self.size(); ++i)
{
//double ir = i / (self.size() - 1.);
if (self[i].m > 0.0) continue;
ct_rgb lcolor = self[i].color;
ct_vfield_plot_line_para(self, plot, self[i].p, power, astart,
astart_radius, n, lcolor);
std::printf("line:(%u of %llu)\r", i, (unsigned long
long)self.size() - 1);
}
}








/* Spiral Fractal
By: Chris M. Thomasson
______________________________________________________*/

void ct_spiral_fractal(
unsigned int ri,
unsigned int rn,
struct ct_plot* const plot,
ct_vfield_vector& vfield,
ct_complex p0,
ct_complex p1,
unsigned int n
) {
if (ri >= rn) return;

ct_complex dif = p1 - p0;

double sangle = std::arg(dif);
double sradius = std::abs(dif);

double rbase = sradius / n;
double abase = CT_PI / 2 / n;

for (unsigned long i = 0; i < n; ++i)
{
double radius = rbase * i;

double angle = 0 + std::sin(abase * i) * CT_PI2 + sangle;
double anglenx = 0 + std::sin(abase * (i + 1)) * CT_PI2 + sangle;

ct_complex s0_raw = { std::cos(angle), std::sin(angle) };
ct_complex s1_raw = { std::cos(anglenx), std::sin(anglenx) };

ct_complex s0 = p0 + s0_raw * radius;
ct_complex s1 = p0 + s1_raw * (radius + rbase);

ct_rgb color = { 0, 0, 0, { .7, .6, .5 } };

ct_plot_add(plot, s0, &color);
ct_plot_add(plot, s1, &color);

//ct_plot_add(plot, -s0, &color);
// ct_plot_add(plot, -s1, &color);


if (ri == rn - 1)
{
//ct_line(plot, s0, s1, 128);

vfield.push_back({ s0, -1, { 1, 0, 0, { 0, 0, 0 } } });
//vfield.push_back({ s1, -1, { 1, 0, 0, { 0, 0, 0 } } });
}


if (!((i + 1) % 2))
{
ct_spiral_fractal(ri + 1, rn, plot, vfield, s0, s1, n * 2);
}

else
{
ct_spiral_fractal(ri + 1, rn, plot, vfield, s1, s0, n * 2);
}
}
}




/* Entry
______________________________________________________*/
int main(void)
{
std::printf("\n\nChris M. Thomasson's Vector Field based on Spiral
Plot\n\n");


struct ct_canvas canvas;

if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT))
{
ct_complex plane_origin = { 0, 0 };
double plane_radius = 2.0;

struct ct_axes axes = ct_axes_from_point(plane_origin,
plane_radius);

struct ct_plot plot;
ct_plot_init(&plot, &axes, &canvas);

// Our unit circle
ct_circle(&plot, { 0, 0 }, 1.0, 2048);
ct_circle(&plot, { 2, 0 }, 2.0, 2048);
ct_circle(&plot, { -2, 0 }, 2.0, 2048);
ct_circle(&plot, { 0, 2 }, 2.0, 2048);
ct_circle(&plot, { 0, -2 }, 2.0, 2048);

{
ct_vfield_vector vfield;

std::printf("\n\nct_spiral_fractal plotting...\n\n");
ct_spiral_fractal(0, 3, &plot, vfield, { -1, 0 }, { 1, 0 }, 3);

ct_vfield_plot_lines(vfield, &plot, 2., 1.2, 0.04, 6);
}

ct_canvas_save_ppm(&canvas, "ct_plane.ppm");

ct_canvas_destroy(&canvas);

std::printf("\n\nPlot Complete! :^)\n\n");

return EXIT_SUCCESS;
}

return EXIT_FAILURE;
}
__________________________________

[Chris M. Thomasson]

On 1/8/2021 5:49 PM, Chris M. Thomasson wrote:
> Fwiw, here is a highly experimental and very crude test of my vector
> field algorihtm that produces a file called ct_plane.ppm. The code was
> quickly ported from pure C99. However, its going to make a C++
> programmers eyes bleed! Using defines and C io, ect... My port was
> basically just enough to get it to compile in C++11. I need to work on
> it to make proper classes out of ct_canvas, ct_axes, ct_plane and ct_rgb.
>
> I am just wondering if anybody else can get it to compile on their C++11
> compiler of choice? Does it produce the ppm; can you view it? How many
> warnings do you get?
>
>
> Any advise?
>
>
> Here is the code:
> __________________________________

[...]

So far, I see a couple of places where I failed to add in a std:: in
some math functions during the port. fabs, log10, cos and sin.

ct_circle
ct_plane_init
ct_canvas_log_density_post_processing

[Chris M. Thomasson]

On 1/8/2021 5:49 PM, Chris M. Thomasson wrote:

> Fwiw, here is a highly experimental and very crude test of my vector
> field algorihtm that produces a file called ct_plane.ppm. The code was
> quickly ported from pure C99. However, its going to make a C++
> programmers eyes bleed! Using defines and C io, ect... My port was
> basically just enough to get it to compile in C++11. I need to work on
> it to make proper classes out of ct_canvas, ct_axes, ct_plane and ct_rgb.
>
> I am just wondering if anybody else can get it to compile on their C++11
> compiler of choice? Does it produce the ppm; can you view it? How many
> warnings do you get?
>
>
> Any advise?

[...]

An easier way to get the code:

https://pastebin.com/raw/txz3pQ2Y