1 // Hyperbolic Rogue -- hyperbolic graphics
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file hypgraph.cpp
5 * \brief mapping hyperpoints to the screen, and related functions
6 */
7
8 #include "hyper.h"
9 namespace hr {
10
11 hyperpoint ghxy, ghgxy;
12 shiftpoint ghpm = shiftless(C02);
13
14 #if HDR
sphereflipped()15 inline bool sphereflipped() { return sphere && pconf.alpha > 1.1; }
16 #endif
17
ghcheck(hyperpoint & ret,const shiftpoint & H)18 void ghcheck(hyperpoint &ret, const shiftpoint &H) {
19 if(hypot_d(2, ret-ghxy) < hypot_d(2, ghgxy-ghxy)) {
20 ghpm = H; ghgxy = ret;
21 }
22 }
23
camrotate(ld & hx,ld & hy)24 EX void camrotate(ld& hx, ld& hy) {
25 ld cam = pconf.camera_angle * degree;
26 GLfloat cc = cos(cam);
27 GLfloat ss = sin(cam);
28 ld ux = hx, uy = hy * cc + ss, uz = cc - ss * hy;
29 hx = ux / uz, hy = uy / uz;
30 }
31
non_spatial_model()32 EX bool non_spatial_model() {
33 if(among(pmodel, mdRotatedHyperboles, mdJoukowsky, mdJoukowskyInverted, mdPolygonal, mdPolynomial))
34 return true;
35 if(pmodel == mdSpiral && euclid)
36 return true;
37 #if CAP_GL
38 return pmodel && vid.consider_shader_projection && (get_shader_flags() & SF_DIRECT);
39 #else
40 return false;
41 #endif
42 }
43
44 EX hyperpoint perspective_to_space(hyperpoint h, ld alpha IS(pconf.alpha), eGeometryClass gc IS(ginf[geometry].cclass)) {
45 ld hx = h[0], hy = h[1];
46
47 if(gc == gcEuclid)
48 return hpxy(hx * (1 + alpha), hy * (1 + alpha));
49
50 ld hr = hx*hx+hy*hy;
51 if(LDIM == 3) hr += h[2]*h[2];
52
53 if(hr > .9999 && gc == gcHyperbolic) return Hypc;
54
55 ld A, B, C;
56
57 ld curv = gc == gcSphere ? 1 : -1;
58
59 A = 1+curv*hr;
60 B = 2*hr*alpha*-curv;
61 C = 1 - curv*hr*alpha*alpha;
62
63 B /= A; C /= A;
64
65 ld rootsign = 1;
66 if(gc == gcSphere && pconf.alpha > 1) rootsign = -1;
67
68 ld hz = B / 2 + rootsign * sqrt(C + B*B/4);
69
70 hyperpoint H;
71 H[0] = hx * (hz+alpha);
72 H[1] = hy * (hz+alpha);
73 if(LDIM == 3) H[2] = h[2] * (hz + alpha);
74 H[LDIM] = hz;
75
76 return H;
77 }
78
79 EX hyperpoint space_to_perspective(hyperpoint z, ld alpha IS(pconf.alpha)) {
80 ld s = 1 / (alpha + z[LDIM]);
81 z[0] *= s;
82 z[1] *= s;
83 if(GDIM == 3) {
84 z[2] *= s;
85 z[3] = 0;
86 }
87 else
88 z[2] = 0;
89 return z;
90 }
91
pointable()92 EX hyperpoint pointable() {
93 return WDIM == 2 && GDIM == 3 ? zpush0(cgi.FLOOR) : C0;
94 }
95
96 /** find a shiftpoint which minimizes value -- we represent points by matrices to make things a bit simpler */
97
minimize_point_value(shiftmatrix T,function<ld (const shiftmatrix &)> value)98 EX shiftmatrix minimize_point_value(shiftmatrix T, function<ld(const shiftmatrix&)> value) {
99
100 ld best = value(T);
101
102 for(int it=0; it<50; it++)
103 for(int s=0; s<2*WDIM; s++) {
104 shiftmatrix T1 = T * cpush(s/2, (s&1?1:-1) * pow(1.2, -it));
105 ld dist = value(T1);
106 if(dist < best) best = dist, T = T1;
107 if(mdBandAny()) {
108 T1.shift += 2 * M_PI;
109 dist = value(T1);
110 if(dist < best) best = dist, T = T1;
111 T1.shift -= 4 * M_PI;
112 dist = value(T1);
113 if(dist < best) best = dist, T = T1;
114 T1.shift += 2 * M_PI;
115 }
116 }
117
118 return T;
119 }
120
find_on_screen(hyperpoint hxy,const shiftmatrix & T)121 EX shiftpoint find_on_screen(hyperpoint hxy, const shiftmatrix& T) {
122 hyperpoint rel = pointable();
123 auto distance_at = [&] (const shiftmatrix& T1) {
124 hyperpoint h1;
125 applymodel(T1*rel, h1);
126 return sqhypot_d(2, hxy - h1);
127 };
128 return minimize_point_value(T, distance_at) * rel;
129 }
130
gethyper(ld x,ld y)131 EX shiftpoint gethyper(ld x, ld y) {
132
133 ld hx = (x - current_display->xcenter) / current_display->radius;
134 ld hy = (y - current_display->ycenter) / current_display->radius / pconf.stretch;
135 hyperpoint hxy = point3(hx, hy, 0);
136
137 if(WDIM == 2 && GDIM == 3) {
138 return mouseover ? find_on_screen(hxy, ggmatrix(mouseover)): shiftless(Hypc);
139 }
140
141 if(pmodel) {
142 ghxy = hxy;
143 return find_on_screen(hxy, rgpushxto0(ghpm));
144 }
145
146 if(pconf.camera_angle) camrotate(hx, hy);
147
148 return shiftless(perspective_to_space(hpxyz(hx, hy, 0)));
149 }
150
ballmodel(hyperpoint & ret,double alpha,double d,double zl)151 void ballmodel(hyperpoint& ret, double alpha, double d, double zl) {
152 hyperpoint H = ypush(vid.camera) * xpush(d) * ypush(zl) * C0;
153 ld tzh = pconf.ballproj + H[LDIM];
154 ld ax = H[0] / tzh;
155 ld ay = H[1] / tzh;
156
157 ld ca = cos(alpha), sa = sin(alpha);
158
159 ret[0] = ax * ca;
160 ret[1] = ay;
161 ret[2] = ax * sa;
162
163 models::apply_ball(ret[2], ret[1]);
164 }
165
use_z_coordinate()166 bool use_z_coordinate() {
167 #if CAP_VR
168 if(vrhr::rendering()) return true;
169 #endif
170 return current_display->stereo_active();
171 }
172
apply_depth(hyperpoint & f,ld z)173 void apply_depth(hyperpoint &f, ld z) {
174 if(vid.usingGL)
175 f[2] = z * pconf.depth_scaling;
176 else {
177 z = z * current_display->radius * pconf.depth_scaling;
178 ld mul = current_display->radius / (current_display->radius + z);
179 f[0] = f[0] * mul;
180 f[1] = f[1] * mul;
181 f[2] = vid.xres * current_display->eyewidth() / 2 / current_display->radius + vid.ipd * mul / 2;
182 }
183 }
184
hypot_zlev(ld zlev,ld & d,ld & df,ld & zf)185 bool hypot_zlev(ld zlev, ld& d, ld& df, ld& zf) {
186 if(zlev == 1) {
187 df = 1; zf = 0;
188 return false;
189 }
190 else {
191 // (0,0,1) -> (0, sin z, cos z) -> (sin d cos z, sin z, cos d cos z)
192 ld z = geom3::factor_to_lev(zlev);
193
194 ld tz = sin_auto(z);
195 ld td = sin_auto(abs(d)) * cos_auto(z);
196 ld h = hypot(td, tz);
197 zf = tz / h, df = td / h;
198
199 if(d > 0)
200 d = hypot_auto(d, z);
201 else {
202 d = -hypot_auto(-d, z);
203 zf = -zf;
204 }
205 return true;
206 }
207 }
208
209 int twopoint_sphere_flips;
210 bool twopoint_do_flips;
211
find_zlev(hyperpoint & H)212 ld find_zlev(hyperpoint& H) {
213
214 if(spatial_graphics) {
215 ld zlev = zlevel(H);
216 if(zlev > 1-1e-6 && zlev < 1+1e-6) return 1;
217 H /= zlev;
218 return zlev;
219 }
220
221 return 1;
222 }
223
get_tz(hyperpoint H)224 ld get_tz(hyperpoint H) {
225 ld tz = pconf.alpha+H[LDIM];
226 if(tz < BEHIND_LIMIT && tz > -BEHIND_LIMIT) tz = BEHIND_LIMIT;
227 return tz;
228 }
229
atan2(hyperpoint h)230 EX ld atan2(hyperpoint h) {
231 return atan2(h[1], h[0]);
232 }
233
move_z_to_y(hyperpoint & H)234 pair<ld, ld> move_z_to_y(hyperpoint& H) {
235 if(GDIM == 2) return make_pair(0, 0);
236 ld R = hypot(H[1], H[2]);
237 pair<ld, ld> res = { H[1] / R, H[2] / R };
238 H[1] = R; H[2] = 0;
239 return res;
240 }
241
move_y_to_z(hyperpoint & H,pair<ld,ld> coef)242 void move_y_to_z(hyperpoint& H, pair<ld, ld> coef) {
243 if(GDIM == 3) {
244 H[2] = H[1] * coef.second;
245 H[1] = H[1] * coef.first;
246 #if MAXMDIM >= 4
247 H[3] = 1;
248 #endif
249 }
250 }
251
makeband(shiftpoint H,hyperpoint & ret,const T & f)252 template<class T> void makeband(shiftpoint H, hyperpoint& ret, const T& f) {
253 ld zlev = find_zlev(H.h);
254 models::apply_orientation_yz(H[1], H[2]);
255 models::apply_orientation(H[0], H[1]);
256 auto r = move_z_to_y(H.h);
257
258 ld x, y, yf, zf=0;
259 y = asin_auto(H[1]);
260 x = asin_auto_clamp(H[0] / cos_auto(y));
261 if(sphere) {
262 if(H[LDIM] < 0 && x > 0) x = M_PI - x;
263 else if(H[LDIM] < 0 && x <= 0) x = -M_PI - x;
264 }
265 x += H.shift;
266 hypot_zlev(zlev, y, yf, zf);
267
268 f(x, y);
269
270 ld yzf = y * zf; y *= yf;
271 ret = hpxyz(x / M_PI, y / M_PI, 0);
272 move_y_to_z(ret, r);
273 models::apply_orientation(ret[1], ret[0]);
274 models::apply_orientation_yz(ret[2], ret[1]);
275 if(zlev != 1 && use_z_coordinate())
276 apply_depth(ret, yzf / M_PI);
277 return;
278 }
279
band_conformal(ld & x,ld & y)280 void band_conformal(ld& x, ld& y) {
281 switch(cgclass) {
282 case gcSphere:
283 y = atanh(sin(y));
284 x *= 2; y *= 2;
285 break;
286 case gcHyperbolic:
287 y = 2 * atan(tanh(y/2));
288 x *= 2; y *= 2;
289 break;
290 case gcEuclid:
291 default:
292 // y = y;
293 y *= 2; x *= 2;
294 break;
295 }
296 }
297
make_twopoint(ld & x,ld & y)298 void make_twopoint(ld& x, ld& y) {
299 auto p = pconf.twopoint_param;
300 ld dleft = hypot_auto(x-p, y);
301 ld dright = hypot_auto(x+p, y);
302 if(sphere) {
303 int tss = twopoint_sphere_flips;
304 if(tss&1) { tss--;
305 dleft = 2*M_PI - 2*p - dleft;
306 dright = 2*M_PI - 2*p - dright;
307 swap(dleft, dright);
308 y = -y;
309 }
310 while(tss) { tss -= 2;
311 dleft = 2*M_PI - 4*p + dleft;
312 dright = 2*M_PI - 4*p + dright;
313 }
314 }
315 x = (dright*dright-dleft*dleft) / 4 / p;
316 y = (y>0?1:-1) * sqrt(dleft * dleft - (x-p)*(x-p) + 1e-9);
317 }
318
mobius(hyperpoint h,ld angle,ld scale=1)319 hyperpoint mobius(hyperpoint h, ld angle, ld scale = 1) {
320 h = perspective_to_space(h * scale, 1, gcSphere);
321 h = cspin(1, 2, angle * degree) * h;
322 return space_to_perspective(h, 1) / scale;
323 }
324
compute_hybrid(hyperpoint H,int rootid)325 hyperpoint compute_hybrid(hyperpoint H, int rootid) {
326 auto& t = pconf.twopoint_param;
327 hyperpoint Hl = xpush(+t) * H;
328 hyperpoint Hr = xpush(-t) * H;
329 ld g = (Hl[0] + 1e-7) / (Hl[1] + 1e-8);
330 ld d = hdist0(Hr);
331
332 hyperpoint spinned = spintox(Hl) * xpush0(2*t);
333 if(Hl[0] < 0) spinned = pispin * spinned;
334
335 ld y = asin_auto(spinned[1]);
336 ld x = asin_auto_clamp(spinned[0] / cos_auto(y));
337
338 int sign = (Hl[0] > 0 ? 1 : -1) * hdist0(Hl) < x ? -1 : 1;
339
340 switch(rootid & 3) {
341 case 1: sign = -sign; break;
342 case 2: sign = 1; break;
343 case 3: sign = -1; break;
344 }
345
346 // (x + t) / g = y
347 // yy + (x-t)(x-t) = dd
348 // (x+t)*(x+t)/g*g + x*x + t*t - 2*x*t = dd
349
350 // x*x*(1+1/g*g) + t*t*(1+1/g*g) + 2xt (1/gg-1) = dd
351 // xx + 2xt (1/gg-1) / (1+1/gg) = dd / (1+1/gg) - tt
352
353 ld b = t*(1/g/g - 1) / (1+1/g/g);
354 ld c = d*d / (1+1/g/g) - t*t;
355
356 // xx + 2bx = c
357 // xx + 2bx + bb = c + bb
358 // (x+b)^2 = c+bb
359 // x = +/- sqrt(c+bb) - b
360
361 ld a = c+b*b;
362
363 hyperpoint ret;
364 ret[0] = (a > 0 ? sign * sqrt(a) : 0) - b;
365 ret[1] = (ret[0] + t) / g;
366 ret[2] = 0;
367
368 return ret;
369 }
370
signed_sqrt(ld x)371 EX ld signed_sqrt(ld x) { return x > 0 ? sqrt(x) : -sqrt(-x); }
372
373 EX int axial_x, axial_y;
374
apply_perspective(const hyperpoint & H,hyperpoint & ret)375 EX void apply_perspective(const hyperpoint& H, hyperpoint& ret) {
376 if(H[2] == 0) { ret[0] = 1e6; ret[1] = 1e6; ret[2] = 1; return; }
377 ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2;
378 ret[0] = H[0]/H[2] * ratio;
379 ret[1] = H[1]/H[2] * ratio;
380 ret[2] = 1;
381 ret[3] = 1;
382 }
383
apply_nil_rotation(hyperpoint & H)384 EX void apply_nil_rotation(hyperpoint& H) {
385 if(nil) {
386 H[2] -= H[0] * H[1] / 2;
387 models::apply_orientation(H[0], H[1]);
388 H[2] += H[0] * H[1] / 2 * pconf.rotational_nil;
389 models::apply_orientation(H[1], H[0]);
390 }
391 }
392
applymodel(shiftpoint H_orig,hyperpoint & ret)393 EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
394 apply_other_model(H_orig, ret, pmodel);
395 }
396
vr_sphere(hyperpoint & ret,hyperpoint & H,eModel md)397 EX void vr_sphere(hyperpoint& ret, hyperpoint& H, eModel md) {
398 ret = H;
399 int flip = 1;
400 if(md == mdHalfplane) flip = -flip;
401 if(pconf.alpha < 1) flip = -flip;
402 ret *= pow(sqhypot_d(3, H), (flip * pconf.depth_scaling-1) / 2);
403 ret[2] += pconf.alpha;
404 if(md == mdHalfplane) {
405 ld d = sqhypot_d(3, ret);
406 ret /= abs(d);
407 }
408 }
409
vr_disk(hyperpoint & ret,hyperpoint & H)410 void vr_disk(hyperpoint& ret, hyperpoint& H) {
411 if(euclid) {
412 ret = H;
413 ret[2] = vid.depth * (1 - (ret[2] - 1) * pconf.depth_scaling) + pconf.alpha + vid.camera;
414 }
415 else if(sphere) {
416 vr_sphere(ret, H, mdDisk);
417 return;
418 }
419 else {
420 ld zlev = find_zlev(H);
421 ld zl = vid.depth-geom3::factor_to_lev(zlev) * pconf.depth_scaling;
422
423 ld d = hdist0(H);
424 ld dd = hypot_d(2, H);
425
426 hyperpoint H1 = ypush(vid.camera) * xpush(d) * ypush0(zl);
427 ld tzh = pconf.alpha + H1[2];
428 ld ax = H1[0] / tzh;
429 ld ay = H1[1] / tzh;
430
431 ret[0] = ax * H[0] / dd;
432 ret[1] = ax * H[1] / dd;
433 ret[2] = ay;
434 }
435 }
436
437 #if MAXMDIM >= 4
438 /** Compute the three-point projection. Currently only works in isotropic 3D spaces. */
threepoint_projection(const hyperpoint & H,hyperpoint & ret)439 EX void threepoint_projection(const hyperpoint& H, hyperpoint& ret) {
440
441 hyperpoint H1 = H;
442 find_zlev(H1);
443 if(true) {
444 models::apply_orientation_yz(H1[1], H1[2]);
445 models::apply_orientation(H1[0], H1[1]);
446 }
447
448 auto p = pconf.twopoint_param;
449
450 ld dist[3];
451 for(int i=0; i<3; i++) {
452 hyperpoint h1 = xspinpush0(2*M_PI*i/3, p);
453 dist[i] = geo_dist(h1, H1);
454 }
455
456 /* we are looking for the points (x,y,z) such that:
457 (x-xi)^2 + (y-yi)^2 + z^2 = di^2
458
459 which is equivalent to:
460 x^2+y^2+z^2 -2xxi -2yyi = di^2-xi^2-yi^2
461
462 After setting s = x^2+y^2+z^2, we get a system of linear equations for (x,y,s)
463 */
464
465 dynamicval<eGeometry> g(geometry, gEuclid);
466
467 transmatrix T = Id;
468 hyperpoint v = C0;
469 for(int i=0; i<3; i++) {
470 hyperpoint pp = xspinpush0(2*M_PI*i/3, p);
471 v[i] = dist[i]*dist[i] - p*p;
472 T[i][0] = -2 * pp[0];
473 T[i][1] = -2 * pp[1];
474 T[i][2] = 1;
475 }
476
477 transmatrix U = inverse3(T);
478 hyperpoint sxy = U * v;
479
480 // compute the actual z based on s
481 sxy[2] = sxy[2] - sqhypot_d(2, sxy);
482 sxy[2] = sxy[2] > 0 ? sqrt(sxy[2]) : 0;
483
484 if(H1[2] < 0) sxy[2] *= -1;
485
486 sxy[3] = 1;
487
488 geometry = gCubeTiling;
489
490 ret = sxy;
491 models::apply_orientation(ret[1], ret[0]);
492 models::apply_orientation_yz(ret[2], ret[1]);
493 }
494 #endif
495
apply_other_model(shiftpoint H_orig,hyperpoint & ret,eModel md)496 EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) {
497
498 hyperpoint H = H_orig.h;
499
500 if(models::product_model(md)) {
501 ld zlev = zlevel(H_orig.h);
502 H_orig.h /= exp(zlev);
503 hybrid::in_underlying_geometry([&] { applymodel(H_orig, ret); });
504 ret[2] = zlev * pconf.product_z_scale;
505 ret = NLP * ret;
506 return;
507 }
508
509 switch(md) {
510 case mdPerspective: {
511 if(prod) H = product::inverse_exp(H);
512 apply_nil_rotation(H);
513 H = lp_apply(H);
514 apply_perspective(H, ret);
515 return;
516 }
517
518 case mdGeodesic: {
519 auto S = lp_apply(inverse_exp(H_orig, pNORMAL | pfNO_DISTANCE));
520 apply_perspective(S, ret);
521 return;
522 }
523
524 case mdPixel:
525 ret = H / current_display->radius;
526 return;
527
528 case mdBall: {
529 if(vrhr::rendering()) {
530 vr_disk(ret, H);
531 return;
532 }
533 ld zlev = find_zlev(H);
534
535 ld zl = vid.depth-geom3::factor_to_lev(zlev) * pconf.depth_scaling;
536
537 ballmodel(ret, atan2(H), hdist0(H), zl);
538 break;
539 }
540
541 case mdDisk: {
542 if(nonisotropic) {
543 ret = lp_apply(inverse_exp(H_orig, pNORMAL | pfNO_DISTANCE));
544 ld w;
545 if(sn::in()) {
546 // w = 1 / sqrt(1 - sqhypot_d(3, ret));
547 // w = w / (pconf.alpha + w);
548 w = 1 / (sqrt(1 - sqhypot_d(3, ret)) * pconf.alpha + 1);
549 }
550 else {
551 w = hypot_d(3, ret);
552 w = sinh(w) / ((pconf.alpha + cosh(w)) * w);
553 }
554 for(int i=0; i<3; i++) ret[i] *= w;
555 ret[3] = 1;
556 break;
557 }
558 if(vrhr::rendering() && WDIM == 2) {
559 vr_disk(ret, H);
560 return;
561 }
562 ld tz = get_tz(H);
563 if(!pconf.camera_angle) {
564 ret[0] = H[0] / tz;
565 ret[1] = H[1] / tz;
566 if(GDIM == 3) ret[2] = H[2] / tz;
567 else ret[2] = vid.xres * current_display->eyewidth() / 2 / current_display->radius - vid.ipd / tz / 2;
568 if(MAXMDIM == 4) ret[3] = 1;
569 }
570 else {
571 ld tx = H[0];
572 ld ty = H[1];
573 ld cam = pconf.camera_angle * degree;
574 GLfloat cc = cos(cam);
575 GLfloat ss = sin(cam);
576 ld ux = tx, uy = ty * cc - ss * tz, uz = tz * cc + ss * ty;
577 ret[0] = ux / uz;
578 ret[1] = uy / uz;
579 ret[2] = vid.xres * current_display->eyewidth() / 2 / current_display->radius - vid.ipd / uz / 2;
580 }
581 return;
582 }
583
584 case mdCentralInversion: {
585 ld tz = get_tz(H);
586 for(int d=0; d<GDIM; d++) ret[d] = H[d] / tz;
587 for(int d=GDIM; d<MAXMDIM; d++) ret[d] = 1;
588 ld r = 0;
589 for(int d=0; d<GDIM; d++) r += ret[d]*ret[d];
590 for(int d=0; d<GDIM; d++) ret[d] /= r;
591 return;
592 }
593
594 case mdHalfplane: {
595 if(sphere && vrhr::rendering()) {
596 vr_sphere(ret, H, md);
597 return;
598 }
599 // Poincare to half-plane
600
601 ld zlev = find_zlev(H);
602 H = space_to_perspective(H);
603
604 models::apply_orientation_yz(H[1], H[2]);
605 models::apply_orientation(H[0], H[1]);
606
607 H[1] += 1;
608 double rad = sqhypot_d(GDIM, H);
609 H /= -rad;
610 H[1] += .5;
611
612 if(GDIM == 3) {
613 // a bit simpler when we do not care about 3D
614 H *= pconf.halfplane_scale;
615 ret[0] = -H[0];
616 ret[1] = 1 + H[1];
617 ret[2] = H[2];
618 ret[3] = 1;
619 models::apply_orientation(ret[1], ret[0]);
620 models::apply_orientation_yz(ret[2], ret[1]);
621 break;
622 }
623
624 models::apply_orientation(H[0], H[1]);
625
626 H *= pconf.halfplane_scale;
627
628 ret[0] = -models::osin - H[0];
629 ld height = 0;
630 if(zlev != 1) {
631 if(abs(models::ocos) > 1e-5)
632 height += H[1] * (pow(zlev, models::ocos) - 1);
633 if(abs(models::ocos) > 1e-5 && models::osin)
634 height += H[0] * models::osin * (pow(zlev, models::ocos) - 1) / models::ocos;
635 else if(models::osin)
636 height += H[0] * models::osin * log(zlev);
637 }
638 ret[1] = models::ocos + H[1];
639 ret[2] = GDIM == 3 ? H[2] : 0;
640 if(MAXMDIM == 4) ret[3] = 1;
641 if(zlev != 1 && use_z_coordinate())
642 apply_depth(ret, height);
643 else
644 ret[1] += height * pconf.depth_scaling;
645 break;
646 }
647
648 case mdAxial: {
649 models::apply_orientation_yz(H[1], H[2]);
650 models::apply_orientation(H[0], H[1]);
651
652 ld& mt = pconf.model_transition;
653
654 ld z = H[LDIM];
655 if(mt != 1) z += (1-mt) * pconf.alpha;
656
657 ret[0] = H[0] / z;
658 ret[1] = H[1] / z;
659 if(GDIM == 3) ret[2] = H[2] / z;
660 else ret[2] = 0;
661 ret[3] = 1;
662
663 if(mt) for(int i=0; i<LDIM; i++) {
664 if(mt < 1)
665 ret[i] *= mt;
666 ret[i] = atan_auto(ret[i]);
667 if(mt < 1)
668 ret[i] /= mt;
669 }
670
671 if(sphere) ret[0] += axial_x * M_PI, ret[1] += axial_y * M_PI;
672
673 models::apply_orientation(ret[1], ret[0]);
674 models::apply_orientation_yz(ret[2], ret[1]);
675 break;
676 }
677
678 case mdAntiAxial: {
679 models::apply_orientation_yz(H[1], H[2]);
680 models::apply_orientation(H[0], H[1]);
681
682 ret[0] = asin_auto(H[0]);
683 ret[1] = asin_auto(H[1]);
684
685 ret[2] = 0; ret[3] = 1;
686
687 models::apply_orientation(ret[1], ret[0]);
688 models::apply_orientation_yz(ret[2], ret[1]);
689 break;
690 }
691
692 case mdQuadrant: {
693 H = space_to_perspective(H);
694 models::apply_orientation_yz(H[1], H[2]);
695 models::apply_orientation(H[0], H[1]);
696
697 tie(H[0], H[1]) = make_pair((H[0] + H[1]) / sqrt(2), (H[1] - H[0]) / sqrt(2));
698
699 H[1] += 1;
700 double rad = sqhypot_d(GDIM, H);
701 H /= -rad;
702 H[1] += .5;
703
704 H *= 2;
705
706 ld x = exp(-H[0]/2);
707 ret[0] = -H[1] * x - 1;
708 ret[1] = H[1] / x + 1;
709
710 models::apply_orientation(ret[1], ret[0]);
711 models::apply_orientation_yz(ret[2], ret[1]);
712 break;
713 }
714
715 case mdHorocyclic: {
716
717 find_zlev(H);
718
719 apply_nil_rotation(H);
720
721 if(hyperbolic) {
722 models::apply_orientation_yz(H[1], H[2]);
723 models::apply_orientation(H[0], H[1]);
724 }
725
726 ret = hyperbolic ? deparabolic10(H) : H;
727 ret *= .5;
728 ret[LDIM] = 1;
729
730 if(hyperbolic) {
731 models::apply_orientation(ret[1], ret[0]);
732 models::apply_orientation_yz(ret[2], ret[1]);
733 }
734
735 if(nonisotropic && !vrhr::rendering()) ret = lp_apply(ret);
736
737 break;
738 }
739
740 case mdHemisphere: {
741
742 #if CAP_VR
743 ld dir = vrhr::rendering() ? -1:1;
744 #else
745 constexpr ld dir = 1;
746 #endif
747
748 switch(cgclass) {
749 case gcHyperbolic: {
750 ld zl = zlevel(H);
751 ret = H / H[2];
752 ret[2] = sqrt(1 - sqhypot_d(2, ret));
753 // need to reverse in VR
754 ret = ret * (1 + (zl - 1) * ret[2] * pconf.depth_scaling * dir);
755 break;
756 }
757
758 case gcEuclid: default: {
759 // stereographic projection to a sphere
760 auto hd = hdist0(H) / pconf.euclid_to_sphere;
761 if(hd == 0) ret = hpxyz(0, 0, -1);
762 else {
763 ld x = 2 * hd / (1 + hd * hd);
764 ld y = x / hd;
765 ret = H * x / hd / pconf.euclid_to_sphere;
766 ret[2] = (1 - y);
767 ret[2] *= dir;
768 ret = ret * (1 + (H[2]-1) * y * pconf.depth_scaling * dir / pconf.euclid_to_sphere);
769 }
770 break;
771 }
772
773 case gcSphere: {
774 if(vrhr::rendering()) { vr_sphere(ret, H, md); return; }
775 ret = H;
776 if(pconf.depth_scaling != 1) {
777 ld v = intval(H, Hypc);
778 ret *= pow(v, (dir * pconf.depth_scaling-1) / 2);
779 }
780 break;
781 }
782 }
783
784 if(vrhr::rendering()) return;
785
786 swap(ret[1], ret[2]);
787
788 models::apply_ball(ret[2], ret[1]);
789
790 break;
791 }
792
793 case mdHyperboloidFlat:
794 case mdHyperboloid: {
795
796 if(nonisotropic) {
797 // if(nisot::local_perspective_used()) H = NLP * H;
798 ret = lp_apply(H);
799 break;
800 }
801 if(prod) {
802 ret = H;
803 break;
804 }
805
806 #if CAP_VR
807 if(vrhr::rendering()) {
808 if(sphere) { vr_sphere(ret, H, md); return; }
809 ret[0] = H[0] * pconf.hyperboloid_scaling;
810 ret[1] = H[1] * pconf.hyperboloid_scaling;
811 ret[2] = (pconf.alpha + H[2]);
812 if(pconf.depth_scaling != 1) {
813 ld v = intval(H, Hypc);
814 ret *= pow(v, (pconf.depth_scaling-1) / 2);
815 }
816 return;
817 }
818 #endif
819
820 if(pconf.depth_scaling != 1) {
821 ld v = intval(H, Hypc);
822 H *= pow(v, (pconf.depth_scaling-1) / 2);
823 }
824
825 if(pmodel == mdHyperboloid) {
826 ld& topz = pconf.top_z;
827 if(H[2] > topz) {
828 ld scale = sqrt(topz*topz-1) / hypot_d(2, H);
829 H *= scale;
830 H[2] = topz;
831 }
832 }
833 else {
834 H = space_to_perspective(H, pconf.alpha);
835 H[2] = 1 - pconf.alpha;
836 }
837
838 ret[0] = H[0] / 3;
839 ret[1] = (1 - H[2]) / 3;
840 ret[2] = H[1] / 3;
841
842 models::apply_ball(ret[2], ret[1]);
843 break;
844 }
845
846 case mdFisheye: {
847 ld zlev;
848 if(nonisotropic) {
849 H = lp_apply(inverse_exp(H_orig));
850 zlev = 1;
851 }
852 else {
853 zlev = find_zlev(H);
854 H = space_to_perspective(H);
855 }
856 H /= pconf.fisheye_param;
857 H[LDIM] = zlev;
858 ret = H / sqrt(1 + sqhypot_d(GDIM+1, H));
859 if(GDIM == 3) ret[LDIM] = zlev;
860 break;
861 }
862
863 case mdSimulatedPerspective: {
864 models::apply_orientation_yz(H[1], H[2]);
865 models::apply_orientation(H[0], H[1]);
866 auto yz = move_z_to_y(H);
867 hyperpoint Hl = xpush(-pconf.twopoint_param) * H;
868 hyperpoint Hr = xpush(+pconf.twopoint_param) * H;
869 ld lyx = (Hl[1] + 1e-7) / (Hl[0] + 1e-8);
870 ld ryx = (Hr[1] + 1e-7) / (Hr[0] + 1e-8);
871 // (r.x + t) * lyx = (r.x - t) * ryx = r.y
872 // r.x * lyx + t * lyx = r.x * ryx - t * ryx
873 // r.x * (lyx-ryx) = - t * (ryx + lyx)
874 // r.x = -t * (ryx+lyx) / (lyx-ryx)
875 // r.x = - 2 * t * lyx * ryx / lyx / ryx
876
877 ret[0] = -pconf.twopoint_param * (ryx + lyx) / (lyx - ryx);
878 ret[1] = (ret[0] + pconf.twopoint_param) * lyx;
879 ret[2] = 0;
880
881 move_y_to_z(ret, yz);
882 models::apply_orientation(ret[0], ret[1]);
883 models::apply_orientation_yz(ret[2], ret[1]);
884 break;
885 }
886
887 case mdTwoHybrid: {
888 models::apply_orientation_yz(H[1], H[2]);
889 models::apply_orientation(H[0], H[1]);
890 auto yz = move_z_to_y(H);
891
892 ret = compute_hybrid(H, whateveri[0]);
893
894 move_y_to_z(ret, yz);
895 models::apply_orientation(ret[0], ret[1]);
896 models::apply_orientation_yz(ret[2], ret[1]);
897 break;
898 }
899
900 case mdJoukowsky:
901 case mdJoukowskyInverted: {
902 models::apply_orientation_yz(H[1], H[2]);
903 models::apply_orientation(H[0], H[1]);
904 // with equal speed skiprope: models::apply_orientation(H[1], H[0]);
905
906 if(pconf.skiprope) {
907 static ld last_skiprope = 0;
908 static transmatrix lastmatrix;
909 if(pconf.skiprope != last_skiprope) {
910 ret = mobius(C0, -pconf.skiprope, 2);
911 const cld c1(1, 0);
912 const cld c2(2, 0);
913 const cld c4(4, 0);
914 cld w(ret[0], ret[1]);
915 cld z = sqrt(c4*w*w-c1) + c2*w;
916 if(abs(z) > 1) z = c1 / z;
917 hyperpoint zr = hpxyz(real(z), imag(z), 0);
918
919 hyperpoint inhyp = perspective_to_space(zr, 1, gcHyperbolic);
920 last_skiprope = pconf.skiprope;
921 lastmatrix = rgpushxto0(inhyp);
922 }
923 H = lastmatrix * H;
924 }
925
926 H = space_to_perspective(H);
927 auto yz = move_z_to_y(H);
928 ld r = hypot_d(2, H);
929 ld c = H[0] / r;
930 ld s = H[1] / r;
931 ld& mt = pconf.model_transition;
932 ld a = 1 - .5 * mt, b = .5 * mt;
933 swap(a, b);
934
935 ret[0] = (a * r + b/r) * c / 2;
936 ret[1] = (a * r - b/r) * s / 2;
937 ret[2] = 0;
938
939 if(pconf.skiprope)
940 ret = mobius(ret, pconf.skiprope, 2);
941
942 if(pmodel == mdJoukowskyInverted) {
943 ld r2 = sqhypot_d(2, ret);
944 ret[0] = ret[0] / r2;
945 ret[1] = -ret[1] / r2;
946 move_y_to_z(ret, yz);
947 models::apply_orientation(ret[1], ret[0]);
948
949 /*
950
951 ret[0] += 1;
952 ld alpha = atan2(ret[1], ret[0]);
953 ld mod = hypot(ret[0], ret[1]);
954 // ret[0] = cos(alpha/2) * sqrt(mod);
955 // ret[1] = sin(alpha/2) * sqrt(mod);
956 ret[0] = alpha;
957 ret[1] = log(mod); */
958 }
959 else {
960 move_y_to_z(ret, yz);
961 models::apply_orientation(ret[0], ret[1]);
962 }
963 models::apply_orientation_yz(ret[2], ret[1]);
964
965 break;
966 }
967
968 case mdPolygonal: case mdPolynomial: {
969
970 H = space_to_perspective(H);
971
972 models::apply_orientation(H[0], H[1]);
973
974 pair<long double, long double> p = polygonal::compute(H[0], H[1]);
975
976 models::apply_orientation(p.second, p.first);
977 ret[0] = p.first;
978 ret[1] = p.second;
979 ret[2] = 0;
980 break;
981 }
982
983 case mdBand:
984 if(pconf.model_transition != 1) {
985 ld& mt = pconf.model_transition;
986
987 H = space_to_perspective(H);
988
989 models::apply_orientation(H[0], H[1]);
990
991 H[0] += 1;
992 double rad = H[0]*H[0] + H[1]*H[1];
993 H[1] /= rad;
994 H[0] /= rad;
995 H[0] -= .5;
996
997 ld phi = atan2(H);
998 ld r = hypot_d(2, H);
999
1000 r = pow(r, 1 - mt);
1001 phi *= (1 - mt);
1002 ret[0] = r * cos(phi);
1003 ret[1] = r * sin(phi);
1004 ret[2] = 0;
1005
1006 ret[0] -= pow(0.5, 1-mt);
1007 ret[0] /= -(1-mt) * M_PI / 2;
1008 ret[1] /= (1-mt) * M_PI / 2;
1009
1010 models::apply_orientation(ret[1], ret[0]);
1011 }
1012 else
1013 makeband(H_orig, ret, band_conformal);
1014 break;
1015
1016 case mdMiller:
1017 makeband(H_orig, ret, [] (ld& x, ld& y) {
1018 y *= pconf.miller_parameter;
1019 band_conformal(x, y);
1020 y /= pconf.miller_parameter;
1021 });
1022 break;
1023
1024 case mdLoximuthal:
1025 makeband(H_orig, ret, [] (ld&x, ld &y) {
1026 ld orig_y = y;
1027 band_conformal(x, y);
1028 ld x0 = 0, y0 = pconf.loximuthal_parameter; band_conformal(x0, y0);
1029 y -= y0;
1030
1031 orig_y -= pconf.loximuthal_parameter;
1032
1033 if(y) x = x * orig_y / y;
1034 y = orig_y;
1035 });
1036 break;
1037
1038 case mdTwoPoint:
1039 makeband(H_orig, ret, make_twopoint);
1040 break;
1041
1042 case mdThreePoint:
1043 #if MAXMDIM >= 4
1044 threepoint_projection(H, ret);
1045 #else
1046 throw hr_exception();
1047 #endif
1048 break;
1049
1050 case mdMollweide:
1051 makeband(H_orig, ret, [] (ld& x, ld& y) {
1052 ld theta =
1053 hyperbolic ? min(y / 2 + 0.572365, y * 0.78509) :
1054 euclid ? y :
1055 y > 0 ? max(y * 0.012/0.015, M_PI/2 - (M_PI/2-y) * 0.066262/0.015708) :
1056 min(y * 0.012/0.015, -M_PI/2 + (M_PI/2+y) * 0.066262/0.015708);
1057
1058 if(sphere && abs(theta) >= M_PI/2 - 1e-6) ;
1059 else {
1060 for(int it=0; it<4; it++) {
1061 auto a = (sin_auto(2*theta) +2*theta - M_PI * sin_auto(y));
1062 auto b = (2 + 2 * cos_auto(2*theta));
1063 theta = theta - a / b;
1064 } }
1065 y = M_PI * sin_auto(theta) / 2;
1066 x = x * cos_auto(theta);
1067 });
1068 break;
1069
1070 case mdCentralCyl:
1071 makeband(H_orig, ret, [] (ld& x, ld& y) { y = tan_auto(y); ld top = vid.yres * M_PI / current_display->radius; if(y>top) y=top; if(y<-top) y=-top; });
1072 break;
1073
1074 case mdGallStereographic:
1075 makeband(H_orig, ret, [] (ld& x, ld& y) {
1076 y = 2 * sin_auto(y) / (1 + cos_auto(y));
1077 ld top = vid.yres * M_PI / current_display->radius; if(y>top) y=top; if(y<-top) y=-top;
1078 });
1079 break;
1080
1081 case mdAitoff: case mdHammer: case mdWinkelTripel:
1082 makeband(H_orig, ret, [&] (ld& x, ld& y) {
1083 ld ox = x, oy = y;
1084 x *= pconf.aitoff_parameter;
1085
1086 ld x0 = sin_auto(x) * cos_auto(y);
1087 ld y0 = cos_auto(x) * cos_auto(y);
1088 ld z0 = sin_auto(y);
1089
1090 ld d = acos_auto(y0);
1091 ld d0 = hypot(x0, z0);
1092
1093 if(md == mdAitoff || md == mdWinkelTripel) ;
1094 else if(sphere) d = sqrt(2*(1 - cos(d))) * M_PI / 2;
1095 else d = sqrt(2*(cosh(d) - 1)) / 1.5;
1096
1097 x = x0 * d / d0 / pconf.aitoff_parameter, y = z0 * d / d0;
1098 if(md == mdWinkelTripel)
1099 x = lerp(x, ox, pconf.winkel_parameter),
1100 y = lerp(y, oy, pconf.winkel_parameter);
1101
1102 });
1103 break;
1104
1105 case mdWerner: {
1106
1107 models::apply_orientation_yz(H[1], H[2]);
1108 models::apply_orientation(H[0], H[1]);
1109
1110 find_zlev(H); // ignored for now
1111
1112 ld r = hdist0(H);
1113 if(r == 0) { ret = H; return; }
1114 ld angle = atan2(H[0], H[1]);
1115 angle *= sin_auto(r) / r;
1116
1117 ret[0] = sin(angle) * r;
1118 ret[1] = cos(angle) * r;
1119 ret[2] = 0;
1120 ret[3] = 1;
1121
1122 models::apply_orientation(ret[1], ret[0]);
1123 models::apply_orientation_yz(ret[2], ret[1]);
1124 break;
1125 }
1126
1127 case mdCollignon:
1128 find_zlev(H_orig.h);
1129 makeband(H_orig, ret, [] (ld& x, ld& y) {
1130 ld sgn = 1;
1131 if(pconf.collignon_reflected && y > 0) y = -y, sgn = -1;
1132 y = signed_sqrt(sin_auto(y) + pconf.collignon_parameter);
1133 x *= y / 1.2;
1134 y -= signed_sqrt(pconf.collignon_parameter);
1135 y *= sgn;
1136 y *= M_PI;
1137 });
1138 break;
1139
1140 case mdBandEquiarea:
1141 makeband(H_orig, ret, [] (ld& x, ld& y) { y = sin_auto(y); });
1142 break;
1143
1144 case mdBandEquidistant:
1145 makeband(H_orig, ret, [] (ld& x, ld& y) { });
1146 break;
1147
1148 case mdSinusoidal:
1149 makeband(H_orig, ret, [] (ld& x, ld& y) { x *= cos_auto(y); });
1150 break;
1151
1152 case mdEquidistant: case mdEquiarea: case mdEquivolume: {
1153 if(vrhr::rendering() && GDIM == 3 && pmodel == mdEquidistant) {
1154 ret = inverse_exp(H_orig);
1155 ret[3] = 1;
1156 return;
1157 }
1158
1159 if(nonisotropic || prod) {
1160 ret = lp_apply(inverse_exp(H_orig));
1161 ret[3] = 1;
1162 break;
1163 }
1164 ld zlev = find_zlev(H);
1165
1166 ld rad = hypot_d(GDIM, H);
1167 if(rad == 0) rad = 1;
1168 ld d = hdist0(H);
1169 ld df, zf;
1170 hypot_zlev(zlev, d, df, zf);
1171
1172 // 4 pi / 2pi = M_PI
1173
1174 if(pmodel == mdEquivolume)
1175 d = pow(volume_auto(d), 1/3.) * pow(M_PI / 2, 1/3.);
1176 else if(pmodel == mdEquiarea && sphere)
1177 d = sqrt(2*(1 - cos(d))) * M_PI / 2;
1178 else if(pmodel == mdEquiarea && hyperbolic)
1179 d = sqrt(2*(cosh(d) - 1)) / 1.5;
1180
1181 ld factor = d * df / rad;
1182 if(!vrhr::rendering()) factor /= M_PI;
1183
1184 ret = H * factor;
1185 if(GDIM == 2) ret[2] = 0;
1186 if(MAXMDIM == 4) ret[3] = 1;
1187 if(zlev != 1 && use_z_coordinate())
1188 apply_depth(ret, d * zf / M_PI);
1189
1190 break;
1191 }
1192
1193 case mdRotatedHyperboles: {
1194 // ld zlev = <- not implemented
1195 find_zlev(H); // + vid.depth;
1196 models::apply_orientation(H[0], H[1]);
1197
1198 ld y = asin_auto(H[1]);
1199 ld x = asin_auto_clamp(H[0] / cos_auto(y));
1200 // ld z = zlev == 1 ? 0 : geom3::factor_to_lev(zlev);
1201
1202 ld factor = geom3::lev_to_factor(y + vid.depth);
1203
1204 ret[0] = sinh(x) * factor;
1205 ret[1] = cosh(x) * factor;
1206 ret[2] = 0;
1207
1208 if(pconf.use_atan) {
1209 ret[0] = atan(ret[0]);
1210 ret[1] = atan(ret[1]);
1211 }
1212
1213 break;
1214 }
1215
1216 case mdFormula: {
1217 dynamicval<eModel> m(pmodel, pconf.basic_model);
1218 applymodel(H_orig, ret);
1219 exp_parser ep;
1220 ep.extra_params["z"] = cld(ret[0], ret[1]);
1221 ep.extra_params["cx"] = ret[0];
1222 ep.extra_params["cy"] = ret[1];
1223 ep.extra_params["cz"] = ret[2];
1224 ep.extra_params["ux"] = H[0];
1225 ep.extra_params["uy"] = H[1];
1226 ep.extra_params["uz"] = H[2];
1227 ep.s = pconf.formula;
1228 cld res;
1229 try {
1230 res = ep.parse();
1231 }
1232 catch(hr_parse_exception&) {
1233 res = 0;
1234 }
1235 ret[0] = real(res);
1236 ret[1] = imag(res);
1237 ret[2] = 0;
1238 break;
1239 }
1240
1241 case mdSpiral: {
1242 cld z;
1243 if(hyperbolic || sphere) makeband(H_orig, ret, band_conformal);
1244 else ret = H;
1245 z = cld(ret[0], ret[1]) * models::spiral_multiplier;
1246
1247 if(pconf.spiral_cone < 360) {
1248 ld alpha = imag(z) * 360 / pconf.spiral_cone;
1249 ld r = real(z);
1250 r = exp(r);
1251
1252 ret[0] = -sin(alpha) * r;
1253 ret[1] = cos(alpha) * r;
1254 if(euclid) ret = models::euclidean_spin * ret;
1255 ret[2] = (r-1) * sqrt( pow(360/pconf.spiral_cone, 2) - 1);
1256
1257 models::apply_ball(ret[2], ret[1]);
1258 }
1259 else {
1260 z = exp(z);
1261 ret[0] = real(z);
1262 ret[1] = imag(z);
1263 if(euclid) ret = models::euclidean_spin * ret;
1264
1265 if(pconf.skiprope)
1266 ret = mobius(ret, pconf.skiprope, 1);
1267 }
1268 break;
1269 }
1270
1271 case mdRetroCraig: {
1272 makeband(H_orig, ret, [] (ld& x, ld& y) {
1273 if(x)
1274 y = x / sin_auto(x) * (sin_auto(y) * cos_auto(x) - tan_auto(pconf.loximuthal_parameter) * cos_auto(y));
1275 else
1276 y = sin_auto(y) - tan_auto(pconf.loximuthal_parameter) * cos_auto(y);
1277 });
1278 break;
1279 }
1280
1281 case mdRetroLittrow: {
1282 makeband(H_orig, ret, [] (ld& x, ld& y) {
1283 tie(x, y) = make_pair(
1284 sin_auto(x) / cos_auto(y),
1285 cos_auto(x) * tan_auto(y)
1286 );
1287 });
1288 break;
1289 }
1290
1291 case mdRetroHammer: {
1292 ld d = hdist(H, ypush0(pconf.loximuthal_parameter));
1293 makeband(H_orig, ret, [d,H] (ld& x, ld& y) {
1294 if(x == 0 && y == 0) return;
1295
1296 if(x)
1297 y = x / sin_auto(x) * (sin_auto(y) * cos_auto(x) - tan_auto(pconf.loximuthal_parameter) * cos_auto(y));
1298 else
1299 y = sin_auto(y) - tan_auto(pconf.loximuthal_parameter) * cos_auto(y);
1300
1301 ld scale = d / hypot(x, y);
1302 if(H[2] < 0) scale = -scale;
1303 x *= scale;
1304 y *= scale;
1305 });
1306 break;
1307 }
1308
1309 case mdPanini: {
1310 ld proh = sqrt(H[2]*H[2] + curvature() * H[0] * H[0]);
1311 H /= proh;
1312 H /= (H[2] + pconf.alpha);
1313 ret = H;
1314 ret[2] = 0; ret[3] = 1;
1315 break;
1316 }
1317
1318 case mdPoorMan: {
1319 find_zlev(H);
1320 H = space_to_perspective(H);
1321
1322 models::apply_orientation_yz(H[1], H[2]);
1323 models::apply_orientation(H[0], H[1]);
1324
1325 ld u = H[0], v = H[1];
1326 if(abs(u) > 1e-3 && abs(v) > 1e-3) {
1327 ld r2 = u*u+v*v;
1328 ld scale = sqrt((-r2+sqrt(r2*(r2+4*u*u*v*v*(r2-2))))/(2*(r2-2))) / u / v;
1329 if(u*v<0) scale = -scale;
1330 H = scale * H;
1331 }
1332 ret = H;
1333 ret[2] = 0;
1334 ret[3] = 1;
1335
1336 models::apply_orientation(ret[1], ret[0]);
1337 models::apply_orientation_yz(ret[2], ret[1]);
1338 break;
1339 }
1340
1341 case mdGUARD: case mdManual: break;
1342 }
1343
1344 ghcheck(ret,H_orig);
1345 }
1346
1347 // game-related graphics
1348
1349 EX transmatrix sphereflip; // on the sphere, flip
1350 EX bool playerfound; // has player been found in the last drawing?
1351
outofmap(hyperpoint h)1352 EX bool outofmap(hyperpoint h) {
1353 if(GDIM == 3)
1354 return false;
1355 else if(euclid)
1356 return h[2] < .5; // false; // h[0] * h[0] + h[1] * h[1] > 15 * cgi.crossf;
1357 else if(sphere)
1358 return h[2] < .1 && h[2] > -.1 && h[1] > -.1 && h[1] < .1 && h[0] > -.1 && h[0] < .1;
1359 else
1360 return h[2] < .5;
1361 }
1362
mirrorif(const hyperpoint & V,bool b)1363 EX hyperpoint mirrorif(const hyperpoint& V, bool b) {
1364 if(b) return Mirror*V;
1365 else return V;
1366 }
1367
mirrorif(const shiftmatrix & V,bool b)1368 EX shiftmatrix mirrorif(const shiftmatrix& V, bool b) {
1369 if(b) return V*Mirror;
1370 else return V;
1371 }
1372
1373 // -1 if away, 0 if not away
away(const transmatrix & V2)1374 EX int away(const transmatrix& V2) {
1375 return (intval(C0, V2 * xpush0(.1)) > intval(C0, tC0(V2))) ? -1 : 0;
1376 }
1377
1378 /* double zgrad(double f1, double f2, int nom, int den) {
1379 using namespace geom3;
1380 ld fo1 = factor_to_lev(f1);
1381 ld fo2 = factor_to_lev(f2);
1382 return lev_to_factor(fo1 + (fo2-fo1) * nom / den);
1383 } */
1384
zgrad0(double l1,double l2,int nom,int den)1385 EX double zgrad0(double l1, double l2, int nom, int den) {
1386 using namespace geom3;
1387 return lev_to_factor(l1 + (l2-l1) * nom / den);
1388 }
1389
behindsphere(const hyperpoint & h)1390 EX bool behindsphere(const hyperpoint& h) {
1391 if(!sphere) return false;
1392
1393 if(mdBandAny()) return false;
1394
1395 if(pconf.alpha > 1) {
1396 if(h[LDIM] > -1/pconf.alpha) return true;
1397 }
1398
1399 if(pconf.alpha <= 1) {
1400 if(h[LDIM] < .2-pconf.alpha) return true;
1401 }
1402
1403 return false;
1404 }
1405
to01(ld a0,ld a1,ld x)1406 ld to01(ld a0, ld a1, ld x) {
1407 if(x < a0) return 0;
1408 if(x > a1) return 1;
1409 return (x-a0) / (a1-a0);
1410 }
1411
spherity(const hyperpoint & h)1412 EX ld spherity(const hyperpoint& h) {
1413 if(!sphere) return 1;
1414
1415 if(pconf.alpha > 1) {
1416 return to01(1/pconf.alpha, 1, abs(h[2]));
1417 }
1418
1419 if(pconf.alpha <= 1) {
1420 return to01(-1.5, 1, h[2]);
1421 }
1422
1423 return 1;
1424 }
1425
behindsphere(const transmatrix & V)1426 EX bool behindsphere(const transmatrix& V) {
1427 return behindsphere(tC0(V));
1428 }
1429
behindsphere(const shiftmatrix & V)1430 EX bool behindsphere(const shiftmatrix& V) {
1431 return behindsphere(tC0(V.T));
1432 }
1433
spherity(const transmatrix & V)1434 EX ld spherity(const transmatrix& V) {
1435 return spherity(tC0(V));
1436 }
1437
confusingGeometry()1438 EX bool confusingGeometry() {
1439 #if MAXMDIM >= 4
1440 if(reg3::ultra_mirror_in()) return true;
1441 #endif
1442 return quotient || elliptic || (fake::in() && fake::multiple);
1443 }
1444
master_to_c7_angle()1445 EX ld master_to_c7_angle() {
1446 if(hybri) return hybrid::in_underlying_geometry(master_to_c7_angle);
1447 if(WDIM == 3) return 0;
1448 ld alpha = 0;
1449 #if CAP_GP
1450 if(cgi.gpdata) alpha = cgi.gpdata->alpha;
1451 #endif
1452 return (!BITRUNCATED && !bt::in() && !arcm::in()) ? M_PI + alpha : 0;
1453 }
1454
actualV(const heptspin & hs,const transmatrix & V)1455 EX transmatrix actualV(const heptspin& hs, const transmatrix& V) {
1456 if(prod) return PIU(actualV(hs, V));
1457 if(WDIM == 3) return V;
1458 #if CAP_IRR
1459 if(IRREGULAR)
1460 return V * spin(M_PI + 2 * M_PI / S7 * (hs.spin + irr::periodmap[hs.at].base.spin));
1461 #endif
1462 #if CAP_ARCM
1463 if(arcm::in()) return V * spin(-arcm::current.triangles[arcm::id_of(hs.at)][hs.spin].first);
1464 #endif
1465 #if CAP_BT
1466 if(bt::in()) return V;
1467 #endif
1468 if(kite::in()) return V;
1469 return (hs.spin || !BITRUNCATED) ? V * spin(hs.spin*2*M_PI/hs.at->type + master_to_c7_angle()) : V;
1470 }
1471
actualV(const heptspin & hs,const shiftmatrix & V)1472 EX shiftmatrix actualV(const heptspin& hs, const shiftmatrix& V) {
1473 return shiftless(actualV(hs, V.T), V.shift);
1474 }
1475
point_behind(const shiftpoint h)1476 EX bool point_behind(const shiftpoint h) {
1477 if(sphere) return false;
1478 if(!in_perspective()) return false;
1479 hyperpoint h1;
1480 if(pmodel == mdGeodesic) h1 = inverse_exp(h, pQUICK);
1481 if(pmodel == mdPerspective && prod) h1 = product::inverse_exp(h.h);
1482 h1 = lp_apply(h1);
1483 return h1[2] < 1e-8;
1484 }
1485
raise_error()1486 void raise_error() {
1487 println(hlog, "something wrong");
1488 }
1489
invalid_matrix(const transmatrix T)1490 EX bool invalid_matrix(const transmatrix T) {
1491 for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++)
1492 if(std::isnan(T[i][j]) || T[i][j] > 1e8 || T[i][j] < -1e8 || std::isinf(T[i][j]))
1493 return true;
1494 if(prod || (cgflags & qAFFINE)) {
1495 for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++) if(abs(T[i][j]) > 1e-60) return false;
1496 }
1497 else
1498 for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++) if(T[i][j] > .5 || T[i][j] < -.5) return false;
1499 return true;
1500 }
1501
invalid_point(const hyperpoint h)1502 EX bool invalid_point(const hyperpoint h) {
1503 return std::isnan(h[LDIM]) || h[LDIM] > 1e8 || std::isinf(h[LDIM]);
1504 }
1505
invalid_point(const shiftpoint h)1506 EX bool invalid_point(const shiftpoint h) { return invalid_point(h.h); }
1507
in_smart_range(const shiftmatrix & T)1508 EX bool in_smart_range(const shiftmatrix& T) {
1509 shiftpoint h = tC0(T);
1510 if(invalid_point(h)) return false;
1511 if(nil || nih) return true;
1512 #if CAP_SOLV
1513 if(pmodel == mdGeodesic) return sn::in_table_range(h.h);
1514 #endif
1515 hyperpoint h1;
1516 applymodel(h, h1);
1517 if(invalid_point(h1)) return false;
1518 ld x = current_display->xcenter + current_display->radius * h1[0];
1519 ld y = current_display->ycenter + current_display->radius * h1[1] * pconf.stretch;
1520
1521 if(x > current_display->xtop + current_display->xsize * 2) return false;
1522 if(x < current_display->xtop - current_display->xsize * 1) return false;
1523 if(y > current_display->ytop + current_display->ysize * 2) return false;
1524 if(y < current_display->ytop - current_display->ysize * 1) return false;
1525 if(GDIM == 3) {
1526 if(-h1[2] < pconf.clip_min * 2 - pconf.clip_max) return false;
1527 if(-h1[2] > pconf.clip_max * 2 - pconf.clip_min) return false;
1528 }
1529
1530 ld epsilon = 0.01;
1531
1532 transmatrix ar;
1533
1534 ld dx = 0, dy = 0, dz = 0, dh[MAXMDIM];
1535 for(int i=0; i<GDIM; i++) {
1536 hyperpoint h2;
1537 applymodel(T * cpush0(i, epsilon), h2);
1538 ld x1 = current_display->radius * abs(h2[0] - h1[0]) / epsilon;
1539 ld y1 = current_display->radius * abs(h2[1] - h1[1]) * pconf.stretch / epsilon;
1540
1541 for(int j=0; j<GDIM; j++) ar[i][j] = current_display->radius * (h2[j]-h1[j]) / epsilon;
1542
1543 dx = max(dx, x1); dy = max(dy, y1);
1544 if(GDIM == 3) dz = max(dz, abs(h2[2] - h1[2]));
1545 dh[i] = hypot(x1, y1);
1546 }
1547
1548 if(GDIM == 2 && vid.smart_area_based) {
1549 ld area = det2(ar);
1550 ld scale = sqrt(area) * cgi.scalefactor * hcrossf7;
1551 if(scale <= vid.smart_range_detail) return false;
1552 }
1553
1554 else if(GDIM == 3) {
1555 if(-h1[2] + 2 * dz < pconf.clip_min || -h1[2] - 2 * dz > pconf.clip_max) return false;
1556 sort(dh, dh+GDIM);
1557 ld scale = sqrt(dh[1] * dh[2]) * cgi.scalefactor * hcrossf7;
1558 if(scale <= (WDIM == 2 ? vid.smart_range_detail : vid.smart_range_detail_3)) return false;
1559 }
1560 else {
1561 ld scale = sqrt(dh[0] * dh[1]) * cgi.scalefactor * hcrossf7;
1562 if(scale <= vid.smart_range_detail) return false;
1563 }
1564
1565 return
1566 x - 2 * dx < current_display->xtop + current_display->xsize &&
1567 x + 2 * dx > current_display->xtop &&
1568 y - 2 * dy < current_display->ytop + current_display->ysize &&
1569 y + 2 * dy > current_display->ytop;
1570 }
1571
1572 #if CAP_GP
1573 namespace gp {
1574
1575 /*
1576 void drawrec(cell *c, const transmatrix& V) {
1577 if(dodrawcell(c))
1578 drawcell(c, V, 0, false);
1579 for(int i=0; i<c->type; i++) {
1580 cell *c2 = c->move(i);
1581 if(!c2) continue;
1582 if(c2->move(0) != c) continue;
1583 if(c2 == c2->master->c7) continue;
1584 transmatrix V1 = V * ddspin(c, i) * xpush(cgi.crossf) * iddspin(c2, 0) * spin(M_PI);
1585 drawrec(c2, V1);
1586 }
1587 } */
1588
drawrec(cell * c,const shiftmatrix & V,gp::loc at,int dir,int maindir,local_info li)1589 bool drawrec(cell *c, const shiftmatrix& V, gp::loc at, int dir, int maindir, local_info li) {
1590 bool res = false;
1591 shiftmatrix V1 = V * cgi.gpdata->Tf[li.last_dir][at.first&GOLDBERG_MASK][at.second&GOLDBERG_MASK][fixg6(dir)];
1592 if(do_draw(c, V1)) {
1593 /* auto li = get_local_info(c);
1594 if((dir - li.total_dir) % S6) printf("totaldir %d/%d\n", dir, li.total_dir);
1595 if(at != li.relative) printf("at %s/%s\n", disp(at), disp(li.relative));
1596 if(maindir != li.last_dir) printf("ld %d/%d\n", maindir, li.last_dir); */
1597 li.relative = at;
1598 li.total_dir = fixg6(dir);
1599 current_li = li;
1600 li_for = c;
1601 drawcell(c, V1);
1602 res = true;
1603 }
1604 for(int i=0; i<c->type; i++) {
1605 cell *c2 = c->move(i);
1606 if(!c2) continue;
1607 if(c2->move(0) != c) continue;
1608 if(c2 == c2->master->c7) continue;
1609 res |= drawrec(c2, V, at + eudir(dir+i), dir + i + SG3, maindir, li);
1610 }
1611 return res;
1612 }
1613
drawrec(cell * c,const shiftmatrix & V)1614 bool drawrec(cell *c, const shiftmatrix& V) {
1615 local_info li;
1616 li.relative = loc(0,0);
1617 li.total_dir = 0;
1618 li.last_dir = -1;
1619 li.first_dir = -1;
1620 li_for = c;
1621 current_li = li;
1622 bool res = false;
1623 if(do_draw(c, V))
1624 drawcell(c, V), res = true;
1625 for(int i=0; i<c->type; i++) {
1626 cell *c2 = c->move(i);
1627 if(!c2) continue;
1628 if(c2->move(0) != c) continue;
1629 if(c2 == c2->master->c7) continue;
1630 li.last_dir = i;
1631 res |= drawrec(c2, V, gp::loc(1,0), SG3, i, li);
1632 }
1633 return res;
1634 }
1635 }
1636 #endif
1637
1638 vector<tuple<heptspin, hstate, shiftmatrix> > drawn_cells;
1639
drawcell_subs(cell * c,const shiftmatrix & V)1640 EX bool drawcell_subs(cell *c, const shiftmatrix& V) {
1641
1642 #if CAP_GP
1643 if(GOLDBERG) {
1644 return gp::drawrec(c, V);
1645 }
1646 #endif
1647
1648 bool draw = false;
1649
1650 #if CAP_IRR
1651 if(IRREGULAR) {
1652 auto& hi = irr::periodmap[c->master];
1653 auto& vc = irr::cells_of_heptagon[hi.base.at];
1654 for(int i=0; i<isize(vc); i++) {
1655 cell *c = hi.subcells[i];
1656 shiftmatrix V1 = V * irr::cells[vc[i]].pusher;
1657 if(do_draw(c, V1))
1658 draw = true,
1659 drawcell(hi.subcells[i], V * irr::cells[vc[i]].pusher);
1660 }
1661 return draw;
1662 }
1663 #endif
1664
1665 if(do_draw(c, V)) {
1666 draw = true;
1667 drawcell(c, V);
1668 }
1669
1670 if(BITRUNCATED) forCellIdEx(c1, d, c) {
1671 if(c->c.spin(d) == 0) {
1672 shiftmatrix V2 = V * currentmap->adj(c, d);
1673 if(do_draw(c1, V2))
1674 draw = true,
1675 drawcell(c1, V2);
1676 }
1677 }
1678
1679 return draw;
1680 }
1681
draw_all()1682 void hrmap::draw_all() {
1683 if(sphere && pmodel == mdSpiral) {
1684 if(models::ring_not_spiral) {
1685 int qty = ceil(1. / pconf.sphere_spiral_multiplier);
1686 if(qty > 100) qty = 100;
1687 for(int i=-qty; i < qty; i++)
1688 draw_at(centerover, cview(2 * M_PI * i));
1689 }
1690 else {
1691 draw_at(centerover, cview());
1692 if(vid.use_smart_range) for(int i=1;; i++) {
1693 int drawn = cells_drawn;
1694 draw_at(centerover, cview(2 * M_PI * i));
1695 draw_at(centerover, cview(-2 * M_PI * i));
1696 if(drawn == cells_drawn) break;
1697 }
1698 }
1699 }
1700 else
1701 draw_at(centerover, cview());
1702 }
1703
draw_at(cell * at,const shiftmatrix & where)1704 void hrmap::draw_at(cell *at, const shiftmatrix& where) {
1705 dq::clear_all();
1706 auto& enq = confusingGeometry() ? dq::enqueue_by_matrix_c : dq::enqueue_c;
1707
1708 enq(at, where);
1709
1710 while(!dq::drawqueue_c.empty()) {
1711 auto& p = dq::drawqueue_c.front();
1712 cell *c = p.first;
1713 shiftmatrix V = p.second;
1714 dq::drawqueue_c.pop();
1715
1716 if(!do_draw(c, V)) continue;
1717 drawcell(c, V);
1718 if(in_wallopt() && isWall3(c) && isize(dq::drawqueue) > 1000) continue;
1719
1720 #if MAXMDIM >= 4
1721 if(reg3::ultra_mirror_in())
1722 for(auto& T: cgi.ultra_mirrors)
1723 enq(c, optimized_shift(V * T));
1724 #endif
1725
1726 for(int i=0; i<c->type; i++) {
1727 // note: need do cmove before c.spin
1728 cell *c1 = c->cmove(i);
1729 if(c1 == &out_of_bounds) continue;
1730 enq(c1, optimized_shift(V * adj(c, i)));
1731 }
1732 }
1733 }
1734
draw_at(cell * at,const shiftmatrix & where)1735 void hrmap_standard::draw_at(cell *at, const shiftmatrix& where) {
1736 if(S3 > 4) {
1737 hrmap::draw_at(at, where);
1738 return;
1739 }
1740 drawn_cells.clear();
1741 drawn_cells.emplace_back(at->master, hsOrigin, where * master_relative(at, true));
1742 for(int i=0; i<isize(drawn_cells); i++) {
1743 // prevent reallocation due to insertion
1744 if(drawn_cells.capacity() < drawn_cells.size() + 16)
1745 drawn_cells.reserve(max<size_t>(2 * drawn_cells.size(), 128));
1746
1747 const auto& dc = drawn_cells[i];
1748 auto& hs = get<0>(dc);
1749 auto& s = get<1>(dc);
1750 auto& V = get<2>(dc);
1751
1752 cell *c = hs.at->c7;
1753
1754 const shiftmatrix& V1 = hs.mirrored ? V * Mirror : V;
1755
1756 bool draw = drawcell_subs(c, actualV(hs, V1));
1757
1758 if(sphere) draw = true;
1759
1760 if(draw) for(int d=0; d<c->master->type; d++) {
1761 hstate s2 = transition(s, d);
1762 if(s2 == hsError) continue;
1763 heptspin hs2 = hs + d + wstep;
1764 shiftmatrix Vd;
1765 if(inforder::mixed()) {
1766 int d1 = gmod(hs.spin+d, c->type);
1767 Vd = V * spin(-2*M_PI*d/c->type) * xpush(spacedist(c, d1)) * spin(M_PI);
1768 }
1769 else
1770 Vd = V * cgi.heptmove[d];
1771 optimize_shift(Vd);
1772 drawn_cells.emplace_back(hs2, s2, Vd);
1773 }
1774 }
1775 }
1776
keep_vertical()1777 EX bool keep_vertical() {
1778 if(CAP_ORIENTATION) return false;
1779 if((WDIM == 2 || prod) && GDIM == 3 && vid.fixed_yz) return true;
1780 if(downseek.qty) return true;
1781 return false;
1782 }
1783
vertical_vector()1784 EX hyperpoint vertical_vector() {
1785 auto& ds = downseek;
1786 if((WDIM == 2 || prod) && GDIM == 3 && vid.fixed_yz)
1787 return get_view_orientation() * ztangent(1);
1788 else if(ds.qty && prod)
1789 return get_view_orientation() * product::inverse_exp(ds.point);
1790 else if(ds.qty)
1791 return ds.point;
1792 return C0;
1793 }
1794
1795 EX bool down_is_forward;
1796
spinEdge(ld aspd)1797 EX void spinEdge(ld aspd) {
1798
1799 #if CAP_VR
1800 if(vrhr::active() && keep_vertical() && !vrhr::first) {
1801 transmatrix T = vrhr::hmd_ref_at;
1802 T = vrhr::sm * inverse(T);
1803 vrhr::be_33(T);
1804
1805 transmatrix V = T * get_view_orientation();
1806
1807 hyperpoint h = inverse(V) * C0;
1808 if(!prod) {
1809 V = V * rgpushxto0(h);
1810 }
1811
1812 int dir = down_is_forward ? 0 : 1;
1813
1814 V = cspin(2, dir, 90 * degree) * V;
1815
1816 if(1) {
1817 dynamicval<eGeometry> g(geometry, gSphere);
1818 bool b = vid.always3;
1819 vid.always3 = false;
1820 geom3::apply_always3();
1821 V = gpushxto0(V*C0) * V;
1822 fixmatrix(V);
1823 if(b) {
1824 vid.always3 = b;
1825 geom3::apply_always3();
1826 }
1827 }
1828
1829 vrhr::be_33(V);
1830
1831 V = cspin(dir, 2, 90 * degree) * V;
1832 V = inverse(T) * V;
1833 if(!prod) V = V * gpushxto0(h);
1834 get_view_orientation() = V;
1835 return;
1836 }
1837 #endif
1838
1839 ld downspin = 0;
1840 auto& ds = downseek;
1841 if(dual::state == 2 && (dual::one_euclidean ? !euclid : dual::currently_loaded != dual::main_side)) {
1842 transmatrix our = dual::get_orientation();
1843 transmatrix their = dual::player_orientation[dual::main_side];
1844 fixmatrix(our);
1845 fixmatrix(their);
1846 if(GDIM == 2) {
1847 transmatrix T = their * iso_inverse(our);
1848 hyperpoint H = T * xpush0(1);
1849 downspin = -atan2(H[1], H[0]);
1850 }
1851 else rotate_view(their * iso_inverse(our));
1852 }
1853 else if(playerfound && vid.fixed_facing) {
1854 hyperpoint H = gpushxto0(unshift(playerV) * C0) * unshift(playerV) * xpush0(5);
1855 downspin = atan2(H[1], H[0]);
1856 downspin += vid.fixed_facing_dir * degree;
1857 if(flipplayer) downspin += M_PI;
1858 cyclefix(downspin, 0);
1859 aspd = (1 + 2 * abs(downspin)) * aspd;
1860 }
1861 else if(keep_vertical()) {
1862 hyperpoint h = vertical_vector();
1863 downspin = -atan2(h[0], h[1]);
1864 if(ds.qty && GDIM == 2) {
1865 downspin += models::rotation * degree;
1866 }
1867 if(ds.qty) {
1868 cyclefix(downspin, 0);
1869 downspin = downspin * min(ds.speed, (double)1);
1870 }
1871 else aspd = 999999;
1872 }
1873 if(downspin > aspd) downspin = aspd;
1874 if(downspin < -aspd) downspin = -aspd;
1875 rotate_view(spin(downspin));
1876 }
1877
1878 /** \brief convert a shiftmatrix to the coordinate system of View
1879 * usually used to set which_copy
1880 */
back_to_view(const shiftmatrix & V)1881 EX transmatrix back_to_view(const shiftmatrix& V) {
1882 // ortho_inverse does not work in 2.5D, iso_inverse does not work in Nil.
1883 // just use inverse
1884 return inverse(actual_view_transform) * unshift(V);
1885 }
1886
fix_whichcopy(cell * c)1887 EX void fix_whichcopy(cell *c) {
1888 if(!gmatrix.count(cwt.at)) return;
1889 current_display->which_copy = back_to_view(gmatrix[c]);
1890 }
1891
fix_whichcopy_if_near()1892 void fix_whichcopy_if_near() {
1893 if(!gmatrix.count(cwt.at)) return;
1894 transmatrix T = back_to_view(gmatrix[cwt.at]);
1895 if(!eqmatrix(T, current_display->which_copy)) return;
1896 current_display->which_copy = T;
1897 }
1898
centerpc(ld aspd)1899 EX void centerpc(ld aspd) {
1900
1901 if(subscreens::split([=] () {centerpc(aspd);})) return;
1902 if(dual::split([=] () { centerpc(aspd); })) return;
1903
1904 #if CAP_CRYSTAL && CAP_RUG
1905 if(cryst)
1906 crystal::centerrug(aspd);
1907 #endif
1908
1909 #if CAP_RACING
1910 if(racing::on && racing::set_view()) return;
1911 #endif
1912
1913 #if MAXMDIM >= 4
1914 if(shmup::on && vid.sspeed > -5 && GDIM == 3) {
1915 int id = subscreens::in ? subscreens::current_player : 0;
1916 auto& pc = shmup::pc[id];
1917 centerover = pc->base;
1918 transmatrix T = pc->at;
1919 int sl = snakelevel(cwt.at);
1920 if((sl || vid.eye) && WDIM == 2) T = T * zpush(cgi.SLEV[sl] - cgi.FLOOR + vid.eye);
1921 /* in nonisotropic geometries, T is isometry * rotation, so iso_inverse does not work */
1922 if(nonisotropic)
1923 View = inverse(T);
1924 else
1925 View = iso_inverse(T);
1926 if(prod) NLP = ortho_inverse(pc->ori);
1927 if(WDIM == 2) rotate_view( cspin(0, 1, M_PI) * cspin(2, 1, M_PI/2 + shmup::playerturny[id]) * spin(-M_PI/2) );
1928 return;
1929 }
1930 #endif
1931
1932 if(ors::mode == 2 && vid.sspeed < 5) return;
1933 if(vid.sspeed >= 4.99) aspd = 1000;
1934 DEBBI(DF_GRAPH, ("center pc"));
1935
1936 auto& W = current_display->which_copy;
1937 ors::unrotate(W); ors::unrotate(View); ors::unrotate(cwtV.T);
1938
1939 /* what should we center? */
1940 transmatrix T;
1941 if(multi::players > 1)
1942 T = unshift(cwtV); /* do not even try */
1943 else {
1944 T = W;
1945 if(shmup::on)
1946 T = T * shmup::pc[0]->at;
1947 }
1948
1949 if(invalid_matrix(T)) return;
1950
1951 #if MAXMDIM >= 4
1952 if(GDIM == 3 && WDIM == 2) {
1953 geom3::do_auto_eye();
1954 int sl = snakelevel(cwt.at);
1955 if(isWorm(cwt.at->monst) && sl < 3) sl++;
1956 if(sl || vid.eye) T = T * zpush(cgi.SLEV[sl] - cgi.FLOOR + vid.eye);
1957 }
1958 #endif
1959
1960 hyperpoint H = tC0(T);
1961 ld R = (zero_d(GDIM, H) && !prod) ? 0 : hdist0(H);
1962 if(R < 1e-9) {
1963 // either already centered or direction unknown
1964 /* if(playerfoundL && playerfoundR) {
1965
1966 } */
1967 spinEdge(aspd);
1968 fixmatrix(View);
1969 fix_whichcopy(cwt.at);
1970 fixmatrix(current_display->which_copy);
1971 }
1972
1973 else {
1974 aspd *= euclid ? (2+3*R*R) : (1+R+(shmup::on?1:0));
1975
1976 if(R < aspd) fix_whichcopy_if_near();
1977
1978 if(R < aspd)
1979 shift_view_to(shiftless(H));
1980 else
1981 shift_view_towards(shiftless(H), aspd);
1982
1983 fixmatrix(View);
1984 fixmatrix(current_display->which_copy);
1985 spinEdge(aspd);
1986 }
1987
1988 ors::rerotate(W); ors::rerotate(cwtV.T); ors::rerotate(View);
1989 }
1990
1991 EX transmatrix oView;
1992
1993 EX purehookset hooks_preoptimize, hooks_postoptimize;
1994
optimizeview()1995 EX void optimizeview() {
1996
1997 if(subscreens::split(optimizeview)) return;
1998 if(dual::split(optimizeview)) return;
1999
2000 cell *c = centerover;
2001 transmatrix iView = view_inverse(View);
2002 callhooks(hooks_preoptimize);
2003 virtualRebase(centerover, iView);
2004 if(c != centerover && (sphere || sl2)) {
2005 transmatrix T = currentmap->relative_matrix(centerover, c, C0);
2006 T = stretch::itranslate(tC0(T)) * T;
2007 stretch::mstretch_matrix = T * stretch::mstretch_matrix;
2008 }
2009
2010 View = iview_inverse(iView);
2011 fixmatrix(View);
2012 callhooks(hooks_postoptimize);
2013
2014 if(is_boundary(centerover))
2015 centerover = c, View = oView;
2016 else
2017 oView = View;
2018
2019 #if CAP_ANIMATIONS
2020 if(centerover && inmirror(centerover)) {
2021 anims::reflect_view();
2022 }
2023 #endif
2024 }
2025
addball(ld a,ld b,ld c)2026 void addball(ld a, ld b, ld c) {
2027 hyperpoint h;
2028 ballmodel(h, a, b, c);
2029 for(int i=0; i<3; i++) h[i] *= current_display->radius;
2030 curvepoint(h);
2031 }
2032
ballgeometry()2033 void ballgeometry() {
2034 queuereset(mdPixel, PPR::CIRCLE);
2035 for(int i=0; i<60; i++)
2036 addball(i * M_PI/30, 10, 0);
2037 for(double d=10; d>=-10; d-=.2)
2038 addball(0, d, 0);
2039 for(double d=-10; d<=10; d+=.2)
2040 addball(0, d, vid.depth);
2041 addball(0, 0, -vid.camera);
2042 addball(0, 0, vid.depth);
2043 addball(0, 0, -vid.camera);
2044 addball(0, -10, 0);
2045 addball(0, 0, -vid.camera);
2046 queuecurve(shiftless(Id), darkena(0xFF, 0, 0x80), 0, PPR::CIRCLE);
2047 queuereset(pmodel, PPR::CIRCLE);
2048 }
2049
resetview()2050 EX void resetview() {
2051 DEBBI(DF_GRAPH, ("reset view"));
2052 // EUCLIDEAN
2053 NLP = Id;
2054 stretch::mstretch_matrix = Id;
2055 auto lView = View;
2056 if(cwt.at) {
2057 centerover = cwt.at;
2058 View = iddspin(cwt.at, cwt.spin);
2059 if(!flipplayer) View = pispin * View;
2060 if(cwt.mirrored) View = Mirror * View;
2061
2062 if(centering) {
2063 hyperpoint vl = View * get_corner_position(cwt.at, cwt.spin);
2064 hyperpoint vr = View * get_corner_position(cwt.at, (cwt.spin+1) % cwt.at->type);
2065
2066 hyperpoint vm = (centering == eCentering::edge) ? mid(vl, vr) : vl;
2067
2068 transmatrix rm = gpushxto0(vm);
2069
2070 View = spintox(rm*vr) * rm * View;
2071 }
2072
2073 if(GDIM == 2) View = spin(M_PI + vid.fixed_facing_dir * degree) * View;
2074 if(GDIM == 3 && !prod) View = cspin(0, 2, M_PI/2) * View;
2075 if(prod) NLP = cspin(0, 2, M_PI/2);
2076 if(cheater && eqmatrix(View, lView) && !centering) {
2077 View = Id;
2078 static ld cyc = 0;
2079 cyc += 90 * degree;
2080 View = spin(cyc) * View;
2081 if(GDIM == 2) View = spin(M_PI + vid.fixed_facing_dir * degree) * View;
2082 if(GDIM == 3 && !prod) View = cspin(0, 2, M_PI/2) * View;
2083 }
2084 }
2085 else if(currentmap) {
2086 centerover = currentmap->gamestart();
2087 View = Id;
2088 }
2089 cwtV = shiftless(View);
2090 current_display->which_copy =
2091 nonisotropic ? gpushxto0(tC0(view_inverse(View))) :
2092 View;
2093 // SDL_LockSurface(s);
2094 // SDL_UnlockSurface(s);
2095 }
2096
2097
panning(shiftpoint hf0,shiftpoint ht0)2098 EX void panning(shiftpoint hf0, shiftpoint ht0) {
2099 hyperpoint hf = hf0.h;
2100 hyperpoint ht = unshift(ht0, hf0.shift);
2101 View =
2102 rgpushxto0(hf) * rgpushxto0(gpushxto0(hf) * ht) * gpushxto0(hf) * View;
2103 playermoved = false;
2104 }
2105
2106 EX int cells_drawn, cells_generated;
2107
fullcenter()2108 EX void fullcenter() {
2109 if(history::saved_ends == 0)
2110 history::path_for_lineanimation.clear();
2111 if(playerfound && false) centerpc(INF);
2112 else {
2113 bfs();
2114 resetview();
2115 drawthemap();
2116 if(!centering) centerpc(INF);
2117 centerover = cwt.at;
2118 }
2119 playermoved = !centering;
2120 }
2121
screenpos(ld x,ld y)2122 transmatrix screenpos(ld x, ld y) {
2123 transmatrix V = Id;
2124 V[0][2] += (x - current_display->xcenter) / current_display->radius * (1+pconf.alpha);
2125 V[1][2] += (y - current_display->ycenter) / current_display->radius * (1+pconf.alpha);
2126 return V;
2127 }
2128
2129 /**
2130 In 3D, we use the standard translation matrices to place stuff on the screen.
2131 In 2D, this does not work (as HyperRogue reduces matrices to 3x3) so we use the native disk projection
2132 */
2133
2134 EX int flat_on;
2135 eGeometry backup_geometry;
2136 eVariation backup_variation;
2137 videopar backup_vid;
2138
2139 /** \brief enable the 'flat' model for drawing HUD. See hr::flat_model_enabler */
enable_flat_model(int val)2140 EX void enable_flat_model(int val) {
2141 if(flat_on < 1 && flat_on + val >= 1) {
2142 #if CAP_GL
2143 glClear(GL_DEPTH_BUFFER_BIT);
2144 #endif
2145 backup_geometry = geometry;
2146 backup_variation = variation;
2147 backup_vid = vid;
2148 geometry = gNormal;
2149 variation = eVariation::bitruncated;
2150 pmodel = mdDisk;
2151 pconf.alpha = 1;
2152 pconf.scale = 1;
2153 pconf.camera_angle = 0;
2154 pconf.stretch = 1;
2155
2156 vid.always3 = false;
2157 vid.wall_height = .3;
2158 vid.human_wall_ratio = .7;
2159 vid.camera = 1;
2160 vid.depth = 1;
2161 geom3::apply_always3();
2162 check_cgi();
2163 cgi.require_shapes();
2164 calcparam();
2165 }
2166 if(flat_on >= 1 && flat_on + val < 1) {
2167 geometry = backup_geometry;
2168 variation = backup_variation;
2169 vid = backup_vid;
2170 geom3::apply_always3();
2171 calcparam();
2172 check_cgi();
2173 }
2174 flat_on += val;
2175 }
2176
2177 #if HDR
2178 struct flat_model_enabler {
flat_model_enablerhr::flat_model_enabler2179 flat_model_enabler() { enable_flat_model(+1); }
~flat_model_enablerhr::flat_model_enabler2180 ~flat_model_enabler() { enable_flat_model(-1); }
2181 };
2182 #endif
2183
atscreenpos(ld x,ld y,ld size)2184 EX transmatrix atscreenpos(ld x, ld y, ld size) {
2185 transmatrix V = Id;
2186
2187 if(pmodel == mdPixel) {
2188 V[0][3] += (x - current_display->xcenter);
2189 V[1][3] += (y - current_display->ycenter);
2190 V[0][0] = size * 2 * cgi.hcrossf / cgi.crossf;
2191 V[1][1] = size * 2 * cgi.hcrossf / cgi.crossf;
2192 if(WDIM == 3) V[2][2] = -1;
2193 }
2194 else if(pmodel == mdHorocyclic) {
2195 V[0][3] += (x - current_display->xcenter) * 2 / current_display->radius;
2196 V[1][3] += (y - current_display->ycenter) * 2/ current_display->radius;
2197 V[0][0] = size * 2 / current_display->radius;
2198 V[1][1] = size * 2 / current_display->radius;
2199 }
2200 else {
2201 V[0][2] += (x - current_display->xcenter);
2202 V[1][2] += (y - current_display->ycenter);
2203 V[0][0] = size * 2 * cgi.hcrossf / cgi.crossf;
2204 V[1][1] = size * 2 * cgi.hcrossf / cgi.crossf;
2205 V[2][2] = current_display->radius;
2206 if(S3 >= OINF) V[0][0] /= 5, V[1][1] /= 5;
2207 }
2208
2209 return V;
2210 }
2211
circle_around_center(ld radius,color_t linecol,color_t fillcol,PPR prio)2212 void circle_around_center(ld radius, color_t linecol, color_t fillcol, PPR prio) {
2213 #if CAP_QUEUE
2214 if(among(pmodel, mdDisk, mdEquiarea, mdEquidistant, mdFisheye) && !(pmodel == mdDisk && hyperbolic && pconf.alpha <= -1) && pconf.camera_angle == 0) {
2215 hyperpoint ret;
2216 applymodel(shiftless(xpush0(radius)), ret);
2217 ld r = hypot_d(2, ret);
2218 queuecircle(current_display->xcenter, current_display->ycenter, r * current_display->radius, linecol, prio, fillcol);
2219 return;
2220 }
2221 #endif
2222 #if CAP_QUEUE
2223 for(int i=0; i<=360; i++) curvepoint(xspinpush0(i * degree, 10));
2224 auto& c = queuecurve(shiftless(Id), linecol, fillcol, prio);
2225 if(pmodel == mdDisk && hyperbolic && pconf.alpha <= -1)
2226 c.flags |= POLY_FORCE_INVERTED;
2227 if(pmodel == mdJoukowsky)
2228 c.flags |= POLY_FORCE_INVERTED;
2229 c.flags |= POLY_ALWAYS_IN;
2230 #endif
2231 }
2232
2233 EX color_t periodcolor = 0x00FF0080;
2234 EX color_t ringcolor = 0xFFFF;
2235 EX color_t modelcolor = 0;
2236
2237 #if CAP_QUEUE
draw_model_elements()2238 EX void draw_model_elements() {
2239
2240 #if CAP_VR
2241 if(vrhr::active() && pmodel == mdHyperboloid) return;
2242 #endif
2243
2244 dynamicval<ld> lw(vid.linewidth, vid.linewidth * vid.multiplier_ring);
2245 switch(pmodel) {
2246
2247 case mdRotatedHyperboles: {
2248 queuestr(current_display->xcenter, current_display->ycenter + current_display->radius * pconf.alpha, 0, vid.fsize, "X", ringcolor, 1, 8);
2249 return;
2250 }
2251
2252 case mdTwoHybrid: {
2253 queuereset(mdPixel, PPR::CIRCLE);
2254
2255 for(int mode=0; mode<4; mode++) {
2256 for(int s=-200; s<=200; s ++) {
2257 ld p = tanh(s / 40.);
2258 ld a = pconf.twopoint_param * (1+p);
2259 ld b = pconf.twopoint_param * (1-p);
2260 ld h = ((mode & 2) ? -1 : 1) * sqrt(asin_auto(tan_auto(a) * tan_auto(b)));
2261
2262 hyperpoint H = xpush(p * pconf.twopoint_param) * ypush0(h);
2263
2264 hyperpoint res = compute_hybrid(H, 2 | mode);
2265 models::apply_orientation(res[0], res[1]);
2266 models::apply_orientation_yz(res[2], res[1]);
2267 curvepoint(res * current_display->radius);
2268 }
2269 queuecurve(shiftless(Id), ringcolor, 0, PPR::CIRCLE);
2270 }
2271
2272 queuereset(pmodel, PPR::CIRCLE);
2273 goto fallthrough;
2274 }
2275
2276 case mdTwoPoint: case mdSimulatedPerspective: fallthrough: {
2277 ld a = -pconf.model_orientation * degree;
2278 queuestr(shiftless(xspinpush0(a, +pconf.twopoint_param)), vid.xres / 100, "X", ringcolor >> 8);
2279 queuestr(shiftless(xspinpush0(a, -pconf.twopoint_param)), vid.xres / 100, "X", ringcolor >> 8);
2280 return;
2281 }
2282
2283 case mdThreePoint: {
2284 vid.linewidth *= 5;
2285 for(int i=0; i<=3; i++) {
2286 hyperpoint h = xspinpush0(2*M_PI*i/3, pconf.twopoint_param);
2287 models::apply_orientation(h[1], h[0]);
2288 models::apply_orientation_yz(h[2], h[1]);
2289 curvepoint(h);
2290 }
2291
2292 queuecurve(shiftless(Id), ringcolor, 0, PPR::SUPERLINE);
2293 vid.linewidth /= 5;
2294 return;
2295 }
2296
2297 case mdBall: {
2298 queuecircle(current_display->xcenter, current_display->ycenter, current_display->radius, ringcolor, PPR::OUTCIRCLE, modelcolor);
2299 ballgeometry();
2300 return;
2301 }
2302
2303 case mdHyperboloid: {
2304 if(!pconf.show_hyperboloid_flat) return;
2305 if(hyperbolic) {
2306 #if CAP_QUEUE
2307 curvepoint(point3(0,0,1));
2308 curvepoint(point3(0,0,-pconf.alpha));
2309 queuecurve(shiftless(Id), ringcolor, 0, PPR::CIRCLE);
2310
2311 ld& tz = pconf.top_z;
2312 ld z = acosh(tz);
2313
2314 hyperpoint a = xpush0(z);
2315 ld cb = models::cos_ball;
2316 ld sb = models::sin_ball;
2317
2318 a[1] = sb * a[2] / -cb;
2319 a[0] = sqrt(-1 + a[2] * a[2] - a[1] * a[1]);
2320
2321 curvepoint(point3(0,0,-pconf.alpha));
2322 curvepoint(a);
2323 curvepoint(point3(0,0,0));
2324 a[0] = -a[0];
2325 curvepoint(a);
2326 curvepoint(point3(0,0,-pconf.alpha));
2327 queuecurve(shiftless(Id), ringcolor, 0, PPR::CIRCLE);
2328
2329 curvepoint(point3(-1,0,0));
2330 curvepoint(point3(1,0,0));
2331 queuecurve(shiftless(Id), ringcolor, 0, PPR::CIRCLE);
2332
2333 a[1] = sb * tz / -cb;
2334 a[0] = sqrt(tz * tz - a[1] * a[1]);
2335 a[2] = tz - pconf.alpha;
2336
2337 curvepoint(a);
2338 curvepoint(point3(0,0,-pconf.alpha));
2339 a[0] = -a[0];
2340 curvepoint(a);
2341 queuecurve(shiftless(Id), ringcolor, 0, PPR::CIRCLE);
2342 #endif
2343 }
2344 return;
2345 }
2346
2347 default: break;
2348 }
2349 }
2350
queuestraight(hyperpoint X,int style,color_t lc,color_t fc,PPR p)2351 void queuestraight(hyperpoint X, int style, color_t lc, color_t fc, PPR p) {
2352
2353 hyperpoint H0, H1;
2354 applymodel(shiftless(X), H0);
2355 H0 *= current_display->radius;
2356 ld mul0 = hypot(vid.xres, vid.yres) / hypot_d(2, H0);
2357
2358 if(style == 1) {
2359 H1 = H0 * -mul0;
2360 }
2361 else {
2362 applymodel(shiftless(pispin * X), H1);
2363 H1 *= current_display->radius;
2364 }
2365
2366 ld mul1 = hypot(vid.xres, vid.yres) / hypot_d(2, H1);
2367
2368 queuereset(mdPixel, p);
2369 curvepoint(H0 + spin(M_PI/2) * H0 * mul0);
2370 curvepoint(H0 - spin(M_PI/2) * H0 * mul0);
2371 curvepoint(H1 + spin(M_PI/2) * H1 * mul1);
2372 curvepoint(H1 - spin(M_PI/2) * H1 * mul1);
2373 curvepoint(H0 + spin(M_PI/2) * H0 * mul0);
2374
2375 queuecurve(shiftless(Id), lc, fc, p).flags |= POLY_ALWAYS_IN | POLY_FORCEWIDE;
2376 queuereset(pmodel, p);
2377 /*
2378 for(int i=0; i<1; i++) {
2379 hyperpoint h = spin(i * 45 * degree) * X;
2380 hyperpoint res;
2381 applymodel(h, res);
2382 if(hypot2(res) < 1000 && !std::isnan(res[0]) && !std::isnan(res[1]))
2383 queuestr(h, 16, "X", 0xFF0000 + i * 0x20);
2384 } */
2385 }
2386
draw_boundary(int w)2387 EX void draw_boundary(int w) {
2388
2389 if(w == 1) return;
2390 if(nonisotropic || euclid || prod) return;
2391 #if CAP_VR
2392 if(vrhr::active() && pmodel == mdHyperboloid) return;
2393 #endif
2394
2395 dynamicval<ld> lw(vid.linewidth, vid.linewidth * vid.multiplier_ring);
2396
2397 color_t lc = ringcolor;
2398 color_t fc = modelcolor;
2399 PPR p = PPR::OUTCIRCLE;
2400
2401 if(haveaura()) lc = 0;
2402 if(lc == 0 && fc == 0) return;
2403 if(pmodel == mdRotatedHyperboles) return;
2404
2405 ld fakeinf = sphere ? M_PI-1e-5 : hyperbolic ? 10 : exp(10);
2406
2407 #if CAP_SVG
2408 dynamicval<ld> dw(vid.linewidth, vid.linewidth * (svg::in ? svg::divby : 1));
2409 #endif
2410
2411 if(elliptic && !among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal, mdMollweide, mdCollignon))
2412 circle_around_center(M_PI/2, periodcolor, 0, PPR::CIRCLE);
2413
2414 int broken_coord = models::get_broken_coord(pmodel);
2415 if(broken_coord) {
2416 int unbroken_coord = 3 - broken_coord;
2417 const ld eps = 1e-3;
2418 const ld rem = sqrt(1-eps*eps);
2419 for(int s: {-1, 1}) {
2420 for(int a=1; a<180; a++) {
2421 hyperpoint h = Hypc;
2422 h[broken_coord] = -sin_auto(a*degree) * rem;
2423 h[0] = sin_auto(a*degree) * eps * s;
2424 h[unbroken_coord] = cos_auto(a*degree);
2425 models::apply_orientation(h[1], h[0]);
2426 curvepoint(h);
2427 }
2428 queuecurve(shiftless(Id), periodcolor, 0, PPR::CIRCLE).flags |= POLY_FORCEWIDE;
2429 }
2430 }
2431
2432 if(pmodel == mdWerner && hyperbolic) return;
2433
2434 switch(pmodel) {
2435
2436 case mdTwoPoint: {
2437 if(twopoint_do_flips || current_display->stereo_active() || !sphere) return;
2438 queuereset(mdPixel, p);
2439
2440 for(int b=-1; b<=1; b+=2)
2441 for(ld a=-90; a<=90+1e-6; a+=pow(.5, vid.linequality)) {
2442 ld x = sin(a * pconf.twopoint_param * b / 90);
2443 ld y = 0;
2444 ld z = -sqrt(1 - x*x);
2445 models::apply_orientation(y, x);
2446 hyperpoint h1;
2447 applymodel(shiftless(hpxyz(x,y,z)), h1);
2448
2449 models::apply_orientation(h1[0], h1[1]);
2450 h1[1] = abs(h1[1]) * b;
2451 models::apply_orientation(h1[1], h1[0]);
2452 curvepoint(h1);
2453 }
2454
2455 queuecurve(shiftless(Id), lc, fc, p).flags |= POLY_FORCEWIDE;
2456 queuereset(pmodel, p);
2457 return;
2458 }
2459
2460 case mdBand: case mdBandEquidistant: case mdBandEquiarea: case mdSinusoidal: case mdMollweide: case mdCentralCyl: case mdCollignon:
2461 case mdGallStereographic: case mdMiller:
2462 {
2463 if(GDIM == 3) return;
2464 if(pmodel == mdBand && pconf.model_transition != 1) return;
2465 bool bndband = (among(pmodel, mdBand, mdMiller, mdGallStereographic, mdCentralCyl) ? hyperbolic : sphere);
2466 transmatrix T = spin(-pconf.model_orientation * degree);
2467 ld right = M_PI/2 - 1e-5;
2468 if(bndband)
2469 queuestraight(T * ypush0(hyperbolic ? 10 : right), 2, lc, fc, p);
2470 ld xperiod = elliptic ? fakeinf/2 : fakeinf;
2471 if(sphere && !bndband) {
2472 queuestraight(T * xpush0(xperiod), 2, periodcolor, 0, PPR::CIRCLE);
2473 }
2474 if(sphere && bndband) {
2475 ld adegree = degree-1e-6;
2476 for(ld a=-90; a<90+1e-6; a+=pow(.5, vid.linequality)) {
2477 curvepoint(T * xpush(xperiod) * ypush0(a * adegree));
2478 }
2479 for(ld a=-90; a<90+1e-6; a+=pow(.5, vid.linequality)) {
2480 curvepoint(T * xpush(-xperiod) * ypush0(-a * adegree));
2481 }
2482 curvepoint(T * xpush(xperiod) * ypush0(-90 * adegree));
2483 queuecurve(shiftless(Id), periodcolor, 0, PPR::CIRCLE).flags |= POLY_FORCEWIDE;
2484 }
2485 return;
2486 }
2487
2488 case mdHalfplane:
2489 if(hyperbolic && GDIM == 2) {
2490 queuestraight(xspinpush0(-pconf.model_orientation * degree - M_PI/2, fakeinf), 1, lc, fc, p);
2491 return;
2492 }
2493 break;
2494
2495 case mdHemisphere: {
2496 if(hyperbolic) {
2497 queuereset(mdPixel, p);
2498 for(int i=0; i<=360; i++) {
2499 ld s = sin(i * degree);
2500 curvepoint(point3(current_display->radius * cos(i * degree), current_display->radius * s * (models::cos_ball * s >= 0 - 1e-6 ? 1 : abs(models::sin_ball)), 0));
2501 }
2502 queuecurve(shiftless(Id), lc, fc, p);
2503 queuereset(pmodel, p);
2504 p = PPR::CIRCLE; fc = 0;
2505 queuereset(mdPixel, p);
2506
2507 for(int i=0; i<=360; i++) {
2508 ld s = sin(i * degree);
2509 curvepoint(point3(current_display->radius * cos(i * degree), current_display->radius * s * models::sin_ball, 0));
2510 }
2511 queuecurve(shiftless(Id), lc, fc, p);
2512 queuereset(pmodel, p);
2513 }
2514 if(euclid || sphere) {
2515 queuereset(mdPixel, p);
2516 for(int i=0; i<=360; i++) {
2517 curvepoint(point3(current_display->radius * cos(i * degree), current_display->radius * sin(i * degree), 0));
2518 }
2519 queuecurve(shiftless(Id), lc, fc, p);
2520 queuereset(pmodel, p);
2521 }
2522 return;
2523 }
2524
2525 case mdHyperboloid: {
2526 if(hyperbolic) {
2527 ld& tz = pconf.top_z;
2528 ld mz = acosh(tz);
2529 ld cb = models::cos_ball;
2530 ld sb = models::sin_ball;
2531
2532 if(abs(sb) <= abs(cb) + 1e-5) {
2533 ld step = .01 / (1 << vid.linequality);
2534
2535 hyperpoint a;
2536
2537 for(ld t=-1; t<=1; t += step) {
2538
2539 a = xpush0(t * mz);
2540
2541 if(t != 0) {
2542 a[1] = sb * a[2] / -cb;
2543 ld v = -1 + a[2] * a[2] - a[1] * a[1];
2544 if(v < 0) continue;
2545 a[0] = sqrt(v);
2546 if(t < 0) a[0] = -a[0];
2547 }
2548
2549 curvepoint(a);
2550 }
2551
2552 if((sb > 0) ^ (cb < 0)) {
2553 ld alpha = M_PI - atan2(a[0], -a[1]);
2554
2555 for(ld t=-1; t<=1; t += step)
2556 curvepoint(xspinpush0(-M_PI/2 - t * alpha, mz));
2557 }
2558 else {
2559 ld alpha = - atan2(a[0], -a[1]);
2560
2561 for(ld t=-1; t<=1; t += step)
2562 curvepoint(xspinpush0(+M_PI/2 - t * alpha, mz));
2563 }
2564
2565 queuecurve(shiftless(Id), lc, fc, p);
2566 fc = 0; p = PPR::CIRCLE;
2567 }
2568
2569 for(ld t=0; t<=360; t ++)
2570 curvepoint(xspinpush0(t * degree, mz));
2571
2572 queuecurve(shiftless(Id), lc, fc, p);
2573 }
2574 return;
2575 }
2576
2577 case mdSpiral: {
2578 if(euclid) return;
2579 if(models::ring_not_spiral) return;
2580 // if(p == PPR::CIRCLE) p = PPR::OUTCIRCLE;
2581 auto& sm = models::spiral_multiplier;
2582 ld u = hypot(1, imag(sm) / real(sm));
2583 if(real(sm)) {
2584 queuereset(mdPixel, p);
2585 for(ld a=-10; a<=10; a+=0.01 / (1 << vid.linequality) / u) {
2586 cld z = exp(cld(a, a * imag(sm) / real(sm) + M_PI));
2587 hyperpoint ret = point2(real(z), imag(z));
2588 ret = mobius(ret, pconf.skiprope, 1);
2589 ret *= current_display->radius;
2590 curvepoint(ret);
2591 }
2592 queuecurve(shiftless(Id), ringcolor, 0, p).flags |= POLY_ALWAYS_IN;
2593 queuereset(pmodel, p);
2594 }
2595 return;
2596 }
2597
2598 default: break;
2599 }
2600
2601 if(sphere && pmodel == mdDisk && pconf.alpha > 1) {
2602 double rad = current_display->radius / sqrt(pconf.alpha*pconf.alpha - 1);
2603 queuecircle(current_display->xcenter, current_display->ycenter, rad, lc, p, fc);
2604 return;
2605 }
2606
2607 if(sphere && !among(pmodel, mdEquidistant, mdEquiarea)) return;
2608 circle_around_center(fakeinf, lc, fc, p);
2609 }
2610 #endif
2611
change_shift(shiftpoint & h,ld by)2612 EX void change_shift(shiftpoint& h, ld by) {
2613 if(!by) return;
2614 h.shift += by;
2615 if((mdinf[pmodel].flags & mf::uses_bandshift) || (sphere && pmodel == mdSpiral)) {
2616 h.h = spin(pconf.model_orientation * degree) * h.h;
2617 h.h = xpush(-by) * h.h;
2618 h.h = spin(-pconf.model_orientation * degree) * h.h;
2619 }
2620 if(sl2) {
2621 ld ca = cos(by), sa = sin(by);
2622 tie(h[2], h[3]) = make_pair(h[2] * ca - h[3] * sa, h[3] * ca + h[2] * sa);
2623 tie(h[0], h[1]) = make_pair(h[0] * ca - h[1] * sa, h[1] * ca + h[0] * sa);
2624 }
2625 }
2626
change_shift(shiftmatrix & T,ld by)2627 EX void change_shift(shiftmatrix& T, ld by) {
2628 if(!by) return;
2629 T.shift += by;
2630 if((mdinf[pmodel].flags & mf::uses_bandshift) || (sphere && pmodel == mdSpiral)) {
2631 T.T = spin(pconf.model_orientation * degree) * T.T;
2632 T.T = xpush(-by) * T.T;
2633 fixmatrix(T.T);
2634 T.T = spin(-pconf.model_orientation * degree) * T.T;
2635 }
2636 if(sl2) {
2637 ld ca = cos(by), sa = sin(by);
2638 for(int a=0; a<4; a++) {
2639 tie(T[2][a], T[3][a]) = make_pair(T[2][a] * ca - T[3][a] * sa, T[3][a] * ca + T[2][a] * sa);
2640 tie(T[0][a], T[1][a]) = make_pair(T[0][a] * ca - T[1][a] * sa, T[1][a] * ca + T[0][a] * sa);
2641 }
2642 }
2643 }
2644
2645 EX transmatrix unshift(shiftmatrix T, ld to IS(0)) {
2646 change_shift(T, to - T.shift);
2647 return T.T;
2648 }
2649
2650 EX hyperpoint unshift(shiftpoint T, ld to IS(0)) {
2651 change_shift(T, to - T.shift);
2652 return T.h;
2653 }
2654
inverse_shift(const shiftmatrix & T1,const shiftmatrix & T2)2655 EX transmatrix inverse_shift(const shiftmatrix& T1, const shiftmatrix& T2) {
2656 return iso_inverse(T1.T) * unshift(T2, T1.shift);
2657 }
2658
inverse_shift(const shiftmatrix & T1,const shiftpoint & T2)2659 EX hyperpoint inverse_shift(const shiftmatrix& T1, const shiftpoint& T2) {
2660 return iso_inverse(T1.T) * unshift(T2, T1.shift);
2661 }
2662
optimize_shift(shiftmatrix & T)2663 EX void optimize_shift(shiftmatrix& T) {
2664 if(((mdinf[pmodel].flags & mf::uses_bandshift) && T[LDIM][LDIM] > 1e6) || (sphere && pmodel == mdSpiral)) {
2665 T.T = spin(pconf.model_orientation * degree) * T.T;
2666 hyperpoint H = tC0(T.T);
2667 find_zlev(H);
2668
2669 ld y = asin_auto(H[1]);
2670 ld x = asin_auto_clamp(H[0] / cos_auto(y));
2671 if(sphere) {
2672 if(H[LDIM] < 0 && x > 0) x = M_PI - x;
2673 else if(H[LDIM] < 0 && x <= 0) x = -M_PI - x;
2674 }
2675 T.shift += x;
2676 T.T = xpush(-x) * T.T;
2677 fixmatrix(T.T);
2678 T.T = spin(-pconf.model_orientation * degree) * T.T;
2679 }
2680
2681 if(sl2) {
2682 change_shift(T, atan2(T[2][3], T[3][3]));
2683 if(hybrid::csteps) {
2684 auto period = (M_PI * hybrid::csteps) / cgi.psl_steps;
2685 while(T.shift > period*.4999)
2686 T.shift -= period;
2687 while(T.shift < -period*.5001)
2688 T.shift += period;
2689 }
2690 }
2691 }
2692
optimized_shift(const shiftmatrix & T)2693 EX shiftmatrix optimized_shift(const shiftmatrix& T) {
2694 shiftmatrix U = T;
2695 optimize_shift(U);
2696 return U;
2697 }
2698
2699 EX namespace dq {
2700 EX queue<pair<heptagon*, shiftmatrix>> drawqueue;
2701
bucketer(const shiftpoint & T)2702 EX unsigned bucketer(const shiftpoint& T) {
2703 return bucketer(T.h) + unsigned(floor(T.shift*81527+.5));
2704 }
2705
2706 EX set<heptagon*> visited;
enqueue(heptagon * h,const shiftmatrix & T)2707 EX void enqueue(heptagon *h, const shiftmatrix& T) {
2708 if(!h || visited.count(h)) { return; }
2709 visited.insert(h);
2710 drawqueue.emplace(h, T);
2711 }
2712
2713 EX set<unsigned> visited_by_matrix;
enqueue_by_matrix(heptagon * h,const shiftmatrix & T)2714 EX void enqueue_by_matrix(heptagon *h, const shiftmatrix& T) {
2715 if(!h) return;
2716 unsigned b = bucketer(tC0(T));
2717 if(visited_by_matrix.count(b)) { return; }
2718 visited_by_matrix.insert(b);
2719 drawqueue.emplace(h, T);
2720 }
2721
2722 EX queue<pair<cell*, shiftmatrix>> drawqueue_c;
2723 EX set<cell*> visited_c;
2724
enqueue_c(cell * c,const shiftmatrix & T)2725 EX void enqueue_c(cell *c, const shiftmatrix& T) {
2726 if(!c || visited_c.count(c)) { return; }
2727 visited_c.insert(c);
2728 drawqueue_c.emplace(c, T);
2729 }
2730
enqueue_by_matrix_c(cell * c,const shiftmatrix & T)2731 EX void enqueue_by_matrix_c(cell *c, const shiftmatrix& T) {
2732 if(!c) return;
2733 unsigned b = bucketer(tC0(T));
2734 if(visited_by_matrix.count(b)) { return; }
2735 visited_by_matrix.insert(b);
2736 drawqueue_c.emplace(c, T);
2737 }
2738
clear_all()2739 EX void clear_all() {
2740 visited.clear();
2741 visited_by_matrix.clear();
2742 visited_c.clear();
2743 while(!drawqueue_c.empty()) drawqueue_c.pop();
2744 while(!drawqueue.empty()) drawqueue.pop();
2745 }
2746
2747
2748 EX }
2749
do_draw(cell * c)2750 EX bool do_draw(cell *c) {
2751 // do not display out of range cells, unless on torus
2752 if(c->pathdist == PINFD && !(euclid && quotient) && vid.use_smart_range == 0)
2753 return false;
2754 // do not display not fully generated cells, unless changing range allowed
2755 if(c->mpdist > 7 && !allowChangeRange()) return false;
2756 // in the Yendor Challenge, scrolling back is forbidden
2757 if(c->cpdist > 7 && (yendor::on || isHaunted(cwt.at->land)) && !cheater && !autocheat) return false;
2758
2759 return true;
2760 }
2761
2762 EX ld extra_generation_distance = 99;
2763
2764 // returns false if limited
limited_generation(cell * c)2765 bool limited_generation(cell *c) {
2766 if(c->mpdist <= 7) return true;
2767 if(cells_generated > vid.cells_generated_limit) return false;
2768 setdist(c, 7, c);
2769 cells_generated++;
2770 return true;
2771 }
2772
do_draw(cell * c,const shiftmatrix & T)2773 EX bool do_draw(cell *c, const shiftmatrix& T) {
2774
2775 if(WDIM == 3) {
2776 // do not care about cells outside of the track
2777 if(GDIM == 3 && racing::on && c->land == laMemory && cells_drawn >= S7+1) return false;
2778
2779 if(cells_drawn > vid.cells_drawn_limit) return false;
2780 if(cells_drawn < 50) { limited_generation(c); return true; }
2781 #if MAXMDIM >= 4
2782 if(nil && pmodel == mdGeodesic) {
2783 ld dist = hypot_d(3, inverse_exp(tC0(T), pQUICK));
2784 if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : 0.9)) return false;
2785 if(dist <= extra_generation_distance && !limited_generation(c)) return false;
2786 }
2787 else if(pmodel == mdGeodesic && sol) {
2788 if(!nisot::in_table_range(tC0(T.T))) return false;
2789 if(!limited_generation(c)) return false;
2790 }
2791 else if(pmodel == mdGeodesic && nih) {
2792 hyperpoint h = inverse_exp(tC0(T), pQUICK);
2793 ld dist = hypot_d(3, h);
2794 if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : cgi.corner_bonus)) return false;
2795 if(dist <= extra_generation_distance && !limited_generation(c)) return false;
2796 }
2797 else if(pmodel == mdGeodesic && sl2) {
2798 if(hypot(tC0(T.T)[2], tC0(T.T)[3]) > cosh(slr::range_xy)) return false;
2799 if(abs(T.shift * stretch::not_squared()) > sightranges[geometry]) return false;
2800 if(!limited_generation(c)) return false;
2801 }
2802 #endif
2803 else if(vid.use_smart_range) {
2804 if(cells_drawn >= 50 && !in_smart_range(T)) return false;
2805 if(!limited_generation(c)) return false;
2806 }
2807 else {
2808 ld dist = hdist0(tC0(T.T));
2809 if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : cgi.corner_bonus)) return false;
2810 if(dist <= extra_generation_distance && !limited_generation(c)) return false;
2811 }
2812 return true;
2813 }
2814
2815 #if MAXMDIM >= 4
2816 if(rots::drawing_underlying && euclid && hdist0(tC0(T)) > 6) return false;
2817 #endif
2818 if(just_gmatrix && sphere) return true;
2819 if(!do_draw(c)) return false;
2820 if(euclid && pmodel == mdSpiral) {
2821 hyperpoint h = tC0(T.T);
2822 cld z(h[0], h[1]);
2823 z = z * models::spiral_multiplier;
2824 ld iz = imag(z) + 1.14279e-2; // make it never fall exactly on PI
2825 if(iz < -M_PI || iz >= M_PI) return false;
2826 }
2827 if(pmodel == mdSpiral && models::ring_not_spiral) {
2828 cld z;
2829 shiftpoint H = tC0(T);
2830 hyperpoint ret;
2831 makeband(H, ret, band_conformal);
2832 z = cld(ret[0], ret[1]) * models::spiral_multiplier;
2833 if(imag(z) < -models::spiral_cone_rad/2-1e-5 || imag(z) >= models::spiral_cone_rad/2-1e-5) return false;
2834 }
2835 if(cells_drawn > vid.cells_drawn_limit) return false;
2836 bool usr = vid.use_smart_range || quotient;
2837 if(usr && cells_drawn >= 50 && !in_smart_range(T) && !(WDIM == 2 && GDIM == 3 && hdist0(tC0(T)) < 2.5)) return false;
2838 if(vid.use_smart_range == 2 && !limited_generation(c)) return false;
2839 return true;
2840 }
2841
cone_side(const shiftpoint H)2842 EX int cone_side(const shiftpoint H) {
2843 hyperpoint ret;
2844 if(hyperbolic) makeband(H, ret, band_conformal);
2845 else ret = unshift(H);
2846 cld z = cld(ret[0], ret[1]) * models::spiral_multiplier;
2847
2848 auto zth = [&] (cld z) {
2849 ld alpha = imag(z) * 360 / pconf.spiral_cone;
2850 ld r = real(z);
2851 r = exp(r);
2852
2853 hyperpoint ret;
2854
2855 ret[0] = -sin(alpha) * r;
2856 ret[1] = cos(alpha) * r;
2857 ret[2] = (r-1) * sqrt( pow(360/pconf.spiral_cone, 2) - 1);
2858
2859 models::apply_ball(ret[2], ret[1]);
2860 return ret;
2861 };
2862
2863 hyperpoint ret0 = zth(z);
2864 hyperpoint ret1 = zth(z + cld(1e-3, 0));
2865 hyperpoint ret2 = zth(z + cld(0, 1e-3));
2866
2867 return (ret1[1] - ret0[1]) * (ret2[0] - ret0[0]) < (ret2[1] - ret0[1]) * (ret1[0] - ret0[0]) ? 1 : -1;
2868 }
2869
2870 /** get the current orientation of the view */
get_view_orientation()2871 EX transmatrix& get_view_orientation() {
2872 return prod ? NLP : View;
2873 }
2874
2875 EX hookset<bool(const transmatrix&)> hooks_rotate_view;
2876 EX hookset<bool(const hyperpoint&)> hooks_shift_view;
2877
2878 /** rotate the view using the given rotation matrix */
rotate_view(transmatrix T)2879 EX void rotate_view(transmatrix T) {
2880 if(callhandlers(false, hooks_rotate_view, T)) return;
2881 transmatrix& which = get_view_orientation();
2882 which = T * which;
2883 if(!prod && !nonisotropic && !rug::rugged) current_display->which_copy = T * current_display->which_copy;
2884 }
2885
2886 /** shift the view according to the given tangent vector */
get_shift_view_of(const hyperpoint H,const transmatrix V)2887 EX transmatrix get_shift_view_of(const hyperpoint H, const transmatrix V) {
2888 if(!nonisotropic && !stretch::in()) {
2889 return rgpushxto0(direct_exp(lp_iapply(H))) * V;
2890 }
2891 else if(!nisot::geodesic_movement) {
2892 transmatrix IV = view_inverse(V);
2893 nisot::fixmatrix(IV);
2894 const transmatrix V1 = iview_inverse(IV);
2895 return V1 * eupush(IV * eupush(H) * V1 * C0);
2896 }
2897 else {
2898 return iview_inverse(nisot::parallel_transport(view_inverse(V), -H));
2899 }
2900 }
2901
2902 /** shift the view according to the given tangent vector */
shift_view(hyperpoint H)2903 EX void shift_view(hyperpoint H) {
2904 if(callhandlers(false, hooks_shift_view, H)) return;
2905 auto oView = View;
2906 View = get_shift_view_of(H, View);
2907 auto& wc = current_display->which_copy;
2908 if(nonisotropic || stretch::in()) {
2909 transmatrix ioldv = eupush(tC0(view_inverse(oView)));
2910 transmatrix newv = inverse(eupush(tC0(view_inverse(View))));
2911 wc = newv * ioldv * wc;
2912 }
2913 else
2914 wc = get_shift_view_of(H, wc);
2915 }
2916
multiply_view(transmatrix T)2917 void multiply_view(transmatrix T) {
2918 View = T * View;
2919 auto& wc = current_display->which_copy;
2920 wc = T * wc;
2921 }
2922
shift_view_to(shiftpoint H)2923 EX void shift_view_to(shiftpoint H) {
2924 if(!nonisotropic) multiply_view(gpushxto0(unshift(H)));
2925 else shift_view(-inverse_exp(H));
2926 }
2927
shift_view_towards(shiftpoint H,ld l)2928 EX void shift_view_towards(shiftpoint H, ld l) {
2929 if(!nonisotropic && !prod)
2930 multiply_view(rspintox(unshift(H)) * xpush(-l) * spintox(unshift(H)));
2931 else if(nonisotropic && !nisot::geodesic_movement)
2932 shift_view(tangent_length(unshift(H)-C0, -l));
2933 else {
2934 hyperpoint ie = inverse_exp(H, pNORMAL | pfNO_DISTANCE);
2935 if(prod) ie = lp_apply(ie);
2936 shift_view(tangent_length(ie, -l));
2937 }
2938 }
2939
set_view(hyperpoint camera,hyperpoint forward,hyperpoint upward)2940 EX void set_view(hyperpoint camera, hyperpoint forward, hyperpoint upward) {
2941 if(WDIM == 2) {
2942 View = gpushxto0(camera);
2943 View = spin(90*degree) * spintox(View * upward) * View;
2944 return;
2945 }
2946
2947 transmatrix V = gpushxto0(camera);
2948 forward = V * forward;
2949 upward = V * upward;
2950
2951 if(pmodel == mdGeodesic || hyperbolic || sphere) {
2952 forward = inverse_exp(shiftless(forward));
2953 }
2954 else {
2955 // apply_nil_rotation(forward);
2956 forward -= C0;
2957 }
2958
2959 if(hypot_d(3, forward) == 0) forward[0] = 1;
2960
2961 forward /= hypot_d(3, forward);
2962
2963 if(pmodel == mdGeodesic || hyperbolic || sphere)
2964 upward = inverse_exp(shiftless(upward));
2965 else {
2966 // apply_nil_rotation(upward);
2967 upward -= C0;
2968 }
2969
2970 upward -= (forward|upward) * forward;
2971 if(hypot_d(3, upward) == 0) return;
2972
2973 upward /= hypot_d(3, upward);
2974
2975 hyperpoint rightward = (forward ^ upward);
2976
2977 transmatrix rotator = Id;
2978 rotator[2] = forward;
2979 rotator[1] = -upward;
2980 rotator[0] = -rightward;
2981
2982 if(det(rotator) < 0) rotator[0] = -rotator[0];
2983
2984 View = iso_inverse(rgpushxto0(camera));
2985 if(prod)
2986 NLP = rotator;
2987 else
2988 View = rotator * View;
2989 }
2990
2991 }
2992