Площадь самопересекающегося многоугольника

32

Рассмотрим потенциально самопересекающийся многоугольник, определенный списком вершин в двумерном пространстве. Например

{{0, 0}, {5, 0}, {5, 4}, {1, 4}, {1, 2}, {3, 2}, {3, 3}, {2, 3}, {2, 1}, {4, 1}, {4, 5}, {0, 5}}

Есть несколько способов определить площадь такого многоугольника, но наиболее интересным является правило четного нечетного. Взяв любую точку на плоскости, проведите линию от точки до бесконечности (в любом направлении). Если эта линия пересекает многоугольник нечетное число раз, точка является частью области многоугольника, если она пересекает многоугольник четное число раз, точка не является частью многоугольника. Для приведенного выше примера многоугольника здесь представлены как его контур, так и четно-нечетная область:

КонтурПлощадь

Многоугольник в общем случае не будет ортогональным. Я выбрал только такой простой пример, чтобы было легче подсчитать площадь.

Область этого примера 17(нет 24или 33как могут дать другие определения или область).

Обратите внимание, что согласно этому определению площадь многоугольника не зависит от порядка его намотки.

Соревнование

Учитывая список вершин с целочисленными координатами, определяющими многоугольник, определите его площадь по четно-нечетному правилу.

Вы можете написать функцию или программу, используя ввод через STDIN или ближайшую альтернативу, аргумент командной строки или аргумент функции, и либо вернуть результат, либо распечатать его в STDOUT или ближайшую альтернативу.

Вы можете принимать ввод в любом удобном формате списка или строки, если он не был предварительно обработан.

Результатом должно быть либо число с плавающей запятой с точностью до 6 значащих (десятичных) цифр, либо рациональный результат, представление которого с плавающей запятой с точностью до 6 значащих цифр. (Если вы дадите рациональные результаты, они, вероятно, будут точными, но я не могу требовать этого, поскольку у меня нет точных результатов для справки.)

Вы должны быть в состоянии решить каждый из приведенных ниже тестовых случаев в течение 10 секунд на подходящем настольном компьютере. (В этом правиле есть некоторая свобода действий, поэтому примите во внимание свое мнение. Если на моем ноутбуке это займет 20 секунд, я извлеку выгоду из сомнений, если это займет минуту, я не буду.) Я думаю, что это ограничение должно быть очень щедрым, но предполагается исключить подходы, в которых вы просто дискретизируете многоугольник на достаточно тонкой сетке и рассчитываете, или используете вероятностные подходы, такие как Монте-Карло. Будьте хорошим спортсменом и не пытайтесь оптимизировать эти подходы так, чтобы вы все равно могли уложиться в сроки. ;)

Вы не должны использовать какие-либо существующие функции, связанные непосредственно с полигонами.

Это код гольф, поэтому выигрывает самое короткое представление (в байтах).

Предположения

  • Все координаты являются целыми числами в диапазоне 0 ≤ x ≤ 100, 0 ≤ y ≤ 100.
  • Там будет хотя бы 3и не больше 50вершин.
  • Повторных вершин не будет. Ни одна из вершин не будет лежать на другом ребре. ( Хотя в списке могут быть коллинеарные точки.)

Тестовые случаи

{{0, 0}, {5, 0}, {5, 4}, {1, 4}, {1, 2}, {3, 2}, {3, 3}, {2, 3}, {2, 1}, {4, 1}, {4, 5}, {0, 5}}
17.0000

{{22, 87}, {6, 3}, {98, 77}, {20, 56}, {96, 52}, {79, 34}, {46, 78}, {52, 73}, {81, 85}, {90, 43}}
2788.39

{{90, 43}, {81, 85}, {52, 73}, {46, 78}, {79, 34}, {96, 52}, {20, 56}, {98, 77}, {6, 3}, {22, 87}}
2788.39

{{70, 33}, {53, 89}, {76, 35}, {14, 56}, {14, 47}, {59, 49}, {12, 32}, {22, 66}, {85, 2}, {2, 81},
 {61, 39}, {1, 49}, {91, 62}, {67, 7}, {19, 55}, {47, 44}, {8, 24}, {46, 18}, {63, 64}, {23, 30}}
2037.98

{{42, 65}, {14, 59}, {97, 10}, {13, 1}, {2, 8}, {88, 80}, {24, 36}, {95, 94}, {18, 9}, {66, 64},
 {91, 5}, {99, 25}, {6, 66}, {48, 55}, {83, 54}, {15, 65}, {10, 60}, {35, 86}, {44, 19}, {48, 43},
 {47, 86}, {29, 5}, {15, 45}, {75, 41}, {9, 9}, {23, 100}, {22, 82}, {34, 21}, {7, 34}, {54, 83}}
3382.46

{{68, 35}, {43, 63}, {66, 98}, {60, 56}, {57, 44}, {90, 52}, {36, 26}, {23, 55}, {66, 1}, {25, 6},
 {84, 65}, {38, 16}, {47, 31}, {44, 90}, {2, 30}, {87, 40}, {19, 51}, {75, 5}, {31, 94}, {85, 56},
 {95, 81}, {79, 80}, {82, 45}, {95, 10}, {27, 15}, {18, 70}, {24, 6}, {12, 73}, {10, 31}, {4, 29},
 {79, 93}, {45, 85}, {12, 10}, {89, 70}, {46, 5}, {56, 67}, {58, 59}, {92, 19}, {83, 49}, {22,77}}
3337.62

{{15, 22}, {71, 65}, {12, 35}, {30, 92}, {12, 92}, {97, 31}, {4, 32}, {39, 43}, {11, 40}, 
 {20, 15}, {71, 100}, {84, 76}, {51, 98}, {35, 94}, {46, 54}, {89, 49}, {28, 35}, {65, 42}, 
 {31, 41}, {48, 34}, {57, 46}, {14, 20}, {45, 28}, {82, 65}, {88, 78}, {55, 30}, {30, 27}, 
 {26, 47}, {51, 93}, {9, 95}, {56, 82}, {86, 56}, {46, 28}, {62, 70}, {98, 10}, {3, 39}, 
 {11, 34}, {17, 64}, {36, 42}, {52, 100}, {38, 11}, {83, 14}, {5, 17}, {72, 70}, {3, 97}, 
 {8, 94}, {64, 60}, {47, 25}, {99, 26}, {99, 69}}
3514.46
Мартин Эндер
источник
1
В частности, я хотел бы заменить разделители таким образом, чтобы сделать список допустимым PostScript User Path, чтобы я мог проанализировать все это с помощью одного upathоператора. (На самом деле это чрезвычайно простое преобразование 1: 1 между разделителями. }, {Просто становится lineto, и запятая между x и y удаляется, а открывающие и закрывающие скобки заменяются статическим
верхним
1
@AJMansfield Я обычно не возражаю против использования удобных родных представлений списков, но использование upathи linetoзвучит так, будто вы на самом деле предварительно обрабатываете ввод. Т.е. вы берете не список координат, а фактический многоугольник.
Мартин Эндер
1
@MattNoonan О, это хороший момент. Да, вы можете предположить, что.
Мартин Эндер
2
@Ray Хотя направление может влиять на количество пересечений, оно будет только увеличиваться или уменьшаться на 2, сохраняя паритет. Я постараюсь найти ссылку. Для начала SVG использует то же определение.
Мартин Эндер
1
Mathematica 12,0 имеет новый встроенный в функции для этого: CrossingPolygon.
алефальфа

Ответы:

14

Mathematica, 247 225 222

p=Partition[#,2,1,1]&;{a_,b_}~r~{c_,d_}=Det/@{{a-c,c-d},{a,c}-b}/Det@{a-b,c-d};f=Abs@Tr@MapIndexed[Det@#(-1)^Tr@#2&,p[Join@@MapThread[{1-#,#}&/@#.#2&,{Sort/@Cases[{s_,t_}/;0<=s<=1&&0<=t<=1:>s]/@Outer[r,#,#,1],#}]&@p@#]]/2&

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

введите описание изображения здесь

Пример:

In[2]:= f[{{15, 22}, {71, 65}, {12, 35}, {30, 92}, {12, 92}, {97, 31}, {4, 32}, {39, 43}, {11, 40}, 
 {20, 15}, {71, 100}, {84, 76}, {51, 98}, {35, 94}, {46, 54}, {89, 49}, {28, 35}, {65, 42}, 
 {31, 41}, {48, 34}, {57, 46}, {14, 20}, {45, 28}, {82, 65}, {88, 78}, {55, 30}, {30, 27}, 
 {26, 47}, {51, 93}, {9, 95}, {56, 82}, {86, 56}, {46, 28}, {62, 70}, {98, 10}, {3, 39}, 
 {11, 34}, {17, 64}, {36, 42}, {52, 100}, {38, 11}, {83, 14}, {5, 17}, {72, 70}, {3, 97}, 
 {8, 94}, {64, 60}, {47, 25}, {99, 26}, {99, 69}}]

Out[2]= 3387239559852305316061173112486233884246606945138074528363622677708164\
 6419838924305735780894917246019722157041758816629529815853144003636562\
 9161985438389053702901286180223793349646170997160308182712593965484705\
 3835036745220226127640955614326918918917441670126958689133216326862597\
 0109115619/\
 9638019709367685232385259132839493819254557312303005906194701440047547\
 1858644412915045826470099500628074171987058850811809594585138874868123\
 9385516082170539979030155851141050766098510400285425157652696115518756\
 3100504682294718279622934291498595327654955812053471272558217892957057\
 556160

In[3]:= N[%] (*The numerical value of the last output*)

Out[3]= 3514.46
alephalpha
источник
К сожалению, я не уверен, что эта логика будет работать для всех ситуаций. Ты можешь попробовать {1,2},{4,4},{4,2},{2,4},{2,1},{5,3}? Вы должны выйти с 3.433333333333309. Я посмотрел на использование аналогичной логики.
MickyT
@ MickyT Да, это работает. Возвращено 103/30, и числовое значение равно 3.43333.
алефальфа
Прости за это. Хорошее решение
MickyT
44

Python 2, 323 319 байт

exec u"def I(s,a,b=1j):c,d=s;d-=c;c-=a;e=(d*bX;return e*(0<=(b*cX*e<=e*e)and[a+(d*cX*b/e]or[]\nE=lambda p:zip(p,p[1:]+p);S=sorted;P=E(input());print sum((t-b)*(r-l)/2Fl,r@E(S(i.realFa,b@PFe@PFi@I(e,a,b-a)))[:-1]Fb,t@E(S(((i+j)XFe@PFi@I(e,l)Fj@I(e,r)))[::2])".translate({70:u" for ",64:u" in ",88:u".conjugate()).imag"})

Принимает список вершин через STDIN как комплексные числа в следующем виде

[  X + Yj,  X + Yj,  ...  ]

и записывает результат в STDOUT.

Тот же код после замены строки и некоторый интервал:

def I(s, a, b = 1j):
    c, d = s; d -= c; c -= a;
    e = (d*b.conjugate()).imag;
    return e * (0 <= (b*c.conjugate()).imag * e <= e*e) and \
           [a + (d*c.conjugate()).imag * b/e] or []

E = lambda p: zip(p, p[1:] + p);
S = sorted;

P = E(input());

print sum(
    (t - b) * (r - l) / 2

    for l, r in E(S(
        i.real for a, b in P for e in P for i in I(e, a, b - a)
    ))[:-1]

    for b, t in E(S(
        ((i + j).conjugate()).imag for e in P for i in I(e, l) for j in I(e, r)
    ))[::2]
)

объяснение

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

Рисунок 1

(На самом деле, из-за игры в гольф программа пропускает еще несколько линий; это не имеет значения, если мы проходим хотя бы эти линии.) Тело многоугольника между любыми двумя последовательными линиями состоит из вертикальных трапеций ( и треугольники, и отрезки, как особые случаи). Это должно быть так, поскольку, если бы у какой-либо из этих фигур была дополнительная вершина между двумя основаниями, через эту точку была бы другая вертикальная линия между двумя рассматриваемыми линиями. Сумма площадей всех таких трапеций представляет собой площадь многоугольника.

Вот как мы находим эти трапеции: для каждой пары последовательных вертикальных линий мы находим сегменты каждой стороны многоугольника, которые (должным образом) лежат между этими двумя линиями (которые могут не существовать для некоторых сторон). На приведенном выше рисунке это шесть красных сегментов при рассмотрении двух красных вертикальных линий. Обратите внимание, что эти сегменты не пересекаются должным образом (т. Е. Они могут встречаться только в конечных точках, полностью совпадать или вообще не пересекаться, поскольку, опять же, если они правильно пересекаются, между ними будет другая вертикальная линия;) и поэтому имеет смысл поговорить о порядке их сверху вниз, что мы и делаем. Согласно четно-нечетному правилу, когда мы пересекаем первый сегмент, мы находимся внутри многоугольника; как только мы пересекаем второй, мы уходим; третий снова; четвертый, вне; и так далее...

В целом, это алгоритм O ( n 3 log n ).

флигель
источник
4
Это великолепно! Я знал, что могу рассчитывать на тебя в этом. ;) (Возможно, вы захотите ответить на этот вопрос более подробно о переполнении стека.)
Мартин Эндер
@ MartinBüttner Держите их в пути :)
Ell
7
Отличная работа и отличное объяснение
MickyT
1
Это впечатляющий ответ. Вы сами разработали алгоритм или уже есть работа над этой проблемой? Если есть существующая работа, я был бы признателен за указатель, где я могу ее найти. Я понятия не имел, как с этим справиться.
Логический рыцарь
5
@CarpetPython Я разработал его сам, но я был бы очень удивлен, если бы это не было сделано раньше.
Ell
9

Хаскелл, 549

Не похоже, что я могу сыграть в эту игру достаточно далеко, но концепция отличалась от двух других ответов, поэтому я решил, что все равно поделюсь этим. Он выполняет O (N ^ 2) рациональных операций для вычисления площади.

import Data.List
_%0=2;x%y=x/y
h=sort
z f w@(x:y)=zipWith f(y++[x])w
a=(%2).sum.z(#);(a,b)#(c,d)=b*c-a*d
(r,p)?(s,q)=[(0,p)|p==q]++[(t,v t p r)|u t,u$f r]where f x=(d q p#x)%(r#s);t=f s;u x=x^2<x
v t(x,y)(a,b)=(x+t*a,y+t*b);d=v(-1)
s x=zip(z d x)x
i y=h.(=<<y).(?)=<<y
[]!x=[x];x!_=x
e n(a@(x,p):y)|x>0=(n!y,a):(e(n!y)$tail$dropWhile((/=p).snd)y)|0<1=(n,a):e n y
c[p]k=w x[]where((_,q):x)=e[]p;w((n,y):z)b|q==y=(k,map snd(q:b)):c n(-k)|0<1=w z(y:b);c[]_=[]
b(s,p)=s*a p
u(_,x)(_,y)=h x==h y
f p=abs$sum$map b$nubBy u$take(length p^2)$c[cycle$i$s p]1

Пример:

λ> f test''
33872395598523053160611731124862338842466069451380745283636226777081646419838924305735780894917246019722157041758816629529815853144003636562916198543838905370290128618022379334964617099716030818271259396548470538350367452202261276409556143269189189174416701269586891332163268625970109115619 % 9638019709367685232385259132839493819254557312303005906194701440047547185864441291504582647009950062807417198705885081180959458513887486812393855160821705399790301558511410507660985104002854251576526961155187563100504682294718279622934291498595327654955812053471272558217892957057556160
λ> fromRational (f test'')
3514.4559380388832

Идея состоит в том, чтобы заново соединить многоугольник при каждом пересечении, что приведет к объединению многоугольников без пересекающихся ребер. Затем мы можем вычислить (подписанную) область каждого из многоугольников, используя формулу Гаусса для шнурков ( http://en.wikipedia.org/wiki/Shoelace_formula ). Четно-нечетное правило требует, чтобы при пересечении пересечения площадь нового многоугольника считалась отрицательно относительно старого многоугольника.

Например, рассмотрим многоугольник в исходном вопросе. Пересечение в верхнем левом углу преобразуется в два пути, которые встречаются только в точке; оба пути ориентированы по часовой стрелке, поэтому области каждого будут положительными, если мы не объявим, что внутренний путь взвешен на -1 относительно внешнего пути. Это эквивалентно изменению пути альфаальфы.

Полигоны, полученные из оригинального примера

В качестве другого примера рассмотрим многоугольник из комментария MickyT:

Полигоны, полученные из комментария MickyT

Здесь некоторые полигоны ориентированы по часовой стрелке, а некоторые против часовой стрелки. Правило «знак переворачивания» гарантирует, что области, ориентированные по часовой стрелке, получают дополнительный коэффициент -1, что приводит к тому, что они вносят положительный вклад в площадь.

Вот как работает программа:

import Data.List  -- for sort and nubBy

-- Rational division, with the unusual convention that x/0 = 2
_%0=2;x%y=x/y

-- Golf
h=sort

-- Define a "cyclic zipWith" operation. Given a list [a,b,c,...z] and a binary
-- operation (@), z (@) [a,b,c,..z] computes the list [b@a, c@b, ..., z@y, a@z]
z f w@(x:y)=zipWith f(y++[x])w

-- The shoelace formula for the signed area of a polygon
a=(%2).sum.z(#)

-- The "cross-product" of two 2d vectors, resulting in a scalar.
(a,b)#(c,d)=b*c-a*d

-- Determine if the line segment from p to p+r intersects the segment from
-- q to q+s.  Evaluates to the singleton list [(t,x)] where p + tr = x is the
-- point of intersection, or the empty list if there is no intersection. 
(r,p)?(s,q)=[(0,p)|p==q]++[(t,v t p r)|u t,u$f r]where f x=(d q p#x)%(r#s);t=f s;u x=x^2<x

-- v computes an affine combination of two vectors; d computes the difference
-- of two vectors.
v t(x,y)(a,b)=(x+t*a,y+t*b);d=v(-1)

-- If x is a list of points describing a polygon, s x will be the list of
-- (displacement, point) pairs describing the edges.
s x=zip(z d x)x

-- Given a list of (displacement, point) pairs describing a polygon's edges,
-- create a new polygon which also has a vertex at every point of intersection.
-- Mercilessly golfed.
i y=h.(=<<y).(?)=<<y


-- Extract a simple polygon; when an intersection point is reached, fast-forward
-- through the polygon until we return to the same point, then continue.  This
-- implements the edge rewiring operation. Also keep track of the first
-- intersection point we saw, so that we can process that polygon next and with
-- opposite sign.
[]!x=[x];x!_=x
e n(a@(x,p):y)|x>0=(n!y,a):(e(n!y)$tail$dropWhile((/=p).snd)y)|0<1=(n,a):e n y

-- Traverse the polygon from some arbitrary starting point, using e to extract
-- simple polygons marked with +/-1 weights.
c[p]k=w x[]where((_,q):x)=e[]p;w((n,y):z)b|q==y=(k,map snd(q:b)):c n(-k)|0<1=w z(y:b);c[]_=[]

-- If the original polygon had N vertices, there could (very conservatively)
-- be up to N^2 points of intersection.  So extract N^2 polygons using c,
-- throwing away duplicates, and add up the weighted areas of each polygon.
b(s,p)=s*a p
u(_,x)(_,y)=h x==h y
f p=abs$sum$map b$nubBy u$take(length p^2)$c[cycle$i$s p]1
Мэтт Нунан
источник