From 1b7f4b869ed828993b275340e1e6b8cde3a4b08c Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Sun, 19 Mar 2023 12:21:05 +0100 Subject: [PATCH] new Conformal Square projection (with transition and shader and Euclidean form) --- classes.cpp | 7 +++-- drawing.cpp | 56 ++++++++++++++++++++++++++++++---- hypgraph.cpp | 85 ++++++++++++++++++++++++++++++++++++++++++++++++++++ shaders.cpp | 68 ++++++++++++++++++++++++++++++++++++++++- 4 files changed, 206 insertions(+), 10 deletions(-) diff --git a/classes.cpp b/classes.cpp index 1c2af09d..79458e1c 100644 --- a/classes.cpp +++ b/classes.cpp @@ -1014,9 +1014,9 @@ enum eModel : int { mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdMiller, mdGallStereographic, mdWinkelTripel, // 39..48 mdPoorMan, mdPanini, mdRetroCraig, mdRetroLittrow, mdRetroHammer, mdThreePoint, mdLiePerspective, mdLieOrthogonal, mdRelPerspective, mdRelOrthogonal, - // 49 - mdHorocyclicEqa, - // 50.. + // 49..50 + mdHorocyclicEqa, mdConformalSquare, + // 51.. mdGUARD, mdPixel, mdHyperboloidFlat, mdPolynomial, mdManual }; #endif @@ -1078,6 +1078,7 @@ EX vector mdinf = { {X3("relativistic perspective"), mf::euc_boring | mf::perspective}, {X3("relativistic orthogonal"), mf::euc_boring}, {X3("horocyclic equal-area"), mf::euc_boring | mf::equiarea | mf::orientation | mf::horocyclic}, + {X3("conformal square"), mf::orientation | mf::broken | mf::transition}, {X3("guard"), mf::technical}, {X3("pixel"), mf::technical}, {X3("hypflat"), mf::technical}, diff --git a/drawing.cpp b/drawing.cpp index 687b2ac4..e1523394 100644 --- a/drawing.cpp +++ b/drawing.cpp @@ -1622,10 +1622,14 @@ EX namespace ods { bool broken_projection(dqi_poly& p0) { int broken_coord = models::get_broken_coord(pmodel); static bool in_broken = false; + bool both_broken = pmodel == mdConformalSquare; if(broken_coord && !in_broken) { int zcoord = broken_coord; int ycoord = 3 - zcoord; + int xcoord = 0; + + zcoord = 2; vector all; for(int i=0; i ib(in_broken, true); - ld part = ilerp(all[last_fail][0], all[last_fail+1][0], 0); + + ld part = ilerp(all[last_fail][xcoord], all[last_fail+1][xcoord], 0); + if(both_broken && all[last_fail][ycoord] * all[last_fail+1][ycoord] < 0) { + ld part2 = ilerp(all[last_fail][ycoord], all[last_fail+1][ycoord], 0); + if(part2 > part) part = part2, swap(xcoord, ycoord); + } hyperpoint initial = normalize(lerp(all[last_fail], all[last_fail+1], 1 - (1-part) * .99)); bool have_initial = true; v.push_back(glhr::pointtogl(initial)); @@ -1663,19 +1680,46 @@ bool broken_projection(dqi_poly& p0) { int at = last_fail; do { v.push_back(glhr::pointtogl(all[at])); - if(at == p0.cnt-1 && all[at] != all[0]) { + if(at == p0.cnt-1 && sqhypot_d(2, all[at] - all[0]) > 1e-6) { p.cnt = isize(v); p.draw(); v.clear(); at = 0; have_initial = false; } int next = at+1; if(next == p0.cnt) next = 0; if(break_in(all[at], all[next])) { - ld part = ilerp(all[at][0], all[next][0], 0); + ld part = ilerp(all[at][xcoord], all[next][xcoord], 0); + if(both_broken && all[at][ycoord] * all[next][ycoord] < 0) { + ld part2 = ilerp(all[at][ycoord], all[next][ycoord], 0); + if(part2 < part) part = part2, swap(xcoord, ycoord); + } + hyperpoint final = normalize(lerp(all[at], all[next], part * .99)); v.push_back(glhr::pointtogl(final)); if(have_initial) { int max = 4 << vid.linequality; - if(final[0] * initial[0] > 0) { + if(both_broken) { + auto square_close_corner = [] (hyperpoint h) { + hyperpoint end = -C0; + if(abs(h[0]) > abs(h[1])) + end[0] = 0.01 * signum(h[0]), end[1] = 0.001 * signum(h[1]); + else + end[1] = 0.01 * signum(h[1]), end[0] = 0.001 * signum(h[0]); + return normalize(end); + }; + hyperpoint endf = square_close_corner(final); + hyperpoint endi = square_close_corner(initial); + if(endf != endi) { + for(int i=1; i<=max; i++) + v.push_back(glhr::pointtogl(lerp(final, endf, i * 1. / max))); + for(int i=0; i<=max; i++) + v.push_back(glhr::pointtogl(lerp(endi, initial, i * 1. / max))); + } + else { + for(int i=1; i<=max; i++) + v.push_back(glhr::pointtogl(lerp(final, initial, i * 1. / max))); + } + } + else if(final[xcoord] * initial[xcoord] > 0) { for(int i=1; i<=max; i++) v.push_back(glhr::pointtogl(lerp(final, initial, i * 1. / max))); } diff --git a/hypgraph.cpp b/hypgraph.cpp index bc891857..07a6ef95 100644 --- a/hypgraph.cpp +++ b/hypgraph.cpp @@ -528,6 +528,81 @@ EX void make_axial(hyperpoint H, hyperpoint& ret, const hr::function eps || abs(dely) > eps || abs(delz) > eps) { + ld sx = sqrt(x); + ld sy = sqrt(y); + ld sz = sqrt(z); + ld len = sx * (sy+sz) + sy * sz; + x = .25 * (x+len); + y = .25 * (y+len); + z = .25 * (z+len); + mean = (x+y+z)/3; + delx = (mean-x) / mean; + dely = (mean-y) / mean; + delz = (mean-z) / mean; + } + ld e2 = delx * dely - delz * delz; + ld e3 = delx * dely * delz; + return ((1.0 + (e2 / 24.0 - 0.1 - 3.0 * e3 / 44.0) * e2+ e3 / 14) / sqrt(mean)); + } + +ld ellFaux(ld cos_phi, ld sin_phi, ld k) { + ld x = cos_phi * cos_phi; + ld y = 1 - k * k * sin_phi * sin_phi; + return sin_phi * ellRF(x, y, 1); + } + +ld sqrt_clamp(ld x) { if(x<0) return 0; return sqrt(x); } + +hyperpoint to_square(hyperpoint H) { + + ld d = hypot_d(2, H); + ld x = d / (H[2] + 1); + x *= pconf.model_transition; + + ld cos_phiosqrt2 = sqrt(2) / (x + 1/x); + ld cos_lambda = -H[1] / d; + ld sin_lambda = H[0] / d; + ld cos_a = cos_phiosqrt2 * (sin_lambda + cos_lambda); + ld cos_b = cos_phiosqrt2 * (sin_lambda - cos_lambda); + ld sin_a = sqrt(1 - cos_a * cos_a); + ld sin_b = sqrt(1 - cos_b * cos_b); + ld cos_a_cos_b = cos_a * cos_b; + ld sin_a_sin_b = sin_a * sin_b; + ld sin2_m = 1.0 + cos_a_cos_b - sin_a_sin_b; + ld sin2_n = 1.0 - cos_a_cos_b - sin_a_sin_b; + ld sin_m = sqrt_clamp(sin2_m); + ld cos_m = sqrt_clamp(1 - sin2_m); + if(sin_lambda < 0) sin_m = -sin_m; + ld sin_n = sqrt_clamp(sin2_n); + ld cos_n = sqrt_clamp(1.0 - sin2_n); + if(cos_lambda > 0.0) sin_n = -sin_n; + + hyperpoint res; + ld divby = 0.53935260118837935472; + res[0] = ellFaux(cos_m,sin_m,sqrt(2)/2.) * divby; + res[1] = ellFaux(cos_n,sin_n,sqrt(2)/2.) * divby; + res[2] = 0; res[3] = 1; + + if(x > 1) { + if(abs(res[0]) > abs(res[1])) { + if(res[0] > 0) res[0] = 2 - res[0]; else res[0] = -2 - res[0]; + } + else { + if(res[1] > 0) res[1] = 2 - res[1]; else res[1] = -2 - res[1]; + } + } + + res /= pconf.model_transition; + return res; + } + EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) { hyperpoint H = H_orig.h; @@ -810,6 +885,16 @@ EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) { break; } + case mdConformalSquare: { + find_zlev(H); + models::apply_orientation_yz(H[1], H[2]); + models::apply_orientation(H[0], H[1]); + ret = to_square(H); + models::apply_orientation(ret[1], ret[0]); + models::apply_orientation_yz(ret[2], ret[1]); + break; + } + case mdLieOrthogonal: { ret = lie_log_correct(H_orig, H); diff --git a/shaders.cpp b/shaders.cpp index d9bdfad4..4d7427c8 100644 --- a/shaders.cpp +++ b/shaders.cpp @@ -276,10 +276,15 @@ shared_ptr write_shader(flagtype shader_flags) { } else if(pmodel == mdDisk && GDIM == 3 && !spherespecial && !nonisotropic && !gproduct) { coordinator += "t /= (t[3] + uAlpha);\n"; - vsh += "uniform mediump float uAlpha;"; + vsh += "uniform mediump float uAlpha;\n"; shader_flags |= SF_DIRECT | SF_BOX | SF_ZFOG; treset = true; } + else if(pmodel == mdConformalSquare && pconf.model_transition == 1) { + shader_flags |= SF_ORIENT | SF_DIRECT; + coordinator += "t = uPP * t;", vsh += "uniform mediump mat4 uPP;"; + coordinator += "t = to_square(t);"; + } else if(pmodel == mdBand && hyperbolic) { shader_flags |= SF_BAND | SF_ORIENT | SF_BOX | SF_DIRECT; coordinator += "t = uPP * t;", vsh += "uniform mediump mat4 uPP;"; @@ -848,6 +853,67 @@ EX void add_if(string& shader, const string& seek, const string& function) { EX void add_fixed_functions(string& shader) { /* from the most complex to the simplest */ + add_if(shader, "to_square", + "mediump vec4 to_square(mediump vec4 h) {\n" + "float d = length(h.xy);\n" + "float x = d / (h.z + 1.);\n" + + "float cos_phiosqrt2 = sqrt(2.) / (x + 1./x);\n" + "float cos_lambda = -h.y / d;\n" + "float sin_lambda = h.x / d;\n" + "float cos_a = cos_phiosqrt2 * (sin_lambda + cos_lambda);\n" + "float cos_b = cos_phiosqrt2 * (sin_lambda - cos_lambda);\n" + "float sin_a = sqrt(1. - cos_a * cos_a);\n" + "float sin_b = sqrt(1. - cos_b * cos_b);\n" + "float cos_a_cos_b = cos_a * cos_b;\n" + "float sin_a_sin_b = sin_a * sin_b;\n" + "float sin2_m = 1.0 + cos_a_cos_b - sin_a_sin_b;\n" + "float sin2_n = 1.0 - cos_a_cos_b - sin_a_sin_b;\n" + "float sin_m = sqrt_clamp(sin2_m);\n" + "float cos_m = sqrt_clamp(1. - sin2_m);\n" + "if(sin_lambda < 0.) sin_m = -sin_m;\n" + "float sin_n = sqrt_clamp(sin2_n);\n" + "float cos_n = sqrt_clamp(1.0 - sin2_n);\n" + "if(cos_lambda > 0.0) sin_n = -sin_n;\n" + "#define divby 0.53935260118837935472\n" + "vec4 res = vec4(ellFaux(cos_m,sin_m,sqrt(2.)/2.) * divby, ellFaux(cos_n,sin_n,sqrt(2.)/2.) * divby, 0, 1);\n" + "if(x > 1.) {\n" + " if(abs(res[0]) > abs(res[1])) {\n" + " if(res[0] > 0.) res[0] = 2. - res[0]; else res[0] = -2. - res[0];\n" + " }\n" + " else {\n" + " if(res[1] > 0.) res[1] = 2. - res[1]; else res[1] = -2. - res[1];\n" + " }\n" + " }\n" + "return res;\n" + "}\n"); + + add_if(shader, "sqrt_clamp", "mediump float sqrt_clamp(mediump float x) { return x >= 0. ? sqrt(x) : 0.; }\n"); + add_if(shader, "ellFaux", "mediump float ellFaux(mediump float cos_phi, mediump float sin_phi, mediump float k) {\n" + "return sin_phi * ellRF(cos_phi * cos_phi, 1. - k * k * sin_phi * sin_phi, 1.);\n" + "}\n"); + add_if(shader, "ellRF", "mediump float ellRF(mediump float x, mediump float y, mediump float z) {\n" + "float delx = 1., dely = 1., delz = 1.;\n" + "const float eps = 0.0025;\n" + "float mean;\n" + "while(abs(delx) > eps || abs(dely) > eps || abs(delz) > eps) {\n" + " float sx = sqrt(x);\n" + " float sy = sqrt(y);\n" + " float sz = sqrt(z);\n" + " float len = sx * (sy+sz) + sy * sz;\n" + " float x = .25 * (x+len);\n" + " float y = .25 * (y+len);\n" + " float z = .25 * (z+len);\n" + " mean = (x+y+z)/3.;\n" + " delx = (mean-x) / mean;\n" + " dely = (mean-y) / mean;\n" + " delz = (mean-z) / mean;\n" + " }\n" + "float e2 = delx * dely - delz * delz;\n" + "float e3 = delx * dely * delz;\n" + "return ((1.0 + (e2 / 24.0 - 0.1 - 3.0 * e3 / 44.0) * e2+ e3 / 14.) / sqrt(mean));\n" + "}\n"); + add_if(shader, "tanh", "mediump float tanh(mediump float x) { return sinh(x) / cosh(x); }\n"); add_if(shader, "sinh", "mediump float sinh(mediump float x) { return (exp(x) - exp(-x)) / 2.0; }\n"); add_if(shader, "asin_clamp", "mediump float asin_clamp(mediump float x) { return x > 1. ? PI/2. : x < -1. ? -PI/2. : asin(x); }\n");