Find closest point between two sets of points

Hello I have a matrix A=[x,y,c] in which there are around 10^5 points stored. There are 10^5 rows and 3 columns. Columns x,y define the coordinates of a point and column c has a value stored.

The points stored in the matrix define a numerical mesh and the value C is the value in that mesh point. The matrix A is not ordered, but the grid consists of concentric discrete circles of different radius. Transforming Matrix A into polar coordinates makes a lot of sense and i have done so using the cart2pol function. Resulting in A2=[theta,radius,c]. Now i have a second solution matrix B2=[theta,radius,c] from a different solution time and i want to do an average between both solution steps A2 and B2. The mesh [x,y] is the same in both cases, My problem is that the solution in Matrix B has been rotated by Phi degrees. So first i want to rotate B by Phi radians, and then i want to match the new positions of the mesh points from B to the positions of the mesh points in A. and then do the average between both solutions

for example i have

A=[theta,radius,c]=[0,0,1; 0,1,2; 0,2,3] and B = [theta,radius,c]=[0,0,7; 0,1,8; 0,2,9; 1,0,1.1; 1,1,2.1; 1,2,3.1]

Now i rotate B by -1 degree back so that its mesh overlaps with the mesh from Solution A

A=[theta,radius,c]=[0,0,1; 0,1,2; 0,2,3] and B = [theta,radius,c]=[-1,0,7; -1,1,8; -1,2,9; 0,0,1.1; 0,1,2.1; 0,2,3.1]

Now elements 4,5 and 6 of B overlap with elements 1, 2 and 3 of A, but of course my code doesnt know that

I want to find which elements(indices) of B correspond to which elements(indices) of A in order to do a simple average between A and B

and get my correct result

 A_averaged=[0,0,1.05; 0,1,2.05; 0,2,3.05]

I think the first step would be to order A2. How can i do this? I am thinking about something like A2= [theta1,radius1,c;theta1,radius2,c;theta1,radiusn,c;theta2,radius1,c;theta2,radius2,c;.....]

Rotating B is easy, i would just use a rotation matrix

  R = [cos(angle) -sin(angle); sin(angle) cos(angle)];

to rotate the x,y values stored in B. But this operations results in A(1) no longer having the same x,y coordinates as B(1) so i would need to sort B again? or find which point in B is nearest to A(1)? Any help is appreciated, i dont know how to sort a matrix nor what solution path is more effcient and surely there are paths i havent thought about.

[dpb]

I the x,y mesh is fixed, then simply averaging them on an element-by-element basis would seem to be all that would need to be done?

If not, explain why not...and like always, it would help immensely to visualize the problem and provide actual code if you would attach a (small) subset of the data files that illustrates the problem and provides folks something to work with rather than trying to make up something...

[Ángel Brito Gadeschi]

The system is rotating but the mesh is defined in absolute coordinates(nonrotating) so first i need to either rotate the first solution on top of the second or the second on top of the first. So before rotation the first mesh element of A is the same as the first mesh element of B, but after rotation i need to find which index of B has the x,y coordinates of the first element in A and do this for all mesh points. I have tried sorting both matrixes after the radius and angle, but rounding errors lead to some points not completely overlapping. I.e if A has a point at r=4mm and phi=0° sometimes it happens that the corresponding point in B after rotation has r=3.99 and phi=0.001° which after sorting both A and B matrix leads to them having different indexes and a simple (A(:,3)+B(:,3))/2 doesnt cut it anymore.

I have 350 different radii and 400 different angles. Rounding the x,y coordinates to 1 decimal place in radius and angle leads to the sorting function to sort A and B into 450 radii and 440 angles. So its not a perfect solution

I am unsure about the rumber of radii. but i am sure there are only 400 points at each radius, i.e one mesh element every 0.9°. Can I tell the sorting function to always match 400 points to the samish ~radius and skip the rounding step altogether?

I have done plots with polarscatter, which is results in it plotting a lot of discrete points. Is it possible to do a continous plot and have the rotation and average operation happen on that continuous plane?

The plane is a circle with a hole in the middle i.e circle with origin at (x,y)=0, inner radius of r=1 and outer radius of r=2

Using scattered interpolation has worked pretty good i think

[Matt J]

Using scattered interpolation has worked pretty good i think

If you are content with it, please Accept-click the answer.

Answer 1

[Matt J]

One way is to use scattered interpolation.

 F=scatteredInterpolant(xA,yA,cA);
 
 result=0.5*( F(xB,yB)+cB); %average the interpolated points with cB.

Here, I assume the (xB,yB) coordinates for B are already rotation corrected.

[Ángel Brito Gadeschi]

Hello, i have tried your suggestion and it has worked pretty well. If i understand scatteredInterpolant correctly, there is no need to sort my matrices so that A(1) and B(1) have the same coordinates?

What happens if the points of B dont completely overlap the ones in A after rotation? I.e. a point in A has values (x,y)=(1,1) and the corresponding point of B after rotation has (x,y)=(0.98,0.98)?

How would i scale your solution for more than 2 solutions? lets say i have 100 matrices and i have already rotated them all

[Matt J]

What happens if the points of B dont completely overlap the ones in A after rotation? I.e. a point in A has values (x,y)=(1,1) and the corresponding point of B after rotation has (x,y)=(0.98,0.98)?

That's why we are interpolating. F(0.98,0.98) estimates, through interpolation, what A(0.98,0.98) would be.

How would i scale your solution for more than 2 solutions? lets say i have 100 matrices and i have already rotated them all

You would loop.

[Ángel Brito Gadeschi]

OK, thanks a lot for the help until now. My only question remaining, the values stored in "result" are using the mesh defined with xA and yA correct?

[Matt J]

The mesh being interpolated is A, but the result is ordered according to B.