Encuentra el centroide de un polígono

16

De Wikipedia :

El centroide de un polígono cerrado sin intersección propia definido por n vértices ( x 0 , y 0 ), ( x 1 , y 1 ), ..., ( x n - 1 , y n − 1 ) es el punto ( C x , C y ), donde

Fórmula para el centroide

y donde A es el área firmada del polígono,

Fórmula para el área del polígono

En estas fórmulas, se supone que los vértices están numerados en orden de ocurrencia a lo largo del perímetro del polígono. Además, se supone que el vértice ( x n , y n ) es el mismo que ( x 0 , y 0 ), lo que significa que i + 1 en el último caso debe girar alrededor de i = 0 . Tenga en cuenta que si los puntos están numerados en el sentido de las agujas del reloj, el área A , calculada como se indica arriba, tendrá un signo negativo; pero las coordenadas del centroide serán correctas incluso en este caso.


  • Dada una lista de vértices en orden (ya sea en sentido horario o en sentido antihorario), encuentre el centroide del polígono cerrado que no se cruza automáticamente representado por los vértices.
    • Si ayuda, puede suponer que la entrada es solo CW o solo CCW. Dígalo en su respuesta si necesita esto.
  • No se requiere que las coordenadas sean números enteros, y pueden contener números negativos.
  • La entrada siempre será válida y contendrá al menos tres vértices.
  • Solo es necesario manejar las entradas que se ajusten al tipo de datos de punto flotante nativo de su idioma.
  • Puede suponer que los números de entrada siempre contendrán un punto decimal.
  • Puede suponer que los enteros de entrada terminan en .o .0.
  • Puede usar números complejos para la entrada.
  • La salida debe ser precisa a la milésima más cercana.

Ejemplos

[(0.,0.), (1.,0.), (1.,1.), (0.,1.)]        -> (0.5, 0.5)
[(-15.21,0.8), (10.1,-0.3), (-0.07,23.55)]  -> -1.727 8.017
[(-39.00,-55.94), (-56.08,-4.73), (-72.64,12.12), (-31.04,53.58), (-30.36,28.29), (17.96,59.17), (0.00,0.00), (10.00,0.00), (20.00,0.00), (148.63,114.32), (8.06,-41.04), (-41.25,34.43)]   -> 5.80104769975, 15.0673812762

Para ver cada polígono en un plano de coordenadas, pegue las coordenadas sin los corchetes en el menú "Editar" de esta página .

Confirmé mis resultados usando esta calculadora de puntos de centroide poligonal , que es horrible. No pude encontrar uno en el que pueda ingresar todos los vértices a la vez, o que no intentara borrar su -signo cuando lo escribe primero. Publicaré mi solución Python para su uso después de que las personas hayan tenido la oportunidad de responder.

mbomb007
fuente
La técnica mucho más simple de promediar todas las x e y funciona para los dos primeros conjuntos, pero no para el tercero. Me pregunto qué hace la diferencia ...
ETHproductions
1
@ETHproductions El tercer polígono no es convexo.
JungHwan Min
1
@ETHproductions Si aproxima un círculo con un polígono, puede mover el punto promedio arbitrario cerca de un punto en el círculo utilizando más puntos cerca de ese punto, mientras casi no afecta el centroide y mantiene el polígono convexo.
Christian Sievers
2
@ETHproductions En realidad, la convexidad no es la razón. Un promedio de todos los xs y ys pone todo el peso en los vértices en lugar de distribuida sobre el cuerpo. El primero funciona porque es regular, por lo que ambos métodos terminan en el centro de simetría. El segundo funciona porque para los triángulos ambos métodos conducen al mismo punto.
Ton Hospel
1
¿Podemos usar números complejos para E / S?
xnor

Respuestas:

16

Jalea , 25 24 22 21 18 bytes

S×3÷@×"
ṙ-żµÆḊçS€S

Aplica la fórmula que se muestra en el problema.

Guardado 3 bytes con ayuda de @ Jonathan Allan.

Pruébalo en línea! o Verificar todos los casos de prueba.

Explicación

S×3÷@×"  Helper link. Input: determinants on LHS, sum of pairs on RHS
S        Sum the determinants
 ×3      Multiply by 3
     ×"  Vectorized multiply between determinants and sums
   ÷@    Divide that by the determinant sum multipled by 3 and return

ṙ-żµÆḊçS€S  Main link. Input: 2d list of points
ṙ-          Rotate the list of points by 1 to the right
  ż         Interleave those with the original points
            This creates all overlapping slices of length 2
   µ        Start new monadic chain
    ÆḊ      Get the determinant of each slice
       S€   Get the sum of each slice (sum of pairs of points)
      ç     Call the helper link
         S  Sum and return
millas
fuente
Puede reemplazar ṁL‘$ṡ2con ṙ1ż@ożṙ1$
Jonathan Allan
@JonathanAllan Gracias, también puedo girar ṙ-żpara evitar el intercambio y guardar otro byte
millas
¡Oh si por supuesto!
Jonathan Allan
17

Mathematica, 23 bytes

RegionCentroid@*Polygon

Toma eso , Jelly!

Editar: Uno no simplemente vence a Jelly ...

Explicación

Polygon

Genere un polígono con vértices en los puntos especificados.

RegionCentroid

Encuentra el centroide del polígono.

JungHwan Min
fuente
2
Bueno, me ganaste, pero probablemente haya un camino más corto que el que tengo, todavía no tengo una comprensión completa de Jelly
millas del
3
@miles aw ... :(
JungHwan Min
4

J, 29 bytes

2+/@(+/\(*%3*1#.])-/ .*\)],{.

Aplica la fórmula que se muestra en el problema.

Uso

   f =: 2+/@(+/\(*%3*1#.])-/ .*\)],{.
   f 0 0 , 1 0 , 1 1 ,: 0 1
0.5 0.5
   f _15.21 0.8 , 10.1 _0.3 ,: _0.07 23.55
_1.72667 8.01667
   f _39 _55.94 , _56.08 _4.73 , _72.64 12.12 , _31.04 53.58 , _30.36 28.29 , 17.96 59.17 , 0 0 , 10 0 , 20 0 , 148.63 114.32 , 8.06 _41.04 ,: _41.25 34.43
5.80105 15.0674

Explicación

2+/@(+/\(*%3*1#.])-/ .*\)],{.  Input: 2d array of points P [[x1 y1] [x2 y2] ...]
                           {.  Head of P
                         ]     Get P
                          ,    Join, makes the end cycle back to the front
2                              The constant 2
2                      \       For each pair of points
                  -/ .*        Take the determinant
2    +/\                       Sum each pair of points
         *                     Multiply the sum of each pair by its determinant
          %                    Divide each by
             1#.]              The sum of the determinants
           3*                  Multiplied by 3
 +/@                           Sum and return
millas
fuente
4

Máxima, 124 118 116 112 106 byte

f(l):=(l:endcons(l[1],l),l:sum([3,l[i-1]+l[i]]*determinant(matrix(l[i-1],l[i])),i,2,length(l)),l[2]/l[1]);

No tengo experiencia con Maxima, por lo que cualquier sugerencia es bienvenida.

Uso:

(%i6) f([[-15.21,0.8], [10.1,-0.3], [-0.07,23.55]]);
(%o6)              [- 1.726666666666668, 8.016666666666668]
Christian Sievers
fuente
3

Raqueta 420 bytes

(let*((lr list-ref)(getx(lambda(i)(lr(lr l i)0)))(gety(lambda(i)(lr(lr l i)1)))(n(length l))(j(λ(i)(if(= i(sub1 n))0(add1 i))))
(A(/(for/sum((i n))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i))))2))
(cx(/(for/sum((i n))(*(+(getx i)(getx(j i)))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i)))))(* 6 A)))
(cy(/(for/sum((i n))(*(+(gety i)(gety(j i)))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i)))))(* 6 A))))
(list cx cy))

Sin golf:

(define(f l)
  (let* ((lr list-ref)
         (getx (lambda(i)(lr (lr l i)0)))
         (gety (lambda(i)(lr (lr l i)1)))
         (n (length l))
         (j (lambda(i) (if (= i (sub1 n)) 0 (add1 i))))
         (A (/(for/sum ((i n))
                (-(* (getx i) (gety (j i)))
                  (* (getx (j i)) (gety i))))
              2))
         (cx (/(for/sum ((i n))
                 (*(+(getx i)(getx (j i)))
                   (-(*(getx i)(gety (j i)))
                     (*(getx (j i))(gety i)))))
               (* 6 A)))
         (cy (/(for/sum ((i n))
                 (*(+(gety i)(gety (j i)))
                   (-(*(getx i)(gety (j i)))
                     (*(getx (j i))(gety i)))))
               (* 6 A))))
    (list cx cy)))

Pruebas:

(f '[(-15.21 0.8)  (10.1 -0.3)  (-0.07 23.55)] ) 
(f '[(-39.00 -55.94)  (-56.08 -4.73)  (-72.64 12.12)  (-31.04 53.58) 
     (-30.36 28.29)  (17.96 59.17)  (0.00 0.00)  (10.00 0.00)  
     (20.00 0.00) (148.63 114.32)  (8.06 -41.04)  (-41.25 34.43)])

Salida:

'(-1.7266666666666677 8.01666666666667)
'(5.8010476997538465 15.067381276150996)
rnso
fuente
3

R, 129127 bytes

function(l){s=sapply;x=s(l,`[`,1);y=s(l,`[`,2);X=c(x[-1],x[1]);Y=c(y[-1],y[1]);p=x*Y-X*y;c(sum((x+X)*p),sum((y+Y)*p))/sum(p)/3}

Función sin nombre que toma una lista R de tuplas como entrada. El equivalente nombrado se puede llamar usando, por ejemplo:

f(list(c(-15.21,0.8),c(10.1,-0.3),c(-0.07,23.55)))

Desengañado y explicado

f=function(l){s=sapply;                           # Alias for sapply
              x=s(l,`[`,1);                       # Split list of tuples into vector of first elements
              y=s(l,`[`,2);                       # =||= but for second element 
              X=c(x[-1],x[1]);                    # Generate a vector for x(i+1)
              Y=c(y[-1],y[1]);                    # Generate a vector for y(i+1)
              p=x*Y-X*y;                          # Calculate the outer product used in both A, Cx and Cy
              c(sum((x+X)*p),sum((y+Y)*p))/sum(p)/3    # See post for explanation
}

El paso final ( c(sum((x+X)*p),sum((y+Y)*p))/sum(p)*2/6) es una forma vectorizada de calcular ambos Cxy Cy. La suma en las fórmulas para Cxy Cyse almacena en un vector y, en consecuencia, se divide por la "suma enA " *2/6. P.ej:

(SUMinCx, SUMinCy) / SUMinA / 3

, y luego implícitamente impreso.

Pruébalo en R-Fiddle

Billywob
fuente
*2/6probablemente podría ser /3?
mbomb007
@ mbomb007 Es tan minuciosamente obvio, supongo que me quedé atrapado en el golf la otra parte. / encogimiento de hombros
Billywob
Elegante, me gusta tu uso sapplypara lidiar con esas listas. Podría haber margen para jugar golf aquí, no estoy seguro de cuán flexible es la entrada permitida. Si se le permite ingresar solo una secuencia de coordenadas, como c(-15.21,0.8,10.1,-0.3,-0.07,23.55), entonces puede guardar 17 bytes reemplazando las primeras líneas de su función con y=l[s<-seq(2,sum(1|l),2)];x=l[-s];. Es decir, establecer ycada elemento indexado par de l, y xser cada elemento indexado impar.
rturnbull
Sin embargo, sería aún mejor si pudiéramos ingresar una matriz (o matriz), como matrix(c(-15.21,0.8,10.1,-0.3,-0.07,23.55),2), como puede ser el comienzo de su función x=l[1,];y=l[2,];, que ahorra 35 bytes. (La matriz de entrada podría transponerse, en cuyo caso x=l[,1];y=l[,2];). Por supuesto, la solución más fácil de todas es si los puntos xy yson simplemente ingresados ​​como vectores separados function(x,y), pero no creo que eso esté permitido ...
rturnbull
@rturnbull Le pregunté a OP en los comentarios y él quería específicamente una lista de tuplas (muy inconveniente en R, por supuesto), así que no creo que el enfoque matricial esté permitido. E incluso si lo fuera, la entrada tendría que ser la parte del vector (es decir c(...)) y la conversión de la matriz tendría que hacerse dentro de la función.
Billywob
2

Python, 156127 bytes

def f(p):n=len(p);p=p+p[:1];i=s=0;exec'd=(p[i].conjugate()*p[i+1]).imag;s+=d;p[i]=(p[i]+p[i+1])*d;i+=1;'*n;print sum(p[:n])/s/3

Sin golf:

def f(points):
  n = len(points)
  points = points + [points[0]]
  determinantSum = 0
  for i in range(n):
    determinant = (points[i].conjugate() * points[i+1]).imag
    determinantSum += determinant
    points[i] = (points[i] + points[i+1]) * determinant
  print sum(points[:n]) / determinantSum / 3

Ideónalo.

Esto toma cada par de puntos [x, y]como un número complejox + y*j y genera el centroide resultante como un número complejo en el mismo formato.

Para el par de puntos [a, b]y [c, d], el valor a*d - b*cque se necesita para cada par de puntos puede calcularse a partir del determinante de la matriz

| a b |
| c d |

Usando aritmética compleja, los valores complejos a + b*jy c + d*jse pueden usar como

conjugate(a + b*j) * (c + d*j)
(a - b*j) * (c + d*j)
(a*c + b*d) + (a*d - b*c)*j

Observe que la parte imaginaria es equivalente al determinante. Además, el uso de valores complejos permite que los puntos se sumen fácilmente por componentes en las otras operaciones.

millas
fuente
2

R + sp (46 bytes)

Asume que el sppaquete está instalado ( https://cran.r-project.org/web/packages/sp/ )

Toma una lista de vértices (por ejemplo list(c(0.,0.), c(1.,0.), c(1.,1.), c(0.,1.)) )

Aprovecha el hecho de que el "labpt" de un polígono es el centroide.

function(l)sp::Polygon(do.call(rbind,l))@labpt
mnel
fuente
2

JavaScript (ES6), 102

Implementación directa de la fórmula

l=>[...l,l[0]].map(([x,y],i)=>(i?(a+=w=t*y-x*u,X+=(t+x)*w,Y+=(u+y)*w):X=Y=a=0,t=x,u=y))&&[X/3/a,Y/3/a]

Prueba

f=
l=>[...l,l[0]].map(([x,y],i)=>(i?(a+=w=t*y-x*u,X+=(t+x)*w,Y+=(u+y)*w):X=Y=a=0,t=x,u=y))&&[X/3/a,Y/3/a]

function go()
{
  var c=[],cx,cy;
  // build coordinates array
  I.value.match(/-?[\d.]+/g).map((v,i)=>i&1?t[1]=+v:c.push(t=[+v]));
  console.log(c+''),
  [cx,cy]=f(c);
  O.textContent='CX:'+cx+' CY:'+cy;
  // try to display the polygon
  var mx=Math.max(...c.map(v=>v[0])),
    nx=Math.min(...c.map(v=>v[0])),
    my=Math.max(...c.map(v=>v[1])),
    ny=Math.min(...c.map(v=>v[1])),  
    dx=mx-nx, dy=my-ny,
    ctx=C.getContext("2d"),
    cw=C.width, ch=C.height,
    fx=(mx-nx)/cw, fy=(my-ny)/ch, fs=Math.max(fx,fy)
  C.width=cw
  ctx.setTransform(1,0,0,1,0,0);
  ctx.beginPath();
  c.forEach(([x,y],i)=>ctx.lineTo((x-nx)/fs,(y-ny)/fs));
  ctx.closePath();
  ctx.stroke();
  ctx.fillStyle='#ff0000';
  ctx.fillRect((cx-nx)/fs-2,(cy-ny)/fs-2,5,5);
}
go()
#I { width:90% }
#C { width:90%; height:200px;}
<input id=I value='[[-15.21,0.8], [10.1,-0.3], [-0.07,23.55]]'>
<button onclick='go()'>GO</button>
<pre id=O></pre>
<canvas id=C></canvas>

edc65
fuente
1

Python 2, 153 bytes

No usa números complejos.

P=input()
A=x=y=0;n=len(P)
for i in range(n):m=-~i%n;a=P[i][0];b=P[i][1];c=P[m][0];d=P[m][1];t=a*d-b*c;A+=t;x+=t*(a+c);y+=t*(b+d)
k=1/(3*A);print x*k,y*k

Pruébalo en línea

Sin golf:

def centroid(P):
    A=x=y=0
    n=len(P)
    for i in range(n):
        m=-~i%n
        x0=P[i][0];y0=P[i][1]
        x1=P[m][0];y1=P[m][1]
        t = x0*y1 - y0*x1
        A += t/2.
        x += t * (x0 + x1)
        y += t * (y0 + y1)
    k = 1/(6*A)
    x *= k
    y *= k
    return x,y
mbomb007
fuente
1

En realidad, 45 40 39 bytes

Esto utiliza un algoritmo similar a la respuesta Jelly de miles . Hay una forma más corta de calcular los determinantes usando un producto de puntos, pero actualmente hay un error con el producto de puntos de Actually donde no funcionará con listas de flotadores. Sugerencias de golf bienvenidas. Pruébalo en línea!

;\Z♂#;`i¥`M@`i│N@F*)F@N*-`M;Σ3*)♀*┬♂Σ♀/

Ungolfing

         Implicit input pts.
;\       Duplicate pts, rotate right.
Z        Zip rot_pts and pts together.
♂#       Convert the iterables inside the zip to lists
         (currently necessary due to a bug with duplicate)
;        Duplicate the zip.
`...`M   Get the sum each pair of points in the zip.
  i        Flatten the pair to the stack.
  ¥        Pairwise add the two coordinate vectors.
@        Swap with the other zip.
`...`M   Get the determinants of the zip.
  i│       Flatten to stack and duplicate entire stack.
           Stack: [a,b], [c,d], [a,b], [c,d]
  N@F*)    Push b*c and move it to BOS.
  F@N*     Push a*d.
  -        Get a*d-b*c.
;Σ3*)    Push 3 * sum(determinants) and move it to BOS.
♀*       Vector multiply the determinants and the sums.
┬        Transpose the coordinate pairs in the vector.
♂Σ       Sum the x's, then the y's.
♀/       Divide the x and y of this last coordinate pair by 3*sum(determinants).
         Implicit return.

Una versión más corta y no competitiva.

Esta es otra versión de 24 bytes que usa números complejos. No es competitivo porque se basa en correcciones de errores posteriores a este desafío. Pruébalo en línea!

;\│¥)Z`iá*╫@X`M;Σ3*)♀*Σ/

Ungolfing

         Implicit input a list of complex numbers, pts.
;\       Duplicate pts, rotate right.
│        Duplicate stack. Stack: rot_pts, pts, rot_pts, pts.
¥)       Pairwise sum the two lists of points together and rotate to BOS.
Z        Zip rot_pts and pts together.
`...`M   Map the following function over the zipped points to get our determinants.
  i        Flatten the list of [a+b*i, c+d*i].
  á        Push the complex conjugate of a+bi, i.e. a-b*i.
  *        Multiply a-b*i by c+d*i, getting (a*c+b*d)+(a*d-b*c)*i.
           Our determinant is the imaginary part of this result.
  ╫@X      Push Re(z), Im(z) to the stack, and immediately discard Re(z).
           This map returns a list of these determinants.
;        Duplicate list_determinants.
Σ3*)     Push 3 * sum(list_determinants) and rotate that to BOS.
♀*Σ      Pairwise multiply the sums of pairs of points and the determinants and sum.
/        Divide that sum by 3*sum(list_determinants).
         Implicit return.
Sherlock9
fuente
1

C ++ 14, 241 bytes

struct P{float x;float y;};
#define S(N,T)auto N(P){return 0;}auto N(P a,P b,auto...V){return(T)*(a.x*b.y-b.x*a.y)+N(b,V...);}
S(A,1)S(X,a.x+b.x)S(Y,a.y+b.y)auto f(auto q,auto...p){auto a=A(q,p...,q)*3;return P{X(q,p...,q)/a,Y(q,p...,q)/a};}

La salida es la estructura auxiliar P,

Sin golf:

 //helper struct
struct P{float x;float y;};

//Area, Cx and Cy are quite similar
#define S(N,T)\  //N is the function name, T is the term in the sum
auto N(P){return 0;} \   //end of recursion for only 1 element
auto N(P a,P b,auto...V){ \ //extract the first two elements
  return (T)*(a.x*b.y-b.x*a.y) //compute with a and b
         + N(b,V...); \        //recursion without first element
}

//instantiate the 3 formulas
S(A,1)
S(X,a.x+b.x)
S(Y,a.y+b.y)


auto f(auto q,auto...p){
  auto a=A(q,p...,q)*3; //q,p...,q appends the first element to the end
  return P{X(q,p...,q)/a,Y(q,p...,q)/a};
}

Uso:

f(P{0.,0.}, P{1.,0.}, P{1.,1.}, P{0.,1.})
f(P{-15.21,0.8}, P{10.1,-0.3}, P{-0.07,23.55})
Karl Napf
fuente
1

Clojure, 177 156 143 bytes

Actualización: en lugar de una devolución de llamada, estoy usando [a b c d 1]una función y el argumento es solo una lista de índices para este vector. 1se usa como valor centinela al calcular A.

Actualización 2: No se calcula previamente Aen let, utilizando (rest(cycle %))para obtener los vectores de entrada compensados ​​por uno.

#(let[F(fn[I](apply +(map(fn[[a b][c d]](*(apply +(map[a b c d 1]I))(-(* a d)(* c b))))%(rest(cycle %)))))](for[i[[0 2][1 3]]](/(F i)(F[4])3)))

Versión original:

#(let[F(fn[L](apply +(map(fn[[a b][c d]](*(L[a b c d])(-(* a d)(* c b))))%(conj(subvec % 1)(% 0)))))A(*(F(fn[& l]1))3)](map F[(fn[v](/(+(v 0)(v 2))A))(fn[v](/(+(v 1)(v 3))A))]))

En la etapa menos golfizada:

(def f (fn[v](let[F (fn[l](apply +(map
                                    (fn[[a b][c d]](*(l a b c d)(-(* a d)(* c b))))
                                    v
                                    (conj(subvec v 1)(v 0)))))
                  A (* (F(fn[& l] 1)) 3)]
                [(F (fn[a b c d](/(+ a c)A)))
                 (F (fn[a b c d](/(+ b d)A)))])))

Crea una función auxiliar Fque implementa la suma con cualquier devolución de llamada l. Porque Ala devolución de llamada regresa constantemente, 1mientras que las coordenadas X e Y tienen sus propias funciones. (conj(subvec v 1)(v 0))suelta el primer elemento y se agrega al final, de esta manera es fácil hacer un seguimiento de x_iy x_(i+1). Tal vez todavía quede algo de repetición por eliminar, especialmente al final (map F[....

NikoNyrh
fuente