4-dihedron works

MANIFEST.in

@@ -1 +1 @@
-include data/*.off
+include data/*

README.md

@@ -1,18 +1,32 @@
 # antitile
-Manipulation of tilings using Python. This package is designed to work with [Antiprism](https://github.com/antiprism/antiprism) but can be used on its own.
-
-Under construction.
+Manipulation of polyhedra and tilings using Python. This package is designed to work with [Antiprism](https://github.com/antiprism/antiprism) but can be used on its own.
 
 ## Installation
 
     pip3 install git+git://github.com/brsr/antitile.git
 
 ## Usage
-There are currently 4 scripts included with this package:
+There are currently 4 programs included with this package:
 * `view_off.py`: A viewer for OFF files using matplotlib.
 * `balloon.py`: Balloon tiling of the sphere
 * `sgs.py`: Similar grid subdivision of tilings
 * `sgsstats.py`: Statistics of SGS tilings
+These can be chained together with programs from Antiprism.
+
+OFF files for the regular icosahedron, octahedron, tetrahedron, cube, and 3- and 4-point dihedra are included (although where Python puts them may depend on your system).
+
+## Examples
+Statistics of a geodesic polyhedron (using what geodesic dome people call Method 1):
+
+	sgs.py -a 5 -b 3 icosahedron.off | sgsstats.py
+
+Visualize a Goldberg polyhedron, with color:
+
+    sgs.py -a 5 -b 3 icosahedron.off | off_color -v M -m color_map.txt | pol_recip | view_off.py
+
+Canonical form (no skew faces) of a quadrilateral-faced similar grid subdivision polyhedron:
+
+    sgs.py -a 5 -b 3 cube.off | canonical | view_off.py
 
 ## For Contributors
 This code makes heavy use of vectorized operations on NumPy multidimensional arrays, which are honestly pretty impenetrable until you get familiar with them. (And, uh, even after that.) I use the convention that the last axis of an array specifies the spatial coordinates:

TODO.md

@@ -0,0 +1,14 @@
+#To Do
+* Clean up the code
+* Document
+* Replace the stitching logic in `antitile/sgs.py` so that for Class I and II
+subdivision, the faces do not need to be oriented. In theory this would 
+allow SGS to operate on tilings on unorientable surfaces, although 
+self-intersecting surfaces might have quirks to deal with.
+* Allow "mixed" triangle & quad grids, at least for Class I. (Triangle 
+class I and quad class I work together, but triangle class II works with
+particular quad class III subdivisions, not every triangle class III has a
+compatible quad class III, and quad class II doesn't work with triangle at 
+all. Probably need to work out the theory of what's compatible with what.)
+* Allow arbitrary face figures other than SGS (would probably need whole new 
+program)

antitile/breakdown.py

@@ -75,11 +75,13 @@ def _t(self):
         x = vertices[:, 0]
         y = vertices[:, 1]
 
+
         #featural lines
+        q = a - b        
         line1 = x == 0
         line1_seg = y <= max(a,b)
         line2 = y == b
-        line2_seg = x >= (0 if a > b else a-b)     
+        line2_seg = x >= min(0,q)     
         line3 = y == a - x
         line3_seg = y >= min(a,b)
         group[line1 & line1_seg] = 1
@@ -128,23 +130,24 @@ def _q(self):
         x = vertices[:, 0]
         y = vertices[:, 1]
         #featural lines
+        q = a - b
         line1 = x == 0
-        line1_seg = y <= a
+        line1_seg = y <= max(a,b)
         line2 = y == b
-        line2_seg = x >= 0      
+        line2_seg = x >= min(0,q)      
         line3 = x == a-b
-        line3_seg = y >= b
+        line3_seg = y >= min(a,b)
         line4 = y == a
-        line4_seg = x <= a - b   
+        line4_seg = x <= max(0,q)  
         group[line1 & line1_seg] = 1
         group[line2 & line2_seg] = 2
         group[line3 & line3_seg] = 3
         group[line4 & line4_seg] = 4
-        group[line1 & ~line1_seg] = 11
-        group[line2 & ~line2_seg] = 12
-        group[line3 & ~line3_seg] = 13
-        group[line4 & ~line4_seg] = 14
-
+        if a == b:             
+            group[line1 & ~line1_seg] = 11
+            group[line2 & ~line2_seg] = 12
+            group[line3 & ~line3_seg] = 13
+            group[line4 & ~line4_seg] = 14             
         right = a*y - b*x
         left = a*x + b*y
         anorm = (a**2+b**2)

antitile/projection.py

@@ -50,6 +50,7 @@ def bary_tri(tri, vertices):
     ainv = np.linalg.inv(a)
     return b.dot(ainv)
 
+
 def square_to_circle(xy, rotation=1):#np.exp(1j*np.pi/4)):
     """Transforms square on [0,1]^2 to unit circle"""
     pts = 2*xy - 1
@@ -71,9 +72,14 @@ def tri_to_circle(beta, rotation=1, pts=TRI_C):
     result = r*np.exp(1j*angle)
     return xmath.complex_to_float2d(result)
 
+SQ_C = np.array([1, 1j, -1, -1j])
+SQ_R = xmath.complex_to_float2d(SQ_C)
+
 def _sq_cir(bkdn, abc, freq, tweak):
-    #FIXME this is hella broken
-    return square_to_quad(square_to_circle(bkdn.coord)[:, np.newaxis], abc)
+    sc = square_to_circle(bkdn.coord)
+    sc2 = sc/np.sqrt(2) + 0.5
+    result = square_to_quad(sc2[:, np.newaxis], abc)
+    return result
 
 
 def _tri_cir(bkdn, abc, freq, tweak):
@@ -107,13 +113,8 @@ def equidistant(disk):
 #triangles -> spherical triangle
 
 def tri_naive_slerp(bary, base_pts):
-    omegas = xmath.spherical_distance(base_pts, np.roll(base_pts, 1, axis=0))
-    max_diff = np.abs(omegas - np.roll(omegas, 1)).max()
-    if not np.isclose(max_diff, 0):
-        raise ValueError("naive_slerp used with non-equilateral face. " +
-                         "Difference is " + str(max_diff) + " radians.")
-    omega = omegas[0]
-    b = np.sin(omega * bary) / np.sin(omega)
+    angle = xmath.central_angle_equilateral(base_pts)
+    b = np.sin(angle * bary) / np.sin(angle)
     return b.dot(base_pts)
 
 def tri_areal(beta, triangle):
@@ -141,7 +142,7 @@ def tri_areal(beta, triangle):
            [ 0.        ,  0.92387953,  0.38268343]])
     """
     triangle = xmath.normalize(triangle)
-    area = xmath.spherical_triangle_area(triangle[0], triangle[1], triangle[2])
+    area = xmath.triangle_solid_angle(triangle[0], triangle[1], triangle[2])
     area_i = beta * area
     triangle_iplus1 = np.roll(triangle, -1, axis=0)
     triangle_iplus2 = np.roll(triangle, 1, axis=0)
@@ -211,6 +212,20 @@ def tri_intersections(lindex, base_pts, freq, tweak=False):
 
 
 #squares -> spherical quadrilateral
+def square_naive_slerp(xy, base_pts):
+    angle = xmath.central_angle_equilateral(base_pts)
+    x, y = xy[..., 0], xy[..., 1]
+    sx = np.sin(x*angle)
+    sy = np.sin(y*angle)    
+    scx = np.sin((1-x)*angle)
+    scy = np.sin((1-y)*angle)
+    a = scx * scy
+    b = sx * scy
+    c = sx * sy
+    d = scx * sy
+    mat = np.stack([a, b, c, d], axis=-1) / np.sin(angle)**2
+    return mat.dot(base_pts)
+
 def square_slerp(xy, abcd):
     """Transforms a square in [0,1]^2 to a spherical quadrilateral
     defined by abcd, using spherical linear interpolation
@@ -224,6 +239,9 @@ def square_slerp(xy, abcd):
             (4, ..., 3)
     Returns:
         Coordinates on the sphere."""
+    #FIXME returns coordinates off the sphere.
+    #square_slerp(xy, abcd) == square_slerp(xy[:,::-1], abcd[[0,3,2,1]]) ?
+    #if not, what can we do to make it isotropic?
     a, b, c, d = abcd[0], abcd[1], abcd[2], abcd[3]
     x, y = xy[..., 0], xy[..., 1]
     ab = xmath.slerp(a, b, x)
@@ -326,9 +344,11 @@ def project_sphere(sphcoords, zfunc=np.arcsin, scale=180 / np.pi):
 
 SLERP = {3: lambda bkdn, abc, freq, tweak: tri_naive_slerp(bkdn.coord, abc),
          4: lambda bkdn, abc, freq, tweak:
-             square_slerp(bkdn.coord[:, np.newaxis], abc)}
+             square_naive_slerp(bkdn.coord, abc)}
 
-AREAL = {3: lambda bkdn, abc, freq, tweak: tri_areal(bkdn.coord, abc)}
+OTHER = {3: lambda bkdn, abc, freq, tweak: tri_areal(bkdn.coord, abc),
+         4: lambda bkdn, abc, freq, tweak:
+             square_slerp(bkdn.coord[:, np.newaxis], abc)}
 
 GC = {3: lambda bkdn, abc, freq, tweak:
             tri_intersections(bkdn.lindex, abc, freq, tweak),
@@ -340,7 +360,7 @@ def project_sphere(sphcoords, zfunc=np.arcsin, scale=180 / np.pi):
 
 PROJECTIONS = {'flat':  FLAT,
                'slerp': SLERP,
-               'areal': AREAL,
+               'other': OTHER,
                'gc':    GC,
                'disk':  DISK}
 

antitile/sgs.py

@@ -160,7 +160,8 @@ def subdiv(base, freq=(2, 0), proj='flat', tweak=False):
         index = (rbkdn.base_face == i)
         base_face = base_faces[i]
         vbf = base.vertices[base_face]
-        rbkdn.vertices[index] = proj_fun(rbkdn[index], vbf, freq, tweak)
+        rbkdn_i = rbkdn[index]
+        rbkdn.vertices[index] = proj_fun(rbkdn_i, vbf, freq, tweak)
     result = tiling.Tiling(rbkdn.vertices, faces)
     result.group = group
     result.base_face = bf

antitile/xmath.py

@@ -146,6 +146,15 @@ def normalize(vectors, axis=-1):
     return np.where(n <= 0, 0, vectors / n)
 
 
+def central_angle_equilateral(pts):
+    omegas = central_angle(pts, np.roll(pts, 1, axis=0))
+    max_diff = np.abs(omegas - np.roll(omegas, 1)).max()
+    if not np.isclose(max_diff, 0):
+        raise ValueError("naive_slerp used with non-equilateral face. " +
+                         "Difference is " + str(max_diff) + " radians.")
+    return omegas[0]
+
+
 def slerp(pt1, pt2, intervals, axis=-1):
     """Spherical linear interpolation.
     >>> x = np.array([1,0,0])
@@ -158,8 +167,8 @@ def slerp(pt1, pt2, intervals, axis=-1):
            [ 0.       ,  0.       ,  1.       ]])
     """
     t = intervals
-    a = central_angle(pt1, pt2)[..., np.newaxis]
-    return (np.sin((1 - t)*a)*pt1 + np.sin((t)*a)*pt2)/np.sin(a)
+    angle = central_angle(pt1, pt2)[..., np.newaxis]
+    return (np.sin((1 - t)*angle)*pt1 + np.sin((t)*angle)*pt2)/np.sin(angle)
 
 
 def lerp(pt1, pt2, intervals):
@@ -252,7 +261,7 @@ def bearing(origin, destination, pole=np.array([0, 0, 1])):
     return np.arctan2(x, d)
 
 
-def central_angle(x, y):
+def central_angle(x, y, signed=False):
     """Central angle between vectors with respect to 0. If vectors have norm 
     1, this is the spherical distance between them.
     Args:
@@ -271,7 +280,11 @@ def central_angle(x, y):
     """
     cos = np.sum(x*y, axis=-1)
     sin = norm(np.cross(x, y), axis=-1)
-    return np.arctan(sin/cos)
+    result = np.arctan2(sin, cos)
+    if signed:
+        return result
+    else:
+        return abs(result)
 
 
 def triangle_solid_angle(a, b, c):