Fläche eines sich selbst schneidenden Polygons

32

Betrachten Sie ein sich möglicherweise selbst schneidendes Polygon, das durch eine Liste von Scheitelpunkten im 2D-Raum definiert wird. Z.B

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

Es gibt mehrere Möglichkeiten, die Fläche eines solchen Polygons zu definieren, aber die interessanteste ist die Gerade-Ungerade-Regel. Zeichnen Sie für einen beliebigen Punkt in der Ebene eine Linie vom Punkt bis ins Unendliche (in eine beliebige Richtung). Wenn diese Linie das Polygon ungerade oft schneidet, ist der Punkt Teil der Fläche des Polygons. Wenn sie das Polygon gerade oft schneidet, ist der Punkt nicht Teil des Polygons. Für das obige Beispielpolygon sind hier sowohl der Umriss als auch der gerade und ungerade Bereich:

GliederungBereich

Das Polygon ist im Allgemeinen nicht orthogonal. Ich habe nur ein so einfaches Beispiel gewählt, um das Zählen der Fläche zu erleichtern.

Die Fläche in diesem Beispiel ist 17(nicht 24oder 33wie andere Definitionen oder Flächen ergeben könnten).

Beachten Sie, dass unter dieser Definition die Fläche des Polygons unabhängig von seiner Wicklungsreihenfolge ist.

Die Herausforderung

Bestimmen Sie anhand einer Liste von Eckpunkten mit ganzzahligen Koordinaten, die ein Polygon definieren, dessen Fläche nach der Gerade-Ungerade-Regel.

Sie können eine Funktion oder ein Programm schreiben, Eingaben über STDIN oder die nächstgelegene Alternative, ein Befehlszeilenargument oder ein Funktionsargument vornehmen und das Ergebnis entweder zurückgeben oder an STDOUT oder die nächstgelegene Alternative ausgeben.

Sie können Eingaben in beliebigen Listen- oder Zeichenfolgenformaten vornehmen, sofern diese nicht vorverarbeitet wurden.

Das Ergebnis sollte entweder eine Gleitkommazahl sein, die auf 6 signifikante (Dezimal-) Stellen genau ist, oder ein rationales Ergebnis, dessen Gleitkommadarstellung auf 6 signifikante Stellen genau ist. (Wenn Sie rationale Ergebnisse erzielen, sind diese wahrscheinlich genau, aber ich kann dies nicht verlangen, da ich keine genauen Ergebnisse als Referenz habe.)

Sie müssen in der Lage sein, jeden der folgenden Testfälle innerhalb von 10 Sekunden auf einem vernünftigen Desktop-Computer zu lösen. (In dieser Regel gibt es einen gewissen Spielraum. Wenn es auf meinem Laptop 20 Sekunden dauert, gebe ich Ihnen den Vorteil des Zweifels, wenn es eine Minute dauert, werde ich es nicht tun.) Ich denke, diese Grenze sollte sehr großzügig sein, aber es soll Ansätze ausschließen, bei denen Sie das Polygon nur auf einem ausreichend feinen Raster diskretisieren und zählen, oder probabilistische Ansätze wie Monte Carlo verwenden. Seien Sie ein guter Sportler und versuchen Sie nicht, diese Ansätze so zu optimieren, dass Sie das Zeitlimit trotzdem einhalten können. ;)

Sie dürfen keine vorhandenen Funktionen verwenden, die sich direkt auf Polygone beziehen.

Dies ist Codegolf, daher gewinnt die kürzeste Übermittlung (in Bytes).

Annahmen

  • Alle Koordinaten sind ganze Zahlen im Bereich 0 ≤ x ≤ 100, 0 ≤ y ≤ 100.
  • Es wird mindestens 3und höchstens 50Eckpunkte geben.
  • Es wird keine wiederholten Eckpunkte geben. Auch werden keine Eckpunkte auf einer anderen Kante liegen. (Es kann jedoch kollineare Punkte in der Liste geben.)

Testfälle

{{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
Martin Ender
quelle
1
Insbesondere möchte ich die Trennzeichen so ersetzen, dass die Liste einen gültigen PostScript-Benutzerpfad enthält, sodass ich das Ganze mit einem upathOperator analysieren kann . (Es ist eigentlich eine extrem einfache 1: 1-Konvertierung zwischen den Trennzeichen. Wird }, {einfach entfernt linetound das Komma zwischen x und y wird entfernt und die öffnenden und schließenden Klammern werden durch eine statische Kopf- und Fußzeile ersetzt ...)
AJMansfield
1
@AJMansfield Normalerweise macht es mir nichts aus, praktische, native Listendarstellungen zu verwenden, aber mit upathund linetoklingt es, als würden Sie die Eingabe tatsächlich vorverarbeiten. Dh Sie nehmen keine Koordinatenliste, sondern ein tatsächliches Polygon.
Martin Ender
1
@MattNoonan Oh, das ist ein guter Punkt. Ja, das darfst du annehmen.
Martin Ender
2
@Ray Während die Richtung die Anzahl der Überfahrten beeinflussen kann, wird sie immer nur um 2 erhöht oder verringert, wobei die Parität erhalten bleibt. Ich werde versuchen, eine Referenz zu finden. Zunächst verwendet SVG dieselbe Definition.
Martin Ender
1
12,0 Mathematica hat eine neue integrierte Funktion dafür: CrossingPolygon.
Alephalpha

Antworten:

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&

Fügen Sie zuerst die Schnittpunkte zum Polygon hinzu und kehren Sie dann einige der Kanten um. Anschließend kann die Fläche genau wie bei einem einfachen Polygon berechnet werden.

Bildbeschreibung hier eingeben

Beispiel:

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
quelle
Leider bin ich mir nicht sicher, ob diese Logik in allen Situationen funktioniert. Kannst du es versuchen {1,2},{4,4},{4,2},{2,4},{2,1},{5,3}? Sie sollten mit 3.433333333333309 herauskommen. Ich habe nach einer ähnlichen Logik gesucht.
MickyT
@ MickyT Ja, es funktioniert. Es wurde zurückgegeben 103/30, und der numerische Wert ist 3.43333.
Alephalpha
Das tut mir leid. Gute Lösung
MickyT
44

Python 2, 323 319 Bytes

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"})

Nimmt eine Liste von Eckpunkten durch STDIN als komplexe Zahlen in der folgenden Form

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

und schreibt das Ergebnis nach STDOUT.

Gleicher Code nach dem Ersetzen der Zeichenfolge und etwas Abstand:

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]
)

Erläuterung

Führen Sie für jeden Schnittpunkt zweier Seiten des Eingabepolygons (einschließlich der Scheitelpunkte) eine vertikale Linie durch diesen Punkt.

Abbildung 1

(Tatsächlich übergibt das Programm aufgrund des Golfspiels ein paar Zeilen mehr; es spielt keine Rolle, solange mindestens diese Zeilen übergeben werden.) Der Körper des Polygons zwischen zwei aufeinanderfolgenden Zeilen besteht aus vertikalen Trapezoiden ( und Dreiecke und Liniensegmente, als Sonderfälle davon). Es muss der Fall sein, da, wenn eine dieser Formen einen zusätzlichen Scheitelpunkt zwischen den beiden Basen hätte, es eine weitere vertikale Linie durch diesen Punkt zwischen den beiden fraglichen Linien geben würde. Die Summe der Flächen aller solcher Trapezoide ist die Fläche des Polygons.

So finden wir diese Trapezoide: Für jedes Paar aufeinanderfolgender vertikaler Linien finden wir die Segmente jeder Seite des Polygons, die (ordnungsgemäß) zwischen diesen beiden Linien liegen (die für einige der Seiten möglicherweise nicht vorhanden sind). In der obigen Abbildung sind dies die sechs roten Segmente unter Berücksichtigung der beiden roten vertikalen Linien. Beachten Sie, dass sich diese Segmente nicht richtig kreuzen (dh sie treffen sich möglicherweise nur an ihren Endpunkten, fallen vollständig zusammen oder kreuzen sich überhaupt nicht, da, wenn sie sich richtig kreuzen, eine weitere vertikale Linie dazwischen liegen würde;) und so ist es sinnvoll darüber zu sprechen, sie von oben nach unten zu ordnen, was wir tun. Nach der Gerade-Ungerade-Regel befinden wir uns innerhalb des Polygons, sobald wir das erste Segment überquert haben. Sobald wir den zweiten überqueren, sind wir raus; der dritte ist wieder da; der vierte raus; und so weiter...

Insgesamt ist dies ein O ( n 3 log n ) -Algorithmus.

Ell
quelle
4
Das ist brilliant! Ich wusste, dass ich auf dich zählen kann. ;) (Vielleicht möchten Sie diese Frage bei Stack Overflow beantworten .)
Martin Ender
@ MartinBüttner Keep 'em coming :)
Ell
7
Tolle Arbeit und eine tolle Erklärung
MickyT
1
Dies ist eine beeindruckende Antwort. Haben Sie den Algorithmus selbst entwickelt oder gibt es bereits Arbeiten zu diesem Problem? Wenn es bereits eine Arbeit gibt, würde ich mich über einen Hinweis freuen, wo ich ihn finden kann. Ich hatte keine Ahnung, wie ich das angehen sollte.
Logic Knight
5
@CarpetPython Ich habe es selbst entwickelt, aber ich wäre sehr überrascht, wenn es noch nicht geschehen wäre.
Ell
9

Haskell, 549

Es sieht nicht so aus, als könnte ich dieses Spiel weit genug runter spielen, aber das Konzept war anders als die beiden anderen Antworten, also dachte ich, ich würde es trotzdem teilen. Es führt rationale O (N ^ 2) -Operationen durch, um die Fläche zu berechnen.

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

Beispiel:

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

Die Idee ist, das Polygon an jeder Kreuzung neu zu verdrahten, was zu einer Vereinigung von Polygonen ohne sich überschneidende Kanten führt. Wir können dann den (vorzeichenbehafteten) Bereich jedes der Polygone unter Verwendung der Gauss-Formel für Schnürsenkel berechnen ( http://en.wikipedia.org/wiki/Shoelace_formula ). Die Gerade-Ungerade-Regel verlangt, dass beim Konvertieren einer Kreuzung die Fläche des neuen Polygons im Verhältnis zum alten Polygon negativ gezählt wird.

Betrachten Sie beispielsweise das Polygon in der ursprünglichen Frage. Die Kreuzung oben links wird in zwei Pfade umgewandelt, die sich nur an einem Punkt treffen. Die beiden Pfade sind beide im Uhrzeigersinn ausgerichtet, sodass die Bereiche jedes Pfads positiv wären, wenn wir nicht deklarieren, dass der innere Pfad gegenüber dem äußeren Pfad mit -1 gewichtet ist. Dies entspricht der Pfadumkehr von alphaalpha.

Polygone abgeleitet vom ursprünglichen Beispiel

Als weiteres Beispiel betrachten wir das Polygon aus MickyTs Kommentar:

Polygone abgeleitet von MickyTs Kommentar

Hier sind einige der Polygone im Uhrzeigersinn und einige gegen den Uhrzeigersinn ausgerichtet. Die Vorzeichen-Kipp-Kreuzungsregel stellt sicher, dass die im Uhrzeigersinn ausgerichteten Bereiche einen zusätzlichen Faktor von -1 aufnehmen, wodurch sie einen positiven Beitrag zur Fläche leisten.

So funktioniert das Programm:

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
Matt Noonan
quelle