From f6fa051eabbe6e2bb517ba0031bbaaa8b389a1eb Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Wed, 27 Dec 2017 06:31:47 +0100 Subject: [PATCH] torus movement --- hyperpoint.cpp | 36 ++++++++----- pattern2.cpp | 1 + rug.cpp | 144 ++++++++++++++++++++++++++++--------------------- 3 files changed, 106 insertions(+), 75 deletions(-) diff --git a/hyperpoint.cpp b/hyperpoint.cpp index ac093b29..888e4968 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -86,6 +86,10 @@ ld intvalxyz(const hyperpoint &h1, const hyperpoint &h2) { return squar(h1[0]-h2[0]) + squar(h1[1]-h2[1]) + squar(h1[2]-h2[2]); } +ld hypot3(const hyperpoint& h) { + return sqrt(h[0]*h[0]+h[1]*h[1]+h[2]*h[2]); + } + ld zlevel(const hyperpoint &h) { if(euclid) return h[2]; else if(sphere) return sqrt(intval(h, Hypc)); @@ -394,26 +398,32 @@ double hdist(const hyperpoint& h1, const hyperpoint& h2) { namespace hyperpoint_vec { - hyperpoint operator * (ld d, hyperpoint h) { - return hpxyz(h[0]*d, h[1]*d, h[2]*d); + hyperpoint& operator *= (hyperpoint& h, ld d) { + h[0] *= d; h[1] *= d; h[2] *= d; + return h; } - hyperpoint operator * (hyperpoint h, ld d) { - return hpxyz(h[0]*d, h[1]*d, h[2]*d); + hyperpoint& operator /= (hyperpoint& h, ld d) { + h[0] /= d; h[1] /= d; h[2] /= d; + return h; } - hyperpoint operator / (hyperpoint h, ld d) { - return hpxyz(h[0]/d, h[1]/d, h[2]/d); + hyperpoint operator += (hyperpoint& h, hyperpoint h2) { + for(int i: {0,1,2}) h[i] += h2[i]; + return h; } - - hyperpoint operator + (hyperpoint h, hyperpoint h2) { - return hpxyz(h[0]+h2[0], h[1]+h2[1], h[2]+h2[2]); - } - - hyperpoint operator - (hyperpoint h, hyperpoint h2) { - return hpxyz(h[0]-h2[0], h[1]-h2[1], h[2]-h2[2]); + + hyperpoint operator -= (hyperpoint& h, hyperpoint h2) { + for(int i: {0,1,2}) h[i] -= h2[i]; + return h; } + hyperpoint operator * (ld d, hyperpoint h) { return h *= d; } + hyperpoint operator * (hyperpoint h, ld d) { return h *= d; } + hyperpoint operator / (hyperpoint h, ld d) { return h /= d; } + hyperpoint operator + (hyperpoint h, hyperpoint h2) { return h += h2; } + hyperpoint operator - (hyperpoint h, hyperpoint h2) { return h -= h2; } + // cross product hyperpoint operator ^ (hyperpoint h1, hyperpoint h2) { return hpxyz( diff --git a/pattern2.cpp b/pattern2.cpp index a0c56b4a..a60bec46 100644 --- a/pattern2.cpp +++ b/pattern2.cpp @@ -19,6 +19,7 @@ int eupattern(cell *c) { } int eupattern4(cell *c) { + if(torus) return 0; eucoord x, y; decodeMaster(c->master, x, y); return (x&1) + ((y&1)) * 2; diff --git a/rug.cpp b/rug.cpp index e02e8214..65da64ca 100644 --- a/rug.cpp +++ b/rug.cpp @@ -33,7 +33,8 @@ GLAPI void APIENTRY glDeleteFramebuffers (GLsizei n, const GLuint *framebuffers) namespace rug { -static const bool fast_euclidean = true; +bool fast_euclidean = true; +bool keep_torus = false; int torus_precision = 2; @@ -51,7 +52,7 @@ bool rendernogl = false; int texturesize = 1024; ld scale = 1; -ld err_zero = 1e-3; +ld err_zero = 1e-3, err_zero_current; int queueiter, qvalid, dt; double err; @@ -69,6 +70,19 @@ struct rugpoint { hyperpoint h; hyperpoint flat; vector edges; + // Find-Union algorithm + rugpoint *glue = NULL; + rugpoint *getglue() { + return glue ? (glue = glue->getglue()) : this; + } + hyperpoint& glueflat() { + return glue->flat; + } + void glueto(rugpoint *x) { + x = x->getglue(); + auto y = getglue(); + if(x != y) y->glue = x; + } }; struct triangle { @@ -206,32 +220,40 @@ rugpoint *addRugpoint(hyperpoint h, double dist) { return m; } -rugpoint *findRugpoint(hyperpoint h, double dist) { +rugpoint *findRugpoint(hyperpoint h) { for(int i=0; ih, h) < 1e-5) return points[i]; - return addRugpoint(h, dist); + return NULL; } -void addNewEdge(rugpoint *e1, rugpoint *e2) { +rugpoint *findOrAddRugpoint(hyperpoint h, double dist) { + rugpoint *r = findRugpoint(h); + return r ? r : addRugpoint(h, dist); + } + +void addNewEdge(rugpoint *e1, rugpoint *e2, ld len = 1) { edge e; e.target = e2; e1->edges.push_back(e); e.target = e1; e2->edges.push_back(e); + e.len = len; } -void addEdge(rugpoint *e1, rugpoint *e2) { +void addEdge(rugpoint *e1, rugpoint *e2, ld len = 1) { for(int i=0; iedges); i++) if(e1->edges[i].target == e2) return; - addNewEdge(e1, e2); + addNewEdge(e1, e2, len); } -void addTriangle(rugpoint *t1, rugpoint *t2, rugpoint *t3) { - addEdge(t1, t2); addEdge(t2, t3); addEdge(t3, t1); +void addTriangle(rugpoint *t1, rugpoint *t2, rugpoint *t3, ld len = 1) { + addEdge(t1->getglue(), t2->getglue(), len); + addEdge(t2->getglue(), t3->getglue(), len); + addEdge(t3->getglue(), t1->getglue(), len); triangles.push_back(triangle(t1,t2,t3)); } void addTriangle1(rugpoint *t1, rugpoint *t2, rugpoint *t3) { - rugpoint *t12 = findRugpoint(mid(t1->h, t2->h), (t1->dist+ t2->dist)/2); - rugpoint *t23 = findRugpoint(mid(t2->h, t3->h), (t1->dist+ t2->dist)/2); - rugpoint *t31 = findRugpoint(mid(t3->h, t1->h), (t1->dist+ t2->dist)/2); + rugpoint *t12 = findOrAddRugpoint(mid(t1->h, t2->h), (t1->dist+ t2->dist)/2); + rugpoint *t23 = findOrAddRugpoint(mid(t2->h, t3->h), (t1->dist+ t2->dist)/2); + rugpoint *t31 = findOrAddRugpoint(mid(t3->h, t1->h), (t1->dist+ t2->dist)/2); addTriangle(t1, t12, t31); addTriangle(t12, t2, t23); addTriangle(t23, t3, t31); @@ -336,6 +358,8 @@ void buildTorusRug() { double sc = (factor+1)/4; r->flat = hpxyz((factor+cos(alpha)) * cos(beta) * sc, (factor+cos(alpha)) * sin(beta) * sc, -sin(alpha) * sc); r->valid = true; + rugpoint *r2 = findRugpoint(r->flat); + if(r2 && r2 != r) r->glueto(r2); return r; }; @@ -351,8 +375,8 @@ void buildTorusRug() { for(int yy=0; yyx1 = (vid.xcenter + vid.radius * vid.scale * p->x1)/ vid.xres, p->y1 = (vid.ycenter - vid.radius * vid.scale * p->y1)/ vid.yres; + qvalid = size(points); + printf("qvalid = %d\n", qvalid); + return; } @@ -425,7 +452,7 @@ void force_euclidean(rugpoint& m1, rugpoint& m2, double rd, double d1=1, double t = sqrt(t); // printf("%lf %lf\n", t, rd); err += (t-rd) * (t-rd); - bool nonzero = abs(t-rd) > err_zero; + bool nonzero = abs(t-rd) > err_zero_current; double force = (t - rd) / t / 2; // 20.0; for(int i=0; i<3; i++) { double di = (m2.flat[i] - m1.flat[i]) * force; @@ -452,7 +479,7 @@ void force(rugpoint& m1, rugpoint& m2, double rd, double d1=1, double d2=1) { ld t = hdist(f1, f2); err += (t-rd) * (t-rd); - bool nonzero = abs(t-rd) > err_zero; + bool nonzero = abs(t-rd) > err_zero_current; double forcev = (t - rd) / 2; // 20.0; transmatrix T = gpushxto0(f1); @@ -543,7 +570,11 @@ bool stop = false; void subdivide() { int N = size(points); // if(euclid && gwhere == gEuclid) return; - if(N >= 8000) {stop = true; return; } + if(N >= 8000 || torus) { + err_zero /= 2; + for(auto p: points) enqueue(p); + return; + } printf("subdivide (%d,%d)\n", N, size(triangles)); divides++; vector otriangles = triangles; @@ -602,7 +633,7 @@ void addNewPoints() { void physics() { - if(torus) return; + if(keep_torus && torus) return; int t = SDL_GetTicks(); @@ -617,7 +648,7 @@ void physics() { pqueue.pop(); m->inqueue = false; for(int j=0; jedges); j++) - force(*m, *m->edges[j].target, m->edges[j].len); + force(*m, *m->edges[j].target, m->edges[j].len); } if(!stop) printf("%5d %10.7lf D%d Q%3d Qv%5d\n", queueiter, err, divides, size(pqueue), qvalid); @@ -628,58 +659,43 @@ void physics() { int eyemod; -void getco(rugpoint& m, double& x, double& y, double& z, int &spherepoints) { - x = m.flat[0]; - y = m.flat[1]; - z = m.flat[2]; - if(gwhere == gSphere && z > 0) { - ld rad = sqrt(x*x + y*y + z*z); +void getco(rugpoint *m, hyperpoint& h, int &spherepoints) { + using namespace hyperpoint_vec; + h = m->getglue()->flat; + if(gwhere == gSphere && h[2] > 0) { + ld rad = hypot3(h); // turn M_PI to -M_PI ld rad_to = M_PI + M_PI - rad; ld r = -rad_to / rad; - x *= r; - y *= r; - z *= r; + h *= r; spherepoints++; } - if(eyemod) x += eyemod * z * vid.eye; + if(eyemod) h[0] += eyemod * h[2] * vid.eye; } extern int besti; void drawTriangle(triangle& t) { - rugpoint& m1 = *t.m[0]; - rugpoint& m2 = *t.m[1]; - rugpoint& m3 = *t.m[2]; - if(!m1.valid || !m2.valid || !m3.valid) return; - if(m1.dist >= sightrange+.51 || m2.dist >= sightrange+.51 || m3.dist >= sightrange+.51) - return; + using namespace hyperpoint_vec; + for(int i: {0,1,2}) { + if(!t.m[i]->valid) return; + if(t.m[i]->dist >= sightrange+.51) return; + } dt++; int spherepoints = 0; - double x1, y1, z1; - double x2, y2, z2; - double x3, y3, z3; - getco(m1,x1,y1,z1, spherepoints); - getco(m2,x2,y2,z2, spherepoints); - getco(m3,x3,y3,z3, spherepoints); + array h; + for(int i: {0,1,2}) getco(t.m[i], h[i], spherepoints); if(spherepoints == 1 || spherepoints == 2) return; - double xa = x2-x1, ya = y2-y1, za = z2-z1; - double xb = x3-x1, yb = y3-y1, zb = z3-z1; + hyperpoint hc = (h[1] - h[0]) ^ (h[2] - h[0]); + double hch = hypot3(hc); - double xn = ya * zb - za * yb; - double yn = za * xb - xa * zb; - double zn = xa * yb - ya * xb; - double h = sqrt(xn*xn+yn*yn+zn*zn); + glNormal3f(hc[0]/hch,hc[1]/hch,hc[2]/hch); - glNormal3f(xn/h,yn/h,zn/h); - - glTexCoord2f(m1.x1, m1.y1); - glVertex3f(x1, y1, z1); - glTexCoord2f(m2.x1, m2.y1); - glVertex3f(x2, y2, z2); - glTexCoord2f(m3.x1, m3.y1); - glVertex3f(x3, y3, z3); + for(int i: {0,1,2}) { + glTexCoord2f(t.m[i]->x1, t.m[i]->y1); + glVertex3f(h[i][0], h[i][1], h[i][2]); + } } GLuint FramebufferName = 0; @@ -913,9 +929,13 @@ void init() { drawthemap(); genrug = false; - buildRug(); qvalid = 0; dt = 0; queueiter = 0; + err_zero = err_zero_current; + buildRug(); + if(rug_perspective) + push_all_points(2, -1); + currentrot = Id; } @@ -1006,12 +1026,12 @@ hyperpoint gethyper(ld x, ld y) { auto r0 = triangles[i].m[0]; auto r1 = triangles[i].m[1]; auto r2 = triangles[i].m[2]; - double dx1 = r1->flat[0] - r0->flat[0]; - double dy1 = r1->flat[1] - r0->flat[1]; - double dx2 = r2->flat[0] - r0->flat[0]; - double dy2 = r2->flat[1] - r0->flat[1]; - double dxm = mx - r0->flat[0]; - double dym = my - r0->flat[1]; + double dx1 = r1->getglue()->flat[0] - r0->getglue()->flat[0]; + double dy1 = r1->getglue()->flat[1] - r0->getglue()->flat[1]; + double dx2 = r2->getglue()->flat[0] - r0->getglue()->flat[0]; + double dy2 = r2->getglue()->flat[1] - r0->getglue()->flat[1]; + double dxm = mx - r0->getglue()->flat[0]; + double dym = my - r0->getglue()->flat[1]; // A (dx1,dy1) = (1,0) // B (dx2,dy2) = (0,1) double det = dx1*dy2 - dy1*dx2;