import numpy as np
class Octahedron:
r2 = 2 ** 0.5
r3 = 3 ** 0.5
r6 = 6 ** 0.5 # r2 * r3
i2 = 1.0 / r2
i3 = 1.0 / r3
i6 = 1.0 / r6
i26 = 2. * i6
ov = {
'N': (0.0, 0.0, +1.0), 'S': (0.0, 0.0, -1.0), # North/South 0 1 (z is vertical) +90, -90 lat
'E': (0.0, +1.0, 0.0), 'W': (0.0, -1.0, 0.0), # East/West 2 3 (y is left to right) +90, -90 lon
'A': (+1.0, 0.0, 0.0), 'P': (-1.0, 0.0, 0.0), # Anterior/Posterior 4 5 (x is front to back) 0, 180 lon
}
# These rotate a 3D face to the (Z) plane
matrices = {
'NEA': np.array([[-i2, 0, i2], [i6, -i26, i6], [i3, i3, i3]]),
'NEP': np.array([[0, -i2, i2], [i26, i6, i6], [-i3, i3, i3]]),
'NWA': np.array([[0, i2, i2], [-i26, -i6, i6], [i3, -i3, i3]]),
'NWP': np.array([[i2, 0, i2], [-i6, i26, i6], [-i3, -i3, i3]]),
'SEA': np.array([[0, -i2, -i2], [-i26, i6, -i6], [i3, i3, -i3]]),
'SEP': np.array([[i2, 0, -i2], [-i6, -i26, -i6], [-i3, i3, -i3]]),
'SWA': np.array([[-i2, 0, -i2], [i6, i26, -i6], [i3, -i3, -i3]]),
'SWP': np.array([[0, i2, -i2], [i26, -i6, -i6], [-i3, -i3, -i3]]),
}
def __init__(self):
self.signs = {_face: np.sum([self.ov[v] for v in _face], axis=0) for _face in self.matrices}
self.pt_signs = {tuple(v): k for k, v in self.signs.items()}
def pt_face(self, uvw): # Which face does this point belong to?
return self.pt_signs[tuple(np.sign(uvw).astype(int).tolist())]
def pts_faces(self, uvw_s): # face.keys from np_array of 3d Octahedron/Spherical points
pts = uvw_s if isinstance(uvw_s, np.ndarray) else np.array(uvw_s)
return np.apply_along_axis(self.pt_face, -1, pts)
@classmethod
def pt_valid(cls, pt): # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron)
return np.abs(np.sum(np.abs(pt)) - 1.) < 1e-15
def pts_valid(self, uvw_s): # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron)
pts = uvw_s if isinstance(uvw_s, np.ndarray) else np.array(uvw_s)
return np.all(np.apply_along_axis(self.pt_valid, -1, pts))
def random_octalrandom_octahedral(n):
# return a random points from OctalOctant 0.
uvw = []
r_xy = np.random.uniform(0, 1, [n, 2])
for (a, b) in r_xy:
if a + b > 1:
a, b = 1 - a, 1 - b
uvw.append([a, b, 1 - a - b])
return np.array(uvw)
if __name__ == '__main__':
o = Octahedron()
pts = random_octalrandom_octahedral(100)
all_good = o.pts_valid(pts)
for k, v in o.signs.items():
f_pts = np.copysign(pts, v)
faces = np.unique(o.pts_faces(f_pts))
if len(faces) != 1:
print(f'{k}: mysterious outlier in {faces}')
elif faces[0] != k:
print(f'{k}: mysterious result in {faces}')
else:
print(f'{k}: looks good!')
Became Hot Network Question
