Ambient Dice

AmbientDice.cpp

struct Vertex {
    vec3 value;
    vec3 directionalDerivativeU;
    vec3 directionalDerivativeV;
};

void AmbientDice::hybridCubicBezierWeights(u32 i0, u32 i1, u32 i2, float b0, float b1, float b2, VertexWeights *w0Out, VertexWeights *w1Out, VertexWeights *w2Out) const {
        const float alpha = 0.5f * sqrt(0.5f * (5.0f + sqrt(5.0f))); // 0.9510565163
        const float beta = -0.5f * sqrt(0.1f * (5.0f + sqrt(5.0f))); // -0.4253254042
        
        const float a0 = (sqrt(5.0f) - 5.0f) / 40.0f; // -0.06909830056
        const float a1 = (11.0f * sqrt(5.0f) - 15.0f) / 40.0f; // 0.2399186938
        const float a2 = sqrt(5.0f) / 10.0f; // 0.2236067977
        
        // Project the edges onto the sphere.
        vec3 v0V1 = AmbientDice::vertexPositions[i1] - AmbientDice::vertexPositions[i0];
        v0V1 = v0V1 - normalize(AmbientDice::vertexPositions[i0]) * dot(normalize(AmbientDice::vertexPositions[i0]), v0V1);
        v0V1 = normalize(v0V1);
        
        vec3 v1V0 = AmbientDice::vertexPositions[i0] - AmbientDice::vertexPositions[i1];
        v1V0 = v1V0 - normalize(AmbientDice::vertexPositions[i1]) * dot(normalize(AmbientDice::vertexPositions[i1]), v1V0);
        v1V0 = normalize(v1V0);
        
        vec3 v0V2 = AmbientDice::vertexPositions[i2] - AmbientDice::vertexPositions[i0];
        v0V2 = v0V2 - normalize(AmbientDice::vertexPositions[i0]) * dot(normalize(AmbientDice::vertexPositions[i0]), v0V2);
        v0V2 = normalize(v0V2);
        
        vec3 v2V0 = AmbientDice::vertexPositions[i0] - AmbientDice::vertexPositions[i2];
        v2V0 = v2V0 - normalize(AmbientDice::vertexPositions[i2]) * dot(normalize(AmbientDice::vertexPositions[i2]), v2V0);
        v2V0 = normalize(v2V0);
        
        vec3 v1V2 = AmbientDice::vertexPositions[i2] - AmbientDice::vertexPositions[i1];
        v1V2 = v1V2 - normalize(AmbientDice::vertexPositions[i1]) * dot(normalize(AmbientDice::vertexPositions[i1]), v1V2);
        v1V2 = normalize(v1V2);
        
        vec3 v2V1 = AmbientDice::vertexPositions[i1] - AmbientDice::vertexPositions[i2];
        v2V1 = v2V1 - normalize(AmbientDice::vertexPositions[i2]) * dot(normalize(AmbientDice::vertexPositions[i2]), v2V1);
        v2V1 = normalize(v2V1);
        
        const float fValueFactor = -beta / alpha; // 0.4472135955
        const float fDerivativeFactor = 1.0 / (3.0 * alpha); // 0.3504874081
        
        const float weightDenom = b1 * b2 + b0 * b2 + b0 * b1;
        
        float w0 = (b1 * b2) / weightDenom;
        float w1 = (b0 * b2) / weightDenom;
        float w2 = (b0 * b1) / weightDenom;
        
        if (b0 == 1.0) {
            w0 = 1.0;
            w1 = 0.0;
            w2 = 0.0;
        } else if (b1 == 1.0) {
            w0 = 0.0;
            w1 = 1.0;
            w2 = 0.0;
        } else if (b2 == 1.0) {
            w0 = 0.0;
            w1 = 0.0;
            w2 = 1.0;
        }
        
        // https://en.wikipedia.org/wiki/Bézier_triangle
        // Notation: cxyz means alpha^x, beta^y, gamma^z.
        
        float v0ValueWeight = 0.0;
        float v1ValueWeight = 0.0;
        float v2ValueWeight = 0.0;
        
        float v0DUWeight = 0.0;
        float v1DUWeight = 0.0;
        float v2DUWeight = 0.0;
        
        float v0DVWeight = 0.0;
        float v1DVWeight = 0.0;
        float v2DVWeight = 0.0;
        
        const float b0_2 = b0 * b0;
        const float b1_2 = b1 * b1;
        const float b2_2 = b2 * b2;
        
        // Add c300, c030, and c003
        float c300Weight = b0_2 * b0;
        float c030Weight = b1_2 * b1;
        float c003Weight = b2_2 * b2;
        
        float c120Weight = 3 * b0 * b1_2;
        float c021Weight = 3 * b1_2 * b2;
        float c210Weight = 3 * b0_2 * b1;
        float c012Weight = 3 * b1 * b2_2;
        float c201Weight = 3 * b0_2 * b2;
        float c102Weight = 3 * b0 * b2_2;
        
        const float c111Weight = 6 * b0 * b1 * b2;
        const float c0_111Weight = w0 * c111Weight;
        const float c1_111Weight = w1 * c111Weight;
        const float c2_111Weight = w2 * c111Weight;
        
        v1ValueWeight += a0 * c0_111Weight;
        v2ValueWeight += a0 * c1_111Weight;
        v0ValueWeight += a0 * c2_111Weight;
        
        c021Weight += a1 * c0_111Weight;
        c012Weight += a1 * c0_111Weight;
        c003Weight += a0 * c0_111Weight;
        c120Weight += a2 * c0_111Weight;
        c102Weight += a2 * c0_111Weight;
        
        c102Weight += a1 * c1_111Weight;
        c201Weight += a1 * c1_111Weight;
        c300Weight += a0 * c1_111Weight;
        c012Weight += a2 * c1_111Weight;
        c210Weight += a2 * c1_111Weight;
        
        c210Weight += a1 * c2_111Weight;
        c120Weight += a1 * c2_111Weight;
        c030Weight += a0 * c2_111Weight;
        c201Weight += a2 * c2_111Weight;
        c021Weight += a2 * c2_111Weight;
        
        
        v0ValueWeight += fValueFactor * c210Weight;
        v0DUWeight += fDerivativeFactor * dot(v0V1, AmbientDice::tangents[i0]) * c210Weight;
        v0DVWeight += fDerivativeFactor * dot(v0V1, AmbientDice::bitangents[i0]) * c210Weight;
        
        v0ValueWeight += fValueFactor * c201Weight;
        v0DUWeight += fDerivativeFactor * dot(v0V2, AmbientDice::tangents[i0]) * c201Weight;
        v0DVWeight += fDerivativeFactor * dot(v0V2, AmbientDice::bitangents[i0]) * c201Weight;
        
        v1ValueWeight += fValueFactor * c120Weight;
        v1DUWeight += fDerivativeFactor * dot(v1V0, AmbientDice::tangents[i1]) * c120Weight;
        v1DVWeight += fDerivativeFactor * dot(v1V0, AmbientDice::bitangents[i1]) * c120Weight;
        
        v1ValueWeight += fValueFactor * c021Weight;
        v1DUWeight += fDerivativeFactor * dot(v1V2, AmbientDice::tangents[i1]) * c021Weight;
        v1DVWeight += fDerivativeFactor * dot(v1V2, AmbientDice::bitangents[i1]) * c021Weight;
        
        v2ValueWeight += fValueFactor * c102Weight;
        v2DUWeight += fDerivativeFactor * dot(v2V0, AmbientDice::tangents[i2]) * c102Weight;
        v2DVWeight += fDerivativeFactor * dot(v2V0, AmbientDice::bitangents[i2]) * c102Weight;
        
        v2ValueWeight += fValueFactor * c012Weight;
        v2DUWeight += fDerivativeFactor * dot(v2V1, AmbientDice::tangents[i2]) * c012Weight;
        v2DVWeight += fDerivativeFactor * dot(v2V1, AmbientDice::bitangents[i2]) * c012Weight;
        
        v0ValueWeight += c300Weight;
        v1ValueWeight += c030Weight;
        v2ValueWeight += c003Weight;
        
//        v0ValueWeight = max(0.f, v0ValueWeight);
//        v1ValueWeight = max(0.f, v1ValueWeight);
//        v2ValueWeight = max(0.f, v2ValueWeight);
//
//        v0DUWeight = max(0.f, v0DUWeight);
//        v0DVWeight = max(0.f, v0DVWeight);
//
//        v1DUWeight = max(0.f, v1DUWeight);
//        v1DVWeight = max(0.f, v1DVWeight);
//
//        v2DUWeight = max(0.f, v2DUWeight);
//        v2DVWeight = max(0.f, v2DVWeight);
        
        assert(v0ValueWeight >= 0.f);
        assert(v1ValueWeight >= 0.f);
        assert(v2ValueWeight >= 0.f);
        
//        assert(v0DUWeight >= 0.f);
//        assert(v1DUWeight >= 0.f);
//        assert(v2DUWeight >= 0.f);
//
//        assert(v0DVWeight >= 0.f);
//        assert(v1DVWeight >= 0.f);
//        assert(v2DVWeight >= 0.f);
        
        
        *w0Out = { v0ValueWeight, v0DUWeight, v0DVWeight };
        *w1Out = { v1ValueWeight, v1DUWeight, v1DVWeight };
        *w2Out = { v2ValueWeight, v2DUWeight, v2DVWeight };
    }
    
    void AmbientDice::hybridCubicBezierWeights(vec3 direction, u32 *i0Out, u32 *i1Out, u32 *i2Out, VertexWeights *w0Out, VertexWeights *w1Out, VertexWeights *w2Out) const {
        
        u32 i0, i1, i2;
        float b0, b1, b2;
        this->computeBarycentrics(direction, &i0, &i1, &i2, &b0, &b1, &b2);
        
        this->hybridCubicBezierWeights(i0, i1, i2, b0, b1, b2, w0Out, w1Out, w2Out);
        
        *i0Out = i0;
        *i1Out = i1;
        *i2Out = i2;
    }