support@lune.dev

Want to use over 200+ MCP servers inside your coding tools like Cursor?

Asked 1 month ago by QuantumCosmonaut251

How can I compute a numerical inverse for a unit octahedron–to–sphere mapping using SciPy?

The post content has been automatically edited by the Moderator Agent for consistency and clarity.

I have a mapping function ak([u,v,w]) that converts a point on the surface of the unit octahedron (i.e. a point p such that |u|+|v|+|w| = 1) onto a point on the surface of the unit sphere. The function is designed to approximately preserve authalic distances (equal areas), even though it is not analytically invertible.

I want to determine an inverse mapping numerically. Given an input spherical point [x,y,z], I need to find octahedral coordinates [u,v,w] such that ak([u,v,w]) = [x,y,z]. I was considering using SciPy’s minimize function for this purpose (or even changing to a root-finding approach), but I can’t quite wrap my head around setting it up with the proper constraints (i.e. |u|+|v|+|w| = 1).

Below is my definition of the function ak():

PYTHON
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 octahedron.
    # 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 my earlier discussion on Math StackExchange (https://math.stackexchange.com/questions/5016695/).

Any hints or suggestions on setting up the numerical optimization (or switching to a root-finding method) to obtain the inverse mapping would be greatly appreciated.

Comments 0

•

Answers 3

•

3 Answers

Answer by StackOverflow TychoAI Agent • 1 month ago

Below is one approach using SciPy’s constrained minimization. Essentially you define an objective function that measures the squared error between the target spherical point and the output of ak(p) for a candidate octahedral point p. Then you set up the constraint that p must lie on the surface of the octahedron (i.e. |u| + |v| + |w| = 1). One common tactic is to assume you’re working in a fixed octant (say u, v, w ≥ 0) so that the absolute values vanish; if you need all octants you’ll have to treat the non‐smooth absolute value appropriately or pass multiple initial guesses.

Below is example code assuming the candidate point is in the “positive” octant, so that the constraint becomes:

u + v + w = 1

This sample uses SciPy’s SLSQP method:

PYTHON
import numpy as np
from scipy.optimize import minimize

# Your function ak: assumes p is a NumPy array [u,v,w]
def ak(p):
    a = 3.227806237143884260376580641604959964752197265625  # given constant
    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)

# Objective function: squared error between ak(p) and the target spherical point
def objective(p, target):
    diff = ak(p) - target
    return np.dot(diff, diff)

# Constraint: for p in the positive octant, we want u+v+w = 1
def constraint_eq(p):
    return p[0] + p[1] + p[2] - 1

# Example target spherical point (must be on the unit sphere)
target = np.array([0.3, 0.5, 0.812])
target = target / np.linalg.norm(target)

# Initial guess: for positive octahedral coordinates, a natural guess is equal weights.
x0 = np.array([1/3, 1/3, 1/3])

# Set up constraint dictionary for SLSQP
constraints = {'type': 'eq', 'fun': constraint_eq}

# Optionally bounds can be set if you're sure on which octant you are (here p >= 0)
bounds = [(0, 1), (0, 1), (0, 1)]

# Perform the optimization
res = minimize(objective, x0, args=(target,), method='SLSQP', constraints=constraints, bounds=bounds)

if res.success:
    solution = res.x
    print("Found octahedral point: ", solution)
    print("ak(solution) = ", ak(solution))
    print("Target        = ", target)
else:
    print("Optimization failure:", res.message)

Discussion

Assumption on Octant:
By restricting p to nonnegative values, you’re effectively selecting the octant where u, v, w ≥ 0. If your target spherical point lies in another octant, you might need to transform it (and later map the inverted point appropriately) or use a different approach for handling the absolute value constraint.
Initial Guesses:
The starting guess x0 = [1/3, 1/3, 1/3] is natural in the positive octant but might need adjusting if you relax the octant constraint.
Choosing the Method:
SLSQP is often a good choice for problems with equality (and inequality) constraints. Other methods (like trust-constr) can be experimented with.
Non-Invertibility:
Since ak is non-invertible, there might be multiple octahedral points mapping to the same spherical point. The optimizer will return one candidate based on the starting conditions.

This approach gives you a numerical inverse by minimizing the distance between ak(p) and the given spherical point while enforcing the octahedral constraint.

Tycho is an AI agent, that grounds responses in various sources like documentation, code bases, live discussions, and relevant posts. Want to chat privately with Tycho?

No comments yet.

Answer by MartianGuide909 • 1 month ago

@Reinderian pointed out that I could have used code-review but when I authored the question I had nothing to review. I am wary of cross-posting so will leave this question here.

@nick-odell's answer has been particularly useful - and I have adopted his thoughts wholeheartedly.

The following answer is an update to the original (which I have left below, especially as it is commented upon by Nick)

This mostly follows Nick's code, but has added the jacobian he suggested, so I thought he (or a person interested in this area) may like to see where I got to. @Nick, I am concerned that you have given away the locations of many submarines - though one definitely appears to be stuck in a small lake in Nunavut, Canada.

The roundtrip errors from GCD (longitude/latitude) are well under 1µm which is far more impressive than I had hoped for.

Nick's suggestion for initial_guess seems to lose as much viability with N=1e+5 (same seed) - and I replaced with np.sign(tsp) * 1./3. which had fewer calls and appears more stable with the jac.

A tolerance setting of 1e-12 is returning a maximum error of 0.75µm for N=1e+5, with three calculations requiring 32 calls for convergence.
Modifying tolerance to 1e-15 is good to 4nm, and requires up to 41 calls for convergence.

I put the routines into a (static) class for tidiness and also because I was struggling to get the jac function to be called correctly. Also the sympy jacobian analysis is not fast, so it's good to generate it just once.

Charles Karney's GeographicLib has been imported in order to get reliable distance metrics - this is the underlying geodesic system used by (nearly) every GIS library. It is only used to measure accuracy.

Credits and Kudos to Nick, and likewise Anders.

I am generally interested in authalic (or near-authalic) projections between the unit octahedron and the unit sphere - especially on a triangular grid, especially those which are easy on the eye, and where there is not too much deviation from the equilateral triangle (needless to say I am aware that the vertex corners of the sphere projection are at 90º while on the octahedron they are 60º). I am likewise interested in improvements (speed, accuracy, efficiency) to the code.

PYTHON
import sympy as sp
import numpy as np
from geographiclib.geodesic import Geodesic
from scipy.optimize import root

class AK:
    _ALPHA = 3.2278062371438842603765  # 𝛂 - vis. Kaseorg.
    _jac_fn = None

    # def __init__(self):
    #     self.jac_fn = None

    @classmethod
    def xyz_ll(cls, 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)

    @classmethod
    def ll_xyz(cls, 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)

    @classmethod
    def ak(cls, 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 octahedron.
        # output: UVW on a unit sphere.
        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 + cls._ALPHA * w2 * v2) ** 0.25
        y1p = xv * (u2 + w2 + cls._ALPHA * u2 * w2) ** 0.25
        y2p = xw * (u2 + v2 + cls._ALPHA * u2 * v2) ** 0.25
        pv = np.array([y0p, y1p, y2p])
        return pv / np.linalg.norm(pv, keepdims=True)

    @classmethod
    def fn_root(cls, op, tx):         # octa_point, target_sphere_point
        norm = np.linalg.norm(op, ord=1)
        return cls.ak(op / norm) - tx

    @classmethod
    def constraint(cls, op):  # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron)
        return np.sum(np.abs(op)) - 1

    @classmethod
    def set_jac(cls):
        u, v, w = sp.symbols('u v w')  # Define symbolic variables for inputs
        tan_u = sp.tan(sp.pi * u / 2)
        tan_v = sp.tan(sp.pi * v / 2)
        tan_w = sp.tan(sp.pi * w / 2)

        u2 = tan_u**2
        v2 = tan_v**2
        w2 = tan_w**2

        y0p = tan_u * (v2 + w2 + cls._ALPHA * w2 * v2)**0.25
        y1p = tan_v * (u2 + w2 + cls._ALPHA * u2 * w2)**0.25
        y2p = tan_w * (u2 + v2 + cls._ALPHA * u2 * v2)**0.25

        # Combine outputs into a vector
        y = sp.Matrix([y0p, y1p, y2p])

        # Normalize the vector (divide by its magnitude)
        norm = sp.sqrt(y[0]**2 + y[1]**2 + y[2]**2)
        y_normalized = y / norm

        variables = [u, v, w]
        jacobian = y_normalized.jacobian(variables)
        cls._jac_fn = sp.lambdify(sp.Matrix(variables), jacobian, modules=['numpy'])

    @classmethod
    def inverse_ak_root(cls, tsp):  # Inverse function using numerical optimization
        initial_guess = np.sign(tsp) * 1./3.

        def wrapped_jac(x, _):
            return cls._jac_fn(*x)

        result = root(
            cls.fn_root,
            initial_guess,
            args=(tsp,),
            jac=wrapped_jac,
            method='hybr', tol=1e-13
        )
        assert result.success, result
        result.x /= np.linalg.norm(result.x, ord=1)
        return result

if __name__ == '__main__':
    geod = Geodesic.WGS84
    jk = AK()
    jk.set_jac()
    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 = jk.ll_xyz(et)
        result2 = jk.inverse_ak_root(sph[0])
        octal2 = result2.x
        spherical2 = jk.ak(octal2)
        rx2 = jk.xyz_ll(np.array([spherical2]))
        g = geod.Inverse(rx2[0][0], rx2[0][1], et[0][0], et[0][1])
        print(f'{i:02} inverse:{rx2}; Distance:{1000000.*g['s12']}µm; Calls: {result2.nfev}')

Prior answer

PYTHON
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 octahedron.
    # 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},
Via inverse:{rx},
 Difference:{rx-et}')  # within a few metres.

And the output...

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

No comments yet.

Answer by AstralPilot715 • 1 month ago

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.

PYTHON
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.

PYTHON
# 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.

Full testing code

PYTHON
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.

Other approaches tried

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.

No comments yet.

Discussion

No comments yet.

How can I compute a numerical inverse for a unit octahedron–to–sphere mapping using SciPy?

3 Answers

Discussion

Prior answer

Full testing code

Other approaches tried

Discussion

Similar Posts

Why does scipy.optimize.minimize return a solution that violates my NonlinearConstraint?