Triangle strip texturing (3 стр)

// This code is in the public domain.

// Inverse interpolation
vec2f g(const vec2f &A,const vec2f &B,const vec2f &C,const vec2f &D,const vec2f &R)
{
    // Endpoints of the segment u=const,v=[0..1] are A*(1-u)+B*u & C*(1-u)+D*u.
    // If R is on it, then cross((C-A)*(1-u)+(D-B)*u,R-(A*(1-u)+B*u))=0.
    // cross((C-A)+(A-B-C+D)*u,(R-A)+(A-B)*u)=0.

    float u,v;

    vec2f L0u=C-A;
    vec2f L1u=A-B-C+D;
    vec2f R0u=R-A;
    vec2f R1u=A-B;
    float au=cross(L1u,R1u);
    float bu=cross(L0u,R1u)+cross(L1u,R0u);
    float cu=cross(L0u,R0u);
    // Case |C-A|=0 gives rise to the false solution:
    // cross-product is zero not due to zero distance, but because
    // the segment length is zero.
    if(dot(L0u,L0u)==0.0f) u=(au!=0.0f?-bu/au:0.0f);
    // For b=0 (somewhat arbitrarily) pick the root, that gives u>=0.
    // Absent round-off, c/a<=0 should hold inside convex quadrilateral.
    else if(bu==0.0f) u=(au!=0.0f?std::sqrt(std::fmax(-cu/au,0.0f)):0.0f);
    else
    {
        float du=bu*bu-4.0f*au*cu;
        // Absent round-off, d>=0 should hold inside convex quadrilateral.
        if(du<0.0f) du=0.0f;
        // The root we need is u=(-bu+std::copysign(std::sqrt(du),bu))/(2.0f*au),
        // so that R=A gives u=0.
        // But for |a*c|<<b^2 can amplify error (up to 0/0), due to cancellation.
        // So we use the root with larger magnitude and x0*x1=a/c.
        float eu=-bu-std::copysign(std::sqrt(du),bu);
        // The denominator is only zero, when b, and a*c are both zero,
        // and this branch only executes when b!=0.
        u=2.0f*cu/eu;
    }

    // Similarly for v=const, u=[0..1].
    vec2f L0v=B-A;
    vec2f L1v=A-B-C+D;
    vec2f R0v=R-A;
    vec2f R1v=A-C;
    float av=cross(L1v,R1v);
    float bv=cross(L0v,R1v)+cross(L1v,R0v);
    float cv=cross(L0v,R0v);
    // Case |B-A|=0 gives rise to the false solution:
    // cross-product is zero not due to zero distance, but because
    // the segment length is zero.
    if(dot(L0v,L0v)==0.0f) v=(av!=0.0f?-bv/av:0.0f);
    // For b=0 (somewhat arbitrarily) pick the root, that gives v>=0.
    // Absent round-off, c/a<=0 should hold inside convex quadrilateral.
    else if(bv==0.0f) u=(av!=0.0f?std::sqrt(std::fmax(-cv/av,0.0f)):0.0f);
    else
    {
        float dv=bv*bv-4.0f*av*cv;
        // Absent round-off, d>=0 should hold inside convex quadrilateral.
        if(dv<0.0f) dv=0.0f;
        // The root we need is v=(-bv+std::copysign(std::sqrt(dv),bv))/(2.0f*av),
        // so that R=A gives v=0.
        // But for |a*c|<<b^2 can amplify error (up to 0/0), due to cancellation.
        // So we use the root with larger magnitude and x0*x1=a/c.
        float ev=-bv-std::copysign(std::sqrt(dv),bv);
        // The denominator is only zero, when b, and a*c are both zero,
        // and this branch only executes when b!=0.
        v=2.0f*cv/ev;
    }

    return {u,v};
}