float alpha_I(in vec3 wi, in vec2 alpha) {
    float inverseSinTheta2 = 1.0 / (1.0 - square(wi.z));
    float cosPhi2 = square(wi.x) * inverseSinTheta2;
    float sinPhi2 = square(wi.y) * inverseSinTheta2;
    float alpha_i = sqrt(cosPhi2 * square(alpha.x) + sinPhi2 * square(alpha.y));
    return alpha_i;
}

float D(in vec3 wm, in vec2 alpha) {
    if(wm.z <= 0.0) {
        return 0.0;
    }
    vec2 slope = -vec2(wm.x, wm.y) / wm.z;

    float value = P22(slope, alpha);
    return value;
}

float D_wi(in vec3 wi, in vec3 wm, in vec2 alpha) {
    if(wm.z <= 0.0) {
        return 0.0;
    }
    float pa = projectedArea(wi, alpha);
    if(pa == 0.0) {
        return 0.0;
    }
    float c = 1.0 / pa;

    float value = c * max(0.0, dot(wi, wm)) * D(wm, alpha);
    return value;
}

vec3 sampleD_wi(in vec3 wi, in vec2 U, in vec2 alpha) {
    vec3 wi_11 = normalize(vec3(alpha.x * wi.x, alpha.y * wi.y, wi.z));

    vec2 slope_11 = sampleP22_11(acos(wi_11.z), U);

    float phi = atan(wi_11.y, wi_11.x);
    vec2 slope = vec2(cos(phi) * slope_11.x - sin(phi) * slope_11.y, sin(phi) * slope_11.x + cos(phi) * slope_11.y);
    slope *= alpha;

    if(slope.x != slope.x || isinf(slope.x)) {
        if(wi.z > 0.0) {
            return vec3(0.0, 0.0, 1.0);
        } else {
            return normalize(vec3(wi.x, wi.y, 0.0));
        }
    }

    vec3 wm = normalize(vec3(-slope.x, -slope.y, 1.0));

    return wm;
}