matlab中带有叉积的复杂外壳脚本



在学校里,我必须写一个脚本,在给定x和y数组的情况下得到一个凸包。为了测试我的代码是否有效,我为x和y设置了一些值,并绘制了图。情节显然不是凸的,但我不确定我做错了什么。我看了叉积的公式,甚至用手计算了它,但我对特定向量的叉积得出了与matlab相同的答案。我研究了使用的公式不正确的可能性,这不太可能,因为我在两个独立的来源中发现了相同的公式。我复制了下面使用的代码。提前感谢!https://en.wikipedia.org/wiki/Graham_scan,截面算法下公式的来源。

x =[2 4 5 3 9 10]
y =[-10 8 4 -5 9 10] 
plot(x,y,'*')
[min_y,i] = min(y);
pivot = [x(i),y(i)];
number = i
%getting the angle
x(i) = [];
y(i) = [];
delta_x = x - ones(1,length(x))*pivot(1);
delta_y = y - ones(1,length(y))*pivot(2);
theta = atan2d(delta_y,delta_x);
%sorting without losing the position
c = [];
x_new= [];
y_new = [];
for n = 1:length(theta)
[min_theta, d] = min(theta);
c = [c d];
x_new = [x_new x(d)];
y_new = [y_new y(d)];
theta(d) = [1000]; %random value need to be higher than 180 degrees
end
x_new = [pivot(1) x_new pivot(1) x_new(1)];
y_new = [pivot(2) y_new pivot(2) y_new(1)];
%using cross product in this exercise none of the lines are collinear
%If the result is 0, the points are collinear;
%if it is positive, the three points constitute a "left turn" or counter-clockwise orientation,
%otherwise a "right turn" or clockwise orientation (for counter-clockwise numbered points)
%source wiki
n_bad = [];
n_good = [];
x_new
y_new
for n = 1:(length(x_new)-3) %this caused by how the array x_new and y_new contain duplicates
x_new(n+1)
cross_product = (x_new(n+1) - x_new(n))*(y_new(n+2)-y_new(n))-((y_new(n+1)-y_new(n))*(x_new(n+2)-x_new(n)))
if cross_product < 0
n_bad = [n_bad n+1];
end
if cross_product > 0 %this is unnecassaray can be used for skipping making the program faster
n_good = [n_good n+1]
end
end
%in case I did something wrong the graph can be plotted here:
figure
x_good = [pivot(1) x_new(n_good) pivot(1)];
y_good = [pivot(2) y_new(n_good) pivot(2)];
plot(x_good,y_good)
%to get the same index as the input
if number>1
n_good = n_good + 1;
n = sort([n_good number]);
else 
n = [number n_good];
end
n

正如注释中所建议的,您的循环没有正确的嵌套while循环。这是一个用你的修改过的脚本,对我有用。最后我在MATLAB自己的凸包函数上测试了它。

首先我添加了一个MATLAB函数

function [ value ] = ccw( p1,p2,p3,x,y )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
value=(x(p2)-x(p1))*(y(p3)-y(p1))-(y(p2)-y(p1))*(x(p3)-x(p1));
end

接下来是修改后的脚本

clear all
x =[2 4 5 3 9 10]
y =[-10 8 4 -5 9 10] 

x=[4,-3,8,12,4,-4,4,1,2,8,6] %additional test
y=[5,8,-6,1,4,7,8,-4,4,2,1]
[min_y,i] = min(y);
pivot = [x(i),y(i)];
number = i
%getting the angle
% x(i) = [];
% y(i) = [];
delta_x = x - ones(1,length(x))*pivot(1);
delta_y = y - ones(1,length(y))*pivot(2);
theta = atan2d(delta_y,delta_x);
%sorting without losing the position
c = [];
x_new= [];
y_new = [];
for n = 1:length(theta)
[min_theta, d] = min(theta);
c = [c d];
x_new = [x_new x(d)];
y_new = [y_new y(d)];
theta(d) = [1000]; %random value need to be higher than 180 degrees
end
x_new = [x_new x_new(1)];
y_new = [y_new y_new(1)];
%using cross product in this exercise none of the lines are collinear
%If the result is 0, the points are collinear;
%if it is positive, the three points constitute a "left turn" or counter-clockwise orientation,
%otherwise a "right turn" or clockwise orientation (for counter-clockwise numbered points)
%source wiki
n_bad = [];
n_good = [1,2];
figure(1)
plot(x_new,y_new,'-o')
for n = 3:(length(x_new)) %this caused by how the array x_new and y_new contain duplicates
%cross_product = (x_new(n+1) - x_new(n))*(y_new(n+2)-y_new(n))-((y_new(n+1)-y_new(n))*(x_new(n+2)-x_new(n)))
%     if cross_product < 0
%         n_bad = [n_bad n+1];
%     end
len=length(n_good);
while (ccw(n_good(len-1),n_good(len),n,x_new,y_new) < 0 && length(n_good)>=2) %this is unnecassaray can be used for skipping making the program faster
n_good(length(n_good))=[]
len=length(n_good);
end
n_good=[n_good n]
end
%in case I did something wrong the graph can be plotted here:
figure(2)
x_good = [pivot(1) x_new(n_good) pivot(1)];
y_good = [pivot(2) y_new(n_good) pivot(2)];
plot(x_good,y_good,'-o',x,y,'b*')
%to get the same index as the input
if number>1
n_good = n_good + 1;
n = sort([n_good number]);
else 
n = [number n_good];
end
figure(3)
k = convhull(x,y);
plot(x(k),y(k),'r-',x,y,'b*')

最新更新