Giter Site home page Giter Site logo

siavashk / pycpd Goto Github PK

View Code? Open in Web Editor NEW
487.0 16.0 112.0 14.72 MB

Pure Numpy Implementation of the Coherent Point Drift Algorithm

License: MIT License

Python 86.60% TeX 12.40% Makefile 1.00%
point-cloud registration expectation-maximization python

pycpd's Introduction

Python-CPD

Paper Citation joss_paper_badge

Build

Status
master master_badge
development development_badge

Documentation

Pure Numpy Implementation of the Coherent Point Drift Algorithm.

MIT License.

Introduction

This is a pure numpy implementation of the coherent point drift CPD algorithm by Myronenko and Song for use by the python community. It provides three registration methods for point clouds: 1) Scale and rigid registration; 2) Affine registration; and 3) Gaussian regularized non-rigid registration.

The CPD algorithm is a registration method for aligning two point clouds. In this method, the moving point cloud is modelled as a Gaussian Mixture Model (GMM) and the fixed point cloud are treated as observations from the GMM. The optimal transformation parameters maximze the Maximum A Posteriori (MAP) estimation that the observed point cloud is drawn from the GMM.

The registration methods work for arbitrary MxN 2D arrays where M is the number of "points" and N is the number of dimensions. A typical point cloud would be Mx2 or Mx3 for 2D and 3D points clouds respectively. For more information, please refer to my blog.

Installation

Install from PyPI

pip install pycpd

Installation from Source

Clone the repository to a location, referred to as the root folder. For example:

git clone https://github.com/siavashk/pycpd.git $HOME/pycpd

Install the package:

pip install .

or

make requirements
make build

Install Matplotlib for Visualization

For running sample registration examples under /examples, you will need Matplotlib to visualize the registration. This can be downloaded by running:

pip install matplotlib

or

make visualize

Usage

Each registration method is contained within a single class inside the pycpd subfolder. To try out the registration, you can simply run:

python examples/fish{Transform}_{Dimension}.py

where Transform is either rigid, affine or deformable and Dimension is either 2D or 3D. Note that examples are meant to be run from the root folder.

Example

Basic Usage

Basic usage includes providing any of the registration methods with 2 arrays that are MxN & BxN. E.g., they can have different numbers of points (M & B) but must have the same number of dimensions per point (N).

from pycpd import RigidRegistration
import numpy as np

# create 2D target points (you can get these from any source you desire)
# creating a square w/ 2 additional points. 
target = np.array([[0, 0], [0, 1], [1, 0], [1, 1], [0.5, 0], [0, 0.5]])
print('Target Points: \n', target)

# create a translation to apply to the target for testing the registration
translation = [1, 0]

# create a fake source by adding a translation to the target.
# in a real use, you would load the source points from a file or other source. 
# the only requirement is that this array also be 2-dimensional and that the 
# second dimension be the same length as the second dimension of the target array.
source = target + translation
print('Source Points: \n', source)

# create a RigidRegistration object
reg = RigidRegistration(X=target, Y=source)
# run the registration & collect the results
TY, (s_reg, R_reg, t_reg) = reg.register()

# TY is the transformed source points
# the values in () are the registration parameters.
# In this case of rigid registration they are:
#     s_reg the scale of the registration
#     R_reg the rotation matrix of the registration
#     t_reg the translation of the registration

The affine and deformable registration methods are used in the same way, but provide their respective transformation parameters.

Apply Transform to Another Point Cloud

Sometimes you may want to apply the transformation parameters to another point cloud. For example, if you have a very large point cloud it is sometimes appropriate to randomly sample some of the points for registration and then apply the transformation to the entire point cloud.

To do this, after fitting the above registration, you would run reg.transform_point_cloud(Y=points_to_transform). This will apply the learned registration parameters to the point cloud points_to_transform and return the transformed point cloud.

Tuning Registration parameters

For rigid and affine registrations the main parameter you can tweak is w. The w parameter is an indication of the amount of noise in the point clouds [0,1], by default it is set to 0 assuming no noise, but can be set at any value 0 <= w <1 with higher values indicating more noise.

For deformable registration, you can also tune alpha, beta, and use low_rank.

The alpha parameter (lambda in the original paper) identifies a tradeoff between making points align & regularization of the deformation. A higher value makes the deformation more rigid, a lower value makes the deformation more flexible.

The beta is the width of the Gaussian kernel used to regularize the deformation and thus identifies how far apart points should be to move them together (coherently). beta depends on the scale/size of your points cloud. Tuning beta can be simplified by normalizing the point cloud to a unit sphere distance.

The low_rank parameter is a boolean that indicates whether to use a regularized form of the deformation field. This further constrains the deformation, while vastly speeding up the optimization. num_eig is the number of eigenvalues to use in the low rank approximation. num_eig should be less than the number of points in the point cloud, the lower the smoother the deformation and the faster the optimization.

Testing

Tests can be run using pytest:

pip install pytest
pytest

or

make dev
make test

Documentation

The documentation can be built using pydoc3

make dev
make doc

Contributing

Contributions are welcome. Please see the guidelines outlined in the document: CONTRIBUTING.

Code of Conduct

We have adopted the code of conduct defined by the Contributor Covenant to clarify expected behavior in our community. For more information see the Code of Conduct.

pycpd's People

Contributors

agporto avatar brandondube avatar ckolluru avatar gattia avatar hugoledoux avatar johanforslund avatar khmariem avatar mdifranc avatar mhr avatar siavashk avatar vincentme avatar wckao avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

pycpd's Issues

repeat self.intialize() fails due to change in shape of self.t

registrations cannot be repeated because self.t changes shape during processing.
self.t.shape == (1,3) after __init__, but self.t.shape == (3,) after self.register(),

so repeat calls to self.register() gives
ValueError: operands could not be broadcast together with shapes
in the line:
self.TY = self.s * np.dot(self.Y, np.transpose(self.R)) + np.repeat(self.t, self.M, axis=0)
(or the corresponding line in affine_registration)

deformation registration not working with different scale points values

Hi,
Thanks for your CPD implementation.

When i try the nonrigid deformation, i follow the deformation registration example,this is ok.

when i change the points scale,such as 10,100,i find the deformation registration is not working wiht different scale value.

i also try affine registration and rigid registration,both that are working with different scale values.

i don't understand what is problem?

I'm attaching my code and the results in both cases.

Thanks,
junqiangchen

here is my code

`from functools import partial
from mpl_toolkits.mplot3d import Axes3D
from pycpd import deformable_registration, affine_registration
import time
import matplotlib.pyplot as plt
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)

factor = 100
target = np.loadtxt("Cpd/data/face.txt") * factor
X1 = np.zeros((target.shape[0], target.shape[1] + 1))
X1[:, :-1] = target
X2 = np.ones((target.shape[0], target.shape[1] + 1))
X2[:, :-1] = target
X = np.vstack((X1, X2))

source = np.loadtxt("Cpd/data/face_distorted.txt") * factor
Y1 = np.zeros((source.shape[0], source.shape[1] + 1))
Y1[:, :-1] = source
Y2 = np.ones((source.shape[0], source.shape[1] + 1))
Y2[:, :-1] = source
Y = np.vstack((Y1, Y2))

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

reg = deformable_registration(**{'X': X, 'Y': Y})
reg.register(callback)
plt.show()`

非刚性变换
100倍非刚性变换

Coherent Point Drift with STL files

You wouldn't have any advice for how to use your code with stl files to perform the rigid and non-rigid registration for statistical shape modelling would you?

Number of Iterations

Hi, is there any way I can increase the number of iterations? I noticed for me, it always stop at 100 iterations.

3D Deformation Registration not working

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

Invert transforms

Hi! First off: Thanks for this great library!

I'm working with data where I frequently also need to apply the inverse transform, i.e target->source instead of source->target. I was wondering if you would consider adding functionality for this? E.g. by implementing a __neg__ method for the transforms.

Best,
Philipp

Division by zero

Hi,

the parameter w is not checked for range, so this line:
c = c * self.w / (1 - self.w)

can divide by zero.
How about an additional line in init:
assert self.w<=0<1

BTW: I find it hard to understand the meaning of the documentation of w without reading the original paper. I suppose w=0 (default) is the case for Gaussian distributions only and w=1 for equal distribution (in the presence of noise), right?

how this package could be used for aligning two 2D data (or two 2d cloud points )

I want to align two TSNE data using this package could you tell me how to use this package ?

from functools import partial
from scipy.io import loadmat
import matplotlib.pyplot as plt
from pycpd import affine_registration,rigid_registration(X, Y)
import numpy as np
import time

def visualize(iteration, error, X, Y, ax):
    plt.cla()
    ax.scatter(X[:,0] ,  X[:,1], color='red')
    ax.scatter(Y[:,0] ,  Y[:,1], color='blue')
    plt.draw()
    print("iteration %d, error %.5f" % (iteration, error))
    plt.pause(0.001)

def main():
    X = np.random.random((10, 2))  # 10 points in 2 dimensions
    Y = np.random.random((10, 2))  # 10 points in 2 dimensions
    fig = plt.figure()
    fig.add_axes([0, 0, 1, 1])
    callback = partial(visualize, ax=fig.axes[0])

    reg = rigid_registration(X, Y)
    reg.register(callback)
    plt.show()

main()

MemoryError for relatively big point clouds

I'm getting the following stack trace on an attempt to register 3D point clouds with about 40K points each:

  File "/home/theuser/anaconda3/lib/python3.5/site-packages/pycpd/affine_registration.py", line 24, in register
    self.initialize()
  File "/home/theuser/anaconda3/lib/python3.5/site-packages/pycpd/affine_registration.py", line 86, in initialize
    XX = np.tile(XX, (self.M, 1, 1))
  File "/home/theuser/anaconda3/lib/python3.5/site-packages/numpy/lib/shape_base.py", line 912, in tile
    c = c.reshape(-1, n).repeat(nrep, 0)
MemoryError

This happens both with rigid_registration and affine_registration. Should this algorithm be so memory consuming? How can I estimate the memory usage in 3D from the number of points?
(I have 16 Gb RAM if it matters)

Incorrect Integer Division (?)

Quick question on this line:

c = (2 * np.pi * self.sigma2) ** (self.D / 2)

In python 3 this is floating point division, but in python 2 this looks like it will perform integer (floor) division, was this supposed to be floating point division under both versions?

[BUG] Missing transpose in AffineRegistration

Describe the bug
Hi! Thanks for the awsome code, it is very efficient!

I noticed that when I use the class AffineRegistration on my dataset, the result is occasionally wrong. I think the issue comes from "update_variance" method, in the calculation of "trBYPYP", "np.trace(np.dot(np.dot(self.B, self.YPY), self.B))" should be "np.trace(np.dot(np.dot(self.B.T, self.YPY), self.B))", which also equals to "np.trace(np.dot(self.A, self.B))". This is verified in another project probreg.

To Reproduce
Here are my source and target points which cause the code to fail using AffineRegistration:
source.txt
target.txt

Expected behavior
AffineRegistration should work for the given dataset.

Screenshots
Current matching results:
Screen Shot 2022-11-10 at 3 33 40 PM
Expected matching results:
Screen Shot 2022-11-10 at 3 35 23 PM

How to export transformed point cloud - non-rigid registration

Hi,

First of all, thank you for the great work. I'm currently trying to use your library to perform data fusion of CMM and 3D scanner data (for metrological application). I'm trying to perform non-rigid registration between two point clouds, but I cannot find how to export the transformed source point cloud.

For rigid transformation, I understand the example:

# create a RigidRegistration object
reg = RigidRegistration(X=target, Y=source)
# run the registration & collect the results
TY, (s_reg, R_reg, t_reg) = reg.register()

where the transformed point cloud and the transformation matrices are in TY, (s_reg, R_reg, t_reg) = reg.register()

But how do you get those for DeformableRegistration?

Thank you!

PS: I know this is not a feature request, but I don't know where to post my question.

Merge Development

Hey,

Just leaving a small note about maybe merging development with master. There are lots of good documentation additions that you added that the end-user is not getting. Also, I know you merged some of my changes with development where there were small errors in the implementation and would be useful for anyone blindly using this implementation to get those.

Does it work for higher dimensional point? Like 4D or more?

Thank you for sharing the code. Appreciate your work.
Just a quick question: does the code work for higher dimensional point sets?
Currently I have 2 datasets containing 4d data (may be the dimension would be higher in the future). Can I use the code?
Thank you so much!

Documentation + self.w unused

Hi Siavash

Just a few details:

  • would be great if there was some bit of code documentation (mainly constructor arguments and returned objects)
  • the self.w does not have any effect on the result

In any case great work! Thanks a lot for your efforts.
Cheers, Norman

Registration of overlapping point clouds

Hello everyone, I am currently working on my master thesis and would like to use this library for it.
To be able to tackle my actual task, I first need to merge two point clouds. These two point clouds represent the ground in front of a vehicle and are obtained with two 3D cameras (see picture).

image

I know the rough orientation of the two cameras, so that I can roughly pre-transform the point clouds. As can be seen, the point clouds overlap in a certain area. The goal now is to create an overall cloud from the two individual point clouds.

image

The point clouds are very large (23232 points per cloud), so I use only the right or left partial image for the registration. This halves the amount of data. Afterwards I always summarize 9 points by averaging.
Thus I receive two times 1276 points. Each point has an X, Y and Z coordinate.

image

image

This amount of data can be processed by my RaspberryPi.

If I now perform a registration (affine or rigid), the result deteriorates very strongly.

image

The code used looks like this:

# self.__teilbilderReduziert[i] contains the reduced point clouds. As the clouds are sorted by pixels, the shape will be adjusted
self.__links = np.reshape(self.__teilbilderReduziert[1], (1276,3))
self.__rechts = np.reshape(self.__teilbilderReduziert[0], (1276,3))  

self.__reg = RigidRegistration(**{'X': self.__links, 'Y': self.__rechts})
self.__reg.register()  

# self.__bildVorausgereichtet[i] includes the overall images. These are also sorted by pixels and are adjusted in shape 
self.__bildVorausgereichtet[1] = np.reshape(self.__bildVorausgereichtet[1], (23232,3))
self.__bildVorausgereichtet[1] = self.__reg.transform_point_cloud(self.__bildVorausgereichtet[1])

# Finally, the untransformed and the transformed image are saved
self.__bildVorausgereichtet[0] = np.reshape(self.__bildVorausgereichtet[0], (23232,3))
np.savetxt("Bild_rechts.txt", self.__bildVorausgereichtet[1])
np.savetxt("Bild_links.txt", self.__bildVorausgereichtet[0])

What can I do to solve my problem?

I would be very pleased about answers and help.
With best regards

René

Transformation with small number of points and fixed scale of 1.0

A point list with about 10 points is to be transformed. However, of these 10 points, only 3 points match the target point set. All others are to be ignored.
The scale factor is fixed at 1.0.
The tests with 10 points did not give a usable result (s_reg= 0.1591229754).
Only with the 3 matching points the result was correct, but the scale factor was 0.9998774093 although s=1.0 was passed.
Is it possible to solve this problem with PYCPD?

Example data sets:

Source Data
120.0000 120.0000
116.8168 91.1093
130.5683 141.9590
98.6609 82.2620
124.3115 124.2831
123.2244 165.7592
163.1087 123.1000
166.3159 121.9226
150.1782 89.1302
160.1123 102.0151

Source Data with only matching points
98.6609 82.2620
123.2244 165.7592
160.1123 102.0151

Target original data
666498.735 5158548.881
666450.520 5158621.326
666523.059 5158608.670

translation by
-666300 -5158400
target test data set
198.7350 148.8810
150.5200 221.3260
223.0590 208.6700

Test_with_matching_points
Test_with_all_points

Nonrigid deformation not working for 3D to 2D projection

Hi,
Thanks for your CPD implementation.

I'm having trouble matching a 3D model to points extracted from a 2D picture. It works fine with the affine deformation, but the result is not very precise (I'm guessing because the pinhole camera model needs nonrigid deformation).

When I try the nonrigid deformation, the source data is not modified and the error is unbelievably small. I don't understand, am I doing anything wrong or is there a problem with the deformable_registration class? I'm attaching my code (very close to the example code) and the results in both cases.

Thanks,

Gaetan

from functools import partial
from mpl_toolkits.mplot3d import Axes3D
from pycpd import deformable_registration, affine_registration
import time

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)


target = np.loadtxt("PointClouds/" + testp[:-4] + ".pcd", skiprows = 11)
X1 = np.zeros((target.shape[0], target.shape[1] + 1))
X1[:,:-1] = target
X2 = np.ones((target.shape[0], target.shape[1] + 1))
X2[:,:-1] = target
X = np.vstack((X1, X2))

source = np.loadtxt('SmallFrame.pcd', skiprows = 11)
Y1 = np.zeros((source.shape[0], source.shape[1] + 1))
Y1[:,:-1] = source
Y2 = np.ones((source.shape[0], source.shape[1] + 1))
Y2[:,:-1] = source
Y = np.vstack((Y1, Y2))

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

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

Affine:
Untitled

Nonrigid:
Screenshot from 2019-10-25 16-41-39

Attempting to perform 3D deformable registration

I have tried running your code using a trimesh import of a 3D mesh object. However, the code breaks saying that this only works on 2D numpy arrays? How do you perform 3D registrations if it's only able to work with 2D arrays?

build from source and can't find "expectation_maximization_registration" module

Hi,
I choose to build from source and I encounter the error message like this:
image

Generally, it says it fails with this line when running __init__.py:
from expectation_maximization_registration import expectation_maximization_registration in affine_registration.py.

I found that change the above line to
from .expectation_maximization_registration import expectation_maximization_registration
in every file in __init__.py solved this import problem.

Transposed rotation matrix

For rigid registration in function transform_point_cloud() the rotation matrix is used incorrectly:

def transform_point_cloud(self, Y=None):
        if Y is None:
            self.TY = self.s * np.dot(self.Y, self.R) + self.t
            return
        else:
            return self.s * np.dot(Y, self.R) + self.t

In these expressions the rotation matrix should be used transposed:

def transform_point_cloud(self, Y=None):
        if Y is None:
            self.TY = self.s * np.dot(self.Y, self.R.T) + self.t
            return
        else:
            return self.s * np.dot(Y, self.R.T) + self.t

Therefore, as it is now, the solution iterates around transposed rotation matrix, and the resulting rotation matrix is also transposed.

Some questions.

If I have two sets of points, one set is the feature points on another set of points, and the number of feature points is less than the target set, how should I adjust the parameters? It seems that it cannot match the target set well. Also, is the number of iterations fixed or can it be controlled when iterating to a certain degree? I look forward to your answer!

Affine3D is broken

The logic for Affine3D registration is broken. Updating the transform results in solving a singular linear system.

[BUG]

Thanks for this fantastic and well-coded library. One problem: I found that DeformableRegistration().transform_point_cloud is raising an exception when the number of points to be transformed is different from those used to compute the registration. one of the examples triggers this:

% python examples/fish_deformable_3D_register_with_subset_of_points.py
Traceback (most recent call last):
File "/fish_deformable_3D_register_with_subset_of_points.py", line 55, in
main()
File "pycpd-master/examples/fish_deformable_3D_register_with_subset_of_points.py", line 43, in main
YT = reg.transform_point_cloud(Y=Y)
File "/lib/python3.9/site-packages/pycpd/deformable_registration.py", line 68, in transform_point_cloud
return Y + np.dot(self.G, self.W)
ValueError: operands could not be broadcast together with shapes (182,3) (91,3)

rigid transformation

Thanks for the code.
Now I have two sets of points, one is a template, and the other is the same as the template but has some noisy points and is rotated 90 degrees. When I used the rigid method to register these two sets, I found that the effect was not ideal. The rigid transformation include translation, rotation, and scaling. What is the reason for this?
Looking forward to your reply, thanks again.
16535_fusept.txt
FormTax_00013_fusept.txt

Running Examples - Plot not showing

Hi,

When running examples only the first iteration and their respective point placements appear on the plot. The following iterations do not appear and I am unsure why.

Thank you in advance!

UserWarning: This figure includes Axes that are not compatible with tight_layout

Hi, when I try to run your examples in Pycharm (fish_deformable_2D.py), I get the following error:
/snap/pycharm-professional/151/helpers/pycharm_matplotlib_backend/backend_interagg.py:64: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. self.figure.tight_layout()
and only the initial figure is displayed correctly with all the iteration figures empty
myplot1
myplot2

Lack of documentation or descriptive coding conventions

This issue was raised in #12 and then closed, despite the actual issue persisting. Think it would be valuable to leave this open so that potential contributors know that it is an open problem.

Would also be useful if the variable names were descriptive, so that the burden of explicitly documenting the code is lessened. Variables like B and t do nothing to explain the functionality and hinder understanding.

3D Deformation Registration: large poinclouds

Hi!
I'm exploring your excellent library, tried to apply 3D Deformation Registration to my data.
The data are the point clouds that I generated from the Snoopy dataset, they can be downloaded here: https://drive.google.com/drive/folders/1_jVk7gNw1p34lU1d_AKBvrEoYTKwDipG?usp=sharing
snoopies
I read and downsample the pointclouds with Open3d:

voxel_size  = 0.02   # 2cm
sourcepld = o3d.io.read_point_cloud("data/Snoopy__000010.pcd")
targetpld = o3d.io.read_point_cloud("data/Snoopy__000055.pcd")

sourcepld = sourcepld.voxel_down_sample(voxel_size)
targetpld = targetpld.voxel_down_sample(voxel_size)

static    = np.asarray(sourcepld.points)
moving = np.asarray(targetpld.points)

The number of points is the following: static - (329, 3); moving - (331, 3).

Then I follow your recommendations from another issue and first perform Affine registration, which matches the point clouds nicely. The following Deform Registration doesn't seem to enhance the matching.
(The code copied from: #41 )

# 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.1, '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])

pypcd_pcd

Changing the alpha and beta parameters does not help: the smaller any of them - the more points tend to shrink to the center.

Could you possibly tell is this method of non-rigid registration applicable to such point clouds, how could the result be improved?

Python3.6 Error + Hack

I got the following errors when working with Python3. The package as-is has no problems with Python2 on my system though.

Traceback (most recent call last):
  File "src/label.py", line 21, in <module>
    from cpd import cpd
  File "/Users/alvinwan/Documents/Aspire/pcmatch/src/cpd.py", line 1, in <module>
    from pycpd import rigid_registration
  File "/Users/alvinwan/anaconda3/envs/matching/lib/python3.6/site-packages/pycpd/__init__.py", line 1, in <module>
    from affine_registration import affine_registration
ModuleNotFoundError: No module named 'affine_registration'

I fixed it by using relative imports:

from .affine_registration import affine_registration
from .rigid_registration import rigid_registration
from .deformable_registration import deformable_registration

Food for thought. I considered making another fork and pushing pycpd3 to pip, but I wanted your @siavashk thoughts/permission first. I think the change would break Python2 compatibility.

The source (moving) point cloud should be immutable

I see the s, t, and R attributes in the rigid_registration object, which is iteratively used to correct the Y point cloud… but as the registration proceeds, those values seem to converge on 1, [0 0 0], and the identity matrix (respectively), and I can’t seem to figure out how to rebuild the “complete” translation vector or rotation matrix that could be used to transform the original Y point cloud in a single step, which I could then use to register channels in my non-point-cloud experimental images

Deformable registration - alpha and beta

Hi, Thanks for this great library. I have a question that I was wondering if you would be able to answer. If I understand correctly, parameters alpha and beta in the deformable registration step should be scale dependent. In other words, if I were doing two separate registrations (one of small objects, other of large objects), the ideal values for those parameters should be different. Is that your understanding? If so, any chance you would have an idea of how alpha and beta change with the size of the structure?

pycpd package citation

Hello,

I am using pycpd in my work and want to acknowledge pycpd package by citing it.
Apart from the original CPD paper by Myronenko and Song, is there any publication associated with the implementation of CPD (pycpd package)?

Thank you,
Bramsh

[JOSS Review] Minor issues

Hi @gattia and @siavashk,

This is part of the JOSS review process. I found two minor issues after testing your PyCPD package. Could you please address them?

  • Found an error when running python examples/fish_deformable_3D_register_with_subset_of_points.py:
Traceback (most recent call last):
  File "examples/fish_deformable_3D_register_with_subset_of_points.py", line 55, in <module>
    main()
  File "examples/fish_deformable_3D_register_with_subset_of_points.py", line 43, in main
    YT = reg.transform_point_cloud(Y=Y)
  File "/home/geo3d/miniconda3/envs/joss-pycpd/lib/python3.8/site-packages/pycpd/deformable_registration.py", line 68, in transform_point_cloud
    return Y + np.dot(self.G, self.W)
ValueError: operands could not be broadcast together with shapes (182,3) (91,3)

This error also failed the pytest:

============================= test session starts ==============================
platform linux -- Python 3.8.13, pytest-7.1.2, pluggy-1.0.0
rootdir: /home/geo3d/Zexin/project/pycpd
collected 7 items                                                              

testing/affine_test.py ..                                                [ 28%]
testing/deformable_test.py ..F                                           [ 71%]
testing/rigid_test.py ..                                                 [100%]

=================================== FAILURES ===================================
_______________________________ test_3D_low_rank _______________________________

    def test_3D_low_rank():
        fish_target = np.loadtxt('data/fish_target.txt')
        X1 = np.zeros((fish_target.shape[0], fish_target.shape[1] + 1))
        X1[:, :-1] = fish_target
        X2 = np.ones((fish_target.shape[0], fish_target.shape[1] + 1))
        X2[:, :-1] = fish_target
        X = np.vstack((X1, X2))
    
        fish_source = np.loadtxt('data/fish_source.txt')
        Y1 = np.zeros((fish_source.shape[0], fish_source.shape[1] + 1))
        Y1[:, :-1] = fish_source
        Y2 = np.ones((fish_source.shape[0], fish_source.shape[1] + 1))
        Y2[:, :-1] = fish_source
        Y = np.vstack((Y1, Y2))
    
        reg = DeformableRegistration(**{'X': X, 'Y': Y, 'low_rank': True})
        TY, _ = reg.register()
        assert_array_almost_equal(TY, X, decimal=0)
    
        rand_pts = np.random.randint(Y.shape[0], size=int(Y.shape[0]/2))
>       TY2 = reg.transform_point_cloud(Y=Y[rand_pts, :])

testing/deformable_test.py:56: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <pycpd.deformable_registration.DeformableRegistration object at 0x7f5ff28b62b0>
Y = array([[-0.28181282, -0.79432332,  1.        ],
       [-0.37639549, -0.58040108,  1.        ],
       [-0.91641226,  ...1006222,  1.        ],
       [-0.10490431, -0.60374743,  1.        ],
       [-0.50945755,  0.00261754,  1.        ]])

    def transform_point_cloud(self, Y=None):
        """
        Update a point cloud using the new estimate of the deformable transformation.
    
        """
        if Y is None:
            self.TY = self.Y + np.dot(self.G, self.W)
            return
        else:
>           return Y + np.dot(self.G, self.W)
E           ValueError: operands could not be broadcast together with shapes (91,3) (182,3)

../../../miniconda3/envs/joss-pycpd/lib/python3.8/site-packages/pycpd/deformable_registration.py:68: ValueError
=========================== short test summary info ============================
FAILED testing/deformable_test.py::test_3D_low_rank - ValueError: operands co...
========================= 1 failed, 6 passed in 0.43s ==========================
  • A few warnings can be fixed by replacing is not with !=:
/home/geo3d/Zexin/project/pycpd/pycpd/rigid_registration.py:45: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if R is not None and (R.ndim is not 2 or R.shape[0] is not self.D or R.shape[1] is not self.D or not is_positive_semi_definite(R)):
/home/geo3d/Zexin/project/pycpd/pycpd/rigid_registration.py:49: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if t is not None and (t.ndim is not 2 or t.shape[0] is not 1 or t.shape[1] is not self.D):
/home/geo3d/Zexin/project/pycpd/pycpd/rigid_registration.py:49: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if t is not None and (t.ndim is not 2 or t.shape[0] is not 1 or t.shape[1] is not self.D):
/home/geo3d/Zexin/project/pycpd/pycpd/affine_registration.py:30: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if B is not None and (B.ndim is not 2 or B.shape[0] is not self.D or B.shape[1] is not self.D or not is_positive_semi_definite(B)):
/home/geo3d/Zexin/project/pycpd/pycpd/affine_registration.py:34: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if t is not None and (t.ndim is not 2 or t.shape[0] is not 1 or t.shape[1] is not self.D):
/home/geo3d/Zexin/project/pycpd/pycpd/affine_registration.py:34: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if t is not None and (t.ndim is not 2 or t.shape[0] is not 1 or t.shape[1] is not self.D):

Numba decorators

I was thinking of trying to implement pycpd using numba decorators, and I was wondering if you had tried it yourself and whether you think it might speed up calculations significantly. Since pycpd is straight numpy, it sounds to me like it would, but I wanted to make sure I wouldn't waste time on something you had tried. Any thoughts?

"Error" sounds misleading...

Hi,
In the bunny example, I found "Error" in the visualization is misleading, because it is increasing all the time, and a the last iteration, error suddenly becomes zero.
for example, iteration 20:
screenshot from 2018-08-10 11-30-56

iteration 39
screenshot from 2018-08-10 11-33-19

when I check the source code and the paper, it sounds like error means the difference between Q (objective function) in current iteration and previous iteration. For example, in rigid_registration.py:

self.q = (xPx - 2 * self.s * trAR + self.s * self.s * self.YPY) / (2 * self.sigma2) + self.D * self.Np/2 * np.log(self.sigma2)
self.err = np.abs(self.q - qprev)

Given the np.abs(), we even can't figure out whether the Q is actually increasing?? or decreasing.
Q might converge through EM algorithm, so error (actually the difference) < tolerance might be a metric to see whether it converges. But I hope the name of the variable could be clearer.

Second question is, while optimizing Q, is the sigma2 also decreasing? and thus could be a metric to see whether it converges?

Question about tolerance

I've been playing with speeding up PyCPD and have been using Cython to write a few functions and have gone back to the original paper / matlab code to implement some of the other methods (low rank in particular).

Anyways, in that process, I've noticed that you don't always use the same error functions (tolerance) as they used in the matlab scripts.

For example, your deformable registration uses the change in sigma between two iterations where as the CPD original paper calculates an objective function that I believe is the negative log likelihood and uses a change in this to determine when the algorithm has converged.

Just curious if you remember if there was a reason you altered the tolerance/error? Im finding that using the original implementation that the algorithm is "converging" too early and therefore isn't actually getting the right transformation - I don't pass the tests. Im using the tests you setup as a means of determining if my changes are breaking anything.

Thanks!

How to fix scale factor and reduce initial displacement

Hi,
I am trying to align two small point clouds which do not have to be scaled.
How do I set the initial scale closer to 1? Is it possible to fix the final scale to 1?
Currently the rigid transform aligns correctly the point clouds. However, it takes too much time since it starts with a very small scale.

Regards,
Michele

Add Stanford Bunny

I want to add the stanford bunny so I can better test my code for 3D point clouds.

bug in pycpd/rigid_registration.py

I believe that line 41 should be
return self.s * np.dot(Y, self.R) + np.tile(self.t, (Y.shape[0], 1))
instead of the current
return self.s * np.dot(Y, self.R) + np.tile(np.transpose(self.t), (Y.shape[0], 1))

[Documentation] Add contribution file with guidelines

Part of reviewing process for JOSS acceptance (openjournals/joss-reviews/issues/4681)

I couldn't find clear instructions on how to participate in this project. Ideally you should provide a CONTRIBUTE file with some informations about:

  • How to contribute to the software (e.g. how to create a pull request or minimum requirements for accept a PR )
  • How to report issues: your procedure of choice for handle bug reports (is there a template to follow? is required to provide a minimum reproducible example? is ok just a simple issue opened on GitHub?)
  • How to reach support: alternate channels where to get help (email, forum, ecc..)

Example fish_affine_2D diverges instead of converges

Hi,

Thank you for making this implementation of CPD.

I am looking to align two point clouds by affine transformation and later add/draw additional points related to the source points with the achieved transformation matrices (B, t) after succesful registration .

Starting out with the provided example -- and using the fish data -- fish_affine_2D.py, I get a divergence terminating in error (Q) = -550 after 100 iterations. I did not make any modifications to the code. What would be your best guess for what is going on?
Installed with pip install pycpd (2.0.0)

Thanks!

Axis out of bounds error

When I tried testing the data containing 3d point clouds an error occurred in the _makeKernel module. The error is as mentioned below

Traceback (most recent call last): File "fishDeformable.py", line 39, in <module> main() File "fishDeformable.py", line 35, in main reg.register(callback) File "/home/ashok/libraries/Python-CPD-master/core/DeformableRegistration.py", line 26, in register self.initialize() File "/home/ashok/libraries/Python-CPD-master/core/DeformableRegistration.py", line 82, in initialize self._makeKernel() File "/home/ashok/libraries/Python-CPD-master/core/DeformableRegistration.py", line 112, in _makeKernel diff = np.sum(np.multiply(diff, diff), self.D); File "/usr/local/lib/python3.4/dist-packages/numpy/core/fromnumeric.py", line 1840, in sum out=out, keepdims=keepdims) File "/usr/local/lib/python3.4/dist-packages/numpy/core/_methods.py", line 32, in _sum return umr_sum(a, axis, dtype, out, keepdims) ValueError: 'axis' entry is out of bounds

Looks like the axis value is not supportive for 3D yet. Otherwise the algorithm works fine for 2D cases. Since I am not that much familiar with this algorithm, is it possible for the author to make known the possible fix?

Regards,

#0K

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.