MATLAB Answers

Find matching points from two coordinate systems

122 views (last 30 days)
Yangfan Peng
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.

  7 Comments

Show 4 older comments
Matt J
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
darova on 8 Oct 2019
What if just find closest pairs of points? Using pdist2() for example?
Yangfan Peng
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?

Sign in to comment.

Accepted Answer

Matt J
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

  8 Comments

Show 5 older comments
Matt J
Matt J on 30 Apr 2020 at 16:13
Do you have the Optimization Toolbox?
Matt J
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

Sign in to comment.

More Answers (1)

Bruno Luong
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.

  2 Comments

Yangfan Peng
Yangfan Peng on 8 Oct 2019
Thank you for the suggestions. I used the fitgeotrans function and got the transformation object. But how do I go from there to determine which point matches to which other point? See above for detailed description of the problem.

Sign in to comment.

Sign in to answer this question.