连接一组点得到一个非自相交的非凸多边形
Connecting a set of points to get a non-self-intersecting non-convex polygon
我有一组无序的二维点，代表建筑物的角。我需要将它们连接起来才能得到建筑物的轮廓。
这些点是通过组合不同个人收集的不同多边形而获得的。我的想法是使用这些多边形来按顺序获取点（例如，取最大和最小多边形之间的区域并连接这些点，使其进入该区域）。
我尝试使用最小距离标准并根据角度连接点。但不幸的是它不起作用。我拥有的一件有用的东西是许多多边形的原始数据，其中的点顺序是正确的。那么是否有可能与那些连接这些点的多边形进行比较？正如我上面提到的，我的教授提出了采用最大和最小多边形并将中间区域用作缓冲区的想法。所有的点都会落在这个缓冲区内。但是我不确定如何实现它。
X = [364.533 372.267 397.067 408.133 382.471 379.533 329.250 257.200 199.412 195.267 184.385 168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 267.714 306.929 312.143 357.733 421.333 431.000 371.867 364.533];
Y = [192.027 233.360 228.627 286.693 314.541 292.960 327.450 340.500 348.671 326.693 269.308 330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 64.643 56.857 77.786 69.493 133.293 180.427 142.160 192.027];
预期结果是代表建筑物平面图的闭合多边形。我有 15 个建筑样本，代码需要适用于所有人。一些建筑物不保留拐角之间的直角标准。我附上了我拥有的数据。我拥有的点是通过整合多边形获得的。那么有什么方法可以使用这个多边形（其中的点是有序的）actual data before integration
编辑
所以，我可以使用下面提到的想法找到解决方案。
备注：我手动补了一个漏点。而且，我去掉了底部的两个小角。要么，这些一共必须是四个角，要么就可以当成没有角。
我明确声明，因为我的想法结合了假设，角通常有 90 度角。
一般做法
- 通过下面提到的方法找到点的顺序。
- 对于所有点，根据找到的顺序在一定范围内确定潜在的“邻居”。
- 对于每两个邻居，确定邻居#1 - 当前点 - 邻居#2 之间的角度。理想情况下，这个角度应该是90度。
- 对于所有候选组合，找到总距离最小的组合，即距离（邻居#1 - 当前点）+距离（当前点 - 邻居#2）。
我意识到通过在所有点上使用 for 循环，导致绘制所有线两次。此外，许多计算可能会被矢量化并从循环中移出。优化不是我现在的意图。 ;-)
% Point data of building corners; modified!
X = [285.400 372.267 397.067 408.133 382.471 379.533 199.412 195.267 184.385 168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 267.714 306.929 312.143 357.733 421.333 431.000 371.867 364.533];
Y = [130.150 233.360 228.627 286.693 314.541 292.960 348.671 326.693 269.308 330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 64.643 56.857 77.786 69.493 133.293 180.427 142.160 192.027];
% Place approximative center of building at (0, 0)
X = X - mean(X);
Y = Y - mean(Y);
C = [mean(X), mean(Y)];
% Sort points by angle with respect to center
[~, idx] = sort(atan2(X, Y));
% Rearrange points with respect to sorted angles
X = X(idx);
Y = Y(idx);
% Number of data points
n = numel(X);
% Calculate direction vectors for X and Y coordinates
dvX = repmat(X.', 1, n);
dvX = dvX - dvX.';
dvY = repmat(Y.', 1, n);
dvY = dvY - dvY.';
% Calculate distances
dst = sqrt(dvX.^2 + dvY.^2);
% Number of "neighbouring" points to be considered with respect to the order
nn = 8;
figure(1);
hold on;
% Center
plot(C(1), C(2), 'kx', 'MarkerSize', 15);
% Plain points
plot(X, Y, '.', 'MarkerSize', 15);
for k = 1:n
% Index
text(X(k) + 0.05, Y(k) + 0.05, num2str(k), 'FontSize', 12);
% Set up neighbourhood
nbh = mod([k-nn/2:k-1 k+1:k+nn/2], n);
nbh(nbh == 0) = n;
% Calculate angles and total distance arrays
ang = Inf(nn);
len = Inf(nn);
for ii = 1:nn
l = nbh(ii);
d1 = [dvX(k, l) dvY(k, l)];
for jj = ii+1:nn
m = nbh(jj);
d2 = [dvX(k, m) dvY(k, m)];
len(ii, jj) = dst(k, l) + dst(k, m);
ang(ii, jj) = abs(pi/2 - acos(dot(d1, d2) / (norm(d1) * norm(d2))));
end
end
% Find candidates with angle difference < 10 degree
cand = find(ang < pi/18);
% For these candidates, find the one with the shortest total distance
[~, I] = min(len(cand));
% Get corresponding indices
[I, J] = ind2sub([nn nn], cand(I));
cand = nbh([I J]);
% Lines
plot([X(k) X(cand(1))], [Y(k) Y(cand(1))], 'b', 'LineWidth', 1);
plot([X(k) X(cand(2))], [Y(k) Y(cand(2))], 'b', 'LineWidth', 1);
end
hold off;
输出图像：
近似（！）解决方案是确定所找到的点描述的轮廓的中心，并使用相对于中心的 atan2 按角度对点进行排序。请参阅以下代码片段以进行可视化：
% Points
X = 2 * rand(1, 15) - 1;
Y = 2 * rand(1, 15) - 1;
% Center
C = [0, 0];
% Determine indices
[~, idx] = sort(atan2(X, Y));
figure(1);
hold on;
% Center
plot(C(1), C(2), 'kx', 'MarkerSize', 15);
% Plain points
plot(X, Y, '.', 'MarkerSize', 15);
% Indices and lines
for k = 1:numel(X)
text(X(idx(k)) + 0.05, Y(idx(k)) + 0.05, num2str(k), 'FontSize', 12);
if (k == numel(X))
plot([X(idx(k)) X(idx(1))], [Y(idx(k)) Y(idx(1))], 'b');
else
plot([X(idx(k)) X(idx(k+1))], [Y(idx(k)) Y(idx(k+1))], 'b');
end
end
hold off;
给出以下输出：
虽然我确信一定数量的凹面会被正确处理，但恐怕对于给定的示例（尤其是上半部分）它会失败。这是因为图像不是完美的俯视图，因此角度有点“扭曲”。
不过，也许排序可以提高您的最小距离方法。
这是一个解决方案，适用于具有由垂直* 线（如您示例中的线）构成的轮廓的形状。思路如下：
- 我们旋转点以将*它们对齐到 XY 网格。
- 我们将具有相同* X 或 Y 坐标的点分组。
- 对于每个点，我们计算两个点：水平方向最近的点和垂直方向最近的点，来自允许的族。
- 构建连接矩阵并转换回来。
就像在 中一样，我必须更改数据集：添加一个缺失的角点（第 30 点），删除两个非角点（第 7-8 点），删除重复的最后一个点。
* - 大约.
function A = q55511236
%% Initialization:
% Define points:
X = [364.533 372.267 397.067 408.133 382.471 379.533 329.250 257.200 199.412 195.267 184.385 ...
168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 ...
267.714 306.929 312.143 357.733 421.333 431.000 371.867];
Y = [192.027 233.360 228.627 286.693 314.541 292.960 327.450 340.500 348.671 326.693 269.308 ...
330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 ...
64.643 56.857 77.786 69.493 133.293 180.427 142.160];
%% Preprocessing:
% Centering:
XY = [X;Y] - [mean(X); mean(Y)];
% Rotation:
[U,~,~] = svd(XY,'econ');
rXY = (U.' * XY).';
% Fixing problems w/ some points:
rXY = vertcat(rXY, [-21.8, 66]); % add missing point
rXY(7:8, :) = NaN; % remove non-corners
% figure(); scatter(rXY(:,1),rXY(:,2));
%% Processing:
% Group points according to same-X and same-Y
CLOSE_ENOUGH_DISTANCE = 10; % found using trial and error
[~,~,sameXpts] = uniquetol(rXY(:,1), CLOSE_ENOUGH_DISTANCE, 'DataScale', 1);
[~,~,sameYpts] = uniquetol(rXY(:,2), CLOSE_ENOUGH_DISTANCE, 'DataScale', 1);
% Create masks for distance evaluations:
nP = size(rXY,1);
[maskX,maskY] = deal(zeros(nP));
maskX(sameXpts == sameXpts.') = Inf;
maskY(sameYpts == sameYpts.') = Inf;
% Compute X and Y distances separately (we can do this in the rotated space)
dX = abs(rXY(:,1) - rXY(:,1).') + maskX + 1./maskY;
dY = abs(rXY(:,2) - rXY(:,2).') + maskY + 1./maskX;
[~,nX] = min(dX);
[~,nY] = min(dY);
% Construct connectivity matrix:
A = false(nP);
idxTrue = sub2ind(size(A), repmat(1:nP, [1,2]), [nX(:).', nY(:).']);
A(idxTrue) = true;
%% Plot result:
% Rotated coordinates:
figure(); gplot(A, rXY, '-o'); text(rXY(:,1), rXY(:,2), string(1:nP));
uXY = (U*rXY.').';
% Original coordinates:
figure(); gplot(A, uXY, '-o'); text(uXY(:,1), uXY(:,2), string(1:nP)); axis ij;
导致：
答案使用的概念是'Travelling salesman problem'。在这些点周围创建了一个缓冲区，这个缓冲区作为一个额外的标准包含在内。
a=[141 188 178 217 229 282 267 307 313 357 372 422 434 365 372 398 411 382 382 233 229 191 185 166 156 183 173 114 97 149 139 139];
b=[109 103 79 76 140 132 64 56 78 72 141 133 180 192 234 228 287 293 315 348 343 348 329 332 270 268 225 240 194 184 108 108];
X=[364.5333 232.1333 397.0667 157.5333 431 421.3333 306.9286 184.3846 357.7333 199.4118 168.6429 179.4737 408.1333 382.4706 150.7333 372.2667 184.8 138.1053 312.1429 108.5333 174.5 195.2667 257.2 99.33333 379.5333 371.8667 329.25 280.7059 267.7143 218.2778];
Y=[192.0267 140.36 228.6267 274.4933 180.4267 133.2933 56.85714 269.3077 69.49333 348.6706 330.8571 80.44211 286.6933 314.5412 182.76 233.36 101.8933 111 77.78571 239.2 226.7857 326.6933 340.5 193.4667 292.96 142.16 327.45 130.5529 64.64286 74.35556];
R = [a' b'];
d = 12;
polyout = polybuffer(R,'lines',d)
figure
%imshow(I2);
hold on
%plot(R(:,1),R(:,2),'r.','MarkerSize',10)
plot(X,Y,'r.', 'MarkerSize', 15)
plot(polyout)
axis equal
hold off
[s,t] = boundary(polyout); %%this is the boundary polygon of the buffer
numPoints = length(clustersCentroids);
x = X; %these are the points to be connected
y = Y;
x([1 2],:)=x([2 1],:);
y([1 2],:)=y([2 1],:);
figure
plot(x, y, 'bo', 'LineWidth', 2, 'MarkerSize', 17);
grid on;
imshow(I2);
xlabel('X', 'FontSize', 10);
ylabel('Y', 'FontSize', 10);
% Make a list of which points have been visited
beenVisited = false(1, numPoints);
% Make an array to store the order in which we visit the points.
visitationOrder = ones(1, numPoints);
% Define a filasafe
maxIterations = numPoints + 1;
iterationCount = 1;
% Define a current index. currentIndex will be 1 to start and then will vary.
currentIndex = 1;
while sum(beenVisited) < numPoints
visitationOrder(iterationCount) = currentIndex;
beenVisited(currentIndex) = true;
% Get the x and y of the current point.
thisX = x(currentIndex);
thisY = y(currentIndex);
%text(thisX + 0.01, thisY, num2str(currentIndex), 'FontSize', 35, 'Color', 'r');
% Compute distances to all other points
distances = sqrt((thisX - x) .^ 2 + (thisY - y) .^ 2);
distances(beenVisited)=inf;
distances(currentIndex) = inf;
% Don't consider visited points by setting their distance to infinity.
[out,idx] = sort(distances);
xx=[x(currentIndex) x(idx(1))]
yy=[y(currentIndex) y(idx(1))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(1);
else
xx=[x(currentIndex) x(idx(2))]
yy=[y(currentIndex) y(idx(2))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(2);
else
xx=[x(currentIndex) x(idx(3))]
yy=[y(currentIndex) y(idx(3))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(3);
else
xx=[x(currentIndex) x(idx(4))]
yy=[y(currentIndex) y(idx(4))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(4);
end
end
end
end
end
% Plot lines in that order.
hold on;
orderedX = [x(visitationOrder); x(1)];
orderedY = [y(visitationOrder) ;y(1)];
plot(orderedX,orderedY, 'm-', 'LineWidth', 2);
title('Result', 'FontSize', 10);
我有一组无序的二维点，代表建筑物的角。我需要将它们连接起来才能得到建筑物的轮廓。
这些点是通过组合不同个人收集的不同多边形而获得的。我的想法是使用这些多边形来按顺序获取点（例如，取最大和最小多边形之间的区域并连接这些点，使其进入该区域）。
我尝试使用最小距离标准并根据角度连接点。但不幸的是它不起作用。我拥有的一件有用的东西是许多多边形的原始数据，其中的点顺序是正确的。那么是否有可能与那些连接这些点的多边形进行比较？正如我上面提到的，我的教授提出了采用最大和最小多边形并将中间区域用作缓冲区的想法。所有的点都会落在这个缓冲区内。但是我不确定如何实现它。
X = [364.533 372.267 397.067 408.133 382.471 379.533 329.250 257.200 199.412 195.267 184.385 168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 267.714 306.929 312.143 357.733 421.333 431.000 371.867 364.533];
Y = [192.027 233.360 228.627 286.693 314.541 292.960 327.450 340.500 348.671 326.693 269.308 330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 64.643 56.857 77.786 69.493 133.293 180.427 142.160 192.027];
预期结果是代表建筑物平面图的闭合多边形。我有 15 个建筑样本，代码需要适用于所有人。一些建筑物不保留拐角之间的直角标准。我附上了我拥有的数据。我拥有的点是通过整合多边形获得的。那么有什么方法可以使用这个多边形（其中的点是有序的）actual data before integration
编辑
所以，我可以使用下面提到的想法找到解决方案。
备注：我手动补了一个漏点。而且，我去掉了底部的两个小角。要么，这些一共必须是四个角，要么就可以当成没有角。
我明确声明，因为我的想法结合了假设，角通常有 90 度角。
一般做法
- 通过下面提到的方法找到点的顺序。
- 对于所有点，根据找到的顺序在一定范围内确定潜在的“邻居”。
- 对于每两个邻居，确定邻居#1 - 当前点 - 邻居#2 之间的角度。理想情况下，这个角度应该是90度。
- 对于所有候选组合，找到总距离最小的组合，即距离（邻居#1 - 当前点）+距离（当前点 - 邻居#2）。
我意识到通过在所有点上使用 for 循环，导致绘制所有线两次。此外，许多计算可能会被矢量化并从循环中移出。优化不是我现在的意图。 ;-)
% Point data of building corners; modified!
X = [285.400 372.267 397.067 408.133 382.471 379.533 199.412 195.267 184.385 168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 267.714 306.929 312.143 357.733 421.333 431.000 371.867 364.533];
Y = [130.150 233.360 228.627 286.693 314.541 292.960 348.671 326.693 269.308 330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 64.643 56.857 77.786 69.493 133.293 180.427 142.160 192.027];
% Place approximative center of building at (0, 0)
X = X - mean(X);
Y = Y - mean(Y);
C = [mean(X), mean(Y)];
% Sort points by angle with respect to center
[~, idx] = sort(atan2(X, Y));
% Rearrange points with respect to sorted angles
X = X(idx);
Y = Y(idx);
% Number of data points
n = numel(X);
% Calculate direction vectors for X and Y coordinates
dvX = repmat(X.', 1, n);
dvX = dvX - dvX.';
dvY = repmat(Y.', 1, n);
dvY = dvY - dvY.';
% Calculate distances
dst = sqrt(dvX.^2 + dvY.^2);
% Number of "neighbouring" points to be considered with respect to the order
nn = 8;
figure(1);
hold on;
% Center
plot(C(1), C(2), 'kx', 'MarkerSize', 15);
% Plain points
plot(X, Y, '.', 'MarkerSize', 15);
for k = 1:n
% Index
text(X(k) + 0.05, Y(k) + 0.05, num2str(k), 'FontSize', 12);
% Set up neighbourhood
nbh = mod([k-nn/2:k-1 k+1:k+nn/2], n);
nbh(nbh == 0) = n;
% Calculate angles and total distance arrays
ang = Inf(nn);
len = Inf(nn);
for ii = 1:nn
l = nbh(ii);
d1 = [dvX(k, l) dvY(k, l)];
for jj = ii+1:nn
m = nbh(jj);
d2 = [dvX(k, m) dvY(k, m)];
len(ii, jj) = dst(k, l) + dst(k, m);
ang(ii, jj) = abs(pi/2 - acos(dot(d1, d2) / (norm(d1) * norm(d2))));
end
end
% Find candidates with angle difference < 10 degree
cand = find(ang < pi/18);
% For these candidates, find the one with the shortest total distance
[~, I] = min(len(cand));
% Get corresponding indices
[I, J] = ind2sub([nn nn], cand(I));
cand = nbh([I J]);
% Lines
plot([X(k) X(cand(1))], [Y(k) Y(cand(1))], 'b', 'LineWidth', 1);
plot([X(k) X(cand(2))], [Y(k) Y(cand(2))], 'b', 'LineWidth', 1);
end
hold off;
输出图像：
近似（！）解决方案是确定所找到的点描述的轮廓的中心，并使用相对于中心的 atan2 按角度对点进行排序。请参阅以下代码片段以进行可视化：
% Points
X = 2 * rand(1, 15) - 1;
Y = 2 * rand(1, 15) - 1;
% Center
C = [0, 0];
% Determine indices
[~, idx] = sort(atan2(X, Y));
figure(1);
hold on;
% Center
plot(C(1), C(2), 'kx', 'MarkerSize', 15);
% Plain points
plot(X, Y, '.', 'MarkerSize', 15);
% Indices and lines
for k = 1:numel(X)
text(X(idx(k)) + 0.05, Y(idx(k)) + 0.05, num2str(k), 'FontSize', 12);
if (k == numel(X))
plot([X(idx(k)) X(idx(1))], [Y(idx(k)) Y(idx(1))], 'b');
else
plot([X(idx(k)) X(idx(k+1))], [Y(idx(k)) Y(idx(k+1))], 'b');
end
end
hold off;
给出以下输出：
虽然我确信一定数量的凹面会被正确处理，但恐怕对于给定的示例（尤其是上半部分）它会失败。这是因为图像不是完美的俯视图，因此角度有点“扭曲”。
不过，也许排序可以提高您的最小距离方法。
这是一个解决方案，适用于具有由垂直* 线（如您示例中的线）构成的轮廓的形状。思路如下：
- 我们旋转点以将*它们对齐到 XY 网格。
- 我们将具有相同* X 或 Y 坐标的点分组。
- 对于每个点，我们计算两个点：水平方向最近的点和垂直方向最近的点，来自允许的族。
- 构建连接矩阵并转换回来。
就像在
* - 大约.
function A = q55511236
%% Initialization:
% Define points:
X = [364.533 372.267 397.067 408.133 382.471 379.533 329.250 257.200 199.412 195.267 184.385 ...
168.643 157.533 174.500 108.533 99.333 150.733 184.800 138.105 179.474 218.278 232.133 ...
267.714 306.929 312.143 357.733 421.333 431.000 371.867];
Y = [192.027 233.360 228.627 286.693 314.541 292.960 327.450 340.500 348.671 326.693 269.308 ...
330.857 274.493 226.786 239.200 193.467 182.760 101.893 111.000 80.442 74.356 140.360 ...
64.643 56.857 77.786 69.493 133.293 180.427 142.160];
%% Preprocessing:
% Centering:
XY = [X;Y] - [mean(X); mean(Y)];
% Rotation:
[U,~,~] = svd(XY,'econ');
rXY = (U.' * XY).';
% Fixing problems w/ some points:
rXY = vertcat(rXY, [-21.8, 66]); % add missing point
rXY(7:8, :) = NaN; % remove non-corners
% figure(); scatter(rXY(:,1),rXY(:,2));
%% Processing:
% Group points according to same-X and same-Y
CLOSE_ENOUGH_DISTANCE = 10; % found using trial and error
[~,~,sameXpts] = uniquetol(rXY(:,1), CLOSE_ENOUGH_DISTANCE, 'DataScale', 1);
[~,~,sameYpts] = uniquetol(rXY(:,2), CLOSE_ENOUGH_DISTANCE, 'DataScale', 1);
% Create masks for distance evaluations:
nP = size(rXY,1);
[maskX,maskY] = deal(zeros(nP));
maskX(sameXpts == sameXpts.') = Inf;
maskY(sameYpts == sameYpts.') = Inf;
% Compute X and Y distances separately (we can do this in the rotated space)
dX = abs(rXY(:,1) - rXY(:,1).') + maskX + 1./maskY;
dY = abs(rXY(:,2) - rXY(:,2).') + maskY + 1./maskX;
[~,nX] = min(dX);
[~,nY] = min(dY);
% Construct connectivity matrix:
A = false(nP);
idxTrue = sub2ind(size(A), repmat(1:nP, [1,2]), [nX(:).', nY(:).']);
A(idxTrue) = true;
%% Plot result:
% Rotated coordinates:
figure(); gplot(A, rXY, '-o'); text(rXY(:,1), rXY(:,2), string(1:nP));
uXY = (U*rXY.').';
% Original coordinates:
figure(); gplot(A, uXY, '-o'); text(uXY(:,1), uXY(:,2), string(1:nP)); axis ij;
导致：
答案使用的概念是'Travelling salesman problem'。在这些点周围创建了一个缓冲区，这个缓冲区作为一个额外的标准包含在内。
a=[141 188 178 217 229 282 267 307 313 357 372 422 434 365 372 398 411 382 382 233 229 191 185 166 156 183 173 114 97 149 139 139];
b=[109 103 79 76 140 132 64 56 78 72 141 133 180 192 234 228 287 293 315 348 343 348 329 332 270 268 225 240 194 184 108 108];
X=[364.5333 232.1333 397.0667 157.5333 431 421.3333 306.9286 184.3846 357.7333 199.4118 168.6429 179.4737 408.1333 382.4706 150.7333 372.2667 184.8 138.1053 312.1429 108.5333 174.5 195.2667 257.2 99.33333 379.5333 371.8667 329.25 280.7059 267.7143 218.2778];
Y=[192.0267 140.36 228.6267 274.4933 180.4267 133.2933 56.85714 269.3077 69.49333 348.6706 330.8571 80.44211 286.6933 314.5412 182.76 233.36 101.8933 111 77.78571 239.2 226.7857 326.6933 340.5 193.4667 292.96 142.16 327.45 130.5529 64.64286 74.35556];
R = [a' b'];
d = 12;
polyout = polybuffer(R,'lines',d)
figure
%imshow(I2);
hold on
%plot(R(:,1),R(:,2),'r.','MarkerSize',10)
plot(X,Y,'r.', 'MarkerSize', 15)
plot(polyout)
axis equal
hold off
[s,t] = boundary(polyout); %%this is the boundary polygon of the buffer
numPoints = length(clustersCentroids);
x = X; %these are the points to be connected
y = Y;
x([1 2],:)=x([2 1],:);
y([1 2],:)=y([2 1],:);
figure
plot(x, y, 'bo', 'LineWidth', 2, 'MarkerSize', 17);
grid on;
imshow(I2);
xlabel('X', 'FontSize', 10);
ylabel('Y', 'FontSize', 10);
% Make a list of which points have been visited
beenVisited = false(1, numPoints);
% Make an array to store the order in which we visit the points.
visitationOrder = ones(1, numPoints);
% Define a filasafe
maxIterations = numPoints + 1;
iterationCount = 1;
% Define a current index. currentIndex will be 1 to start and then will vary.
currentIndex = 1;
while sum(beenVisited) < numPoints
visitationOrder(iterationCount) = currentIndex;
beenVisited(currentIndex) = true;
% Get the x and y of the current point.
thisX = x(currentIndex);
thisY = y(currentIndex);
%text(thisX + 0.01, thisY, num2str(currentIndex), 'FontSize', 35, 'Color', 'r');
% Compute distances to all other points
distances = sqrt((thisX - x) .^ 2 + (thisY - y) .^ 2);
distances(beenVisited)=inf;
distances(currentIndex) = inf;
% Don't consider visited points by setting their distance to infinity.
[out,idx] = sort(distances);
xx=[x(currentIndex) x(idx(1))]
yy=[y(currentIndex) y(idx(1))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(1);
else
xx=[x(currentIndex) x(idx(2))]
yy=[y(currentIndex) y(idx(2))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(2);
else
xx=[x(currentIndex) x(idx(3))]
yy=[y(currentIndex) y(idx(3))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(3);
else
xx=[x(currentIndex) x(idx(4))]
yy=[y(currentIndex) y(idx(4))]
if isempty(polyxpoly(xx,yy,s,t))
iterationCount = iterationCount + 1;
currentIndex =idx(4);
end
end
end
end
end
% Plot lines in that order.
hold on;
orderedX = [x(visitationOrder); x(1)];
orderedY = [y(visitationOrder) ;y(1)];
plot(orderedX,orderedY, 'm-', 'LineWidth', 2);
title('Result', 'FontSize', 10);