improvements in the two-point model on the sphere

This commit is contained in:
Zeno Rogue 2018-04-21 12:18:33 +02:00
parent abbf536466
commit 323893094f
4 changed files with 104 additions and 39 deletions

View File

@ -5026,11 +5026,6 @@ void queuecircleat(cell *c, double rad, int col) {
void drawMarkers() {
if(pmodel == mdTwoPoint) {
queuechr(xpush(+vid.twopoint_param) * C0, 8, 'X', 0xFFFF00);
queuechr(xpush(-vid.twopoint_param) * C0, 8, 'X', 0xFFFF00);
}
if(!(cmode & sm::NORMAL)) return;
for(cell *c1: crush_now)
@ -5500,18 +5495,24 @@ void drawfullmap() {
ptds.clear();
if(!stereo::active() && sphere && pmodel == mdTwoPoint) {
if(pmodel == mdTwoPoint) {
queuechr(xpush(+vid.twopoint_param) * C0, vid.xres / 100, 'X', 0xFF0000);
queuechr(xpush(-vid.twopoint_param) * C0, vid.xres / 100, 'X', 0xFF0000);
}
if(!twopoint_do_flips && !stereo::active() && sphere && pmodel == mdTwoPoint) {
queuereset(vid.usingGL ? mdDisk : mdUnchanged, PPR_CIRCLE);
for(int a=0; a<=180; a++) {
for(int b=-1; b<=1; b+=2)
for(int a=-90; a<=90; a++) {
using namespace hyperpoint_vec;
ld x = cos(a * M_PI / 90);
ld x = sin(a * vid.twopoint_param * b / 90);
ld y = 0;
ld z = -sqrt(1 - x*x);
hyperpoint h1;
applymodel(hpxyz(x,y,z), h1);
if(h1[1] < 0) h1[1] = -h1[1];
if(a >= 90) h1[1] = -h1[1];
h1[1] = abs(h1[1]) * b;
curvepoint(h1 * vid.radius);
}

View File

@ -66,6 +66,15 @@ ld asin_auto(ld x) {
}
}
ld asin_auto_clamp(ld x) {
switch(cgclass) {
case gcEuclid: return x;
case gcHyperbolic: return asinh(x);
case gcSphere: return x>1 ? M_PI/2 : x<-1 ? -M_PI/2 : isnan(x) ? 0 : asin(x);
default: return x;
}
}
ld cos_auto(ld x) {
switch(cgclass) {
case gcEuclid: return 1;

View File

@ -123,6 +123,8 @@ bool hypot_zlev(bool zlev_used, ld& d, ld zlev, ld& df, ld& zf) {
return hypot_zlev(zlev_used, d, zlev, df, zf, z);
}
int twopoint_sphere_flips;
bool twopoint_do_flips;
void applymodel(hyperpoint H, hyperpoint& ret) {
@ -266,7 +268,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
else {
ld x, y, yf, zf;
y = asin_auto(H[1]);
x = asin_auto(H[0] / cos_auto(y));
x = asin_auto_clamp(H[0] / cos_auto(y));
if(sphere) {
if(H[2] < 0 && x > 0) x = M_PI - x;
else if(H[2] < 0 && x <= 0) x = -M_PI - x;
@ -276,6 +278,19 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
auto p = vid.twopoint_param;
ld dleft = hypot_auto(x-p, y);
ld dright = hypot_auto(x+p, y);
if(sphere) {
int tss = twopoint_sphere_flips;
if(tss&1) { tss--;
dleft = 2*M_PI - 2*p - dleft;
dright = 2*M_PI - 2*p - dright;
swap(dleft, dright);
y = -y;
}
while(tss) { tss -= 2;
dleft = 2*M_PI - 4*p + dleft;
dright = 2*M_PI - 4*p + dright;
}
}
x = (dright*dright-dleft*dleft) / 4 / p;
y = (y>0?1:-1) * sqrt(dleft * dleft - (x-p)*(x-p) + 1e-9);
}

View File

@ -469,6 +469,10 @@ unsigned char& part(int& col, int i) {
bool in_twopoint = false;
ld glhypot2(glvertex a, glvertex b) {
return (a[0]-b[0]) * (a[0]-b[0]) + (a[1]-b[1]) * (a[1]-b[1]) + (a[2]-b[2]) * (a[2]-b[2]);
}
void drawpolyline(polytodraw& p) {
auto& pp = p.u.poly;
@ -485,43 +489,79 @@ void drawpolyline(polytodraw& p) {
}
if(sphere && pmodel == mdTwoPoint && !in_twopoint) {
static vector<glvertex> v0, v1;
v0.clear(); v1.clear();
for(int i=0; i<pp.cnt; i++) {
glvertex h = glhr::pointtogl(pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]));
if(h[2] >= 0)
v0.push_back(h), v1.push_back(h);
else if(h[1] < 0)
v0.push_back(h);
else if(h[1] > 0)
v1.push_back(h);
#define MAX_PHASE 4
vector<glvertex> phases[MAX_PHASE];
extern int twopoint_sphere_flips;
extern bool twopoint_do_flips;
int pha;
if(twopoint_do_flips) {
for(int i=0; i<pp.cnt; i++) {
hyperpoint h1 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]);
for(int j=0; j<MAX_PHASE; j++) {
twopoint_sphere_flips = j;
hyperpoint h2; applymodel(h1, h2);
using namespace hyperpoint_vec;
glvertex h = glhr::pointtogl(h2 * vid.radius);
if(i == 0)
phases[j].push_back(h);
else {
int best = -1;
ld bhypot = 1e60;
for(int j0=0; j0<MAX_PHASE; j0++)
if(size(phases[j0]) == i) {
ld chypot = glhypot2(phases[j0].back(), h);
if(chypot < bhypot || best == -1) bhypot = chypot, best = j0;
}
phases[best].push_back(h);
}
}
}
twopoint_sphere_flips = 0;
pha = MAX_PHASE-1;
}
else {
pha = 1;
if(true) {
// a
// b
// lin(a,b) is of form (x, 0, z)
int cpha = 0;
for(int i=0; i<pp.cnt; i++) {
using namespace hyperpoint_vec;
hyperpoint h1 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]);
hyperpoint mh1; applymodel(h1, mh1);
phases[cpha].push_back(glhr::pointtogl(mh1 * vid.radius));
// check if the i-th edge intersects the boundary of the ellipse
// (which corresponds to the segment between the antipodes of foci)
// if yes, switch cpha to the opposite
hyperpoint h2 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+(i+1)%pp.cnt]);
if(h1[1] * h2[1] > 0) continue;
ld c1 = h1[1], c2 = -h2[1];
if(c1 < 0) c1 = -c1, c2 = -c2;
hyperpoint h = h1 * c1 + h2 * c2;
h /= hypot3(h);
if(h[2] < 0 && abs(h[0]) < sin(vid.twopoint_param)) cpha = 1-cpha, pha = 2;
}
if(cpha == 1) pha = 0;
}
}
auto tV = pp.V;
auto toffset = pp.offset;
auto ttab = pp.tab;
auto tcnt = pp.cnt;
vector<glvertex> tv;
if(pp.tinf) {
for(int i=0; i<pp.cnt; i++)
tv.push_back(pp.tinf->tvertices[pp.offset+i]);
swap(pp.tinf->tvertices, tv);
}
pp.V = Id; pp.offset = 0;
in_twopoint = true;
if(size(v0) == pp.cnt) {
pp.tab = &v0;
drawpolyline(p);
}
else if(size(v1) == pp.cnt) {
pp.tab = &v1;
dynamicval<eModel> d1(pmodel, mdUnchanged);
dynamicval<transmatrix> d2(pp.V, Id);
dynamicval<int> d3(pp.offset, 0);
dynamicval<decltype(pp.tab)> d4(pp.tab, pp.tab);
for(int j=0; j<pha; j++) {
dynamicval<int> d5(pp.cnt, size(phases[j]));
pp.tab = &phases[j];
drawpolyline(p);
}
else {
if(size(v0) >= 3) { pp.tab = &v0; pp.cnt = size(v0); drawpolyline(p); }
if(size(v1) >= 3) { pp.tab = &v1; pp.cnt = size(v1); drawpolyline(p); }
}
in_twopoint = false;
pp.V = tV; pp.offset = toffset; pp.tab = ttab; pp.cnt = tcnt;
if(pp.tinf) swap(pp.tinf->tvertices, tv);
return;
}