# Find matching points from two coordinate systems

122 views (last 30 days)
Yangfan Peng on 7 Oct 2019
Edited: Matt J on 9 May 2020 at 18:45
Hello,
I have 10 points with x and y coordinate (cell locations imaged with microscope). I have used two different methods to image these cells, so I have two sets of coordinates that could be rotated, sheared, shifted, mirrored to each other. I would like to assign each point from one coordinate system to the matching point in the other coordinate system. How can I do that? In theory, some kind of transformation optimizing minimum distance between the points would be desirable.
In the end I would like to assign properties I obtained from the cells in one method to properties obtained with the other method.
Thank you.

Matt J on 8 Oct 2019
Generally speaking, it is an impossible problem. For example, suppose there were 3 points in each image with the A-points forming an equilateral triangle and the B-points an equilateral triangle in a different position. Then there is no unique matching between the points that can be determined purely from their relative spacing - the symmetry among the points is too perfect. There would have to be some a priori known asymmetry among the points in order to match them.
darova on 8 Oct 2019
What if just find closest pairs of points? Using pdist2() for example?
Yangfan Peng on 8 Oct 2019
I always have 7-10 points and they are mostly asymmetric. I would like to rotate/mirror the coordinates in a way that the points can be closest to each other. How can I combine pdist2() with the transformation part?

Matt J on 8 Oct 2019
Edited: Matt J on 8 Oct 2019
If the shear component of the deformation isn't too strong, the matchpoints function defined below might work. It will sort the rows of B to correspond with A and also return the corresponding permutation indices. It relies on absor, mentioned by Bruno, which you will have to download
function [permIndices,Bsorted]=matchpoints(A,B)
%
%IN:
%
% A: an Nx2 matrix of points
% B: an Nx2 matrix of points from A transformed and unordered
%
%OUT:
%
% permIndices: permutation indices of rows of B thought to match A
% Bsorted: the Nx2 permuted version of B
La=landmarks(A);
Lb=landmarks(B);
B3=B.'; B3(3,:)=0;
reg=absor( Lb,La,'doScale',1);
C3=(reg.s*reg.R)*B3+reg.t;
C=C3(1:2,:).';
[~,permIndices]=pdist2(C,A,'euclidean','Smallest',1);
if nargout>1
Bsorted=B(permIndices,:);
end
function L=landmarks(P)
G=pdist2(P,P); G(~G)=nan;
[i,j]=find( G==min(G(:)) ,1);
I=P(i,:);
J=P(j,:);
K=mean(P,1);
if norm(I-K)<norm(J-K)
[I,J]=deal(J,I);
end
L=[I;J;K].';
L(3,:)=0;
end
end
Applying it to your example data, I obtain,
A = [376 455;421 489;465 537;353 512;355 535;329 571;377 593;417 598;482 575;355 634];
B = [168 88;107 138;69 194;126 229;163 232;199 267;272 239;228 210;235 155;215 68];
[permIndices,Bsorted]=matchpoints(A,B)
plot(graph(1:10,1:10),'EdgeColor','none','XData',A(:,1),'YData',A(:,2));
hold on
plot(graph(1:10,1:10),'EdgeColor','none','XData',Bsorted(:,1),'YData',Bsorted(:,2));
hold off

Matt J on 30 Apr 2020 at 16:13
Do you have the Optimization Toolbox?
Yangfan Peng on 1 May 2020 at 6:40
Yes I do.
Matt J on 2 May 2020 at 4:48
Here is a modification of the code that works on your example. You have a very large outlier in your coordinates, however, so it can't always be guaranteed that such cases would be robustly handled. Note that the code uses the FEX file minL1intlin.m (Download).
function [permIndices,Bsorted]=matchpoints(A,B)
%
%IN:
%
% A: an Nx2 matrix of points
% B: an Nx2 matrix of points from A transformed and unordered
%
%OUT:
%
% permIndices: permutation indices of rows of B thought to match A
% Bsorted: the Nx2 permuted version of B
La=landmarks(A);
Lb=landmarks(B);
B3=B.'; B3(3,:)=0;
reg=absor( Lb,La,'doScale',1);
C3=(reg.s*reg.R)*B3+reg.t;
C=C3(1:2,:).';
N=size(A,1);
e=ones(1,N);
E=speye(N);
Aeq=[ kron(E,e) ; kron(e,E) ]; beq=[e,e].';
Q=kron(E,C.');
d=reshape(A.',[],1);
lb=zeros(1,N^2);
ub=ones(1,N^2);
intcon=1:N^2;
P=minL1intlin(Q,d,intcon,[],[],Aeq,beq,lb,ub);
P=round(reshape(P,N,[]));
permIndices=(1:N)*P;
if nargout>1
Bsorted=B(permIndices,:);
end
function L=landmarks(P)
G=pdist2(P,P); G(~G)=nan;
[i,j]=find( G==min(G(:)) ,1);
I=P(i,:);
J=P(j,:);
K=mean(P,1);
if norm(I-K)<norm(J-K)
[I,J]=deal(J,I);
end
L=[I;J;K].';
L(3,:)=0;
end
end

Bruno Luong on 8 Oct 2019
Look for the literature of image registration.
For simple rotation/scaling/translation you can use Matt J's submission
For more complex deformation, you need to apply spline deformation
There are a bunch of intermediate method for camera which take into account for camera cushion distortion or higher order. Pick one that is suitable for your need.