intersection.cpp

vec4 Intersection::newtonStep(const glm::vec4& params){
    // The last point found on the intersection curve
    const TracePoint& lastPoint = getLastPoint();
    vec4 lastParam = lastPoint.params;
    
    // Determines the direction
    int r = 1;
    if(currentTraceStatus == TraceStatus::BACKWARDS)
        r = -1;
    
    vec4 result;
    mat4 J;
    vec4 F;

    vec3 Np;
    vec3 Nq;
    vec3 t;

    vec3 P0 = surface1->compute(lastParam.x, lastParam.y);
    vec3 Pu0 = surface1->computeDu(lastParam.x, lastParam.y);
    vec3 Pv0 = surface1->computeDv(lastParam.x, lastParam.y);
    vec3 Qu0  = surface2->computeDu(lastParam.z, lastParam.w);
    vec3 Qv0  = surface2->computeDv(lastParam.z, lastParam.w);

    vec3 P = surface1->compute(params.x, params.y);
    vec3 Q = surface2->compute(params.z, params.w);
    vec3 Pu  = surface1->computeDu(params.x, params.y);
    vec3 Pv  = surface1->computeDv(params.x, params.y);
    vec3 Qu  = surface2->computeDu(params.z, params.w);
    vec3 Qv  = surface2->computeDv(params.z, params.w);

    Np = cross(Pu0, Pv0);
    Nq = cross(Qu0, Qv0);

    t = normalize(cross(Np, Nq));
    t *= r;

    vec3 dP = P-P0;

    auto f4 = [this, P0, t](float x, float y,
                            float z, float w){
        vec3 _P = surface1->compute(x, y);
        vec3 _dP = _P-P0;

        float dotValue = ifc::dot(_dP, t);
        auto value = dotValue - distance;

        return value;
    };

    auto F4du = ifc::derivative(f4, params.x, params.y, params.z, params.w,
                                ifc::DerivativeTypes::DX);
    auto F4dv = ifc::derivative(f4, params.x, params.y, params.z, params.w,
                                ifc::DerivativeTypes::DY);
    auto F4ds = ifc::derivative(f4, params.x, params.y, params.z, params.w,
                                ifc::DerivativeTypes::DZ);
    auto F4dt = ifc::derivative(f4, params.x, params.y, params.z, params.w,
                                ifc::DerivativeTypes::DW);
    
    F.x = r*(P.x - Q.x);
    F.y = r*(P.y - Q.y);
    F.z = r*(P.z - Q.z);
    F.w = ifc::dot(dP, t) - distance;

    // zeros left for clarity
    J[0].x = r*Pu.x;
    J[0].y = r*Pu.y;
    J[0].z = r*Pu.z;
    J[0].w = F4du;

    J[1].x = r*Pv.x;
    J[1].y = r*Pv.y;
    J[1].z = r*Pv.z;
    J[1].w = F4dv;

    J[2].x = r*(-Qu.x);
    J[2].y = r*(-Qu.y);
    J[2].z = r*(-Qu.z);
    J[2].w = F4ds;

    J[3].x = r*(-Qv.x);
    J[3].y = r*(-Qv.y);
    J[3].z = r*(-Qv.z);
    J[3].w = F4dt;

    J = glm::inverse(J);
    result = params - J*F;
    
    return result;
}