new projections based on cartography

classes.cpp

@@ -940,6 +940,7 @@ namespace mf {
   static const flagtype equivolume = 1024;
   static const flagtype twopoint = 2048;
   static const flagtype uses_bandshift = 4096;
+  static const flagtype broken = 8192; /* in spherical case, these are broken along the meridian 180 deg */
   
   static const flagtype band = (cylindrical | pseudocylindrical | uses_bandshift);
   static const flagtype pseudoband = (pseudocylindrical | uses_bandshift);
@@ -972,9 +973,13 @@ enum eModel : int {
   mdRotatedHyperboles, mdSpiral, mdPerspective,
   // 20..24
   mdEquivolume, mdCentralInversion, mdSimulatedPerspective, mdTwoHybrid, mdGeodesic,
-  // 25
+  // 25..27
-  mdMollweide, mdCentralCyl, mdCollignon, mdHorocyclic, mdQuadrant, mdAxial, mdAntiAxial,
+  mdMollweide, mdCentralCyl, mdCollignon, 
-  // 32..
+  // 28..31 
+  mdHorocyclic, mdQuadrant, mdAxial, mdAntiAxial,
+  // 32..38
+  mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdMiller, mdGallStereographic, mdWinkelTripel,
+  // 39..
   mdGUARD, mdPixel, mdHyperboloidFlat, mdPolynomial, mdManual
   };
 #endif
@@ -1019,6 +1024,13 @@ EX vector<modelinfo> mdinf = {
   {X3("quadrant coordinates"), mf::euc_boring, DEFAULTS},
   {X3("axial coordinates"), mf::euc_boring, DEFAULTS},
   {X3("anti-axial coordinates"), mf::euc_boring, DEFAULTS},
+  {X3("Werner projection"), mf::euc_boring | mf::broken, DEFAULTS}, // keep distances from pole, and distances along parallels
+  {X3("Aitoff projection"), mf::euc_boring | mf::broken, DEFAULTS}, // halve longitudes, do azequid, double x
+  {X3("Hammer projection"), mf::euc_boring | mf::broken, DEFAULTS}, // halve longitudes, do azequia, double x
+  {X3("loximuthal projection"), mf::euc_boring | mf::broken, DEFAULTS}, // map loxodromes azimuthally and equidistantly
+  {X3("Miller projection"), mf::euc_boring | mf::band, DEFAULTS}, // scale latitude 4/5 -> Mercator -> 5/4
+  {X3("Gall stereographic"), mf::euc_boring | mf::band, DEFAULTS}, // like central cylindrical but stereographic
+  {X3("Winkel tripel"), mf::euc_boring | mf::broken, DEFAULTS}, // mean of equirec and Aitoff
   {X3("guard"), 0, DEFAULTS},
   {X3("polynomial"), mf::conformal, DEFAULTS},
   };

drawing.cpp

@@ -777,6 +777,7 @@ ld period_at(ld y) {
   
   switch(pmodel) {
     case mdBand:
+    case mdMiller:
       return m * 4;
     case mdSinusoidal:
       return m * 2 * cos(y * M_PI);
@@ -1544,7 +1545,75 @@ void dqi_poly::draw() {
     cnt = i;
     return;
     }
+  
+  int broken_coord = models::get_broken_coord(pmodel);
+  static bool in_broken = false;
+  if(broken_coord && !in_broken) {
 
+    int zcoord = broken_coord;
+    int ycoord = 3 - zcoord;
+    
+    vector<hyperpoint> all;
+    for(int i=0; i<cnt; i++) 
+      all.push_back(V.T * glhr::gltopoint((*tab)[offset+i]));
+    int fail = 0;
+    int last_fail;
+    for(int i=0; i<cnt-1; i++) 
+      if(all[i][0] * all[i+1][0] <= 0 && (all[i][0] * all[i+1][zcoord] - all[i+1][0] * all[i][zcoord]) * (all[i][0] - all[i+1][0]) < 0)
+        last_fail = i, fail++;
+    vector<glvertex> v;
+    dqi_poly p = *this;
+    p.tab = &v;
+    p.offset = 0;
+    p.V.T = Id;
+    if(fail) {
+      dynamicval<bool> ib(in_broken, true);
+      ld part = ilerp(all[last_fail][0], all[last_fail+1][0], 0);
+      hyperpoint initial = normalize(lerp(all[last_fail], all[last_fail+1], 1 - (1-part) * .99));
+      bool have_initial = true;
+      v.push_back(glhr::pointtogl(initial));
+      last_fail++;
+      int at = last_fail;
+      do {
+        v.push_back(glhr::pointtogl(all[at]));
+        if(at == cnt-1 && all[at] != all[0]) {
+          p.cnt = isize(v); p.draw(); v.clear(); at = 0;
+          have_initial = false;
+          }
+        int next = at+1;
+        if(next == cnt) next = 0;
+        if(all[at][0] * all[next][0] <= 0 && (all[at][0] * all[next][zcoord] - all[next][0] * all[at][zcoord]) * (all[at][0] - all[next][0]) < 0) {
+          ld part = ilerp(all[at][0], all[next][0], 0);
+          if(part > .999 || part < .001) { in_broken = false; return; }
+          hyperpoint final = normalize(lerp(all[at], all[next], part * .99));
+          v.push_back(glhr::pointtogl(final));
+          if(have_initial) {
+            int max = 4 << vid.linequality;
+            if(final[0] * initial[0] > 0) {
+              for(int i=1; i<=max; i++)
+                v.push_back(glhr::pointtogl(lerp(final, initial, i * 1. / max)));
+              }
+            else {
+              hyperpoint end = Hypc;
+              end[ycoord] = final[ycoord] > 0 ? 1 : -1;
+              for(int i=1; i<=max; i++)
+                v.push_back(glhr::pointtogl(lerp(final, end, i * 1. / max)));
+              for(int i=1; i<=max; i++)
+                v.push_back(glhr::pointtogl(lerp(end, initial, i * 1. / max)));
+              }
+            }
+          p.cnt = isize(v); p.draw(); v.clear();
+          initial = normalize(lerp(all[at], all[next], 1 - (1-part) * .99));
+          have_initial = true;
+          v.push_back(glhr::pointtogl(initial));
+          }
+        at = next;
+        }
+      while(at != last_fail);
+      return;
+      }
+    }
+  
   if(sphere && pmodel == mdTwoPoint && !in_twopoint) {
     #define MAX_PHASE 4
     vector<glvertex> phases[MAX_PHASE];

hyper.h

@@ -241,6 +241,7 @@ struct projection_configuration {
   ld clip_min, clip_max;
   ld model_orientation, halfplane_scale, model_orientation_yz;  
   ld collignon_parameter; 
+  ld aitoff_parameter, miller_parameter, loximuthal_parameter, winkel_parameter;
   bool collignon_reflected;
   string formula;
   eModel basic_model;
@@ -267,6 +268,10 @@ struct projection_configuration {
     spiral_cone = 360;
     use_atan = false;
     product_z_scale = 1;
+    aitoff_parameter = .5;
+    miller_parameter = .8;
+    loximuthal_parameter = 0;
+    winkel_parameter = .5;
     }
   };
 

hypgraph.cpp

@@ -364,10 +364,14 @@ EX void apply_nil_rotation(hyperpoint& H) {
   }
 
 EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
+  apply_other_model(H_orig, ret, pmodel);
+  }
+
+EX void apply_other_model(shiftpoint H_orig, hyperpoint& ret, eModel md) {
 
   hyperpoint H = H_orig.h;
   
-  if(models::product_model(pmodel)) {
+  if(models::product_model(md)) {
     ld zlev = zlevel(H_orig.h);
     H_orig.h /= exp(zlev);
     hybrid::in_underlying_geometry([&] { applymodel(H_orig, ret); });
@@ -376,7 +380,7 @@ EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
     return;    
     }
   
-  switch(pmodel) {
+  switch(md) {
     case mdPerspective: {
       if(prod) H = product::inverse_exp(H);
       apply_nil_rotation(H);
@@ -832,6 +836,25 @@ EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
       else 
         makeband(H_orig, ret, band_conformal);
       break;
+    
+    case mdMiller:
+      makeband(H_orig, ret, [] (ld& x, ld& y) {
+        // y *= pconf.miller_parameter;
+        band_conformal(x, y);
+        // y /= pconf.miller_parameter;
+        });        
+      break;
+      
+    case mdLoximuthal:
+      makeband(H_orig, ret, [] (ld&x, ld &y) {
+        ld orig_y = y;
+        band_conformal(x, y);
+        y -= pconf.loximuthal_parameter;
+        orig_y -= pconf.loximuthal_parameter;
+        if(y) x = x * orig_y / y;
+        y = orig_y;
+        });
+      break;
       
     case mdTwoPoint: 
       makeband(H_orig, ret, make_twopoint);
@@ -861,6 +884,59 @@ EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
       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; });
       break;
 
+    case mdGallStereographic: 
+      makeband(H_orig, ret, [] (ld& x, ld& y) { 
+        y = 2 * sin_auto(y) / (1 + cos_auto(y));
+        ld top = vid.yres * M_PI / current_display->radius; if(y>top) y=top; if(y<-top) y=-top;
+        });
+      break;
+    
+    case mdAitoff: case mdHammer: case mdWinkelTripel:
+      makeband(H_orig, ret, [&] (ld& x, ld& y) {
+        ld ox = x, oy = y;
+        x *= pconf.aitoff_parameter;
+
+        ld x0 = sin_auto(x) * cos_auto(y);
+        ld y0 = cos_auto(x) * cos_auto(y);
+        ld z0 = sin_auto(y);
+        
+        ld d = acos_auto(y0);
+        ld d0 = hypot(x0, z0);
+        
+        if(md == mdAitoff || md == mdWinkelTripel) ;
+        else if(sphere) d = sqrt(2*(1 - cos(d))) * M_PI / 2;
+        else d = sqrt(2*(cosh(d) - 1)) / 1.5;
+
+        x = x0 * d / d0 / pconf.aitoff_parameter, y = z0 * d / d0;
+        if(md == mdWinkelTripel) 
+          x = lerp(x, ox, pconf.model_transition),
+          y = lerp(y, oy, pconf.model_transition);
+          
+        });
+      break;
+    
+    case mdWerner: {
+
+      models::apply_orientation_yz(H[1], H[2]);
+      models::apply_orientation(H[0], H[1]);
+
+      find_zlev(H); // ignored for now
+      
+      ld r = hdist0(H);
+      if(r == 0) { ret = H; return; }
+      ld angle = atan2(H[0], H[1]);
+      angle *= sin_auto(r) / r;
+      
+      ret[0] = sin(angle) * r;
+      ret[1] = cos(angle) * r;
+      ret[2] = 0;
+      ret[3] = 1;
+
+      models::apply_orientation(ret[1], ret[0]);
+      models::apply_orientation_yz(ret[2], ret[1]);
+      break;
+      }
+
     case mdCollignon: 
       find_zlev(H_orig.h);
       makeband(H_orig, ret, [] (ld& x, ld& y) { 
@@ -1953,6 +2029,23 @@ EX void draw_boundary(int w) {
   if(elliptic && !among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal, mdMollweide, mdCollignon))
     circle_around_center(M_PI/2, periodcolor, 0, PPR::CIRCLE);
   
+  int broken_coord = models::get_broken_coord(pmodel);
+  if(broken_coord) {
+    int unbroken_coord = 3 - broken_coord;
+    const ld eps = 1e-3;
+    const ld rem = sqrt(1-eps*eps);
+    for(int s: {-1, 1}) {
+      for(int a=1; a<180; a++) {
+        hyperpoint h = Hypc;
+        h[broken_coord] = -sin_auto(a*degree) * rem;
+        h[0] = sin_auto(a*degree) * eps * s;
+        h[unbroken_coord] = cos_auto(a*degree);
+        curvepoint(h);
+        }
+      queuecurve(shiftless(Id), periodcolor, 0, PPR::CIRCLE).flags |= POLY_FORCEWIDE;
+      }
+    }
+  
   switch(pmodel) {
   
     case mdTwoPoint: {
@@ -1979,10 +2072,12 @@ EX void draw_boundary(int w) {
       return;
       }
     
-    case mdBand: case mdBandEquidistant: case mdBandEquiarea: case mdSinusoidal: case mdMollweide: case mdCentralCyl: case mdCollignon: {
+    case mdBand: case mdBandEquidistant: case mdBandEquiarea: case mdSinusoidal: case mdMollweide: case mdCentralCyl: case mdCollignon: 
+    case mdGallStereographic: case mdMiller:
+    {
       if(GDIM == 3) return;
       if(pmodel == mdBand && pconf.model_transition != 1) return;
-      bool bndband = (among(pmodel, mdBand, mdCentralCyl) ? hyperbolic : sphere);
+      bool bndband = (among(pmodel, mdBand, mdMiller, mdGallStereographic, mdCentralCyl) ? hyperbolic : sphere);
       transmatrix T = spin(-pconf.model_orientation * degree);
       ld right = M_PI/2 - 1e-5;
       if(bndband) 

models.cpp

@@ -207,7 +207,15 @@ EX namespace models {
     if(m == mdHorocyclic)
       return hyperbolic;
     return
-      among(m, mdHalfplane, mdPolynomial, mdPolygonal, mdTwoPoint, mdJoukowsky, mdJoukowskyInverted, mdSpiral, mdSimulatedPerspective, mdTwoHybrid, mdHorocyclic, mdAxial, mdAntiAxial, mdQuadrant) || mdBandAny();
+      among(m, mdHalfplane, mdPolynomial, mdPolygonal, mdTwoPoint, mdJoukowsky, mdJoukowskyInverted, mdSpiral, mdSimulatedPerspective, mdTwoHybrid, mdHorocyclic, mdAxial, mdAntiAxial, mdQuadrant,
+        mdWerner, mdAitoff, mdHammer, mdLoximuthal, mdWinkelTripel) || mdBandAny();
+    }
+
+  /** @brief returns the broken coordinate, or zero */
+  EX int get_broken_coord(eModel m) {
+    if(m == mdWerner) return 1;
+    if(sphere) return (mdinf[m].flags & mf::broken) ? 2 : 0;
+    return 0;
     }
   
   EX bool is_perspective(eModel m) {