Área de un polígono autoinsectable

32

Considere un polígono potencialmente auto-intersectado, definido por una lista de vértices en el espacio 2D. P.ej

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

Hay varias formas de definir el área de dicho polígono, pero la más interesante es la regla par-impar. Tomando cualquier punto en el plano, dibuja una línea desde el punto hasta el infinito (en cualquier dirección). Si esa línea cruza el polígono un número impar de veces, el punto es parte del área del polígono, si cruza el polígono un número par de veces, el punto no es parte del polígono. Para el polígono de ejemplo anterior, aquí están tanto su contorno como su área par-impar:

contornoZona

El polígono no será en general ortogonal. Solo he elegido un ejemplo tan simple para que sea más fácil contar el área.

El área de este ejemplo es 17(no 24o 33como podrían dar otras definiciones o áreas).

Tenga en cuenta que, según esta definición, el área del polígono es independiente de su orden de bobinado.

El reto

Dada una lista de vértices con coordenadas enteras que definen un polígono, determine su área bajo la regla par-impar.

Puede escribir una función o programa, tomando la entrada a través de STDIN o la alternativa más cercana, el argumento de la línea de comandos o el argumento de la función y devolver el resultado o imprimirlo en STDOUT o la alternativa más cercana.

Puede tomar la entrada en cualquier lista conveniente o formato de cadena, siempre que no esté preprocesada.

El resultado debe ser un número de coma flotante, con una precisión de 6 dígitos significativos (decimales), o un resultado racional cuya representación en coma flotante tenga una precisión de 6 dígitos significativos. (Si produce resultados racionales, es probable que sean exactos, pero no puedo requerir esto, ya que no tengo resultados exactos para referencia).

Debe poder resolver cada uno de los casos de prueba a continuación en 10 segundos en una máquina de escritorio razonable. (Hay un margen de maniobra en esta regla, así que use su mejor criterio. Si tarda 20 segundos en mi computadora portátil, le daré el beneficio de la duda, si tarda un minuto, no lo haré). Creo que este límite debería ser muy generoso, pero se supone que debe descartar enfoques en los que simplemente se discretiza el polígono en una cuadrícula lo suficientemente fina y se cuenta, o se utilizan enfoques probabilísticos como Monte Carlo. Sea un buen deportista y no intente optimizar estos enfoques de modo que pueda cumplir con el límite de tiempo de todos modos. ;)

No debe utilizar ninguna función existente relacionada directamente con polígonos.

Este es el código de golf, por lo que gana el envío más corto (en bytes).

Supuestos

  • Todas las coordenadas son enteros en el rango 0 ≤ x ≤ 100, 0 ≤ y ≤ 100.
  • Habrá al menos 3y a lo sumo 50vértices.
  • No habrá vértices repetidos. Tampoco los vértices se encuentran en otro borde. (Sin embargo, puede haber puntos colineales en la lista).

Casos de prueba

{{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
fuente
1
Específicamente, me gustaría reemplazar los delimitadores de una manera que haga que la lista sea una ruta de usuario PostScript válida, de modo que pueda analizar todo con un upathoperador. (En realidad, es una conversión extremadamente simple de 1: 1 entre los separadores. }, {Simplemente se convierte lineto, y la coma entre x e y se elimina, y las llaves de apertura y cierre se reemplazan con un encabezado y pie de página estáticos ...)
AJMansfield
1
@AJMansfield Por lo general, no me importa usar representaciones de listas nativas convenientes, pero usar upathy linetosuena como si realmente estuvieras preprocesando la entrada. Es decir, no estás tomando una lista de coordenadas sino un polígono real.
Martin Ender
1
@ MattNoonan Oh, ese es un buen punto. Sí, puedes asumir eso.
Martin Ender
2
@Ray Si bien la dirección puede afectar el número de cruces, solo aumentará o disminuirá en 2, conservando la paridad. Trataré de encontrar una referencia. Para empezar, SVG usa la misma definición.
Martin Ender
1
Mathematica 12.0 tiene una nueva función incorporada para esto: CrossingPolygon.
alephalpha

Respuestas:

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&

Primero agregue los puntos de intersección al polígono, luego invierta algunos de los bordes, luego su área puede calcularse como un polígono simple.

ingrese la descripción de la imagen aquí

Ejemplo:

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
fuente
Lamentablemente, no estoy seguro de que esta lógica funcione para todas las situaciones. Se puede tratar {1,2},{4,4},{4,2},{2,4},{2,1},{5,3}? Deberías salir con 3.433333333333309. Lo miré usando una lógica similar.
MickyT
@MickyT Sí, funciona. Regresó 103/30, y el valor numérico es 3.43333.
alephalpha
Lo siento por eso. Buena solución
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"})

Toma una lista de vértices a través de STDIN como números complejos, en el siguiente formulario

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

y escribe el resultado en STDOUT.

Mismo código después del reemplazo de la cadena y algo de espacio:

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

Explicación

Para cada punto de intersección de dos lados del polígono de entrada (incluidos los vértices), pase una línea vertical a través de ese punto.

Figura 1

(De hecho, debido al golf, el programa pasa unas pocas líneas más; en realidad no importa, siempre y cuando pasemos al menos estas líneas). El cuerpo del polígono entre dos líneas consecutivas está compuesto de trapecios verticales ( y triángulos, y segmentos de línea, como casos especiales de esos). Tiene que ser el caso, ya que si alguna de estas formas tuviera un vértice adicional entre las dos bases, habría otra línea vertical a través de ese punto, entre las dos líneas en cuestión. La suma de las áreas de todos esos trapecios es el área del polígono.

Así es como encontramos estos trapecios: para cada par de líneas verticales consecutivas, encontramos los segmentos de cada lado del polígono que (correctamente) se encuentran entre estas dos líneas (que podrían no existir para algunos de los lados). En la ilustración anterior, estos son los seis segmentos rojos, cuando se consideran las dos líneas verticales rojas. Tenga en cuenta que estos segmentos no se cruzan correctamente entre sí (es decir, solo pueden encontrarse en sus puntos finales, coincidir completamente o no cruzarse en absoluto, ya que, una vez más, si se cruzaran correctamente, habría otra línea vertical en el medio;) y tiene sentido hablar de ordenarlos de arriba a abajo, lo cual hacemos. De acuerdo con la regla par-impar, una vez que cruzamos el primer segmento, estamos dentro del polígono; una vez que cruzamos el segundo, estamos fuera; el tercero, de nuevo; el cuarto, fuera; y así...

En general, este es un algoritmo O ( n 3 log n ).

Ana
fuente
44
¡Esto es brillante! Sabía que podía contar contigo para esto. ;) (Es posible que desee responder esta pregunta en Stack Overflow.)
Martin Ender
@ MartinBüttner Sigue viniendo :)
Ell
77
Gran trabajo y una gran explicación
MickyT
1
Esta es una respuesta impresionante. ¿Desarrolló el algoritmo usted mismo o existe algún trabajo sobre este problema? Si hay trabajo existente, agradecería un puntero a donde puedo encontrarlo. No tenía idea de cómo abordar esto.
Logic Knight
55
@CarpetPython Lo desarrollé yo mismo, pero me sorprendería mucho si no se ha hecho antes.
Ell
9

Haskell, 549

No parece que pueda jugar golf lo suficiente, pero el concepto era diferente a las otras dos respuestas, así que pensé que lo compartiría de todos modos. Realiza operaciones racionales O (N ^ 2) para calcular el área.

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

Ejemplo:

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

La idea es volver a cablear el polígono en cada cruce, lo que resulta en una unión de polígonos sin bordes de intersección. Luego podemos calcular el área (firmada) de cada uno de los polígonos usando la fórmula de cordones de Gauss ( http://en.wikipedia.org/wiki/Shoelace_formula ). La regla par-impar exige que cuando se convierte un cruce, el área del nuevo polígono se cuenta negativamente en relación con el antiguo polígono.

Por ejemplo, considere el polígono en la pregunta original. El cruce en la esquina superior izquierda se convierte en dos caminos que solo se encuentran en un punto; las dos rutas están orientadas en sentido horario, por lo que las áreas de cada una serían positivas a menos que declaremos que la ruta interna está ponderada por -1 en relación con la ruta externa. Esto es equivalente a la inversión de ruta de alfaalfa.

Polígonos derivados del ejemplo original.

Como otro ejemplo, considere el polígono del comentario de MickyT:

Polígonos derivados del comentario de MickyT

Aquí, algunos de los polígonos están orientados en sentido horario y otros en sentido antihorario. La regla de signo de volteo al cruzar asegura que las regiones orientadas en sentido horario recojan un factor adicional de -1, lo que hace que contribuyan con una cantidad positiva al área.

Así es como funciona el programa:

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
fuente