Доработал код, теперь считает и вертикальные стороны, убрал цикл - сделал все по-Матлабовски векторно. Кому интересно - у меня решается система с 2 неизвестными (координаты пересечения диагоналей, остальное выражается через них). Матлабовский синтаксис некорректно отображает мой код, поэтому в теге "code".
Код:
function Main()
clear all; close all;
P = [-11, -7; -20, 13; 4, 10; 1, -8];
% P = [-2, 0; 0, 1; 2, 0; 0, -1];
% P = [-1, -2; -1, 2; 1, 2; 1, -2];
function r = diagonal_conditions(p0)
a = p0(1) - P(:, 1);
b = p0(2) - P(:, 2);
c = a.*P(:, 1) + b.*P(:, 2);
ac = circshift(a, -1);
bc = circshift(b, -1);
cc = circshift(c, -1);
z = a.*bc - ac.*b;
pv = [(bc.*c - b.*cc)./z, (cc.*a - c.*ac)./z];
r = [p0(2)*(pv(3,1) - pv(1,1)) - p0(1)*(pv(3,2) - pv(1,2)) - pv(1,2)*pv(3,1) + pv(1,1)*pv(3,2);
p0(2)*(pv(4,1) - pv(2,1)) - p0(1)*(pv(4,2) - pv(2,2)) - pv(2,2)*pv(4,1) + pv(2,1)*pv(4,2)];
end
pv = zeros(4, 2);
P0 = fsolve(@diagonal_conditions, [mean(P(:, 1)), mean(P(:, 2))]);
figure(1); axis on, axis equal, grid off, hold on
plot(P(:, 1), P(:, 2), 'ok', 'LineWidth', 2);
plot(P0(1), P0(2), 'or', 'LineWidth', 4);
pv = [pv; pv(1, :)];
plot(pv(:, 1), pv(:, 2), '-ob');
plot([pv(1, 1), pv(3, 1)], [pv(1, 2), pv(3, 2)], '-g');
plot([pv(2, 1), pv(4, 1)], [pv(2, 2), pv(4, 2)], '-g');
plot([P0(1), P(1, 1)], [P0(2), P(1, 2)], '-k');
plot([P0(1), P(2, 1)], [P0(2), P(2, 2)], '-k');
plot([P0(1), P(3, 1)], [P0(2), P(3, 2)], '-k');
plot([P0(1), P(4, 1)], [P0(2), P(4, 2)], '-k');
end