Сложный сценарий корпуса в Matlab с перекрестным произведением - PullRequest
0 голосов
/ 11 декабря 2018

Для школы мне нужно написать скрипт, чтобы получить выпуклую оболочку, заданную массивом 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

1 Ответ

0 голосов
/ 11 декабря 2018

Как указано в комментариях, ваш цикл не имеет правильного вложенного цикла 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*')
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...