norlab-ulaval / libpointmatcher

An Iterative Closest Point (ICP) library for 2D and 3D mapping in Robotics
BSD 3-Clause "New" or "Revised" License
1.58k stars 542 forks source link

Residual errors - Weighting each point #443

Open maximecharriere opened 3 years ago

maximecharriere commented 3 years ago

Hello, this isn't a issue but i have two questions.

Q1) Best residual error getter

I want to get the exact same value used by the error minimizer to do the gradient descent, so taking into account the effect of the outliners weight and the Weight property spoken in previous question.

I saw two way of getting the residual error. With getResidualError() as @Ellon has shown in #193 or with getWeightedPointUsedRatio() as shown in icp_advance_api.cpp.

getResidualError() seems to be what I want, but are we sure that the outliners weight and the Weight property will be taken into account?

Q2) Weight each point

Each points of my reading have a Weight property. Points with Weight ↦ 1 must be very close of the reference point cloud. Points with Weight ↦ 0 can be further away.

I have already chosen a pseudo-sinus outlier filter, so not all matching pairs between the reading and reference have the save weight in the error minimizer.
But now I want, as said before, that the Weight property influence the error minimizer so points with bigger weight are more important.

Do you know where I must modify the error minimizer code to use my Weight property ?

maximecharriere commented 3 years ago

[UPDATE]

I have begun to study the two questions, so here are my progress. I would need an expertise for the first question.

Q1) Reliable getResidualError()

First of all, getWeightedPointUsedRatio() is definitly not what I want. I need getResidualError().

Current implementation

For now, getResidualError() function compute the error like this: -> Link to code

N = numPoint
p = readingPointCloud
p' = readingPointCloud after fitting
q = referencePointCloud

R = rotation matrix to minimize the error
t = translation matrix to minimize the error
err=\sum_{i=1}^N\|p_i'-q_i\|

The sum of all Euclidean distances between the two point clouds is done to calculate the error.
Weights determined by the outliners are not taken into account.

My implementation

In the report _"Least-Squares Rigid Motion Using SVD by Olga Sorkine and Michael Rabinovich, 2017, ETHZ"_ used in libpointmatcher to compute the error minimization, it is said that the goal of the algorithm is to minimize the following equation:

err=\sum_{i=1}^N\omega_i\|(Rp_i+t)-q_i\|^2

So shouldn't the correct implementation for calculating the residual error be this ? :

err=\sum_{i=1}^N\omega_i\|p_i'-q_i\|^2

Which could be coded like this:

template<typename T>
T PointToPointErrorMinimizer<T>::computeResidualError(const ErrorElements& mPts)
{
    return  (mPts.weights*mPts.matches.dists.transpose())(0,0);
}

The only thing that seems strange to me is that we take the squared distance. The error thus becomes much greater. Do you think it would be better to take the distance and not the squared distance?

Also, do you think that the result should be normalized? So if the reading point cloud has many points, the error is comparable to a point cloud with few points.

err=\frac{\sum\omega_i\|p_i'-q_i\|^2}{\sum\omega_i}

Q2) Weight each point individually

I think I will not modify the ErrorMinimizer to take into account the weight of each point, but I will create a new OutlierFilter that modify the weight of each matches depending on the weight of each point.

pomerlef commented 3 years ago

To give a bit of history, the function computeResidualError was implemented as an afterthought. It was implemented to support other applications requested by users.

So, I think the question of what is a residual is hill posed. It can unfortunately be a lot of things and should most probably follows the error minimizer implementation.

Maybe we should have multiple functions with explicit names (e.g., mean distance, mean weighted distance, squared or not). Otherwise, user won't know what they receive, which is currently the case.