added model: three-point equidistant

classes.cpp

@@ -982,8 +982,8 @@ enum eModel : int {
   // 32..38
   mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdMiller, mdGallStereographic, mdWinkelTripel,
   // 39..
-  mdPoorMan, mdPanini, mdRetroCraig, mdRetroLittrow, mdRetroHammer,
+  mdPoorMan, mdPanini, mdRetroCraig, mdRetroLittrow, mdRetroHammer, mdThreePoint,
-  // 44..
+  // 45..
   mdGUARD, mdPixel, mdHyperboloidFlat, mdPolynomial, mdManual
   };
 #endif
@@ -1040,6 +1040,7 @@ EX vector<modelinfo> mdinf = {
   {X3("Craig retroazimuthal"), mf::euc_boring | mf::broken, DEFAULTS}, // retroazimuthal cylindrical
   {X3("Littrow retroazimuthal"), mf::euc_boring | mf::broken, DEFAULTS}, // retroazimuthal conformal
   {X3("Hammer retroazimuthal"), mf::euc_boring, DEFAULTS}, // retroazimuthal equidistant
+  {X3("three-point equidistant"), mf::euc_boring, DEFAULTS},
   {X3("guard"), 0, DEFAULTS},
   {X3("polynomial"), mf::conformal, DEFAULTS},
   };

hypgraph.cpp

@@ -432,6 +432,63 @@ void vr_disk(hyperpoint& ret, hyperpoint& H) {
     }
   }
 
+/** Compute the three-point projection. Currently only works in isotropic 3D spaces. */
+EX void threepoint_projection(const hyperpoint& H, hyperpoint& ret) {
+  find_zlev(H);
+
+  hyperpoint H1 = H;
+  if(true) {
+    models::apply_orientation_yz(H1[1], H1[2]);
+    models::apply_orientation(H1[0], H1[1]);
+    }
+
+  auto p = pconf.twopoint_param;
+
+  ld dist[3];
+  for(int i=0; i<3; i++) {
+    hyperpoint h1 = xspinpush0(2*M_PI*i/3, p);    
+    dist[i] = geo_dist(h1, H1);
+    }
+  
+  /* we are looking for the points (x,y,z) such that:
+     (x-xi)^2 + (y-yi)^2 + z^2 = di^2
+
+     which is equivalent to:
+     x^2+y^2+z^2 -2xxi -2yyi = di^2-xi^2-yi^2
+     
+     After setting s = x^2+y^2+z^2, we get a system of linear equations for (x,y,s)
+  */
+    
+  dynamicval<eGeometry> g(geometry, gEuclid);
+  
+  transmatrix T = Id;
+  hyperpoint v = C0;
+  for(int i=0; i<3; i++) {
+    hyperpoint pp = xspinpush0(2*M_PI*i/3, p);
+    v[i] = dist[i]*dist[i] - p*p;
+    T[i][0] = -2 * pp[0];
+    T[i][1] = -2 * pp[1];
+    T[i][2] = 1;
+    }
+  
+  transmatrix U = inverse3(T);
+  hyperpoint sxy = U * v;
+  
+  // compute the actual z based on s
+  sxy[2] = sxy[2] - sqhypot_d(2, sxy);
+  sxy[2] = sxy[2] > 0 ? sqrt(sxy[2]) : 0;
+
+  if(H1[2] < 0) sxy[2] *= -1;
+  
+  sxy[3] = 1;
+  
+  geometry = gCubeTiling;
+  
+  ret = sxy;
+  models::apply_orientation(ret[1], ret[0]);
+  models::apply_orientation_yz(ret[2], ret[1]);
+  }
+
 EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) {
 
   hyperpoint H = H_orig.h;
@@ -978,6 +1035,10 @@ EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) {
       makeband(H_orig, ret, make_twopoint);
       break;
     
+    case mdThreePoint: 
+      threepoint_projection(H, ret);
+      break;
+    
     case mdMollweide: 
       makeband(H_orig, ret, [] (ld& x, ld& y) { 
         ld theta = 
@@ -2170,6 +2231,20 @@ EX void draw_model_elements() {
       return;
       }
     
+    case mdThreePoint: {
+      vid.linewidth *= 5;
+      for(int i=0; i<=3; i++) {
+        hyperpoint h = xspinpush0(2*M_PI*i/3, pconf.twopoint_param);
+        models::apply_orientation(h[1], h[0]);
+        models::apply_orientation_yz(h[2], h[1]);
+        curvepoint(h);
+        }
+      
+      queuecurve(shiftless(Id), ringcolor, 0, PPR::SUPERLINE);
+      vid.linewidth /= 5;
+      return;
+      }
+    
     case mdBall: {
       queuecircle(current_display->xcenter, current_display->ycenter, current_display->radius, ringcolor, PPR::OUTCIRCLE, modelcolor);
       ballgeometry();

models.cpp

@@ -198,6 +198,7 @@ EX namespace models {
       return false;
     if(GDIM == 2 && pm == mdPerspective) return false;
     if(GDIM == 2 && pm == mdEquivolume) return false;
+    if(pm == mdThreePoint && !(GDIM == 3 && !nonisotropic && !prod)) return false;
     if(GDIM == 3 && among(pm, mdBall, mdHyperboloid, mdFormula, mdPolygonal, mdRotatedHyperboles, mdSpiral, mdHemisphere)) return false;
     if(pm == mdCentralInversion && !euclid) return false;
     if(pm == mdPoorMan) return hyperbolic;
@@ -211,7 +212,7 @@ EX namespace models {
     if((m == mdPerspective || m == mdGeodesic) && panini_alpha) return true;
     return
       among(m, mdHalfplane, mdPolynomial, mdPolygonal, mdTwoPoint, mdJoukowsky, mdJoukowskyInverted, mdSpiral, mdSimulatedPerspective, mdTwoHybrid, mdHorocyclic, mdAxial, mdAntiAxial, mdQuadrant,
-        mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdWinkelTripel) || mdBandAny();
+        mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdWinkelTripel, mdThreePoint) || mdBandAny();
     }
 
   /** @brief returns the broken coordinate, or zero */
@@ -236,7 +237,7 @@ EX namespace models {
   
   EX bool product_model(eModel m) {
     if(!prod) return false;
-    if(among(m, mdPerspective, mdHyperboloid, mdEquidistant)) return false;
+    if(among(m, mdPerspective, mdHyperboloid, mdEquidistant, mdThreePoint)) return false;
     return true;
     }
   
@@ -602,7 +603,7 @@ EX namespace models {
     if(vpmodel == mdHemisphere && euclid) 
       add_edit(vpconf.euclid_to_sphere);
       
-    if(among(vpmodel, mdTwoPoint, mdSimulatedPerspective, mdTwoHybrid)) 
+    if(among(vpmodel, mdTwoPoint, mdSimulatedPerspective, mdTwoHybrid, mdThreePoint)) 
       add_edit(vpconf.twopoint_param);
     
     if(vpmodel == mdFisheye)