commit 4ceb2dbc7067b26a0e16e74d827fda5ed3e2d72e Author: Joshua Moerman Date: Wed Jan 3 17:08:35 2018 +0100 Initial commit diff --git a/README.md b/README.md new file mode 100644 index 0000000..85c5960 --- /dev/null +++ b/README.md @@ -0,0 +1,7 @@ +Tilings of the plane with triangles ... +======================================= + +... of equal area, bounded perimenter and all non-congruent. Here, I try to +implement the procedure from a paper by Andrey Kupavskii, János Pach and Gábor +Tardos. It is not correct. I interpret generic as random, is probably fine for +this gimmick. [ArXiv article](https://arxiv.org/pdf/1712.03118.pdf) diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..3a6958b --- /dev/null +++ b/main.cpp @@ -0,0 +1,214 @@ +#include +#include +#include +#include + +using namespace std; +using ld = long double; + +struct pt { + ld x = 0, y = 0; +}; + +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 dist(0.8, 1.6); + + uniform_real_distribution 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 lines in the construction + vector> rays; + rays.emplace_back(zero, c1); + rays.emplace_back(zero, c2); + rays.emplace_back(zero, c3); + + vector> lines; + lines.emplace_back(p12a, p12b); + lines.emplace_back(p23a, p23b); + lines.emplace_back(p31a, p31b); + + // I use a priority queue to pick the closest cone + priority_queue work; + work.push(c12); + work.push(c23); + work.push(c31); + + // we pick the closest cone to expand + while(!work.empty() && lines.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.5) * 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.2); + const auto dm = r2 * c.d1 + (1.0 - r2) * c.d2; + + rays.emplace_back(pm, dm); + + 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) + 20*unif01(gen) < lensqr(p1b, c.p2) + 20*unif01(gen)) { + cone c2{c.p1, c.d1, p2b, c.d2}; + work.push(c2); + lines.emplace_back(c.p1, p2b); + } else { + cone c2{p1b, c.d1, c.p2, c.d2}; + work.push(c2); + lines.emplace_back(p1b, c.p2); + } + } + } + + // the print an svg attribute + const auto p = [](auto str, auto v) { + cout << ' ' << str << "=\"" << v << "\""; + }; + + // I roughly estimate the maximal distance from zero + ld max_l = 0; + + shuffle(lines.begin(), lines.end(), gen); + + cout << "\n"; + for(auto && l : lines) { + max_l = 0.5 * max_l + 0.5 * max(max_l, max(norm(l.first), norm(l.second))); + cout << "\n"; + } + for(auto && l : rays) { + const auto p1 = l.first; + const auto l_left = max(ld(1.0), max_l - norm(p1)); + const auto p2 = p1 + (l_left / norm(l.second)) * l.second; + cout << "\n"; + } + cout << "" << endl; +} diff --git a/render1.png b/render1.png new file mode 100644 index 0000000..7b5f24f Binary files /dev/null and b/render1.png differ diff --git a/test.svg b/test.svg new file mode 100644 index 0000000..f30235f --- /dev/null +++ b/test.svg @@ -0,0 +1,2073 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +