import numpy as np

def slopes_to_n(s):
    return np.array([-s[0], -s[1], 1]) / np.sqrt(s[0]**2 + s[1]**2 + 1);

def edge_integral(v1, v2):
    theta = np.arccos(np.dot(v1, v2));
    factor = 1 if theta < 0.0001 else theta/np.sin(theta);
    res = np.cross(v1, v2)[2] * factor;
    return res;

def triangle_integral_ggx(p0, p1, p2, alpha_x, alpha_y):

    # Linear transformation of the space of slopes which
    # reduces the integral to the case alpha_u = alpha_v = 1
    scale = np.array([1/alpha_x, 1/alpha_y]);

    # Coordinates of the corresponding spherical triangle
    n0 = slopes_to_n(p0 * scale);
    n1 = slopes_to_n(p1 * scale);
    n2 = slopes_to_n(p2 * scale);

    # Integral of a cosine weighted spherical triangle
    integral = edge_integral(n0, n1) \
             + edge_integral(n1, n2) \
             + edge_integral(n2, n0);

    return np.abs(integral) / (2 * np.pi);

def ggx_distribution(p, alpha_x , alpha_y):
    return 1/(np.pi * alpha_x * alpha_y * ((p[0]/alpha_x)**2 + (p[1]/alpha_y)**2 + 1)**2)

def triangle_area(p0, p1, p2):
    return np.linalg.norm(np.cross(p1-p0, p2-p0))*0.5

def triangle_integral_ggx_mc(p0, p1, p2, alpha_x , alpha_y):
    integral = 0;
    NB_SAMPLES = 100000
    for i in range(NB_SAMPLES):
        uv = np.random.rand(2)
        if sum(uv) > 1:
            uv = 1-uv
        pi = p0 + uv[0]*(p1-p0) + uv[1]*(p2-p0)
        integral += ggx_distribution(pi, alpha_x , alpha_y) / NB_SAMPLES
    return integral * triangle_area(p0, p1, p2)

# Test: compare results with Monte Carlo integration

# Define a triangle in slope space
p0 = np.array([  0.0,  0.1])
p1 = np.array([ -0.1,  0.5])
p2 = np.array([  0.2, -0.6])

# GGX roughness in X and Y axes
alpha_x = 0.01
alpha_y = 0.5

print("Integral, analytic    : ", triangle_integral_ggx(p0, p1, p2, alpha_x , alpha_y))
print("Estimating the integral with Monte Carlo integration...")
print("Integral, Monte Carlo : ", triangle_integral_ggx_mc(p0, p1, p2, alpha_x , alpha_y))