Find matching points from two coordinate systems

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.

[darova]

Can you please attach some sample data? How final iamges should look like?

[Yangfan Peng]

They are just coordinates, there is no image. Each coordinate represents a cell location in the coordinate system of the respective microscope.

x y

167.54 87.75

106.60 138.29

69.06 193.54

125.94 229.45

162.99 231.56

198.90 267.31

271.54 239.20

228.31 209.63

235.30 154.54

214.66 67.76

And then there is another set of x and y coordinates of the same cells but on a different microscope with a different coordinate system that might be rotated, scaled, distorted...

[darova]

And what do you want to do with these (x,y)?

[Yangfan Peng]

coordinate_a = [376 455;421 489;465 537;353 512;355 535;329 571;377 593;417 598;482 575;355 634]
coordinate_b = [168 88;107 138;69 194;126 229;163 232;199 267;272 239;228 210;235 155;215 68]

This is the code for the two different set of coordinates. They represent the same cells in different microscope images. However, the points are mixed up, so point 1 in coordinate_a is not point 1 in coordinate_b (see image below). I would like to match the points to know which points belong to which in the different coordinate systems, e.g. point 1 in coordinate_a is point 7 in coordinate_b

[Matt J]

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]

What if just find closest pairs of points? Using pdist2() for example?

[Yangfan Peng]

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?

Answer 1

[Matt J]

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

https://www.mathworks.com/matlabcentral/fileexchange/26186-absolute-orientation-horn-s-method

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

[darova]

Why so many times edited? Do you have some problems?

[Matt J]

Just making improvements.

[Yangfan Peng]

Dear Matt,

thanks you very much for the code and the help! I am incredibly sorry that I am just getting back to you now. I have been distracted and busy with other projects, but am continuing on this problem now. Your function works very well for "easy" cases (great help and thanks for that!!). However, if two points get very close together, these two points seem to overlap and match to one coordinate in Bsorted (see attached image). Is there a way to restrict this overlap? Can we add some margin of minimal distance between the Bsorted points while minimizing the distances?

I know that this answer comes late, but I would be grateful for your help.

Thanks. Best, Yangfan

[Matt J]

I cannot copy/paste the A and B data for this example.

[Yangfan Peng]

A = [474 451;440 475;457 491;409 489;522 513;477 540;378 554;481 604;470 586; 375 618];
B = [150 58;122 133;121 195;189 229;264 213;276 168;255 115;187 112;237 131;76 113];

Here are the coordinates for A and B. Thanks! So overall, I am wondering whether there is a way to uniquely assign each coordinate to another without overlap if points get too close...? Thanks.

[Matt J]

Do you have the Optimization Toolbox?

[Yangfan Peng]

Yes I do.

[Matt J]

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

Answer 2

[Bruno Luong]

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.

[Matt J]

See also, fitgeotrans

https://www.mathworks.com/help/images/ref/fitgeotrans.html

[Yangfan Peng]

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.