Error for calculating point edge distance?

[Yanksi]

In MeshCollisionUtils.cpp, the point edge distance was calculated as following.

inline void d_PE(const Eigen::RowVector3d& v0,
    const Eigen::RowVector3d& v1,
    const Eigen::RowVector3d& v2,
    double& d)
{
    d = (v1 - v0).cross(v2 - v0).squaredNorm() / (v2 - v1).squaredNorm();
}

I think this is actually calculating the squared distance between a point and edge.

[liminchen]

Yeah, for numerical robustness all distances for the barrier input in IPC are squared distances so that square root is avoided. When computing normal forces for friction when unit matters, there are cares taken there to obtain the correct value.

[Yanksi]

Then how did the gradient and hessian of the barrier energy get computed? Compute it in the form of the squared distance and apply the chain rule after it? When I am trying to derive the gradient and hessian of the distance myself, Matlab constantly giving me a dirac function in the C code generated by it that I don't know how to get rid of.

[Yanksi]

And also, I want to provide a slightly optimized way of calculating the edge-edge distance. I noticed that you want to make the current code more efficient for that in compute_graphics_edge_edge_constraint in this. I wrote a little benchmark between these two methods.

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <vector>
#include <ctime>

void mysolver(Eigen::Ref<Eigen::Vector2d> result, const Eigen::Matrix<double, 3,4> a) {
    Eigen::Matrix<double, 3, 2> mat;
    mat << a.col(1) - a.col(0), a.col(2) - a.col(3);
    Eigen::Vector3d b = a.col(2) - a.col(0);
    result = (mat.transpose() * mat).ldlt().solve(mat.transpose() * b);
}

void ipcsolver(Eigen::Ref<Eigen::Vector2d> result, const Eigen::Matrix<double, 3,4> a) {
    Eigen::Matrix3d mat;
    Eigen::Vector3d dir0, dir1, dir2;
    dir0 = a.col(1) - a.col(0);
    dir1 = a.col(3) - a.col(2);
    dir2 = dir1.cross(dir0);
    mat << dir0, -dir1, dir2;
    Eigen::Vector3d b = a.col(2) - a.col(0);
    result = mat.lu().solve(b).head<2>();
}

int main() {
    using namespace Eigen;
    int size = 100000;
    int round = 100;
    std::vector<Eigen::Matrix<double, 3, 4>> as;
    Eigen::Array2Xd myresult(2, size);
    Eigen::Array2Xd ipcresult(2, size);
    as.reserve(size);

    for (int i = 0; i < size; ++i) {
        as.emplace_back(Eigen::Matrix<double, 3, 4>::Random());
    }
    std::clock_t start;
    double duration;
    start = std::clock();
    for (int r = 0; r < round; ++r) {
        for (int i = 0; i < size; ++i) {
            mysolver(myresult.col(i), as[i]);
        }
    }
    duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC / round;
    std::cout << "my:" << duration <<'\n';

    start = std::clock();
    for (int r = 0; r < round; ++r) {
        for (int i = 0; i < size; ++i) {
            ipcsolver(ipcresult.col(i), as[i]);
        }
    }
    duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC / round;
    std::cout << "ipc:" << duration <<'\n';

    std::cout << "err:" << ((myresult - ipcresult) / ipcresult).abs().rowwise().maxCoeff();
}

Similarly, the point-triangle distance can also be derived in a similar form, but I am not sure whether the performance will go up as I haven't written benchmark for that.

[liminchen]

Gradient and hessian are also just computed of the squared distance. Yes unsquared distance' gradient and hessian might have dirac() from matlab because of the square root. The optimization looks great! But compute_graphics_edge_edge_constraint is actually only used for comparison to other methods, not in IPC method.

[zfergus]

The gradient and hessian are also computed for squared distances to avoid the division by distance (i.e., avoid d/dx sqrt(distance^2(x)) = 1 / sqrt(distance^2(x)) * distance(x) * distance'(x)).

The code inside src/CollisionObject/CollisionConstraints.cpp is used for comparison with traditional SQP methods. Our edge-edge distance computation can be found in src/CollisionObject/MeshCollisionUtils.hpp.

[Yanksi]

So the barrier energy b(q) is also computed using the squared distance? I tried to derive the formula for dd/dq and d2d/dq2 when we only have dd2/dq and d2d2/dq2 on hand. But I am not sure whether I get them correct. Could you please check the formula for me?

dd2/dq = dd2/dd * dd/dq 
=> dd/dq = (dd2/dq)/2d

d2d/dq2 = (d2d2/dq2)/2d

[zfergus]

The barrier is computed using squared distance. You might also want to take a look at section 10 in our technical supplementary material (https://ipc-sim.github.io/file/IPC-supplement-A-technical.pdf). As for your derivation, I find it hard to read. By dd/dq do you mean d/dq distance(q)?

[Yanksi]

I think the barrier energy you presented in the paper is actually using the euclidean distance? .

dd/dq indeed means d/dq distance(q), while d2d/dq2 means the hessian of distance(q). dd2/dq means the gradient of distance(q)^2

[zfergus]

You are correct, we present it as standard euclidian distance but note in Section 7 of our paper:

While for clarity in the preceding we derive IPC with unmodified distance evaluations, for numerical accuracy and efficiency our implementation applies squared distances for evaluations of the barrier, we use b(d^2,dhat^2), and related computations, thus avoiding squared roots. In turn expressions for contact forces,λk, and related terms must be modified, from our direct exposition and derivations above. To do so we rescale for consistent dimensions and units in our implementation; see our Supplemental for details.

Let me check your math and get back to you.

[zfergus]

I think the math should be

d(x) = sqrt(d2(x))
∇ d(x) = 1 / (2 * d(x)) * ∇ d2(x)
∇² d(x) = -1 / (4 * d2(x)^{3/2}) * ∇d2(x) * ∇d2(x)ᵀ + 1/(2*d(x)) * ∇²d2(x) 

[Yanksi]

Aha, thanks. I still got a little confusion here, as the gradient and hessian of d(x) can be easily derived by the equation you gave while only need to calculate the square root operation once. There should not be too much numerical issue and efficiency problem caused by it?

[james-bern]

This is a good question, and I'm following this repo so I'll try my best :P

I don't think the issue is really the square root operation per se. It's that if you use d(x) = sqrt(d2(x)) as your energy, then as Zachary points out in order to compute ∇ d(x) you have to divide by distance (if distance is very small, then this factor will be huge -> numerical problems). Also as a general rule I think squared distance is preferable because it's C^2 continuous (continuous first and second partials), whereas |x| is just C^0 (not smooth at x = 0). Note: The second derivative of |x| is discontinuous at the origin. Note: I definitely haven't tested d(x) out for this particular problem, but I personally like distance squared and I believe it's the standard choice for stuff like this 👍

[liminchen]

Yes, the issue for square root is not the efficiency or accuracy, but the robustness when distances are tiny. We didn't test unsquared distance either, just implemented all the things with squared distances at the beginning. If one really cares about constitutive modeling, e.g. what contact force magnitude there will be at a certain distance and what's the trend of their relation like in hyperelasticity, luckily log() functions can easily convert exponent to multiplications and it should be easy to customize. In our code we don't necessarily need \nabla d/\nabla x for the optimization, we directly compute \nabla b/\nabla x = \nabla b/\nabla d^2 * \nabla d^2/\nabla x