triangles/main.cpp

#include <random>
#include <cmath>
#include <queue>
#include <iostream>

using namespace std;
using ld = long double;

struct pt {
	ld x = 0, y = 0;
};

bool operator<(pt lhs, pt rhs) {
	return make_pair(lhs.x, lhs.y) < make_pair(rhs.x, rhs.y);
}

ostream& operator<<(ostream& out, pt p) {
	return out << p.x << ',' << p.y;
}

constexpr pt zero = pt{0, 0};

pt operator+(pt a, pt b) {
	a.x += b.x;
	a.y += b.y;
	return a;
}

pt operator-(pt a, pt b) {
	a.x -= b.x;
	a.y -= b.y;
	return a;
}

pt operator*(ld s, pt p) {
	p.x *= s;
	p.y *= s;
	return p;
}

ld dot(pt a, pt b) {
	return a.x * b.x + a.y * b.y;
}

ld normsqr(pt a) {
	return dot(a, a);
}

ld lensqr(pt a, pt b) {
	return normsqr(b - a);
}

ld norm(pt a) {
	return sqrt(normsqr(a));
}

struct cone {
	pt p1, d1;
	pt p2, d2;
};

struct cone_d {
	ld d;
	cone c;

	cone_d(cone const & co) : c(co) {
		// TODO: proper point line segment distance
		d = min(normsqr(c.p1), normsqr(c.p2));
	}
};

bool operator<(cone_d const & lhs, cone_d const & rhs) {
	return lhs.d > rhs.d;
}

// two points and a direction from p2
// returns x such that 
ld unit_triangle(pt p1, pt p2, pt d2) {
	const auto d1 = p1 - p2;
	// from https://en.wikipedia.org/wiki/Triangle#Using_coordinates
	// 1 = area = 1/2 * abs(x * d2.x * d1.y - x * d1.x * d2.y)
	// 2 = x * abs(d2.x * d1.y - d1.x * d2.y)
	// 2 / abs(d2.x * d1.y - d1.x * d2.y) = x
	const auto x = 2.0 / abs(d2.x * d1.y - d1.x * d2.y);
	return x;
}

int main() {
	mt19937 gen(random_device{}());
	
	// This is used for the starting triangles
	uniform_real_distribution<ld> dist(0.8, 1.6);

	uniform_real_distribution<ld> unif01(0, 1);

	// the three starting rays
	pt c1 = {0, 1};
	pt c2 = { sqrt(3) / 2.0, -0.5};
	pt c3 = {-sqrt(3) / 2.0, -0.5};

	// and the random triangles attached to them
	pt p12a = dist(gen) * c1;
	pt p12b = unit_triangle(p12a, {0, 0}, c2) * c2;
	cone c12{p12a, c1, p12b, c2};

	pt p23a = dist(gen) * c2;
	pt p23b = unit_triangle(p23a, {0, 0}, c3) * c3;
	cone c23{p23a, c2, p23b, c3};

	pt p31a = dist(gen) * c3;
	pt p31b = unit_triangle(p31a, {0, 0}, c1) * c1;
	cone c31{p31a, c3, p31b, c1};

	// we keep track of all triangles in the construction
	vector<tuple<pt, pt, pt>> triangles;
	triangles.emplace_back(p12a, p12b, zero);
	triangles.emplace_back(p23a, p23b, zero);
	triangles.emplace_back(p31a, p31b, zero);

	// I use a priority queue to pick the closest cone
	priority_queue<cone_d> work;
	work.push(c12);
	work.push(c23);
	work.push(c31);

	// we pick the closest cone to expand
	while(!work.empty() && triangles.size() < 2000) {
		auto c = work.top().c;
		work.pop();

		// one quadrant is enough for me
		if((c.p1.x < -4 && c.p2.x < -4) || (c.p1.y < -4 && c.p2.y < -4))
			continue;

		// split or not?
		if (lensqr(c.p1, c.p2) > 16) {
			// choose midpoint on c.p2 - c.p1
			// such that both are length >= 2
			// choose direction such that we get two cones
			// The two pow(-) functions skew the distribution,
			// you can artistically choose something here.
			const auto l = sqrt(lensqr(c.p1, c.p2));
			const auto wiggle = l - 4;
			const auto r0 = unif01(gen);
			const auto r = 2.0 / l + pow(r0, 1.1) * wiggle / l;
			const auto pm = (1.0 - r) * c.p1 + r * c.p2;

			const auto r20 = unif01(gen);
			const auto r2 = 0.001 + 0.998 * pow(r20, 1.1);
			const auto dm = r2 * c.d1 + (1.0 - r2) * c.d2;

			cone c1{c.p1, c.d1, pm, dm};
			cone c2{pm, dm, c.p2, c.d2};
			work.push(c1);
			work.push(c2);
		} else {
			// try the two triangles and see which one fits the bill		
			const auto x1 = unit_triangle(c.p2, c.p1, c.d1);
			const auto x2 = unit_triangle(c.p1, c.p2, c.d2);

			const auto p1b = c.p1 + x1 * c.d1;
			const auto p2b = c.p2 + x2 * c.d2;
			
			// Here I do not follow the paper correctly. I should pick the
			// new cone with angles at most bla bla bla. Instead, I choose
			// the one with a shorter edge, sounds fine to me. I reintroduce
			// some randomness for artistic reasons. Increase 20 for wild
			// results.
			if (lensqr(c.p1, p2b) + 34*unif01(gen) < lensqr(p1b, c.p2) + 34*unif01(gen)) {
				cone c2{c.p1, c.d1, p2b, c.d2};
				work.push(c2);
				triangles.emplace_back(c.p1, p2b, c.p2);
			} else {
				cone c2{p1b, c.d1, c.p2, c.d2};
				work.push(c2);
				triangles.emplace_back(c.p1, p1b, c.p2);
			}
		}
	}

	cout << "<svg>\n";
	for(auto && t : triangles) {
		auto p1 = get<0>(t), p2 = get<1>(t), p3 = get<2>(t);
		vector<ld> lengths{{norm(p2 - p1), norm(p3 - p2), norm(p1 - p3)}};
		sort(lengths.begin(), lengths.end());
		auto r = max(0.0l, min(255.0l, 70 * lengths[2]));
		auto g = max(0.0l, min(255.0l, 70 * lengths[1]));
		auto b = max(0.0l, min(255.0l, 256 - 100 * lengths[0]));
		cout << "<polygon";
		cout << " points=\"" << get<0>(t) << ' ' << get<1>(t) << ' ' << get<2>(t) << "\"";
		cout << " style=\"" << "stroke-width:0.01;"
		     << "fill:rgb(" << r << ", " << g << ", " << b << ");"
		     << "stroke:rgb(" << max(0.0l, r-100) << ", " << max(0.0l, g-100) << ", " << max(0.0l, b-100) << ");" << "\"";
		cout << "/>\n";
	}
	cout << "</svg>" << endl;
}