Вычисления площади произвольного многоугольника (стр. 2 из 2)
Из геометрии известно, что если две точки находятся по одну сторону от прямой, то при подстановке их координат в полином левой части получим значения одного знака. Таким образом, если точки B1 и B2 располагаются по разные стороны от прямой, то
(AX1B+BY1B+C)(AX2B+BY2B+C)<0
Но пересечение прямых не является достаточным для пересечения отрезков, например:
 |
Эти отрезки не пересекаются.
Для достаточного доказательства пересечения отрезков необходимо произвести все вышеприведенные операции над точками B1 и B2, т.е. провести через них прямую и выяснить расположение точек A1 и A2 относительно ее.
Таким образом определяем возможность отсечения вершины многоугольника с количеством вершин, больше четырех.
Для некоторых видов четырехугольника это условие не несправедливо, например:
 | ||
 | ||
Здесь вершину A отсечь нельзя. Для четырехугольника определяем расположение отсекаемой вершины и вершины, несмежной с ней (через оставшиеся проходит "линия отреза"). Если они расположены по одну сторону, как на рисунке, то отсекать нельзя. Приведенный алгоритм контроля реализуем в следующей функции:
function cross(c: integer): boolean;
var a, b, i: integer;
AA, BB, CC,
AA1, BB1, CC1: real;
function Mline(x,y: real): real;
begin
Mline:=AA*x+BB*y+CC
end;
function Sline(x,y: real): real;
begin
Sline:=AA1*x+BB1*y+CC1
end;
begin
if c=1 then
begin
a:=n;
b:=2;
end
else if c=n then
begin
a:=n-1;
b:=1;
end
else
begin
a:=c-1;
b:=c+1;
end;
cross:=true;
AA:=sd[b].y-sd[a].y;
BB:=-(sd[b].x-sd[a].x);
CC:=sd[a].y*(sd[b].x-sd[a].x)-sd[a].x*(sd[b].y-sd[a].y);
if n=4 then
begin
for i:=1 to n do
if (Mline(sd[i].x, sd[i].y)*Mline(sd[c].x, sd[c].y)>0) and (i<>c)
then exit;
cross:=false;
exit
end;
for i:=1 to n-1 do
begin
AA1:=sd[i+1].y-sd[i].y;
BB1:=-(sd[i+1].x-sd[i].x);
CC1:=sd[i].y*(sd[i+1].x-sd[i].x)-sd[i].x*(sd[i+1].y-sd[i].y);
if Mline(sd[i].x, sd[i].y)*Mline(sd[i+1].x,sd[i+1].y)<0 then
if Sline(sd[a].x, sd[a].y)*Sline(sd[b].x,sd[b].y)<0 then exit;
end;
AA1:=sd[1].y-sd[n].y;
BB1:=-(sd[1].x-sd[n].x);
CC1:=sd[n].y*(sd[1].x-sd[n].x)-sd[n].x*(sd[1].y-sd[n].y);
if Mline(sd[n].x, sd[n].y)*Mline(sd[1].x,sd[1].y)<0 then
if Sline(sd[a].x, sd[a].y)*Sline(sd[b].x,sd[b].y)<0 then exit;
cross:=false;
end;
4) Вычисление площадей отсеченных треугольников будем по формуле Герона
где
.function St(x1,y1, x2,y2, x3,y3: real): real;
var a, b, c, p: real;
begin
a:=sqrt(sqr(x1-x2)+sqr(y1-y2));
b:=sqrt(sqr(x1-x3)+sqr(y1-y3));
c:=sqrt(sqr(x3-x2)+sqr(y3-y2));
p:=(a+b+c)/2;
St:=sqrt(p*(p-a)*(p-b)*(p-c));
end;
5) Отсечение i-той вершины
dec(n);
for j:=i to n do sd[j]:=sd[j+1];
После отсечения какой-либо вершины необходимо заново рассчитать внутренние углы многоугольника, т.е. вызвать процедуру Angles.
После вычисления площади, выводим ее на экран и ожидаем нажатия любой клавиши.
Writeln('Площадь фигуры: ', S:3:3);
Readkey
Полный текст программы приведен в приложении 2.
Проверка на контрольных примерах.
Проверим работу программы на фигурах, площади которых можно вычислить по формулам.
1) Треугольник
Содержимое файла points.dat
3
0 0
5 0
2 3
Площадь треугольника по формуле:
Результат работы программы:
Площадь фигуры: 7.500
2) Прямоугольник
Содержимое файла points.dat
4
0 0
5 0
5 3
0 3
Площадь прямоугольника по формуле
Результат работы программы
Площадь фигуры: 15.000
 |
3) Невыпуклая фигура
Содержимое файла points.dat
4
0 0
3 2
2 0
3 -2
Эта фигура симметрична относительно оси OX. Ее площадь будет равна
Результат работы программы:
Площадь фигуры: 4.000
Заключение.
Рассмотренный алгоритм является комбинацией аналитического и численного методов. Поэтому он включает в себя достоинства обоих. Т.к. основной операцией вычисления площади многоугольника является вычисление площади элементарного треугольника, то на результат вычисления не будут влиять методические погрешности, т.е. погрешности вносимые самим алгоритмом. Этим приведенный алгоритм отличается от метода Монте-Карло, точность которого зависит от количества точек. Погрешность будет вноситься лишь на этапе вычислений и будет определяться конкретной ЭВМ, на которой ведется расчет. Точность зависит от вещественного типа Real, в котором представляются основные переменные. Данный тип представим в компьютере лишь с определенной точностью, зависящей от внутреннего формата числа. Для персонального компьютера типа IBM PC/AT тип Real имеет следующие параметры:
· Длинна, байт................................. 6
· Количество значащих цифр......... 11…12
· Диапазон десятичного порядка... -39…+38
Эта точность вполне удовлетворительна для нашей задачи.
Приложение 1.
Блок-схема
 |
Приложение 2. Текст программы.
Uses Crt;
const max=100;
var
n, i, j: integer;
sd: array[1..100] of
record
x,y: real;
angle: real;
end;
S: real;
procedure Angles;
var
al1,al2,
dx, dy, dxp, dyp,
s_in, s_out, a: real;
i,j: integer;
function ArcCos(a: real): real;
var res: real;
begin
if abs(a)<1.0E-30 then res:=pi/2
else res:=ArcTan(sqrt(1-a*a)/a);
if dx<0 then
if dy>=0 then res:=pi+res
else res:=-pi-res
else
if dy<0 then res:=-res;
ArcCos:=res
end;
begin
dxp:=sd[1].x-sd[n].x;
dyp:=sd[1].y-sd[n].y;
a:=sqrt(dxp*dxp+dyp*dyp);
dxp:=dxp/a;
dyp:=dyp/a;
i:=1;
while i<=(n-1) do
begin
dx:=sd[i+1].x-sd[i].x;
dy:=sd[i+1].y-sd[i].y;
a:=sqrt(dx*dx+dy*dy);
dx:=dx/a;
dy:=dy/a;
if (dx=dxp) and (dy=dyp) then
begin
dec(n);
for j:=i to n do sd[j]:=sd[j+1];
end;
dxp:=dx; dyp:=dy;
inc(i)
end;
dx:=sd[1].x-sd[n].x;
dy:=sd[1].y-sd[n].y;
al1:=ArcCos(dx/sqrt(dx*dx+dy*dy));
for i:=1 to n-1 do
begin
dx:=sd[i+1].x-sd[i].x;
dy:=sd[i+1].y-sd[i].y;
al2:=ArcCos(dx/sqrt(dx*dx+dy*dy));
sd[i].angle:=pi-al1+al2;
if sd[i].angle>2*pi then sd[i].angle:=sd[i].angle-2*pi
else
if sd[i].angle<0 then sd[i].angle:=2*pi+sd[i].angle;
al1:=al2
end;
dx:=sd[1].x-sd[n].x;
dy:=sd[1].y-sd[n].y;
al2:=ArcCos(dx/sqrt(dx*dx+dy*dy));
sd[n].angle:=pi-al1+al2;
if sd[n].angle>2*pi then sd[n].angle:=sd[n].angle-2*pi
else
if sd[n].angle<0 then sd[n].angle:=2*pi+sd[n].angle;
s_in:=0;
s_out:=0;
for i:=1 to n do
begin
if sd[i].angle<0 then sd[i].angle:=2*pi+sd[i].angle;
S_in:=S_in+sd[i].angle;
S_out:=S_out+(2*pi-sd[i].angle);
end;
if S_out<S_in then
for i:=1 to n do sd[i].angle:=2*pi-sd[i].angle;
end;
procedure input;
var f: text;
i: integer;
begin
Assign(f,'points.dat');
reset(f);
readln(f, n);
for i:=1 to n do readln(f, sd[i].x, sd[i].y);
end;
function St(x1,y1, x2,y2, x3,y3: real): real;
var a, b, c, p: real;
begin
a:=sqrt(sqr(x1-x2)+sqr(y1-y2));
b:=sqrt(sqr(x1-x3)+sqr(y1-y3));
c:=sqrt(sqr(x3-x2)+sqr(y3-y2));
p:=(a+b+c)/2;
St:=sqrt(p*(p-a)*(p-b)*(p-c));
end;
function cross(c: integer): boolean;
var a, b, i: integer;
AA, BB, CC,
AA1, BB1, CC1: real;
function Mline(x,y: real): real;
begin
Mline:=AA*x+BB*y+CC
end;
function Sline(x,y: real): real;
begin
Sline:=AA1*x+BB1*y+CC1
end;
begin
if c=1 then
begin
a:=n;
b:=2;
end
else if c=n then
begin
a:=n-1;
b:=1;
end
else
begin
a:=c-1;
b:=c+1;
end;
cross:=true;
AA:=sd[b].y-sd[a].y;
BB:=-(sd[b].x-sd[a].x);
CC:=sd[a].y*(sd[b].x-sd[a].x)-sd[a].x*(sd[b].y-sd[a].y);
if n=4 then
begin
for i:=1 to n do
if (Mline(sd[i].x, sd[i].y)*Mline(sd[c].x, sd[c].y)>0) and (i<>c) then exit;
cross:=false;
exit
end;
for i:=1 to n-1 do
begin
AA1:=sd[i+1].y-sd[i].y;
BB1:=-(sd[i+1].x-sd[i].x);
CC1:=sd[i].y*(sd[i+1].x-sd[i].x)-sd[i].x*(sd[i+1].y-sd[i].y);
if Mline(sd[i].x, sd[i].y)*Mline(sd[i+1].x,sd[i+1].y)<0 then
if Sline(sd[a].x, sd[a].y)*Sline(sd[b].x,sd[b].y)<0 then exit;
end;
AA1:=sd[1].y-sd[n].y;
BB1:=-(sd[1].x-sd[n].x);
CC1:=sd[n].y*(sd[1].x-sd[n].x)-sd[n].x*(sd[1].y-sd[n].y);
if Mline(sd[n].x, sd[n].y)*Mline(sd[1].x,sd[1].y)<0 then
if Sline(sd[a].x, sd[a].y)*Sline(sd[b].x,sd[b].y)<0 then exit;
cross:=false;
end;
begin
ClrScr;
input;
Angles;
S:=0;
while n>3 do
begin
i:=1;
while (sd[i].angle>pi) or (cross(i)) do
inc(i);
if i=1 then
S:=S+St(sd[1].x,sd[1].y, sd[2].x,sd[2].y, sd[n].x,sd[n].y)
else
if i=n then
S:=S+St(sd[n].x,sd[n].y, sd[1].x,sd[1].y, sd[n-1].x,sd[n-1].y)
else S:=S+St(sd[i].x,sd[i].y, sd[i-1].x,sd[i-1].y, sd[i+1].x,sd[i+1].y);
dec(n);
for j:=i to n do sd[j]:=sd[j+1];
Angles
end;
S:=S+St(sd[1].x,sd[1].y, sd[2].x,sd[2].y, sd[3].x,sd[3].y);
Writeln('Площадь фигуры: ', S:3:3);
Readkey
end.