Вычислительная геометрия (стр. 1 )

11. Вычислительная геометрия

Вычислительная геометрия – один из разделов информатики. Существует большое количество научных и прикладных задач, связанных с вычислительной геометрией.

Рассмотрим сначала типы данных, применяемых для вычислений. Точка на плоскости описывается парой вещественных чисел. При использовании вещественного типа операции сравнения лучше всего оформить специальными функциями, то есть сравнения вида A=B, где A и B вещественные числа, использовать не стоит. Приведем пример [4].

Пусть для хранения вещественного типа используются два десятичных разряда для порядка и шесть разрядов для хранения мантиссы. Есть два числа: A=0.34567´104 и B=0.98765´10-4. При сложении и вычитании происходит выравнивание порядков, сложение (вычитание) мантисс и нормализация результата. После выравнивания порядков получим B=0.0000000098765´104, а после сложения мантисс и будет результат:

A+B = 0. 3456700098765

При шести десятичных разрядах мантиссы после округления окажется, что A+B=A при B>0! В реальных системах тоже есть подобные ограничения на разрядность представления вещественного типа данных. Например, для типа Real в ПАСКАЛе мантиссе отводится 11-12 цифр.

Операцию проверки равенства вещественных чисел можно реализовать следующей функцией:

Function IsEqual(a, b: double): boolean;

Begin

IsEqual:=Abs(a-b)<Eps

{Eps – константа, определяющая точность,

Например, Eps=1e-10}

End;

Cравнения вида «больше» или «меньше» тоже дают ошибки. Операцию проверки «больше» можно реализовать так:

Function IsMore(a, b: double): boolean;

Begin

if IsEqual(a, b) then IsMore:=false

else

if a > b then IsMore:=true

else IsMore:=false;

End;

Остальные функции сравнения пишутся аналогично. Эти функции можно использовать, например, для определения максимального из двух значений.

Во многих задачах используются векторное и скалярное произведения векторов. Пусть, например, треугольник ABC задан координатами точек A(X1,Y1), B(X2,Y2), C(X3,Y3). Требуется:

·  найти площадь треугольника;

·  определить значение Cos A.

Вектора сторон имеют координаты AB (X2 - X1, Y2 - Y1) и AC (X3 - X1, Y3 - Y1). Ориентированная площадь треугольника проще всего находится через векторное произведение, обозначенное символом ‘´ ’:

S = (AB ´ AC) / 2 = ((X2 - X1 )* (Y3 - Y1 ) - (X3 - X1 )* (Y2 - Y1 )) / 2

Ориентация определяется знаком величины S.

Значение Cos A определяется из скалярного произведения векторов AB и AC:

Cos A =((X2 - X1 )* (Y3 - X1 )+ (X2 - X1 )* (Y3 - X1 )) / ( | AB | * | AC | )

Угол B прямой, если числитель (скалярное произведение) обращается в ноль.

Рассмотрим следующую задачу.

Четырехугольник. Заданы координаты четырех вершин четырехугольника в порядке обхода по часовой стрелке. Совпадающих вершин нет. Требуется определить, является ли четырехугольник выпуклым.

Ввод из файла INPUT. TXT. Первая строка содержит число тестовых блоков L (L ≤ 10). Каждый тестовый блок состоит из четырех строк, задающих координаты вершин четырехугольника в виде двух целых чисел X и Y (-1000 ≤ X, Y ≤1000), разделенных пробелом.

Вывод в файл OUTPUT. TXT. Для каждого тестового блока выводится единственная строка со значением Yes или No – является четырехугольник выпуклым или нет.

Пример

Ввод

2

0 0

-1 2

3 5

7 -1

1 0

0 5

5 0

2 1

Вывод

Yes

No

В выпуклом четырехугольнике каждые две противоположные вершины лежат по разные стороны относительно диагонали, соединяющей две другие вершины. Если же четырехугольник не является выпуклым, то относительно одной из диагоналей две вершины находятся по одну сторону. Если вывести уравнение диагонали по двум точкам в виде AX + BY +C = 0 (в частных случаях A = 0 или B = 0), то выражение Z = AX + BY +C обращается в 0 на диагонали и имеет разные знаки в точках по разные стороны диагонали. Остается для каждого четырехугольника вывести уравнения диагоналей и проверить для обеих диагоналей указанные условия, подставляя координаты вершин в выражение для Z.

Приведем текст программы для решения этой несложной задачи.

Program FourAng;

Const Max=100;

Var

L, i,j: integer;

X: array[1..4] of integer;

Y: array[1..4] of integer;

A, B,C: longint;

Fin, Fout: text;

Znak1,Znak2:integer;

Znak3,Znak4:integer;

Function WhatZnak(X1,Y1,X2,Y2,X, Y: integer):integer;

{возвращает знак при подстановке (X, Y) в уравнение прямой,

проходящей через точки (X1,Y1) и (X2,Y2):

1-больше 0, -1-меньше, 0-равно;

уравнение прямой через точки (X1,Y1) и (X2,Y2):

(Y2-Y1)*X+(X1-X2)*Y+X1(Y1-Y2)+Y1*(X2-X1)=0 }

Var

A, B,C, V: longint;

Begin

A:=Y2-Y1; B:=X1-X2;

C:=X1*(Y1-Y2)+Y1*(X2-X1);

V:=A*X+B*Y+C;

if V>0 then WhatZnak:=1

else if V<0 then WhatZnak:=-1

else WhatZnak:=0

End;

Begin

Assign(Fin,'input. txt');

Reset(Fin);

Assign(Fout,'output. txt');

Rewrite(Fout);

ReadLn(Fin, L); {число тестовых блоков}

For i:=1 to L do

begin

For j:=1 to 4 do ReadLn(Fin, X[j],Y[j]);

Znak2:=WhatZnak(X[1],Y[1],X[3],Y[3],X[2],Y[2]);

{знак второй вершины отоносительно первой диагонали}

Znak4:=WhatZnak(X[1],Y[1],X[3],Y[3],X[4],Y[4]);

{знак четвертой вершины отоносительно первой диагонали}

Znak1:=WhatZnak(X[2],Y[2],X[4],Y[4],X[1],Y[1]);

{знак первой вершины отоносительно второй диагонали}

Znak3:=WhatZnak(X[2],Y[2],X[4],Y[4],X[3],Y[3]);

{знак третьей вершины отоносительно второй диагонали}

if (Znak2=Znak4) or (Znak1=Znak3) then

WriteLn(Fout,'No') {невыпуклый}

else

WriteLn(Fout,'Yes') {выпуклый}

end;

Close(Fin);

Close(Fout);

End.

Рассмотрим еще одну базовую задачу вычислительной геометрии [2].

Выпуклая оболочка. На плоскости заданы своими декартовыми координатами N точек. Найти минимальный периметр многоугольника, содержащего все эти точки. Определить площадь этого многоугольника. Гарантируется, что искомый многоугольник имеет ненулевую площадь.

Ввод из файла INPUT. TXT. В первой строке находится число N, далее - N строк с парами координат.

Вывод в файл OUTPUT. TXT. Вывести с одним знаком после запятой: в первой строке - длину периметра, во второй - площадь.

Ограничения: 3 ≤ N ≤ 1000, -1000 ≤ xi, yi ≤ 1000, все числа целые, все точки различны, время 2 с.

Пример

Ввод

5

1 0

0 1

-1 0

0 -1

0 0

Вывод

5.7

2.0

Во-первых, искомый многоугольник должен быть выпуклым. Во-вторых, вершины многоугольника должны совпадать с какими-то из заданных точек. Оба эти свойства доказываются от противного. Данный многоугольник называется выпуклой оболочкой множества точек.

Опишем алгоритм Джарвиса как наиболее известный способ построения выпуклой оболочки [2-4, 6]. Сначала находится первая вершина многоугольника. Например, из множества самых нижних вершин выбирается самая левая. Найдем другие вершины в порядке обхода точек по часовой стрелке, предполагая сначала, что все точки различны.

Возьмем вектор (-1, 0) и обозначим его как (d1, d2). Вторая вершина выбирается так, чтобы угол между вектором (d1, d2) и вектором из первой вершины во вторую был минимален. Можно сказать, что вторая точка получается минимальным поворотом вектора (d1, d2) из первой точки по часовой стрелке. Если таких точек несколько, то возьмем наиболее удаленную из них.

Далее вектором (d1, d2) будем считать вектор из первой точки во вторую и будем находить следующую точку и т. д. Все продолжается, пока снова не придем в первую точку.

Минимальный угол можно найти с помощью скалярного произведения векторов. Ему соответствует максимальное значение косинуса, изменяющегося от -1 до 1. Поскольку некоторые точки могут повторяться, при поиске следующей вершины многоугольника не нужно исследовать точки, совпадающие с текущей вершиной.

Трудоемкость этого алгоритма O(N2), где N – количество точек. Иногда его называют алгоритмом обертки. Действительно, можно представить, что точки на плоскости соответствуют вбитым гвоздикам, вокруг которых натягивают ленту. Есть более быстрый алгоритм Грэхема [3, 4, 6], трудоемкость которого O(N´ logN).

Снова из множества самых нижних вершин выбирается самая левая. Обозначим ее A0. Точка A0 заведомо принадлежит выпуклой оболочке. В нее переносится центр координат. Остальные вершины сортируются одним из методов быстрой сортировки по возрастанию полярного угла. Если несколько точек имеют одинаковый угол, то оставляется одна с максимальным радиус-вектором. Обозначим оставшиеся точки A1, A2, …, AM в соответствии с порядком сортировки. Будем последовательно, т. е. в направлении против часовой стрелки рассматривать тройки точек Ak-1, Ak, Ak+1. С помощью векторного произведения найдем ориентированную площадь треугольника Ak-1AkAk+1. По знаку площади можно определить, находится ли точка Ak внутри треугольника A0Ak-1Ak+1. Если это так, то точка Ak исключается. Это не означает, что предыдущая не исключенная точка Ak-1 входит в выпуклую оболочку. Точка Ak-1 может находиться внутри треугольника A0Ak-2Ak+1, что приведет к ее исключению. Для возврата к предыдущим точкам можно использовать стек. Как только исключение промежуточных точек завершается, можно перейти к следующей еще не рассмотренной точке. По окончанию обхода оставшиеся точки будут углами многоугольника, являющегося выпуклой оболочкой.

Найдем площадь многоугольника. Пусть даны две точки A(X1,Y1) и B(X2,Y2). Площадь S трапеции, образованной отрезком AB, определяется как