diff --git a/rug.cpp b/rug.cpp index cee16bd7..675e46a8 100644 --- a/rug.cpp +++ b/rug.cpp @@ -42,8 +42,8 @@ bool computed = false; bool valid; bool inqueue; double dist; - hyperpoint h; // point in the represented space - hyperpoint flat; // point in the native space, in azeq + hyperpoint h; // point in the represented space + hyperpoint native; // point in the native space hyperpoint precompute; vector edges; vector anticusp_edges; @@ -52,8 +52,8 @@ bool computed = false; rugpoint *getglue() { return glue ? (glue = glue->getglue()) : this; } - hyperpoint& glueflat() { - return glue->flat; + hyperpoint& gluenative() { + return glue->native; } rugpoint() { glue = NULL; } void glueto(rugpoint *x) { @@ -91,9 +91,14 @@ EX void subdivide(); EX ld modelscale = 1; EX ld model_distance = 4; -EX eGeometry gwhere = gEuclid; +const eGeometry rgHyperbolic = gSpace534; +const eGeometry rgEuclid = gCubeTiling; +const eGeometry rgSphere = gCell120; +const eGeometry rgElliptic = gECell120; -#define USING_NATIVE_GEOMETRY dynamicval gw(geometry, gwhere == gElliptic ? gSphere : gwhere) +EX eGeometry gwhere = rgEuclid; + +#define USING_NATIVE_GEOMETRY dynamicval gw(geometry, hr::rug::gwhere) // hypersian rug datatypes and globals //------------------------------------- @@ -121,114 +126,24 @@ EX rugpoint *finger_center; EX ld finger_range = .1; EX ld finger_force = 1; -EX int rugdim; +// EX eModel rug_model = mdEquidistant; -EX bool rug_perspective = ISANDROID; - -// extra geometry functions -//-------------------------- - -// returns a matrix M -// such that inverse(M) * h1 = ( |h1|, 0, 0) and inverse(M) * h2 = ( .., .., 0) - -transmatrix orthonormalize(hyperpoint h1, hyperpoint h2) { - - hyperpoint vec[3] = {h1, h2, h1 ^ h2}; - - for(int i=0; i<3; i++) { - for(int j=0; j1e-4) { - println(hlog, "Error: h[2] = ", h[2]); - rug_failure = true; - } - if(euclid) { - h[2] = 1; - return h; - } - ld d = hypot(h[0], h[1]); - if(d == 0) { - h[2] = 1; - return h; - } - if(sphere) { - ld d0 = d ? d : 1; - h[0] = sin(d) * h[0]/d0; - h[1] = sin(d) * h[1]/d0; - h[2] = cos(d); - } - else { - ld d0 = d ? d : 1; - h[0] = sinh(d) * h[0]/d0; - h[1] = sinh(d) * h[1]/d0; - h[2] = cosh(d); - } - return h; - } - -hyperpoint hyperboloid_to_azeq(hyperpoint h) { - if(euclid) { - h[2] = 0; - return h; - } - else { - ld d = hdist0(h); - if(d == 0) { h[2] = 0; return h; } - ld d2 = hypot_d(2, h); - if(d2 == 0) { h[2] = 0; return h; } - h[0] = d * h[0] / d2; - h[1] = d * h[1] / d2; - h[2] = 0; - return h; - } - } - -struct normalizer { - - transmatrix M, Mi; - - dynamicval gw; - - normalizer (hyperpoint h1, hyperpoint h2) : gw(geometry, gwhere == gElliptic ? gSphere : gwhere) { - M = orthonormalize(h1, h2); - Mi = inverse(M); - } - - hyperpoint operator() (hyperpoint h) { return azeq_to_hyperboloid(Mi*h); } - hyperpoint operator[] (hyperpoint h) { return M*hyperboloid_to_azeq(h); } - }; +EX bool rug_perspective = false; void push_point(hyperpoint& h, int coord, ld val) { - if(fast_euclidean && gwhere == gEuclid) + USING_NATIVE_GEOMETRY; + if(fast_euclidean && euclid) h[coord] += val; else if(!val) return; else { - // if(zero_d(3, h)) { h[0] = 1e-9; h[1] = 1e-10; h[2] = 1e-11; } - normalizer n(hpxyz(coord==0,coord==1,coord==2), h); - hyperpoint f = n(h); - h = n[xpush(val) * f]; + h = cpush(coord, val) * h; } } EX void push_all_points(int coord, ld val) { if(!val) return; else for(int i=0; iflat, coord, val); + push_point(points[i]->native, coord, val); } // construct the graph @@ -236,6 +151,11 @@ EX void push_all_points(int coord, ld val) { int hyprand; +bool rug_euclid() { USING_NATIVE_GEOMETRY; return euclid; } +bool rug_hyperbolic() { USING_NATIVE_GEOMETRY; return hyperbolic; } +bool rug_sphere() { USING_NATIVE_GEOMETRY; return sphere; } +bool rug_elliptic() { USING_NATIVE_GEOMETRY; return elliptic; } + EX rugpoint *addRugpoint(hyperpoint h, double dist) { rugpoint *m = new rugpoint; m->h = h; @@ -269,73 +189,74 @@ EX rugpoint *addRugpoint(hyperpoint h, double dist) { if(z==0) z = 1; hpoint = hpoint * hpdist / z; - m->flat = hpxyz(hpoint[0], hpoint[1] * sin(a*2*M_PI), hpoint[1]*cos(a*2*M_PI)); + m->native = point31(hpoint[0], hpoint[1] * sin(a*2*M_PI), hpoint[1]*cos(a*2*M_PI)); } else if(sphere) { m->valid = good_shape = true; ld scale; - if(gwhere == gEuclid) { + USING_NATIVE_GEOMETRY; + if(euclid) { scale = modelscale; } - else if(gwhere == gNormal) { + else if(hyperbolic) { // sinh(scale) = modelscale scale = asinh(modelscale); } - else /* sphere/elliptic*/ { + else if(sphere) { if(modelscale >= 1) // do as good as we can... scale = M_PI / 2 - 1e-3, good_shape = false, m->valid = false; else scale = asin(modelscale); } - m->flat = h * scale; + else + scale = 1; + m->native = h * scale; + m->native = hpxy3(m->native[0], m->native[1], m->native[2]); } - else if(euclid && gwhere == gEuclid) { - m->flat = h * modelscale; + else if(euclid && rug_euclid()) { + m->native = h * modelscale; + m->native[2] = 0; + m->native[3] = 1; m->valid = good_shape = true; } - else if(gwhere == gNormal && (euclid || (hyperbolic && modelscale >= 1))) { + else if(rug_hyperbolic() && euclid) { + m->valid = good_shape = true; + USING_NATIVE_GEOMETRY; + m->native = tC0(parabolic13(h[0]*modelscale, h[1]*modelscale)); + } + + else if(rug_hyperbolic() && hyperbolic && modelscale >= 1) { m->valid = good_shape = true; - ld d = hdist0(h); - ld d0 = hypot_d(2, h); if(!d0) d0 = 1; + // radius of the equidistant + ld r = acosh(modelscale); + + h[3] = h[2]; h[2] = 0; - hyperpoint hpoint; - bool orig_euclid = euclid; USING_NATIVE_GEOMETRY; - if(orig_euclid) { - d *= modelscale; - // point on a horocycle going through C0, in distance d along the horocycle - hpoint = hpxy(d*d/2, d); - } - else { - // radius of the equidistant - ld r = acosh(modelscale); - // point on an equdistant going through C0 in distance d along the guiding line - // hpoint = hpxy(cosh(r) * sinh(r) * (cosh(d) - 1), sinh(d) * cosh(r)); - hpoint = xpush(r) * ypush(d) * xpush0(-r); - hpoint[0] = -hpoint[0]; - } - - ld hpdist = hdist0(hpoint); - ld z = hypot_d(2, hpoint); - if(z==0) z = 1; - m->flat = hpxyz(hpdist * h[0]/d0 * hpoint[1] / z, hpdist * h[1]/d0 * hpoint[1] / z, -hpdist * hpoint[0] / z); + m->native = rgpushxto0(h) * cpush0(2, r); } else { - m->flat = h; + m->native = h; ld hd = h[LDIM]; - for(int d=GDIM; dflat[d] = (hd - .99) * (rand() % 1000 - rand() % 1000) / 1000; + for(int d=GDIM; d<4; d++) { + m->native[d] = (hd - .99) * (rand() % 1000 - rand() % 1000) / 1000; + } + USING_NATIVE_GEOMETRY; + if(euclid) + m->native[3] = 1; + else + m->native = normalize(m->native); } if(rug_perspective) - push_point(m->flat, 2, -model_distance); + push_point(m->native, 2, -model_distance); - // if(rug_perspective && gwhere == gEuclid) m->flat[2] -= 3; + // if(rug_perspective && gwhere == gEuclid) m->native[2] -= 3; m->inqueue = false; m->dist = dist; points.push_back(m); @@ -343,8 +264,9 @@ EX rugpoint *addRugpoint(hyperpoint h, double dist) { } EX rugpoint *findRugpoint(hyperpoint h) { + USING_NATIVE_GEOMETRY; for(int i=0; ih - h) < 1e-5) return points[i]; + if(geo_dist_q(points[i]->h, h) < 1e-5) return points[i]; return NULL; } @@ -414,7 +336,7 @@ EX void sort_rug_points() { void calcLengths() { for(auto p: points) for(auto& edge: p->edges) - edge.len = hdist(p->h, edge.target->h) * modelscale; + edge.len = geo_dist_q(p->h, edge.target->h) * modelscale; } EX void calcparam_rug() { @@ -480,23 +402,28 @@ EX void buildTorusRug() { ld ax = alpha + 1.124651, bx = beta + 1.214893; ld x = xfactor * sin(ax), y = xfactor * cos(ax), z = yfactor * sin(bx), t = yfactor * cos(bx); - - ld d; - ld hxyz = sqrt(x*x+y*y+z*z); - - if(gwhere == gNormal) { - d = hxyz / (1-t) / mx; - d *= .95; - hyperpoint test = hpxy(d, 0); - test = perspective_to_space(test, 1, gcHyperbolic); - d = acosh(test[2]) / hxyz; - } - else { - d = (gwhere == gSphere) ? acos(t) / hxyz : modelscale * 3 / (1-t) / mx; - } - r->flat = r->h = hpxyz(x*d, y*d, z*d); + if(1) { + hyperpoint hp = hyperpoint(x, y, z, t); + USING_NATIVE_GEOMETRY; + + /* spherical coordinates are already good, otherwise... */ + + if(!sphere) { + /* stereographic projection to get Euclidean conformal torus */ + hp[3] += 1; + hp /= hp[3]; + hp /= mx; + hp[3] = 1; + } + + /* ... in H^3, use inverse Poincare to get hyperbolic conformal torus */ + + if(hyperbolic) + hp = perspective_to_space(hp, 1, gcHyperbolic); + } + r->valid = true; static const int X = 100003; // a prime @@ -562,14 +489,8 @@ EX void verify() { auto m2 = e.target; ld l = e.len; - ld l0; - if(fast_euclidean) l0 = sqhypot_d(rugdim, m->flat - m2->flat); - else { - normalizer n(m->flat, m2->flat); - hyperpoint h1 = n(m->flat); - hyperpoint h2 = n(m2->flat); - l0 = hdist(h1, h2); - } + USING_NATIVE_GEOMETRY; + ld l0 = geo_dist_q(m->native, m2->native); ratios.push_back(l0 / l); } @@ -668,17 +589,18 @@ EX void enqueue(rugpoint *m) { bool force_euclidean(rugpoint& m1, rugpoint& m2, double rd, bool is_anticusp = false, double d1=1, double d2=1) { if(!m1.valid || !m2.valid) return false; - // double rd = hdist(m1.h, m2.h) * xd; - double t = sqhypot_d(rugdim, m1.flat - m2.flat); + USING_NATIVE_GEOMETRY; + // double rd = geo_dist_q(m1.h, m2.h) * xd; + double t = sqhypot_d(3, m1.native - m2.native); if(is_anticusp && t > rd*rd) return false; t = sqrt(t); current_total_error += (t-rd) * (t-rd); bool nonzero = abs(t-rd) > err_zero_current; double force = (t - rd) / t / 2; // 20.0; - for(int i=0; i0) enqueue(&m2); } return nonzero; @@ -686,35 +608,42 @@ bool force_euclidean(rugpoint& m1, rugpoint& m2, double rd, bool is_anticusp = f bool force(rugpoint& m1, rugpoint& m2, double rd, bool is_anticusp=false, double d1=1, double d2=1) { if(!m1.valid || !m2.valid) return false; - if(gwhere == gEuclid && fast_euclidean) { + if(rug_euclid() && fast_euclidean) { return force_euclidean(m1, m2, rd, is_anticusp, d1, d2); } - normalizer n(m1.flat, m2.flat); - hyperpoint f1 = n(m1.flat); - hyperpoint f2 = n(m2.flat); + USING_NATIVE_GEOMETRY; - ld t = hdist(f1, f2); + ld t = geo_dist_q(m1.native, m2.native); if(is_anticusp && t > rd) return false; current_total_error += (t-rd) * (t-rd); bool nonzero = abs(t-rd) > err_zero_current; double forcev = (t - rd) / 2; // 20.0; - transmatrix T = gpushxto0(f1); - transmatrix T1 = spintox(T * f2) * T; + println(hlog, normalize(m1.native) - m1.native); + println(hlog, normalize(m1.native) - m2.native); + + transmatrix T = gpushxto0(m1.native); + transmatrix T1 = spintox(T * m2.native) * T; transmatrix iT1 = inverse(T1); - for(int i=0; i<3; i++) if(std::isnan(m1.flat[i])) { + println(hlog, sqhypot_d(4, iT1 * C0 - m1.native)); + println(hlog, sqhypot_d(4, iT1 * xpush0(t) - m1.native)); + + for(int i=0; i ", tie(m1.native, m2.native)); + if(nonzero && d2>0) enqueue(&m2); return nonzero; } @@ -724,8 +653,8 @@ vector > preset_points; EX void preset(rugpoint *m) { if(GDIM == 3) return; int q = 0; - hyperpoint h; - for(int i=0; i<3; i++) h[i] = 0; + USING_NATIVE_GEOMETRY; + hyperpoint h = Hypc; preset_points.clear(); @@ -757,10 +686,10 @@ EX void preset(rugpoint *m) { double ah = sqrt(a1*a1 - az*az + 1e-10); // c->h = a->h + (b->h-a->h) * cz + ch * ort - hyperpoint ort = (c->flat - a->flat - cz * (b->flat-a->flat)) / ch; + hyperpoint ort = (c->native - a->native - cz * (b->native-a->native)) / ch; // m->h = a->h + (b->h-a->h) * az - ah * ort - hyperpoint res = a->flat + (b->flat-a->flat) * az - ah * ort; + hyperpoint res = a->native + (b->native-a->native) * az - ah * ort; h += res; @@ -769,9 +698,9 @@ EX void preset(rugpoint *m) { } } - if(q>0) m->flat = h/q; - // printf("preset (%d) -> %s\n", q, display(m->flat)); - if(std::isnan(m->flat[0]) || std::isnan(m->flat[1]) || std::isnan(m->flat[2])) + if(q>0) m->native = normalize(h); + // printf("preset (%d) -> %s\n", q, display(m->native)); + if(std::isnan(m->native[0]) || std::isnan(m->native[1]) || std::isnan(m->native[2])) throw rug_exception(); } @@ -779,15 +708,8 @@ ld sse(const hyperpoint& h) { ld sse = 0; for(auto& p: preset_points) { ld l = p.first; - ld l0; - if(fast_euclidean) - l0 = hypot_d(rugdim, h - p.second->flat); - else { - normalizer n(h, p.second->flat); - hyperpoint h1 = n(h); - hyperpoint h2 = n(p.second->flat); - l0 = hdist(h1, h2); - } + USING_NATIVE_GEOMETRY; + ld l0 = geo_dist_q(h, p.second->native); sse += (l0-l) * (l0-l); } @@ -801,16 +723,17 @@ EX void optimize(rugpoint *m, bool do_preset) { // int ed0 = isize(preset_points); for(auto& e: m->edges) if(e.target->valid) preset_points.emplace_back(e.len, e.target); - if(gwhere >= gSphere || GDIM == 3) { - ld cur = sse(m->flat); + if(!rug_euclid()) { + ld cur = sse(m->native); for(int it=0; it<500; it++) { ld ex = exp(-it/60); again: - hyperpoint last = m->flat; - m->flat[(it/2)%rugdim] += ((it&1)?1:-1) * ex; - ld now = sse(m->flat); + hyperpoint last = m->native; + USING_NATIVE_GEOMETRY; + m->native = rgpushxto0(m->native) * cpush0((it/2)%3, (it&1)?ex:-ex); + ld now = sse(m->native); if(now < cur) { cur = now; ex *= 1.2; goto again; } - else m->flat = last; + else m->native = last; } } } @@ -860,10 +783,9 @@ EX void subdivide() { rugpoint *mm = addRugpoint(mid(m->h, m2->h), (m->dist+m2->dist)/2); halves[make_pair(m, m2)] = mm; if(!good_shape) { - normalizer n(m->flat, m2->flat); - hyperpoint h1 = n(m->flat); - hyperpoint h2 = n(m2->flat); - mm->flat = n[mid(h1, h2)]; + USING_NATIVE_GEOMETRY; + mm->native = mid(m->native, m2->native); + println(hlog, "mid of ", m->native, " and ", m2->native, " is ", mm->native); } mm->valid = m->valid && m2->valid; if(mm->valid) qvalid++; @@ -880,33 +802,9 @@ EX void subdivide() { println(hlog, "result ", make_tuple(isize(points), isize(triangles))); } -EX ld slow_modeldist(const hyperpoint& h1, const hyperpoint& h2) { - normalizer n(h1, h2); - hyperpoint f1 = n(h1); - hyperpoint f2 = n(h2); - return hdist(f1, f2); - } - -typedef array hyperpoint4; - -hyperpoint4 azeq_to_4(const hyperpoint& h) { - array res; - ld rad = hypot_d(3, h); - res[3] = cos(rad); - ld sr = sin(rad) / rad; - for(int j=0; j<3; j++) res[j] = h[j] * sr; - return res; - } - EX ld modeldist(const hyperpoint& h1, const hyperpoint& h2) { - if(gwhere == gSphere) { - hyperpoint4 coord[2] = { azeq_to_4(h1), azeq_to_4(h2) }; - ld edist = 0; - for(int j=0; j<4; j++) edist += sqr(coord[0][j] - coord[1][j]); - return 2 * asin(sqrt(edist) / 2); - } - - return slow_modeldist(h1, h2); + USING_NATIVE_GEOMETRY; + return geo_dist_q(h1, h2); } typedef long long bincode; @@ -925,8 +823,7 @@ bincode get_bincode(hyperpoint h) { case gcHyperbolic: return acd_bin(hypot_d(3, h)); case gcSphere: { - auto p = azeq_to_4(h); - return acd_bin(p[0]) + acd_bin(p[1]) * sY + acd_bin(p[2]) * sZ + acd_bin(p[3]) * sT; + return acd_bin(h[0]) + acd_bin(h[1]) * sY + acd_bin(h[2]) * sZ + acd_bin(h[3]) * sT; } } return 0; @@ -946,7 +843,7 @@ void generate_deltas(vector& target, int dim, bincode offset) { int detect_cusp_at(rugpoint *p, rugpoint *q) { if(hdist(p->h, q->h) * modelscale <= anticusp_dist) return 0; - else if(modeldist(p->flat, q->flat) > anticusp_dist - err_zero_current) + else if(modeldist(p->native, q->native) > anticusp_dist - err_zero_current) return 1; else { add_anticusp_edge(p, q); @@ -967,7 +864,7 @@ int detect_cusps() { map > code_to_point; for(auto p: points) if(p->valid) - code_to_point[get_bincode(p->flat)].push_back(p); + code_to_point[get_bincode(p->native)].push_back(p); vector deltas; generate_deltas(deltas, gwhere == gEuclid ? 3 : gwhere == gNormal ? 1 : 4, 0); @@ -1068,35 +965,16 @@ EX void physics() { bool use_precompute; -void getco(rugpoint *m, hyperpoint& h, int &spherepoints) { - h = use_precompute ? m->getglue()->precompute : m->getglue()->flat; - if(rug_perspective && gwhere >= gSphere) { - if(h[2] > 0) { - ld rad = hypot_d(3, h); - // turn M_PI to -M_PI - // the only difference between sphere and elliptic is here: - // in elliptic, we subtract PI from the distance - ld rad_to = (gwhere == gSphere ? M_PI + M_PI : M_PI) - rad; - ld r = -rad_to / rad; - h *= r; - spherepoints++; - } - } - } - extern int besti; #if CAP_ODS /* these functions are for the ODS projection, used in VR videos */ -bool project_ods(hyperpoint azeq, hyperpoint& h1, hyperpoint& h2, bool eye) { +bool project_ods(hyperpoint h, hyperpoint& h1, hyperpoint& h2, bool eye) { USING_NATIVE_GEOMETRY; - ld d = hypot_d(3, azeq); - ld sindbd = sin_auto(d)/d, cosd = cos_auto(d); - - return ods::project(hyperpoint(azeq[0] * sindbd, azeq[1] * sindbd, azeq[2] * sindbd, cosd), h1, h2, eye); - + return ods::project(h, h1, h2, eye); + // printf("%10.5lf %10.5lf %10.5lf ", azeq[0], azeq[1], azeq[2]); // printf(" => %10.5lf %10.5lf %10.5lf %10.5lf", x, y, z, t); @@ -1126,7 +1004,7 @@ void drawTriangle(triangle& t) { if(num == 2) return; for(int i=0; igetglue()->flat; + pts[i] = t.m[i]->getglue()->native; hyperpoint hc = (pts[1] - pts[0]) ^ (pts[2] - pts[0]); double hch = hypot_d(3, hc); @@ -1178,7 +1056,7 @@ void drawTriangle(triangle& t) { if(h[s+1][0] > +M_PI || h[s+2][0] > +M_PI) fst--; for(int x=fst; x<=lst; x++) for(int i=0; i<3; i++) { ct_array.emplace_back( - hpxyz(h[s+i][0] + 2*M_PI*x, h[s+i][1], h[s+i][2]), + point31(h[s+i][0] + 2*M_PI*x, h[s+i][1], h[s+i][2]), t.m[i]->x1, t.m[i]->y1, col); } @@ -1191,7 +1069,10 @@ void drawTriangle(triangle& t) { int spherepoints = 0; array h; - for(int i=0; inative; + // todo if(rug_elliptic() && h[i][3] > 0) { h[i] = -h[i]; spherepoints++; } + } if(spherepoints == 1 || spherepoints == 2) return; ld col = 1; @@ -1311,7 +1192,7 @@ EX void drawRugScene() { else { use_precompute = true; for(auto p: points) { - p->precompute = p->flat; + p->precompute = p->native; push_point(p->precompute, 0, vid.ipd*ed/2); } } @@ -1388,7 +1269,6 @@ EX void ensure_glbuf() { EX void reopen() { if(rugged) return; - rugdim = 2 * GDIM - 1; when_enabled = 0; GLERR("before init"); ensure_glbuf(); @@ -1476,7 +1356,7 @@ EX void apply_rotation(const transmatrix& t) { if(in_crystal()) crystal::apply_rotation(t); else #endif - for(auto p: points) p->flat = t * p->flat; + for(auto p: points) p->native = t * p->native; } EX void move_forward(ld distance) { @@ -1487,6 +1367,7 @@ EX void move_forward(ld distance) { #define CAP_HOLDKEYS (CAP_SDL && !ISWEB) EX bool handlekeys(int sym, int uni) { + USING_NATIVE_GEOMETRY; if(NUMBERKEY == '1') { ld bdist = 1e12; if(finger_center) @@ -1547,7 +1428,7 @@ EX bool handlekeys(int sym, int uni) { EX void finger_on(int coord, ld val) { for(auto p: points) { ld d = hdist(finger_center->h, p->getglue()->h); - push_point(p->flat, coord, val * finger_force * exp( - sqr(d / finger_range))); + push_point(p->native, coord, val * finger_force * exp( - sqr(d / finger_range))); } enqueue(finger_center), good_shape = false; } @@ -1646,6 +1527,7 @@ EX void actDraw() { push_all_points(1, strafey * ruggospeed); } else { + USING_NATIVE_GEOMETRY; if(finger_center) perform_finger(); else { @@ -1676,18 +1558,6 @@ EX void actDraw() { int besti; -void getco_pers(rugpoint *r, hyperpoint& p, int& spherepoints, bool& error) { - getco(r, p, spherepoints); - if(rug_perspective) { - if(p[2] >= 0) - error = true; - else { - p[0] /= p[2]; - p[1] /= p[2]; - } - } - } - static const ld RADAR_INF = 1e12; ld radar_distance = RADAR_INF; @@ -1709,9 +1579,22 @@ EX hyperpoint gethyper(ld x, ld y) { hyperpoint p0, p1, p2; bool error = false; int spherepoints = 0; - getco_pers(r0, p0, spherepoints, error); - getco_pers(r1, p1, spherepoints, error); - getco_pers(r2, p2, spherepoints, error); + if(1) { + USING_NATIVE_GEOMETRY; + // USING_RUG_PMODEL; + // dynamicval m(pmodel, mdEquidistant); + + if(elliptic && pmodel == mdDisk) { + int sp = + (r0->native[3] < 0) + (r1->native[3] < 0) + (r2->native[3] < 0); + if(sp == 1 || sp == 2) continue; + } + + applymodel(r0->native, p0); + applymodel(r1->native, p1); + applymodel(r2->native, p2); + } + if(error || spherepoints == 1 || spherepoints == 2) continue; double dx1 = p1[0] - p0[0]; double dy1 = p1[1] - p0[1]; @@ -1786,6 +1669,27 @@ void change_texturesize() { ld old_distance; +EX void rug_geometry_choice() { + cmode = sm::SIDE | sm::MAYDARK; + gamescreen(0); + dialog::init(XLAT("hypersian rug mode"), iinf[itPalace].color, 150, 100); + + dialog::addBoolItem("Euclidean", rug_euclid(), 'a'); + dialog::add_action([] { gwhere = rgEuclid; popScreen(); }); + + dialog::addBoolItem("hyperbolic", rug_hyperbolic(), 'b'); + dialog::add_action([] { gwhere = rgHyperbolic; popScreen(); }); + + dialog::addBoolItem("spherical", rug_sphere() && !rug_elliptic(), 'c'); + dialog::add_action([] { gwhere = rgSphere; popScreen(); }); + + dialog::addBoolItem("elliptic", rug_elliptic(), 'd'); + dialog::add_action([] { gwhere = rgElliptic; popScreen(); }); + + dialog::addBack(); + dialog::display(); + } + EX void show() { cmode = sm::SIDE | sm::MAYDARK; gamescreen(0); @@ -1808,7 +1712,7 @@ EX void show() { dialog::addSelItem(XLAT("model distance"), fts(model_distance), 'd'); dialog::addBoolItem(XLAT("projection"), rug_perspective, 'p'); dialog::lastItem().value = XLAT(rug_perspective ? "perspective" : - gwhere == gEuclid ? "orthogonal" : "azimuthal equidistant"); + rug_euclid() ? "orthogonal" : "azimuthal equidistant"); if(!rug::rugged) dialog::addSelItem(XLAT("native geometry"), geometry_name(gwhere), 'n'); else @@ -1830,11 +1734,10 @@ EX void show() { edit_levellines('L'); #if CAP_SURFACE - if(hyperbolic) { - if(gwhere == gEuclid) - dialog::addItem(XLAT("smooth surfaces"), 'c'); - else dialog::addBreak(100); - } + if(hyperbolic && rug_euclid()) + dialog::addItem(XLAT("smooth surfaces"), 'c'); + else + dialog::addBreak(100); #endif dialog::addBreak(50); @@ -1914,7 +1817,7 @@ EX void show() { if(!camera_center) push_all_points(2, model_distance); for(auto p:points) { if(adjust_edges) for(auto& e: p->edges) e.len *= modelscale / last; - if(adjust_points) p->flat *= modelscale / last; + if(adjust_points) p->native *= modelscale / last; enqueue(p); } if(adjust_distance) model_distance = model_distance * modelscale / last; @@ -1955,8 +1858,8 @@ EX void show() { } else if(uni == 'f') pushScreen(showStereo); - else if(uni == 'n' && !rug::rugged) - gwhere = eGeometry((gwhere+1) % 4); + else if(uni == 'n' && !rug::rugged) + pushScreen(rug_geometry_choice); else if(uni == 'g' && !rug::rugged && CAP_SDL) rendernogl = !rendernogl; else if(uni == 's') { @@ -1989,6 +1892,10 @@ int rugArgs() { else if(argis("-ruggeo")) { shift(); gwhere = (eGeometry) argi(); + if(gwhere == gEuclid) gwhere = rgEuclid; + if(gwhere == gSphere) gwhere = rgSphere; + if(gwhere == gNormal) gwhere = rgHyperbolic; + if(gwhere == gElliptic) gwhere = rgElliptic; } else if(argis("-rugpers")) { diff --git a/surface.cpp b/surface.cpp index ad2de140..29b5b010 100644 --- a/surface.cpp +++ b/surface.cpp @@ -11,6 +11,8 @@ #if CAP_SURFACE namespace hr { +#define USING_NATIVE_GEOMETRY dynamicval gw(geometry, hr::rug::gwhere) + EX namespace surface { ld sech(ld d) { return 1 / cosh(d); } @@ -24,7 +26,7 @@ string shape_name[] = { "hypersian rug", "tractricoid", "Dini's surface", "Kuen EX eShape sh; -hyperpoint unit_vector[3] = {hpxyz(1,0,0), hpxyz(0,1,0), hpxyz(0,0,1)}; +hyperpoint unit_vector[3] = {point3(1,0,0), point3(0,1,0), point3(0,0,1)}; ld last_int_of = 0, last_int = 0; @@ -67,7 +69,7 @@ hyperpoint coord(hyperpoint h) { ld r = 1 / cosh(t); ld x = t - tanh(t); - return hpxyz( r * sin(v), r * cos(v), x ); + return point31( r * sin(v), r * cos(v), x ); break; } @@ -77,7 +79,7 @@ hyperpoint coord(hyperpoint h) { ld a = sqrt(1-dini_b*dini_b); - return hpxyz( a * sin(v) * sin(t), a * cos(v) * sin(t), a * (cos(t) + log(tan(t/2))) + dini_b * v ); + return point31( a * sin(v) * sin(t), a * cos(v) * sin(t), a * (cos(t) + log(tan(t/2))) + dini_b * v ); break; } @@ -87,7 +89,7 @@ hyperpoint coord(hyperpoint h) { ld deno = 1 / (1 + u * u * sin(v) * sin(v)); - return hpxyz( + return point31( 2 * (cos(u) + u * sin(u)) * sin(v) * deno, 2 * (sin(u) - u * cos(u)) * sin(v) * deno, log(tan(v/2)) + 2 * cos(v) * deno @@ -101,7 +103,7 @@ hyperpoint coord(hyperpoint h) { ld phi = hyper_b * cosh(v); ld psi = integral(v); - return hpxyz( phi * cos(u), phi * sin(u), psi ); + return point31( phi * cos(u), phi * sin(u), psi ); } default: @@ -122,10 +124,10 @@ hyperpoint coord_derivative(hyperpoint h, int cc) { ld v = h[1]; if(cc == 0) { ld phi = hyper_b * cosh(v); - return hpxyz( phi * -sin(u), phi * cos(u), 0 ); + return point3( phi * -sin(u), phi * cos(u), 0 ); } else { - return hpxyz( hyper_b * sinh(v) * cos(u), hyper_b * sinh(v) * sin(u), f(v) ); + return point3( hyper_b * sinh(v) * cos(u), hyper_b * sinh(v) * sin(u), f(v) ); } } case dsKuen: { @@ -134,12 +136,12 @@ hyperpoint coord_derivative(hyperpoint h, int cc) { ld denom = pow(sin(v),2)*(u*u)+1; ld denom2 = denom * denom; if(cc == 1) - return hpxyz ( + return point3( 2*sin(v)/denom*u*cos(u)+-4*(sin(u)*u+cos(u))*pow(sin(v),3)/denom2*u, -4*pow(sin(v),3)*(sin(u)-u*cos(u))/denom2*u+2*sin(u)*sin(v)/denom*u, -4*pow(sin(v),2)/denom2*u*cos(v) ); - else return hpxyz ( + else return point3( 2*(sin(u)*u+cos(u))/denom*cos(v)+-4*(sin(u)*u+cos(u))*pow(sin(v),2)/denom2*(u*u)*cos(v), 2*(sin(u)-u*cos(u))/denom*cos(v)+-4*pow(sin(v),2)*(sin(u)-u*cos(u))/denom2*(u*u)*cos(v), -4*sin(v)/denom2*(u*u)*pow(cos(v),2)+1/tan(v/2)*(pow(tan(v/2),2)+1)/2+-2*sin(v)/denom @@ -166,13 +168,13 @@ ld compute_curvature(hyperpoint at) { hyperpoint shape_origin() { switch(sh) { case dsDini: - return hpxyz(M_PI * .82, 0, 0); + return point31(M_PI * .82, 0, 0); case dsTractricoid: - return hpxyz(1, 0, 0); + return point31(1, 0, 0); case dsKuen: - return hpxyz(M_PI * .500001, M_PI * 1, 0); + return point31(M_PI * .500001, M_PI * 1, 0); case dsHyperlike: - return hpxyz(0,0,0); + return point31(0,0,0); default: return Hypc; } @@ -284,7 +286,7 @@ transmatrix create_M_matrix(hyperpoint zero, hyperpoint v1) { transmatrix T = build_matrix(Te0, Te1, Hypc, C03); v1 = v1 / hypot_d(3, T*v1); - hyperpoint v2 = hpxyz(1e-3, 1e-4, 0); + hyperpoint v2 = point3(1e-3, 1e-4, 0); v2 = v2 - v1 * ((T*v1) | (T*v2)) / hypot_d(3, T*v1); v2 = v2 / hypot_d(3, T*v2); @@ -307,9 +309,9 @@ dexp_origin at_zero(hyperpoint zero, transmatrix start) { println(hlog, "zero = ", zero); println(hlog, "curvature at zero = ", compute_curvature(zero)); - println(hlog, "curvature at X1 = ", compute_curvature(zero + hpxyz(.3, 0, 0))); - println(hlog, "curvature at X2 = ", compute_curvature(zero + hpxyz(0, .3, 0))); - println(hlog, "curvature at X3 = ", compute_curvature(zero + hpxyz(.4, .3, 0))); + println(hlog, "curvature at X1 = ", compute_curvature(zero + point3(.3, 0, 0))); + println(hlog, "curvature at X2 = ", compute_curvature(zero + point3(0, .3, 0))); + println(hlog, "curvature at X3 = ", compute_curvature(zero + point3(.4, .3, 0))); return {start, create_M_matrix(zero, unit_vector[0]), zero}; } @@ -338,14 +340,14 @@ void addTriangleV(rug::rugpoint *t1, rug::rugpoint *t2, rug::rugpoint *t3, ld le } hyperpoint kuen_cross(ld v, ld u) { - auto du = coord_derivative(hpxyz(v,u,0), 0); - auto dv = coord_derivative(hpxyz(v,u,0), 1); + auto du = coord_derivative(point3(v,u,0), 0); + auto dv = coord_derivative(point3(v,u,0), 1); return du^dv; } ld kuen_hypot(ld v, ld u) { - auto du = coord_derivative(hpxyz(v,u,0), 0); - auto dv = coord_derivative(hpxyz(v,u,0), 1); + auto du = coord_derivative(point3(v,u,0), 0); + auto dv = coord_derivative(point3(v,u,0), 1); auto n = hypot_d(3, du^dv); return n; } @@ -385,8 +387,8 @@ void draw_kuen_map() { for(int h=0; h<512; h++) { ld v = M_PI * (r+.5) / 512; ld u = 2 * M_PI * (h+.5) / 512; - auto du = coord_derivative(hpxyz(v,u,0), 0); - auto dv = coord_derivative(hpxyz(v,u,0), 1); + auto du = coord_derivative(point3(v,u,0), 0); + auto dv = coord_derivative(point3(v,u,0), 1); auto n = hypot_d(3, du^dv); if(n > nmax) nmax = n; @@ -459,7 +461,8 @@ void run_hyperlike() { int sgn = y > 0 ? 1 : -1; ld phi = hyper_b * cosh(y); int pt = y * precision * sgn / hyperlike_bound(); - p->flat = hpxyz(phi * cos(x), phi * sin(x), sgn * integral_table[pt]); + USING_NATIVE_GEOMETRY; + p->native = point31(phi * cos(x), phi * sin(x), sgn * integral_table[pt]); p->valid = true; } } @@ -520,7 +523,8 @@ void run_kuen() { np->inqueue = false; np->dist = 0; np->h = p->h; - np->flat = coord(px.params); + USING_NATIVE_GEOMETRY; + np->native = coord(px.params); np->surface_point = px; np->dexp_id = p->dexp_id; coverages[p->dexp_id] |= pid[part]; @@ -534,9 +538,11 @@ void run_kuen() { for(int i=0; i<3; i++) if(!r[i]) looks_good = false; if(!looks_good) continue; - for(int i=0; i<3; i++) - if(hypot_d(3, r[i]->flat - r[(i+1)%3]->flat) > .2) + for(int i=0; i<3; i++) { + USING_NATIVE_GEOMETRY; + if(hypot_d(3, r[i]->native - r[(i+1)%3]->native) > .2) looks_good = false; + } if(looks_good) addTriangleV(r[0], r[1], r[2]); } @@ -557,9 +563,11 @@ void run_kuen() { template void run_function(T f) { full_mesh(); - for(auto p: rug::points) - p->flat = f(p->h), + for(auto p: rug::points) { + USING_NATIVE_GEOMETRY; + p->native = f(p->h), p->valid = true; + } } void run_other() { @@ -572,7 +580,10 @@ void run_other() { auto h = p->h; p->surface_point = map_to_surface(h, dp); - p->flat = coord(p->surface_point.params); + if(1) { + USING_NATIVE_GEOMETRY; + p->native = coord(p->surface_point.params); + } history::progress(XLAT("solving the geodesics on: %1, %2/%3", XLAT(shape_name[sh]), its(it), its(isize(rug::points)))); if(p->surface_point.remaining_distance == 0) coverage.emplace_back(h, rchar(it) + 256 * 7); @@ -641,8 +652,8 @@ EX void run_shape(eShape s) { ld minx = 1e9, maxx = -1e9; for(auto p: rug::points) if(p->valid) { - minx = min(p->flat[2], minx); - maxx = max(p->flat[2], maxx); + minx = min(p->native[2], minx); + maxx = max(p->native[2], maxx); rug::qvalid++; } @@ -650,7 +661,7 @@ EX void run_shape(eShape s) { ld shift = -(minx + maxx) / 2; for(auto p: rug::points) if(p->valid) - p->flat[2] += shift; + p->native[2] += shift; } rug::apply_rotation(crot); @@ -747,7 +758,7 @@ EX void show_surfaces() { } else if(uni == 'x') for(auto p: rug::points) - p->flat = p->surface_point.params; + p->native = p->surface_point.params; else if(uni == '#') dialog::editNumber(dini_b, -1, 1, .05, .15, XLAT("parameter"), XLAT("The larger the number, the more twisted it is.")