Point set registration (Rigid and Non-rigid).

Rigid, Affine and Non-rigid point set registration toolbox. Coherent Point Drift (CPD) algorithm.

Description

Point set registration Matlab toolbox. The Coherent Point Drift (CPD) algorithm. CPD allows to align and to find correspondences between two N-D point sets. You can select to use rigid, affine or non-rigid transformations.

Multiple options and demos are in the 'examples' folder.
Also see the project webpage:
http://www.bme.ogi.edu/~myron/matlab/cpd/

To use:
a) Unzip and add to Matlab path with subfolders.
b) Run cpd_make.m
c) Check for examples in the 'example' folder.

References:
 A. Myronenko, X. Song,"Point set registration: Coherent Point Drift", Technical Report, OHSU, 2009.
 A. Myronenko, X. Song, M. Carreira-Perpinan, "Non-rigid point set registration: Coherent Point Drift", Neural Information Processing Systems (NIPS) 2006.

Enjoy!
Please report any bugs on my email ( see: http://www.bme.ogi.edu/~myron)

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Fast Gaussian Transform mex implementation

• MATLAB release: MATLAB 7.2 (R2006a)

Other Files

cpd/core/cpd_make.m,
cpd/core/cpd_register.m,
cpd/core/FGT/fgt_model.c,
cpd/core/FGT/fgt_predict.c,
cpd/core/mex/cpd_P.c,
cpd/core/mex/cpd_Pappmex.c,
cpd/core/mex/cpd_Pcorrespondence.c,
cpd/core/Nonrigid/cpd_GRBF.m,
cpd/core/Nonrigid/cpd_GRBF_lowrank.m,
cpd/core/Nonrigid/cpd_GRBF_lowrankQS.m,
cpd/core/Rigid/cpd_affine.m,
cpd/core/Rigid/cpd_rigid.m,
cpd/core/utils/cpd_deform.m,
cpd/core/utils/cpd_denormalize.m,
cpd/core/utils/cpd_G.m,
cpd/core/utils/cpd_normalize.m,
cpd/core/utils/cpd_P_FGT.m,
cpd/core/utils/cpd_Pfast.m,
cpd/core/utils/cpd_plot_iter.m,
cpd/core/utils/cpd_R.m,
cpd/core/utils/cpd_transform.m,
cpd/core/utils/disptime.m,
cpd/data/bun_zipper.mat,
cpd/data/bun_zipper_res2.mat,
cpd/data/bun_zipper_res3.mat,
cpd/data/bun_zipper_res4.mat,
cpd/data/cpd_data2D_fish.mat,
cpd/data/cpd_data3D_bunny.mat,
cpd/data/cpd_data3D_face.mat,
cpd/data/cpd_data3D_face50.mat,
cpd/examples/cpd_affine_example1.m,
cpd/examples/cpd_affine_example2.m,
cpd/examples/cpd_Nonrigid_example1.m,
cpd/examples/cpd_Nonrigid_example2.m,
cpd/examples/cpd_Nonrigid_example3.m,
cpd/examples/cpd_Nonrigid_example4.m,
cpd/examples/cpd_Nonrigid_example5.m,
cpd/examples/cpd_Nonrigid_example6.m,
cpd/examples/cpd_rigid_example1.m,
cpd/examples/cpd_rigid_example2.m,
cpd/examples/cpd_rigid_example3.m,
cpd/examples/cpd_rigid_example4.m,
cpd/examples/cpd_rigid_example5.m,
cpd/examples/cpd_rigid_example6.m,
cpd/INSTALL.txt,
cpd/README.txt,
license.txt

[Luigi Rosa]

Excellent work. It could be used for many patter recognition problems.

[Sebastien Paris]

 

[jichao zhao]

why I have the following errors? thanks

??? Undefined function or method 'cpd_normalize' for input arguments of type 'double'.

Error in ==> cpd_register at 129
if opt.normalize, [X,Y,normal]=cpd_normalize(X,Y); end;

Error in ==> cpd_affine_example1 at 13
Transform=cpd_register(X,Y,opt);

[Andriy Myronenko]

2Jichao Zhao:
Looks like you haven't added all CPD subfolders to the Matlab path.
In Matlab, go to File-Set_Path-Add_with_subfolders, select the CPD main folder, and click save.

[jichao zhao]

thanks, Andrily. It is working perfectly, cheers,jichao

[Steffen]

great work, Andriy.

Is there any way to register with translation and rotation only without scaling? I set the .scale option to 0, but the scaling value in the computed transformation is not equal to 1.

[Steffen]

I checked the code and found out that one cannot use normalization and fixed scaling together. So, I set opt.normalize=0. however, the computed transformation seems to be far away from the optimal transformation. here is my test code:

function test_rigidTranform_noScaling( )

X = [121.3102 58.0957 46.3863
    83.3762 82.3011 69.5807
    88.7430 46.7092 84.9432
    51.3302 67.5140 109.3134];

Y = [36.7239 68.1019 52.6539
    61.0168 46.4281 88.7562
    76.5892 82.4971 84.0007
    99.3448 58.1160 121.3272];

% Set the options
opt.method='rigid'; % use rigid registration
opt.viz=0; % show every iteration
opt.outliers=0; % do not assume any noise

opt.normalize=0; % normalize to unit variance and zero mean before registering (default)
opt.scale=0; % fixed global scaling = 1
opt.rot=1; % estimate strictly rotational matrix (default)
opt.corresp=1; % compute the correspondence vector at the end of registration (default)

opt.max_it=100; % max number of iterations
opt.tol=1e-8; % tolerance

% registering Y to X
[Transform] = cpd_register(X,Y,opt);

figure,cpd_plot_iter(X, Y); title('Before');
figure,cpd_plot_iter(X, Transform.Y); title('After registering Y to X');

% compute error as mean minimum distance between points
error_before = computeError(X, Y)
error_after = computeError(X, Transform.Y)

    function [E, D] = computeError(X, Y)
        %compute for each x the euclidean distance to the closest y
        for i = 1 : size(X,1)
            D(i) = min(sqrt(sum((repmat(X(i,:), size(Y,1), 1) - Y).^2,2)));
        end
        E = mean(D);
    end
end

Any idea, what's going on? I noticed that in cpd_rigid.m the value of ntol can be negative: ntol=(L-L_old)/L; I'm not sure if this is intended. I changed it to ntol=abs((L-L_old)/L). The registration result is only slightly better, though.

Another minor thing. In cpd_register.m the line
if ~isfield(opt,'viz') || isempty(opt.tol), opt.viz = 1; end;
should be
if ~isfield(opt,'viz') || isempty(opt.viz), opt.viz = 1; end;

Best,
Steffen.

[Andriy Myronenko]

Hi Steffen,
Thanks for bringing up some bugs in the cpd code.
1) In cpd_rigid, you right, it should be with abs: ntol=abs((L-L_old)/L). It was my mistake. Fix that please. Now your code should work fine when using scaling.
2) Without scaling estimation (rotation and translation only)
(opt.normalize=0; opt.scale=0;), It get stuck into a local minimum of rotation. This is partially because your data set is very sparse and co-planar. Try using a slightly different initialization of one of your point-sets, let's say rotate Y a bit, For instance:

angle=0.1;
ca=cos(angle); sa=sin(angle);
R=[0 ca -sa;
   0 sa ca ;
   1 0 0];
Y=Y*R;

and then run the registration with your settings. You should get an accurate registration result. Perhaps, more systematic solution would be to run the rigid registration with several different rotation initializations.
Please, email my if you find more bugs or have more questions.
I'll fix these bugs and update my code shortly.
-Andriy

[su su]

 

[su su]

Great Work.
I want to register pair of images with your toolbox. But the size of my images are 320 by 240. How can I change them into N by 2 point sets. Thanks.