Spherical triangle sampling according to [Arvo 95]

spherical_tri.py

#!/bin/env python3
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
import math

def vnorm(v):
    return v / np.linalg.norm(v)

def draw_unit_sphere(ax):
    u = np.linspace(0, 2 * np.pi, 20)
    v = np.linspace(0, np.pi, 20)
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones(np.size(u)), np.cos(v))
    ax.plot_surface(x, y, z, rstride=2, cstride=2, alpha=0.3)

fig = plt.figure()
ax = plt.axes(projection='3d')

# Triangle
A0 = np.array([3, 1, 2])
B0 = np.array([2, 3, 1])
C0 = np.array([1, 1, 2])

# Projected points on unit sphere
A = A0 / np.linalg.norm(A0)
B = B0 / np.linalg.norm(B0)
C = C0 / np.linalg.norm(C0)

draw_unit_sphere(ax)

# Draw triangle & projected points
ax.scatter3D([A[0], A0[0]], [A[1], A0[1]], [A[2], A0[2]], label='A')
ax.scatter3D([B[0], B0[0]], [B[1], B0[1]], [B[2], B0[2]], label='B')
ax.scatter3D([C[0], C0[0]], [C[1], C0[1]], [C[2], C0[2]], label='C')
ax.plot([A0[0], 0], [A0[1], 0], [A0[2], 0])
ax.plot([B0[0], 0], [B0[1], 0], [B0[2], 0])
ax.plot([C0[0], 0], [C0[1], 0], [C0[2], 0])

#tri = mplot3d.art3d.Poly3DCollection([A0, B0, C0])
#tri.set_edgecolor('k')
#ax.add_collection3d(tri)

nA = vnorm(np.cross(B0 - A0, A0))
nB = vnorm(np.cross(C0 - B0, B0))
nC = vnorm(np.cross(A0 - C0, C0))

ax.plot(
    [A0[0], A0[0] + nA[0]],
    [A0[1], A0[1] + nA[1]],
    [A0[2], A0[2] + nA[2]])
ax.plot(
    [B0[0], B0[0] + nB[0]],
    [B0[1], B0[1] + nB[1]],
    [B0[2], B0[2] + nB[2]])
ax.plot(
    [C0[0], C0[0] + nC[0]],
    [C0[1], C0[1] + nC[1]],
    [C0[2], C0[2] + nC[2]])

dot_alpha = np.dot(nC, nA)
dot_beta  = np.dot(nA, nB)
dot_gamma = np.dot(nB, nC)

print("dot_alpha = {}".format(dot_alpha))
print("dot_beta  = {}".format(dot_beta))
print("dot_gamma = {}".format(dot_gamma))

cos_alpha = -(dot_alpha)
cos_beta  = -(dot_beta)
cos_gamma = -(dot_gamma)

print("cos(α) = {}".format(cos_alpha))
print("cos(β) = {}".format(cos_beta))
print("cos(γ) = {}".format(cos_gamma))

alpha = math.acos(cos_alpha)
beta  = math.acos(cos_beta)
gamma = math.acos(cos_gamma)

print("α + β + γ = {} + {} + {} = {}".format(alpha, beta, gamma, alpha + beta + gamma))

sin_alpha = math.sin(alpha)
sin_beta  = math.sin(beta)
cos_c = (cos_gamma + cos_beta * cos_alpha) / (sin_beta * sin_alpha)
print("cos_c = {}".format(cos_c))
spherical_area = alpha + beta + gamma - np.pi

print("spherical area: {}".format(spherical_area))

def ortho_component(x, y):
    return vnorm(x - np.dot(x, y) * y);

def sample(e1, e2):
    small_area = e1 * spherical_area
    s = math.sin(small_area - alpha)
    t = math.cos(small_area - alpha)
    u = t - cos_alpha
    v = s + sin_alpha * cos_c
    q = ((v * t - u * s) * cos_alpha - v) / ((v * s + u * t) * sin_alpha)
    c_hat = q * A + math.sqrt(1 - q * q) * ortho_component(C, A)
    z = 1 - e2 * (1 - np.dot(c_hat, B));
    return z * B + math.sqrt(1 - z * z) * ortho_component(c_hat, B)

def plot_samples(ax):
    n = 10
    x = np.empty((n + 1) * (n + 1))
    y = np.empty((n + 1) * (n + 1))
    z = np.empty((n + 1) * (n + 1))
    for i in range(0, n + 1):
        for j in range(0, n + 1):
            d = sample(i / n, j / n)
            x[i * (n + 1) + j] = d[0]
            y[i * (n + 1) + j] = d[1]
            z[i * (n + 1) + j] = d[2]
    ax.scatter3D(x, y, z)

if spherical_area > 0 and spherical_area < 2 * np.pi:
    plot_samples(ax)

plt.legend()
plt.show()