I have a (non-invertible) function ak([u,v,w])
This takes a point on the surface of the unit octahedron (p: such that |u|+|v|+|w| = 1) and returns a point on the surface of the unit sphere. The function isn't perfect but the intention is to keep the distance between points authalic.
I was thinking of using SciPy minimize to provide a numerical inverse, but I cannot wrap my head around it.
input: spherical pt [x,y,z],
output octahedral pts [u,v,w] such that ak([u,v,w])=[x,y,z]
My function ak() is defined like this:
def ak(p):
# Convert point on octahedron onto the sphere.
# Credit to Anders Kaseorg: https://math.stackexchange.com/questions/5016695/
# input: oct_pt is a Euclidean point on the surface of a unit octagon.
# output: UVW on a unit sphere.
a = 3.227806237143884260376580641604959964752197265625 # 𝛂 - vis. Kaseorg.
p1 = (np.pi * p) / 2.0
tp1 = np.tan(p1)
xu, xv, xw = tp1[0], tp1[1], tp1[2]
u2, v2, w2 = xu ** 2, xv ** 2, xw ** 2
y0p = xu * (v2 + w2 + a * w2 * v2) ** 0.25
y1p = xv * (u2 + w2 + a * u2 * w2) ** 0.25
y2p = xw * (u2 + v2 + a * u2 * v2) ** 0.25
pv = np.array([y0p, y1p, y2p])
return pv / np.linalg.norm(pv, keepdims=True)
This function is based on a post that I made on the Math StackExchange.
Any hints?
One idea I found that helped significantly was to convert this from a minimization problem to a root-finding problem.
The root-finders don't support constraints, so you need to change the objective function to force x to be L1 normalized before converting into Cartesian coordinates.
def fn_root(op, tx): # octa_point, target_sphere_point
norm = np.linalg.norm(op, ord=1)
return ak(op / norm) - tx
Once you do that, you can find the point where ak(op / norm) - tx = 0.
# Inverse function using numerical optimization
def inverse_ak_root(tsp):
initial_guess = np.arcsin(tsp) / (np.pi / 2)
result = root(fn_root, initial_guess, args=(tsp,), method='hybr', tol=1e-12)
assert result.success, result
result.x /= np.linalg.norm(result.x, ord=1)
return result
(I also changed the initial guess from the center of the octahedron's face to np.arcsin(tsp) / (np.pi / 2). I found this reduced how long it took to converge.)
I found that this reduced the error of your current solution by about a factor of 10^7, while using fewer calls to ak() than your previous solution.
import numpy as np
from scipy.optimize import minimize, root
def xyz_ll(spt):
x, y, z = spt[:, 0], spt[:, 1], spt[:, 2]
lat = np.degrees(np.arctan2(z, np.sqrt(x**2. + y**2.)))
lon = np.degrees(np.arctan2(y, x))
return np.stack([lat, lon], axis=1)
def ll_xyz(ll): # convert to radians.
phi, theta = np.radians(ll[:, 0]), np.radians(ll[:, 1])
x = np.cos(phi) * np.cos(theta)
y = np.cos(phi) * np.sin(theta)
z = np.sin(phi) # z is 'up'
return np.stack([x, y, z], axis=1)
def ak(p):
# Convert point on octahedron onto the sphere.
# Credit to Anders Kaseorg: https://math.stackexchange.com/questions/5016695/
# input: oct_pt is a Euclidean point on the surface of a unit octagon.
# output: UVW on a unit sphere.
a = 3.227806237143884260376580641604959964752197265625 # 𝛂 - vis. Kaseorg.
p1 = (np.pi * p) / 2.0
tp1 = np.tan(p1)
xu, xv, xw = tp1[0], tp1[1], tp1[2]
u2, v2, w2 = xu ** 2, xv ** 2, xw ** 2
y0p = xu * (v2 + w2 + a * w2 * v2) ** 0.25
y1p = xv * (u2 + w2 + a * u2 * w2) ** 0.25
y2p = xw * (u2 + v2 + a * u2 * v2) ** 0.25
pv = np.array([y0p, y1p, y2p])
return pv / np.linalg.norm(pv, keepdims=True)
def fn(op, tx): # octa_point, target_sphere_point
return np.sum((ak(op) - tx) ** 2.)
def fn_root(op, tx): # octa_point, target_sphere_point
norm = np.linalg.norm(op, ord=1)
return ak(op / norm) - tx
def constraint(op): # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron)
return np.sum(np.abs(op)) - 1
# Inverse function using numerical optimization
def inverse_ak(tsp):
initial_guess = np.sign(tsp) * [0.5, 0.5, 0.5] # the centre of the current side
# initial_guess = np.arcsin(tsp) / (np.pi / 2)
constraints = {'type': 'eq', 'fun': constraint}
result = minimize( # maxiter a bit ott maybe, but works.
fn, initial_guess, args=(tsp,), constraints=constraints,
bounds=[(-1., 1.), (-1., 1.), (-1., 1.)], method='SLSQP',
options={'maxiter': 1000, 'ftol': 1e-15, 'disp': False}
)
# assert result.success, result
return result
# Inverse function using numerical optimization
def inverse_ak_root(tsp):
initial_guess = np.arcsin(tsp) / (np.pi / 2)
result = root(fn_root, initial_guess, args=(tsp,), method='hybr', tol=1e-12)
assert result.success, result
result.x /= np.linalg.norm(result.x, ord=1)
return result
if __name__ == '__main__':
N = 50
np.random.seed(42)
lat = np.random.uniform(-90, 90, N)
lon = np.random.uniform(-180, 180, N)
all_points = np.column_stack([lat, lon])
for i in range(len(all_points)):
et = np.array([all_points[i]])
sph = ll_xyz(et)
result = inverse_ak(sph[0])
octal = result.x
result2 = inverse_ak_root(sph[0])
octal2 = result2.x
spherical = ak(octal)
rx = xyz_ll(np.array([spherical]))
spherical2 = ak(octal2)
rx2 = xyz_ll(np.array([spherical2]))
print(f'Old Approach')
print(f'Original Pt:{et}')
print(f'Via inverse:{rx}')
print(f'Difference:{rx-et}')
print(f'Calls: {result.nfev}')
print(f'New Approach')
print(f'Via inverse:{rx2}')
print(f'Difference:{rx2-et}')
print(f'Ratio: {np.linalg.norm(rx-et) / (np.linalg.norm(rx2-et) + 1e-15):g}')
print(f'Calls: {result2.nfev}')
print()
Note: The parameter 𝛂 has many more digits than necessary. When Python does math on this number, it ignores all digits that do not fit in a double-precision float, which means the number is converted to 3.2278062371438843 internally.
I also tried coming up with an analytic Jacobian for this function using SymPy, but had too much trouble getting automatic differentiation to work. If you were able find this, you'd be able to differentiate the function more quickly and accurately.
I also tried a solution based on scipy.interpolate.RbfInterpolator. It got similar accuracy and speed to the minimize method.