3D Deformation Registration not working

[BramshQamar]

Hello,

Very helpful and much-needed package. Thank you for sharing it with the community.

I am trying to run fish_deformable_3D.py example to perform 3D deformation registration on my data.

I made some minor changes in the code to make it work on my data. It doesn't give me any error but I don't see any changes in the final output. Here's the code I am running:

from functools import partial
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pycpd import DeformableRegistration
import numpy as np
%matplotlib

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}'.format(
        iteration), 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 my_cpd(source, target):

    X = source
    print("x shape = ", X.shape)

    Y = target

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

    reg = DeformableRegistration(**{'X': X, 'Y': Y})
    ty, pr =reg.register(callback)
    plt.show()

    return ty, pr

moved, pr = my_cpd(static, moving)

Here's the sample data I tried:

static object:

array([[-53.186558 ,  14.807759 ,  16.818079 ],
       [-39.329647 ,  12.384467 ,  17.518822 ],
       [-33.80166  ,   2.464514 ,  22.554306 ],
       [-34.9263   , -11.145478 ,  26.075163 ],
       [-33.56905  , -25.149336 ,  27.467941 ],
       [-33.669647 , -38.030262 ,  22.674452 ],
       [-39.25681  , -42.317623 ,  11.026509 ],
       [-44.4608   , -39.757362 ,  -1.8381366],
       [-50.407585 , -30.750744 , -10.419908 ],
       [-59.45718  , -22.393133 , -14.041341 ]], dtype=float32)

moving object:

array([[-39.16961  ,  32.603317 ,  57.101723 ],
       [-29.31882  ,  26.920492 ,  46.608566 ],
       [-20.51905  ,  26.39721  ,  33.882977 ],
       [-28.453857 ,  19.19351  ,  25.273869 ],
       [-31.410952 ,   3.661875 ,  27.189167 ],
       [-32.405537 , -12.103731 ,  24.411377 ],
       [-34.302643 , -26.626593 ,  17.430185 ],
       [-33.778683 , -35.498913 ,   5.182527 ],
       [-45.300465 , -32.50154  ,  -4.5435076],
       [-60.56806  , -30.297039 ,  -7.113083 ]], dtype=float32)

The output I get is the following:

moved object:

array([[-39.16960981,  32.60331622,  57.1017204 ],
       [-29.31892942,  26.92039722,  46.60837435],
       [-20.52011817,  26.39640361,  33.88213857],
       [-28.59517549,  19.13941354,  25.20927783],
       [-31.53889324,   3.71243713,  27.0872747 ],
       [-32.42795935, -12.11223542,  24.42580914],
       [-34.3022238 , -26.6962291 ,  17.50597361],
       [-33.83356926, -35.54964881,   5.23310869],
       [-45.32547539, -32.5305448 ,  -4.55158402],
       [-60.52187787, -30.24801277,  -7.16867977]])

There's an extremely small difference in the output. I am wondering if I am doing something wrong or need to change some parameters to make it work?

I tried changing tolerance value but still no difference.

Thank You, Bramsh

[gattia]

I played around with your data a bit. First, if you want to force the points to try and find a neighbour/match you can set/tweak the alpha:

reg = DeformableRegistration(**{'X': static, 'Y': moving, 'alpha': 0.01, 'beta': 2})

alpha should trade between regularizing the deformation and having points matching very closely. The lower you make alpha the more closely the point sets should match. beta tweaks the smoothness of the deformation.

Lowering alpha will make your points move, and they will end up closer to the target than when they started! So, we have some progress :).

Your points are really sparse and don't have an obvious way they match. At first I wondered if there just wasn't a great way for them to match and that is why CPD was failing. Then, I wondered if it was just a scaling issue - seems that is likely the case. What I did to remedy is I ran an affine registration, and then took the affine registration output and fed it into a deformable. This gave much better results.

Try something like:

# run affine reg
affine_reg = AffineRegistration(**{'X': static, 'Y': moving})
ty_affine, pr_affine = affine_reg.register()
# look how affine reg compares to your static points. 
fig = plt.figure(figsize=(16,9))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(static[:,0], static[:, 1], static[:,2])
ax.scatter(ty_affine[:,0], ty_affine[:,1], ty_affine[:,2])

# now feed points from affine reg into deformable reg
reg = DeformableRegistration(**{'X': static, 'Y': ty_affine, 'alpha': 0.001, 'beta': 5})
ty, pr = reg.register()

# now see how these final points match. 
fig = plt.figure(figsize=(16,9))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(static[:,0], static[:, 1], static[:,2])
ax.scatter(ty[:,0], ty[:,1], ty[:,2])

alpha = 0.001 beta = 5

Gave me an excellent match, with one point in no-mans-land between the two point-clouds. I also don't know what I am looking at, so this could be horrible registration - but its something :). Hope this helps!.

An aside note I will make in hopes that it helps clarify the language/code a bit. I got a bit confused reading your code because you interchange the words moving/target/source/static.

These words should be defined like: X=target=static Y=source=moving

The source points should be moving to get to the target points. I think your current code is:

X = source = static Y = target = moving

You end up with your static and moving inputs (to your function) going to the right location(s) when they got passed to the function, but they do so in a confusing way. I found the language confusing when I first started playing with these, so hopefully this is helpful.

[BramshQamar]

Hello @gattia,

Thank you for the quick response.

Changing alpha and beta parameters worked for me. I think performing affine before deformable registration is a good idea!

Thanks again.

Bramsh

[gattia]

Excellent. Because your issue seems to be resolved I'm closing this issue.

Thanks.