nilv:: Nil shader, and improved geodesics

basegraph.cpp

@@ -233,10 +233,12 @@ void display_data::set_projection(int ed) {
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::halfplane3;
     if(DIM == 3 && hyperbolic && apply_models && pmodel == mdPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardH3, pers3 = true;
-    if(DIM == 3 && (euclid || sol || nil) && apply_models && pmodel == mdPerspective)
+    if(DIM == 3 && translatable && apply_models && pmodel == mdPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true;
-    if(DIM == 3 && apply_models && pmodel == mdGeodesic)
+    if(DIM == 3 && apply_models && pmodel == mdGeodesic && sol)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSolv, pers3 = true;
+    if(DIM == 3 && apply_models && pmodel == mdGeodesic && nil)
+      shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardNil, pers3 = true;
     if(DIM == 3 && sphere && apply_models && pmodel == mdPerspective) {
       shaderside_projection = true; pers3 = true;
       int sp = spherephase & 3;

drawing.cpp

@@ -793,20 +793,6 @@ void debug_this() { }
 void dqi_poly::draw() {
   if(flags & POLY_DEBUG) debug_this();
 
-  if(nil && vid.usingGL && pmodel == mdPerspective && (current_display->set_all(global_projection), shaderside_projection)) {
-    auto npoly = *this;
-    glcoords.clear();
-    for(int i=0; i<cnt; i++)
-      glcoords.push_back(glhr::pointtogl(nisot::inverse_exp(V * glhr::gltopoint( (*tab)[offset+i]), nisot::iTable)));
-      
-    npoly.offset = 0;
-    npoly.tab = &glcoords;
-    npoly.V = Id;
-    set_width(1);
-    npoly.gldraw();
-    return;
-    }
-
   dynamicval<ld> bs(hr::band_shift, band_shift);
   if(!hyperbolic && among(pmodel, mdPolygonal, mdPolynomial)) {
     bool any = false;

graph.cpp

@@ -5057,16 +5057,26 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
       }
     noclipped++;
     }
-  if(pmodel == mdGeodesic) {
+  if(pmodel == mdGeodesic && sol) {
     using namespace hyperpoint_vec;
     hyperpoint H = tC0(V);
-    if(abs(H[0]) <= 2 && abs(H[1]) <= 2 && abs(H[2]) <= 2) ;
+    if(abs(H[0]) <= 3 && abs(H[1]) <= 3 && abs(H[2]) <= 3 ) ;
     else {
       hyperpoint H2 = nisot::inverse_exp(H, nisot::iLazy);
       for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -.2) return;
       }
     noclipped++;
     }
+  if(pmodel == mdGeodesic && nil) {
+    using namespace hyperpoint_vec;
+    hyperpoint H = tC0(V);
+    if(abs(H[0]) <= 3 && abs(H[1]) <= 3 && abs(H[2]) <= 3 ) ;
+    else {
+      hyperpoint H2 = nisot::inverse_exp(H, nisot::iLazy);
+      for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -2) return;
+      }
+    noclipped++;
+    }
 #endif
 
   #if CAP_SHAPES

hyper.h

@@ -4038,7 +4038,7 @@ namespace glhr {
   
   enum class shader_projection { standard, band, halfplane, standardH3, standardR3, 
     standardS30, standardS31, standardS32, standardS33, 
-    ball, halfplane3, band3, flatten, standardSolv,
+    ball, halfplane3, band3, flatten, standardSolv, standardNil,
     MAX 
     };
   
@@ -5663,6 +5663,7 @@ namespace solv {
 
 
 namespace nilv {
+  extern string nilshader;
   extern array<vector<hyperpoint>, 22> facevertices;
   void software_renderer(dqi_poly *p);
   }

nonisotropic.cpp

@@ -341,7 +341,7 @@ namespace nilv {
     // v[1] = c sin alpha
     // v[2] = w
     
-    if(v[0] == 0 && v[1] == 0) return v;
+    if(v[0] == 0 && v[1] == 0) return point31(v[0], v[1], v[2]);
     
     if(v[2] == 0) return point31(v[0], v[1], v[0] * v[1] / 2);
     
@@ -361,7 +361,7 @@ namespace nilv {
     
     ld side = h[2] - h[0] * h[1] / 2;
     
-    if(hypot_d(2, h) < 1e-6) return h;
+    if(hypot_d(2, h) < 1e-6) return point3(h[0], h[1], h[2]);
     else if(side > 1e-6) {
       wmin = 0, wmax = 2 * M_PI;
       }
@@ -370,23 +370,23 @@ namespace nilv {
       }
     else return point3(h[0], h[1], 0);
     
-    ld alpha_total = atan2(h[1], h[0]);
+    ld alpha_total = h[0] ? atan(h[1] / h[0]) : M_PI/2;
   
-    ld cmul;
+    ld b;
     if(abs(h[0]) > abs(h[1]))
-      cmul = h[0] / 2 / cos(alpha_total);
+      b = h[0] / 2 / cos(alpha_total);
     else
-      cmul = h[1] / 2 / sin(alpha_total);  
+      b = h[1] / 2 / sin(alpha_total);  
     
-    ld sa = sin(2 * alpha_total);
+    ld s = sin(2 * alpha_total);
     
     for(int it=0;; it++) {
-      ld w = (wmin + wmax) / 2;    
-      ld c = cmul / sin(w/2);
-      ld z = w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sa));
-      
+      ld w = (wmin + wmax) / 2;
+      ld z = b * b * (s + (sin(w) - w)/(cos(w) - 1)) + w;
+
       if(it == iterations) {
         ld alpha = alpha_total - w/2;
+        ld c = b / sin(w/2);
         return point3(c * w * cos(alpha), c * w * sin(alpha), w);
         }
       if(h[2] > z) wmin = w;
@@ -394,6 +394,36 @@ namespace nilv {
       }
     }
   
+  string nilshader = 
+    "vec4 inverse_exp(vec4 h) {"
+      "float wmin, wmax;"
+      "float side = h[2] - h[0] * h[1] / 2.;"
+      "if(h[0]*h[0] + h[1]*h[1] < 1e-12) return vec4(h[0], h[1], h[2], 1);"
+      "if(side > 1e-6) { wmin = 0.; wmax = 2.*PI; }"
+      "else if(side < -1e-6) { wmin = -2.*PI; wmax = 0.; }"
+      "else return vec4(h[0], h[1], 0., 1.);"
+      "float at = h[0] != 0. ? atan(h[1] / h[0]) : PI/2.;"
+      "float b = abs(h[0]) > abs(h[1]) ? h[0] / 2. / cos(at) : h[1] / 2. / sin(at);"
+      "float s = sin(2. * at);"
+
+      "for(int it=0; it<50; it++) {"
+        "float w = (wmin + wmax) / 2.;"
+        "float z = b * b * (s + (sin(w) - w)/(cos(w) - 1.)) + w;"
+        "if(h[2] > z) wmin = w;"
+        "else wmax = w;"
+        "}"
+
+      "float w = (wmin + wmax) / 2.;"
+      "float alpha = at - w/2.;"
+      "float c = b / sin(w/2.);"
+      "return vec4(c*w*cos(alpha), c*w*sin(alpha), w, 1.);"
+
+      /* "float w = atan(side) * 4.;"
+          "float alpha = at - w/2.;"
+          "float c = b / sin(w/2.);"
+          "return vec4(c*w*cos(alpha), c*w*sin(alpha), w, 1.);" */
+      "}";
+  
   struct mvec : array<int, 3> {
     
     mvec() { }

shaders.cpp

@@ -562,6 +562,7 @@ void init() {
     bool hp = among(sp, shader_projection::halfplane, shader_projection::halfplane3);
     bool sh3 = (sp == shader_projection::standardH3);
     bool ssol = (sp == shader_projection::standardSolv);
+    bool snil = (sp == shader_projection::standardNil);
     bool sr3 = (sp == shader_projection::standardR3);
     bool ss30 = (sp == shader_projection::standardS30);
     bool ss31 = (sp == shader_projection::standardS31);
@@ -569,13 +570,14 @@ void init() {
     bool ss33 = (sp == shader_projection::standardS33);
     bool ss3 = ss30 || ss31 || ss32 || ss33;
     
-    bool s3 = (sh3 || sr3 || ss3 || ssol);
+    bool s3 = (sh3 || sr3 || ss3 || ssol || snil);
     bool dim3 = s3 || among(sp, shader_projection::ball, shader_projection::halfplane3, shader_projection::band3);
     bool dim2 = !dim3;
     bool ball = (sp == shader_projection::ball);
     bool flatten = (sp == shader_projection::flatten);
     
     programs[i][j] = new GLprogram(stringbuilder(
+      1,       "#define PI 3.14159265358979324\n",
 
       1,       "attribute mediump vec4 aPosition;",
       texture, "attribute mediump vec2 aTexture;",
@@ -621,6 +623,7 @@ void init() {
       1,       "  }",
       
       ssol,    solv::solshader,      
+      snil,    nilv::nilshader,
 
       1,       "void main() {",  
       texture,   "vTexCoord = aTexture;",
@@ -659,9 +662,10 @@ void init() {
       ssol,      "float d = sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]);",
       ssol,      "float ad = (d == 0.) ? 0. : (d < 1.) ? min(atanh(d), 10.) : 10.; ",
       ssol,      "float m = ad / d / 11.; t[0] *= m; t[1] *= m; t[2] *= m; ",
+      snil,      "t = inverse_exp(t);",
        
       sh3,       "float fogs = (uFogBase - acosh(t[3]) / uFog);",
-      sr3,       "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / uFog);",
+      sr3||snil, "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / uFog);",
       ssol,      "float fogs = (uFogBase - ad / uFog);",
       
       ss30,      "float fogs = (uFogBase - (6.284 - acos(t[3])) / uFog); t = -t; ",

system.cpp

@@ -1154,9 +1154,9 @@ void set_geometry(eGeometry target) {
     if(binarytiling || WDIM == 3 || penrose) variation = eVariation::pure;
     #endif
     if(DIM == 3 && old_DIM == 2 && pmodel == mdDisk) pmodel = mdPerspective;
-    if(sol && old_DIM == 2) pmodel = mdSolPerspective;
-    if(DIM == 2 && among(pmodel, mdPerspective, mdSolPerspective)) pmodel = mdDisk;
-    if(sol && old_DIM == 2 && vid.texture_step < 4) vid.texture_step = 4;
+    if(nonisotropic && old_DIM == 2) pmodel = mdGeodesic;
+    if(DIM == 2 && among(pmodel, mdPerspective, mdGeodesic)) pmodel = mdDisk;
+    if(nonisotropic && old_DIM == 2 && vid.texture_step < 4) vid.texture_step = 4;
     }
   }