import numpy as np
from scipy.interpolate import RBFInterpolator
# MacAdam (1942) ellipse data, Table III. 25 ellipses at 10-step magnification.
# Source: Wyszecki & Stiles (1982), Color Science, Table 5(5.4.1);
# digitized via the LuxPy colour science library.
# Format: (x, y, semi-major a, semi-minor b, angle in degrees)
# Note: since the Killing equation is homogeneous (L_X g = 0), the 10x
# scaling factor cancels and does not affect whether solutions exist.
macadam_data = [
(0.160, 0.057, 0.0085, 0.0035, 62.5),
(0.187, 0.118, 0.0220, 0.0055, 77.0),
(0.253, 0.125, 0.0250, 0.0050, 55.5),
(0.150, 0.680, 0.0960, 0.0230, 105.0),
(0.131, 0.521, 0.0470, 0.0200, 112.5),
(0.212, 0.550, 0.0580, 0.0230, 100.0),
(0.258, 0.450, 0.0500, 0.0200, 92.0),
(0.152, 0.365, 0.0380, 0.0190, 110.0),
(0.280, 0.385, 0.0400, 0.0150, 75.5),
(0.380, 0.498, 0.0440, 0.0120, 70.0),
(0.160, 0.200, 0.0210, 0.0095, 104.0),
(0.228, 0.250, 0.0310, 0.0090, 72.0),
(0.305, 0.323, 0.0230, 0.0090, 58.0),
(0.385, 0.393, 0.0380, 0.0160, 65.5),
(0.472, 0.399, 0.0320, 0.0140, 51.0),
(0.527, 0.350, 0.0260, 0.0130, 20.0),
(0.475, 0.300, 0.0290, 0.0110, 28.5),
(0.510, 0.236, 0.0240, 0.0120, 29.5),
(0.596, 0.283, 0.0260, 0.0130, 13.0),
(0.344, 0.284, 0.0230, 0.0090, 60.0),
(0.390, 0.237, 0.0250, 0.0100, 47.0),
(0.441, 0.198, 0.0280, 0.0095, 34.5),
(0.278, 0.223, 0.0240, 0.0055, 57.5),
(0.300, 0.163, 0.0290, 0.0060, 54.0),
(0.365, 0.153, 0.0360, 0.0095, 40.0),
]
# Convert each ellipse to metric tensor components
points, g11_vals, g12_vals, g22_vals = [], [], [], []
for (cx, cy, a, b, angle_deg) in macadam_data:
theta = np.radians(angle_deg)
c, s = np.cos(theta), np.sin(theta)
g11 = (c/a)**2 + (s/b)**2
g12 = c * s * (1/a**2 - 1/b**2)
g22 = (s/a)**2 + (c/b)**2
points.append([cx, cy])
g11_vals.append(g11)
g12_vals.append(g12)
g22_vals.append(g22)
points = np.array(points)
# Interpolate each component using thin-plate splines
interp_g11 = RBFInterpolator(points, g11_vals, kernel='thin_plate_spline', smoothing=1.0)
interp_g12 = RBFInterpolator(points, g12_vals, kernel='thin_plate_spline', smoothing=1.0)
interp_g22 = RBFInterpolator(points, g22_vals, kernel='thin_plate_spline', smoothing=1.0)
def get_metric(x, y):
pt = np.array([[x, y]])
return np.array([[interp_g11(pt)[0], interp_g12(pt)[0]],
[interp_g12(pt)[0], interp_g22(pt)[0]]])