another curvature gradient based on trig function

trig_curvature_gradient.h

void computeGradients(const Vector3& pa, const Vector3& pb, const Vector3& pc, const Vector3& pd,
                      Vector3& grad_pa, Vector3& grad_pb, Vector3& grad_pc, Vector3& grad_pd) {
    // Step 1: Compute edge vectors using operator-
    Vector3 eca = pc - pa;
    Vector3 eda = pd - pa;
    Vector3 edb = pd - pb;
    Vector3 ecb = pc - pb;
    Vector3 pd_pc = pd - pc;
    double L = pd_pc.Norm();
    Vector3 ecd = (L < 1e-8) ? Vector3(0.0, 0.0, 0.0) : (pd_pc * (1.0 / L));

    // Step 2: Compute unnormalized normals using Cross
    Vector3 N1_un = Cross(eca, eda);
    Vector3 N2_un = Cross(edb, ecb);

    // Step 3: Compute magnitudes
    double N1_mag = N1_un.Norm();
    double N2_mag = N2_un.Norm();

    // Step 4: Handle degenerate cases
    if (N1_mag < 1e-8 || N2_mag < 1e-8 || L < 1e-8) {
        grad_pa = Vector3(0.0, 0.0, 0.0);
        grad_pb = Vector3(0.0, 0.0, 0.0);
        grad_pc = Vector3(0.0, 0.0, 0.0);
        grad_pd = Vector3(0.0, 0.0, 0.0);
        return;
    }

    // Step 5: Normalize normals using operator*
    Vector3 N1 = N1_un * (1.0 / N1_mag);
    Vector3 N2 = N2_un * (1.0 / N2_mag);

    // Step 6: Compute N_sum using operator+
    Vector3 N_sum = N1 + N2;

    // Step 7: Compute dT using Cross
    Vector3 dT = Cross(N_sum, ecd);

    // Step 8: Compute dT2 with regularization using Dot
    double dT2 = Dot(dT, dT) + 1e-8;

    // Step 9: Compute dN using operator-
    Vector3 dN = N2 - N1;

    // Step 10: Compute kappa using Dot
    double kappa = Dot(dN, dT) / dT2;

    // Step 11: Compute gradients with respect to N1 and N2
    Vector3 cross_ecd_dN = Cross(ecd, dN);
    Vector3 cross_dT_ecd = Cross(dT, ecd);

    Vector3 term1_N1 = dT * (-dT2);
    Vector3 term2_N1 = cross_ecd_dN * dT2;
    Vector3 term3_N1 = cross_dT_ecd * (-2.0 * kappa);
    Vector3 sum_terms_N1 = term1_N1 + term2_N1 + term3_N1;
    Vector3 grad_N1_kappa = sum_terms_N1 * (1.0 / (dT2 * dT2));

    Vector3 term1_N2 = dT * dT2;
    Vector3 term2_N2 = cross_ecd_dN * dT2;
    Vector3 term3_N2 = cross_dT_ecd * (-2.0 * kappa);
    Vector3 sum_terms_N2 = term1_N2 + term2_N2 + term3_N2;
    Vector3 grad_N2_kappa = sum_terms_N2 * (1.0 / (dT2 * dT2));

    // Step 12: Compute gradients with respect to points

    // grad_pa
    Vector3 ecd_cross_grad_N1 = Cross(ecd, grad_N1_kappa);
    Vector3 ecd_cross_N1 = Cross(ecd, N1);
    Vector3 temp_pa = ecd_cross_grad_N1 * (-1.0) + ecd_cross_N1 * Dot(N1, grad_N1_kappa);
    grad_pa = temp_pa * (L / N1_mag);

    // grad_pb
    Vector3 ecd_cross_grad_N2 = Cross(ecd, grad_N2_kappa);
    Vector3 ecd_cross_N2 = Cross(ecd, N2);
    Vector3 temp_pb = ecd_cross_grad_N2 * (-1.0) + ecd_cross_N2 * Dot(N2, grad_N2_kappa);
    grad_pb = temp_pb * (L / N2_mag);

    // For grad_pc and grad_pd, compute partial derivatives
    Vector3 grad_N1_pc = Cross(eda, ecd) * (L / N1_mag); // ∂N1/∂pc
    Vector3 grad_N2_pc = Cross(edb, ecd) * (-L / N2_mag); // ∂N2/∂pc
    Vector3 grad_ecd_pc = (Vector3(0.0, 0.0, 0.0) - ecd) * (1.0 / L); // ∂ecd/∂pc

    Vector3 grad_N1_pd = Cross(ecd, eca) * (L / N1_mag); // ∂N1/∂pd
    Vector3 grad_N2_pd = Cross(ecd, ecb) * (L / N2_mag); // ∂N2/∂pd
    Vector3 grad_ecd_pd = ecd * (1.0 / L); // ∂ecd/∂pd

    // Chain rule for grad_pc and grad_pd (simplified)
    Vector3 dkappa_decd = Cross(N_sum, dN) * (1.0 / dT2);
    grad_pc = grad_N1_pc * Dot(grad_N1_kappa, grad_N1_pc) +
              grad_N2_pc * Dot(grad_N2_kappa, grad_N2_pc) +
              grad_ecd_pc * Dot(dkappa_decd, grad_ecd_pc);

    grad_pd = grad_N1_pd * Dot(grad_N1_kappa, grad_N1_pd) +
              grad_N2_pd * Dot(grad_N2_kappa, grad_N2_pd) +
              grad_ecd_pd * Dot(dkappa_decd, grad_ecd_pd);
}