# Find the perfect overlay of 2 maps of points

16 vues (au cours des 30 derniers jours)
Alexia Bichon le 8 Nov 2019
Modifié(e) : Adam Danz le 30 Déc 2020
I would like to have a perfect overlay of these 2 maps of points. Most of the points have their pair in the other map, some don't.
I already tried to minimize the distance between the points but it doesn't work as the min distance is not the one that overlay perfectly the points.
Do you have a method to do that ?
##### 3 commentairesAfficher 1 commentaire plus ancienMasquer 1 commentaire plus ancien
Alexia Bichon le 8 Nov 2019
There is no 1:1 mapping, some points have a pair some don't.
It's like stars. You get a picture of the sky one day and a second picture the day after. Some stars have appeared, other have vanished but you are still able to see most of the constellations and you want to overlay them.
I tried to minimize the distance between the points red and blue, but the minimum is not the best overlay as the points without a pair are also trying to find their closest neighbor.
Adam Danz le 8 Nov 2019
Modifié(e) : Adam Danz le 8 Nov 2019
That's helpful. How are the data organized? Are there missing values in either dataset? Are there equal number of coordinates between the two data sets? Could you attach a mat file with two sets of coordinates?

Connectez-vous pour commenter.

### Réponse acceptée

Alexia Bichon le 11 Nov 2019
Modifié(e) : Alexia Bichon le 12 Nov 2019
To solve this problem, I used an optimization method in 2 steps :
1) I loop on the lateral shift (X) and stretch of my map to find the greatest number of points with a nearest neighbor at a distance of less than 11 inches : a "pair" of points.
( R = red points, B = blue points )
D_min = zeros(length(R_X),1);
pair_numbers = zeros(length(steps_stretch), length(steps_X));
for l = 1:length(steps_stretch)
for k = 1:length(steps_X)
l
k
B_X_slide_stretch = B_X*steps_stretch(l) + steps_M(k);
B_Y_slide = B_Y;
for i=1:length(R_X)
D_min(i) = min((B_X_slide_stretch - R_X(i)).^2 + (B_Y_slide - R_Y(i)).^2);
end
pair_numbers(l,k) = size(D_min(D_min<121),1);
end
end
[min_per_row, Index_column_local_min_per_row] = max(pair_numbers, [],2);
[min_table, Index_row] = max(min_per_row);
Index_column = Index_column_local_min_per_row(Index_row);
B_Y_pairing = B_Y;
B_X_pairing = B_X * steps_stretch(Index_row) + steps_X(Index_column);
Then, I consider only the R points (red) that I assume have a pair in the B (blue) set of points.
for i=1:length(R_X)
D_min(i) = min(sqrt((B_X_pairing - R_X(i)).^2 + (B_Y_pairing - R_Y(i)).^2));
end
R_X_pairing = R_X(D_min<11);
R_Y_pairing = R_Y(D_min<11);
2) The second part of the optimization is about minimizing the distance between a point and its pair by playing with the 3 parameters : X (lateral shift), Y (vertical shift) and the stretch.
I use the fmincon function from the Optimization_Toolbox to find the minimum of the function @sum_squared_distances.
nterval = bounds_X_translation(2) - bounds_X_translation(1);
nsteps = 5;
x = zeros(nsteps+1,3);
fval = zeros(nsteps+1,1);
exitflag = zeros(nsteps+1,1);
tic
for i = 0:nsteps
i+1
TolX = 10^(-5);
options = optimset('TolX', TolX, 'Tolfun', 1, 'MaxIter', 500, 'MaxFunEvals', 5000, 'Display', 'iter',...
'LargeScale', 'off', 'DiffMaxChange', 500*12, 'DiffMinChange', TolX/1000,...
'LevenbergMarquardt', 'on', 'LineSearchType', 'cubicpoly');
[x(i+1,:),fval(i+1),exitflag(i+1),output] = fmincon(@sum_squared_distances,[bounds_X_translation(1) + i*interval/nsteps 0 1],[],[],[],[],...
[bounds_X_translation(1) bounds_Y_translation(1) bounds_X_stretch(1)],...
[bounds_X_translation(2) bounds_Y_translation(2) bounds_X_stretch(2)],...
[],options);
end
Optimization_time = toc;
save('Results_optimization')
[m, ind] = min(fval);
results = x(ind,:);
B_Y_final = B_Y_pairing + results(1,2);
B_X_final = B_X_pairing * results(1,3) + results(1,1);
@sum_squared_distances : sum of the distances to be minimized
function ssd = sum_squared_distances(inp)
global B_X_pairing....
R_X_pairing_short...
B_Y_pairing...
R_Y_pairing_short
B_X_pairing_slide_stretch = inp(1) + B_X_pairing*inp(3);
B_Y_pairing_slide = inp(2) + B_Y_pairing;
D_min = zeros(length(R_X_pairing_short),1);
for i=1:length(R_X_pairing_short)
D_min(i) = min(sqrt((B_X_pairing_slide_stretch - R_X_pairing_short(i)).^2 + (B_Y_pairing_slide - IR_Y_pairing_short(i)).^2));
end
ssd = sum(D_min);
end
At the end I have a very good overlap of the 2 maps and I can count the number of points with a pair among the blue points.
Now I should write a code to remove the red points that have the same blue point as a closest neigbor.
Another option would be to try to minimize the number of blue pixels.
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Plus de réponses (2)

Thorsten le 8 Nov 2019
Modifié(e) : Thorsten le 8 Nov 2019
You can maximize the number of perfect matches, that is the number of red points with zero distance to a blue point.
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Alexia Bichon le 8 Nov 2019
Thank you for this idea. I think it can work ! It won't be 0 distance but I can put a threshold very low.

Connectez-vous pour commenter.

Adam Danz le 8 Nov 2019
Modifié(e) : Adam Danz le 30 Déc 2020
You can compute the difference between the x values from each set and the y values from each set which will create a matrix of offsets for the x values and the y values. If (and only if) coordinates (x1,y1) are a linear offset from coordinates (x0,y0) plus some extra non-paired coordinates from both sets, the most frequent offset value for x and y will be the amount you need to shift the dataset to match. This can be computed using mode().
% Create a fake dataset (x0,y0)
x0 = randi(100, 1, 50);
y0 = randi(100, 1, 50);
% Create a fake dataset (x1,y1)
x1 = randi(100, 1, 60);
y1 = randi(100, 1, 60);
% Match several of the values between the two sets and offset them a bit.
randSelect = unique(randi(50,1,40));
x1(randSelect) = x0(randSelect) + 0.5; % Offset some of the x values
y1(randSelect) = y0(randSelect) - 3.5; % Offset some of the y values
% Show inputs
clf()
hold on
plot(x0,y0,'bo', 'DisplayName', 'xy0')
plot(x1,y1,'ro', 'DisplayName', 'xy1')
% SOLUTION :
% Get most common x-offset and y-offset
% This approach uses implicit expansion requires matlab r2016b or later.
% It also uses the 'all' flag in mode() which requires r2018b or later.
xDiff = mode(x1(:) - x0(:).','all'); % implicit expansion requires matlab r2016b or later
yDiff = mode(y1(:) - y0(:).','all'); % implicit expansion requires matlab r2016b or later
% Shift (x1,y1) to overlay (x0,y0)
x1shift = x1 - xDiff;
y1shift = y1 - yDiff;
(x1shift, y1shift) are the new coordinate values for (x1,y1) that overlay (x0,y0) when paired.
If x1shift & y1shift does not result in overlap between the two sets, that indiciates that the difference between the two sets may not be linear or that you do not have enough pairs of coordinates that should overlap between the two datasets.
Update plot; recall that some xy1 values do not have a paired xy0 value in this demo.
plot(x1shift, y1shift, 'rx', 'DisplayName', 'xy1 shifted')
legend('Location','BestOutside')
##### 2 commentairesAfficher AucuneMasquer Aucune
Alexia Bichon le 11 Nov 2019
This function seems very powerful. However, I solved my problem with my version of matlab.
Adam Danz le 11 Nov 2019
I didn't realize you had an earlier version. The two lines that require later versions can easily be changed to work with r2007b.

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Surface and Mesh Plots dans Help Center et File Exchange

R2007b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by