From 410cbd7cb310b37adb11d522ea29430cdae59ed8 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Fri, 25 Apr 2025 14:28:07 +0200
Subject: [PATCH] goldberg-sub was missing...

---
 goldberg-sub.cpp | 251 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 251 insertions(+)
 create mode 100644 goldberg-sub.cpp

diff --git a/goldberg-sub.cpp b/goldberg-sub.cpp
new file mode 100644
index 00000000..fa4caa78
--- /dev/null
+++ b/goldberg-sub.cpp
@@ -0,0 +1,251 @@
+#include "hyper.h"
+
+namespace hr {
+
+EX namespace gp {
+
+#if HDR
+enum class subdivision { linear, equalarea, conformal, gnomonic_vertex, gnomonic_face, equaldist };
+#endif
+
+EX subdivision su = subdivision::conformal;
+
+pair<ld, ld> compute_all_areas(bool recolor) {
+  map<ld, int> qty_of;
+  map<ld, int> qty_of_standard;
+  ld totar = 0;
+
+  map<cell*, ld> areas;
+
+  int gt = S3 == 3 ? 6 : 4;
+
+  for(cell *c: currentmap->allcells()) {
+    vector<hyperpoint> hs;
+    int t = c->type;
+    for(int i=0; i<t; i++) hs.push_back(currentmap->get_corner(c, i));
+    hs.push_back(hs[0]);
+
+    ld area = abs(compute_area(hs));
+    totar += area;
+
+    qty_of[area]++;
+    areas[c] = area;
+    if(t == gt) qty_of_standard[area]++;
+    }
+
+  if(qty_of_standard.empty()) qty_of_standard = qty_of;
+  if(qty_of.empty()) return {0, 0};
+
+  auto e = qty_of.end(); e--;
+  ld mina = qty_of.begin()->first;
+  ld maxa = e->first;
+
+  auto e1 = qty_of_standard.end(); e1--;
+  ld mina1 = qty_of_standard.begin()->first;
+  ld maxa1 = e1->first;
+
+  if(recolor) for(auto [c, a]: areas) if(c->land == laCanvas) c->landparam = gradient(0xFF0000, 0x0000FF, mina, a, maxa);
+
+  return {maxa/mina, maxa1/mina1};
+  }
+
+constexpr int qval = 8;
+
+array<ld, qval> vals;
+
+bool cornmul(const transmatrix& corners, const hyperpoint& c, hyperpoint& h) {
+
+  if(su == subdivision::gnomonic_face) {
+    auto co1 = kleinize(get_column(corners, 1));
+    auto co2 = kleinize(get_column(corners, 2));
+    h = c[0] * C0 + c[1] * co1 + c[2] * co2;
+    return true;
+    }
+
+  if(among(su, subdivision::gnomonic_vertex, subdivision::equaldist)) {
+    // to do: for gp_style = false */
+    int sw = c[2] > c[1] ? -1 : 1;
+
+    auto x = 1 - c[1] - c[2];
+    auto y = 1 - (c[1] - c[2]) * sw;
+
+    auto co1 = get_column(corners, 1);
+    auto co2 = get_column(corners, 2);
+    auto& co = sw == 1 ? co1 : co2;
+    // println(hlog, tuple(intval(co1, Hypc), intval(co2, Hypc), hdist0(co1), hdist0(co2)));
+    auto T = gpushxto0(co);
+    auto dx = kleinize(T * mid(co1, co2));
+    auto dxy = kleinize(T * get_column(corners, 0));
+
+    if(su == subdivision::equaldist) {
+      ld lx = hypot_d(2, dx), ly = hypot_d(2, dxy-dx);
+      x = tan_auto(x * atan_auto(lx)) / lx;
+      y = tan_auto(y * atan_auto(ly)) / ly;
+      }
+
+    hyperpoint hres = C0 + (dx-C0) * y + (dxy-dx) * x;
+
+    h = iso_inverse(T) * hres;
+    return true;
+    }
+
+  hyperpoint u;
+  u[0] = c[0] * (vals[5] + vals[0] * (c[2] + c[1]) + vals[1] * c[1] * c[2]);
+  u[1] = c[1] * (vals[6] + vals[2] * c[2] + vals[3] * c[0] + vals[4] * c[0] * c[2]);
+  u[2] = c[2] * (vals[7] + vals[2] * c[1] + vals[3] * c[0] + vals[4] * c[0] * c[1]);
+  u[3] = 0;
+
+  h = corners * u;
+  return true;
+  }
+
+hyperpoint compute(hyperpoint c, const transmatrix& T) {
+  hyperpoint h;
+  cornmul(T, c, h);
+  return normalize(h);
+  }
+
+ld get_ratio(const transmatrix& T) {
+  ld min_area = 1e6, max_area = 0;
+  for(int a=0; a<=10; a++)
+  for(int b=0; b<=10-a; b++) {
+    hyperpoint h;
+    h[1] = a/10.; h[2] = b/10.; h[0] = 1-h[1]-h[2];
+    ld eps = 1e-3;
+    auto r0 = compute(h, T);
+    h[0] += eps; h[2] -= eps;
+    auto r1 = compute(h, T);
+    h[0] -= eps; h[1] += eps;
+    auto r2 = compute(h, T);
+    h[1] -= eps; h[2] += eps;
+    if(hyperbolic) { transmatrix T = gpushxto0(r0); r0 = C0; r1 = T * r1; r2 = T * r2; }
+    ld area = hypot_d(3, (r1-r0) ^ (r2-r0));
+    if(a == 0 && b == 0) area *= 1;
+    if(area < min_area) min_area = area;
+    if(area > max_area) max_area = area;
+    }
+  return max_area / min_area;
+  }
+
+ld get_conf_error(const transmatrix& T) {
+  ld err = 0;
+  int tot = 0;
+  for(int a=0; a<=10; a++)
+  for(int b=0; b<=10-a; b++) if(a+b) {
+    hyperpoint h;
+    h[1] = a/10.; h[2] = b/10.; h[0] = 1-h[1]-h[2];
+    ld eps = 1e-3;
+    auto r0 = compute(h, T);
+    h[1] += eps; h[0] -= eps;
+    auto r1 = compute(h, T);
+    h[1] -= eps; h[0] += eps;
+
+    
+    if(S3 == 4) { h[2] -= eps; h[0] += eps; }
+    else { ld v = 2 / sqrt(3); h[2] += eps * v; h[1] -= eps * v/2; h[0] -= eps * v/2; }
+    auto r2 = compute(h, T);
+
+    transmatrix T = gpushxto0(r0); r0 = C0; r1 = T * r1 - r0; r2 = T * r2 - r0;
+
+    ld ax = r1[0], ay = r1[1];
+    ld bx = r2[0], by = r2[1];
+
+    ld alpha = atan2(-(2*ax*bx+2*ay*by), (-ax*ax-ay*ay+bx*bx+by*by)) / 2;
+
+    ld d1 = hypot(ax * cos(alpha) + bx * sin(alpha), ay * cos(alpha) + by * sin(alpha));
+    ld d2 = hypot(ax * sin(alpha) - bx * cos(alpha), ay * sin(alpha) - by * cos(alpha));
+
+    tot++;
+
+    err += (d1/d2) + (d2/d1) - 2; // exit(1);
+    }
+  return err / tot;
+  }
+
+ld zerofun(const transmatrix& T) { return 0; }
+
+void init_cornmul(const transmatrix& T) {
+
+  for(auto& v: vals) v = 0;
+  if(qval == 8) vals[5] = vals[6] = vals[7] = 1;
+  else vals[0] = vals[1] = vals[2] = 1;
+
+  auto chk = zerofun;
+  if(su == subdivision::conformal) chk = get_conf_error;
+  if(su == subdivision::equalarea) chk = get_ratio;
+
+  auto cur = chk(T);
+
+  for(int it=0; it<500; it++) {
+    array<ld, qval> gra;
+    ld sq = 0;
+    for(int i=0; i<qval; i++) {
+      vals[i] += 1e-4;
+      gra[i] = cur - chk(T);
+      vals[i] -= 1e-4;
+      sq += gra[i] * gra[i];
+      }
+    if(!sq) break;
+    sq = 1e-6 / sqrt(sq);
+
+    while(true) {
+      for(int i=0; i<qval; i++) vals[i] += sq * gra[i];
+      auto cur1 = chk(T);
+      if(cur1 < cur) { cur = cur1; sq *= 2; }
+      else { for(int i=0; i<qval; i++) vals[i] -= sq * gra[i]; break; }
+      }
+    }
+  }
+
+void config_gpsubs() {
+  addHook(gp::hooks_init_cornmul, 100, init_cornmul);
+  // not sure why C++ could not deduce the template parameter
+  using t = bool(const transmatrix&, const hyperpoint&, hyperpoint&);
+  addHook<t,t> (gp::hooks_cornmul, 100, cornmul);
+
+  param_enum(su, parameter_names("gp_sub", "gp_sub"), subdivision::equalarea)
+    ->editable({{"linear", "straight lines remain straight"}, {"equal-area", "try to have tiles of equal area (tiles on corners will still be different)"}, {"conformal", "try to get regularly shaped tiles"}, {"gnomonic (vertex)", "regular in gnomonic projection (centered on pure vertex)"}, {"gnomonic (face)", "regular in gnomonic projection (centered on pure face)"}, {"tangent adjustment", "equal distances on the original edges"}}, "Goldberg mapping", 'S')
+    ->add_extra([] {
+       auto p = compute_all_areas(false);
+       dialog::addSelItem(XLAT("cell area ratio (max/min)"), fts(p.first), 'X');
+       dialog::add_action([] { compute_all_areas(true); });
+       dialog::addSelItem(S3 == 3 ? XLAT("hex only") : XLAT("square only"), fts(p.second), 'Y');
+       dialog::add_action([] { compute_all_areas(true); });
+       })
+    ->reaction = [] { if(game_active) init_cornmul(dir_matrix(0)); };
+  }
+
+auto hookc = addHook(hooks_configfile, 100, config_gpsubs) + arg::add3("-gbs-debug", [] {
+  // should have more configuration, but no need for it at the moment
+  addHook(hooks_frame, 100, [] {
+    vid.linewidth *= 3;
+    for(int a=0; a<S7; a++) {
+      transmatrix T = dir_matrix(a);
+      shiftmatrix S = ggmatrix(cwt.at);
+      int z = 5;
+      for(int ix=0; ix<z; ix++)
+      for(int iy=0; iy<z-ix; iy++) {
+        ld x = ix * 1./z, y = iy*1./z, x1 = (ix+1.)/z, y1 = (iy+1.)/z;
+        queueline(S * normalize(cornmul(T, point3(1-x-y,x,y))), S * normalize(cornmul(T, point3(1-x1-y,x1,y))), 0x80FF80FF);
+        queueline(S * normalize(cornmul(T, point3(1-x-y,x,y))), S * normalize(cornmul(T, point3(1-x-y1,x,y1))), 0x80FF80FF);
+        queueline(S * normalize(cornmul(T, point3(1-x-y,x,y))), S * normalize(cornmul(T, point3(1-x1-y1,x1,y1))), 0xFF80FFFF);
+        queueline(S * normalize(cornmul(T, point3(1-x1-y,x1,y))), S * normalize(cornmul(T, point3(1-x-y1,x,y1))), 0xFF80FFFF);
+
+        hyperpoint c = point3(1-x-y,x,y);
+        hyperpoint h = normalize(cornmul(T, c));
+        
+        if(a == 0) {
+          int sw = c[2] > c[1] ? -1 : 1;
+          auto x = 1 - c[1] - c[2];
+          auto y = 1 - (c[1] - c[2]) * sw;
+          queuestr(S * rgpushxto0(h), 0.4, format("%.2f %.2f", x, y), 0xFFFFFFFF, 2);
+          }
+        }
+      }
+    vid.linewidth /= 3;
+    });
+  });
+    
+
+}
+}
\ No newline at end of file