SciPy 最小化以找到反函数？

Question

我有一个（不可逆）函数 ak([u,v,w]) 这取单位八面体表面上的一个点（p：使得 |u|+|v|+|w| = 1）并返回单位球体表面上的一个点。该功能并不完美，但目的是保持点之间的距离真实。

我正在考虑使用 scipy minus 来提供数值逆，但我无法理解它。

输入：球形 pt [x,y,z]， 输出八面体点 [u,v,w] 使得 ak([u,v,w])=[x,y,z]

有什么提示吗？

Answer 1

我设法解决了这个问题，但我真的对一般正八面体<=>球体映射问题的替代解决方案感兴趣。

这是我的答案，包括

ak

功能。

import numpy as np
from scipy.optimize import minimize


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 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
    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-20, 'disp': False}
    )
    if result.success:
        return result.x
    else:
        return None


if __name__ == '__main__':
    et = np.array([[48.85837183409, 2.294492063847]])  # in the city of lights.
    sph = ll_xyz(et)
    octal = inverse_ak(sph[0])
    spherical = ak(octal)
    rx = xyz_ll(np.array([spherical]))
    print(f'Original Pt:{et},\nVia inverse:{rx},\n Difference:{rx-et}')  # within a few metres.

输出...

Original Pt:[[48.85837183  2.29449206]],
Via inverse:[[48.85837185  2.29449121]],
 Difference:[[ 1.89213694e-08 -8.54396489e-07]]