normalize log probability by max probability before exponentiation to…

[wckao]

… avoid numerical issues

[siavashk]

If you are dividing correspondence probabilities by exp(-p_min / 2sigma^2), then you need to also divide the contribution of the uniform distribution by the same value. However, without knowing why/when this numerical issue arises I do not approve this PR.

[wckao]

The reason is that when log probabilities are all very small ( say -1000), taking exp(-1000) will make all probabilities zero. And then den will be zero too. And P/den will explode.

I didn't get where you are referring to by "you need to also divide the contribution of the uniform distribution by the same value" though.

[siavashk]

The reason is that when log probabilities are all very small ( say -1000), taking exp(-1000) will make all probabilities zero. And then den will be zero too. And P/den will explode.

I am handling division by zero in expectation_maximization_registration.py:74, where I substitute zero denominators with np.finfo(float).eps. Testing it out now in the interactive shell, it seems that it does not do what I thought it does. I probably should change line den[den==0] to den[den < eps] to reflect this.

I didn't get where you are referring to by "you need to also divide the contribution of the uniform distribution by the same value" though.

c, as calculated in lines 67-69 and used in line 75, is the probability that the a sample point from the target point belongs to a uniform distribution. If you are adding Pmin to all exponents, you need to multiply c by exp(-Pmin / (2 * sigma2)) so that the weight of the uniform distribution is not modified. Otherwise you are dynamically modifying this weight.

[wckao]

The reason is that when log probabilities are all very small ( say -1000), taking exp(-1000) will make all probabilities zero. And then den will be zero too. And P/den will explode.

I am handling division by zero in expectation_maximization_registration.py:74, where I substitute zero denominators with np.finfo(float).eps. Testing it out now in the interactive shell, it seems that it does not do what I thought it does. I probably should change line den[den==0] to den[den < eps] to reflect this.

The real issue I encounter is that in deformable_registration.py line 43, the denominator has a Np term, which will be zero since your P will already be all zero.

I didn't get where you are referring to by "you need to also divide the contribution of the uniform distribution by the same value" though.

c, as calculated in lines 67-69 and used in line 75, is the probability that the a sample point from the target point belongs to a uniform distribution. If you are adding Pmin to all exponents, you need to multiply c by exp(-Pmin / (2 * sigma2)) so that the weight of the uniform distribution is not modified. Otherwise you are dynamically modifying this weight.

But if c is also multiplied by by exp(Pmin / (2 * sigma2)) ( no minus sign in front of Pmin? ), the problem will still be there since all P will effectively still be zero?

I haven't read the paper, but should we add c to some part of P as something like a uniform prior?

[siavashk]

The minus in exp(-Pmin / (2 sigma2)) was a typo. It should be exp(Pmin / (2 sigma2)). You are right.

The real issue I encounter is that in deformable_registration.py line 43, the denominator has a Np term, which will be zero since your P will already be all zero.

When does this usually happen? The only way that I see this happening is if the two point clouds are really far, and the initial value for sigma is low. If you are trying to handle this case, there might be a bug in the initial approximation of sigma2.

But if c is also multiplied by by exp(Pmin / (2 * sigma2)) ( no minus sign in front of Pmin? ), the problem will still be there since all P will effectively still be zero? I haven't read the paper, but should we add c to some part of P as something like a uniform prior?

You are on the right track. Instead of hacking the posterior, you should be looking into why the probabilities are near zero. As with c, it is factored in the posterior as part of the denominator. Look at equation 6 of the cpd paper.

[siavashk]

I looked at the data that you sent me through email. The issue is that the source (Y) and target (X) point clouds are initially far from each other.

In this case, you need to perform a rigid registration first to compensate for the gross misalignment of your point clouds. Once the rigid registration has converged, you can take the rigidly transformed source point cloud (TY) and perform a deformable registration between X and TY.

I ran the rigid followed by deformable registration on the data that you sent me. Visually, rigid registration is all you need. I do not see major improvement using deformable registration.

from functools import partial
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pycpd import rigid_registration, deformable_registration
import numpy as np

def visualize(iteration, error, X, Y, ax):
    plt.cla()
    ax.scatter(X[:,0],  X[:,1], X[:,2], color='red', label='Target')
    ax.scatter(Y[:,0],  Y[:,1], Y[:,2], color='blue', label='Source')
    ax.text2D(0.87, 0.92, 'Iteration: {:d}\nError: {:06.4f}'.format(iteration, error), horizontalalignment='center', verticalalignment='center', transform=ax.transAxes, fontsize='x-large')
    ax.legend(loc='upper left', fontsize='x-large')
    plt.draw()
    plt.pause(0.001)

def main():
    data = np.load('data.npz')
    X = data['X']
    Y = data['Y']

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    callback = partial(visualize, ax=ax)

    rigid = rigid_registration(**{ 'X': X, 'Y': Y })
    TY, _ = rigid.register(callback)
    deformable = deformable_registration(**{ 'X': X, 'Y': TY })
    deformable.register(callback)
    plt.show()

if __name__ == '__main__':
    main()