No algebric-ajuste de curvas a lo largo de pointcloud ponderada (si es posible usando python)

Question

Tengo una lista de ponderación de 2D puntos tomados a partir de análisis de simetría de un humano en la superficie de la espalda. Se supone que debo encontrar la "línea media", que representa la trayectoria más probable que describe las vértebras ubicación (en realidad, el proceso espinal de las vértebras).

EDIT: UN representante del conjunto de datos es la siguiente:

[[ -0.7898176   -3.35201728   4.36142086]
 [  2.99221402  -3.35201728   1.11907575]
 [  6.97475149  -3.35201728   2.4320322 ]
 [ -4.82443609  -2.35201728   0.6479064 ]
 [ -1.32418909  -2.35201728   1.88004944]
 [  0.07067882  -2.35201728   1.10982834]
 [  3.09169448  -2.35201728   1.8557436 ]
 [  7.10399403  -2.35201728   2.03906224]
 [ -3.07207606  -1.35201728   0.35500973]
 [  2.63202993  -1.35201728   5.32397834]
 [  5.19884868  -1.35201728   1.63816326]
 [  7.65721835  -1.35201728   1.13843392]
 [  2.48172754  -0.35201728   6.65584512]
 [  6.0905911   -0.35201728   1.15552652]
 [  8.62497546  -0.35201728   0.30407144]
 [ -4.7300089    0.64798272   0.31481496]
 [ -3.03274093   0.64798272   0.95337568]
 [  2.19653614   0.64798272  10.3675204 ]
 [  6.20384058   0.64798272   1.42106077]
 [ -4.08636605   1.64798272   0.28875288]
 [  2.03344989   1.64798272  13.04648211]
 [ -4.11717795   2.64798272   0.39713141]
 [  1.93304283   2.64798272  10.41313242]
 [ -4.37994815   3.64798272   0.84588643]
 [  1.66081408   3.64798272  14.96380955]
 [ -4.19024027   4.64798272   0.73216113]
 [  1.60252433   4.64798272  14.72419286]
 [  6.77837359   4.64798272   0.6186005 ]
 [ -4.14362668   5.64798272   0.93673165]
 [  1.55372968   5.64798272  12.9421123 ]
 [ -4.62223541   6.64798272   0.6510101 ]
 [  1.527865     6.64798272  10.80209351]
 [  6.86820685   6.64798272   0.82550801]
 [ -4.68259732   7.64798272   0.45321369]
 [  1.36167494   7.64798272   6.45338514]
 [ -5.19205787   8.64798272   0.23935013]
 [  1.21003466   8.64798272  10.13528877]
 [  7.6689546    8.64798272   0.32421776]
 [ -5.36436818   9.64798272   0.79809416]
 [  1.26248534   9.64798272   7.67036253]
 [  7.35472418   9.64798272   0.92555691]
 [ -5.61723652  10.64798272   0.4741007 ]
 [  1.23101086  10.64798272   7.97064105]
 [ -7.83024735  11.64798272   0.47557318]
 [  1.20348982  11.64798272   8.20694816]
 [  1.14422758  12.64798272   9.26244889]
 [  9.18164464  12.64798272   0.72428381]
 [  1.0827069   13.64798272  10.08599118]
 [  6.80116007  13.64798272   0.4571425 ]
 [  9.384236    13.64798272   0.42399893]
 [  1.04053491  14.64798272  10.48370805]
 [  9.16197972  14.64798272   0.39930227]
 [ -9.85958581  15.64798272   0.39524976]
 [  0.9942501   15.64798272   8.39992264]
 [  8.07642416  15.64798272   0.61480371]
 [  9.55088151  15.64798272   0.54076473]
 [ -7.13657331  16.64798272   0.32929172]
 [  0.92606211  16.64798272   7.83597033]
 [  8.74291069  16.64798272   0.74246827]
 [ -7.20022443  17.64798272   0.52555351]
 [  0.81344517  17.64798272   6.81654834]
 [  8.52844624  17.64798272   0.70543711]
 [ -6.97465178  18.64798272   1.04527813]
 [  0.61959631  18.64798272  10.33529022]
 [  5.733054    18.64798272   1.2309691 ]
 [  8.14818453  18.64798272   1.37532423]
 [ -6.82823664  19.64798272   2.0314052 ]
 [  0.56391636  19.64798272  13.61447357]
 [  5.79971126  19.64798272   0.30148347]
 [  8.01499476  19.64798272   1.72465327]
 [ -6.78504689  20.64798272   2.88657804]
 [ -4.79580634  20.64798272   0.36201975]
 [  0.548376    20.64798272   7.8414544 ]
 [  7.62258506  20.64798272   1.52817905]
 [-10.50328534  21.64798272   0.90358671]
 [ -6.59976138  21.64798272   2.62980169]
 [ -3.71180255  21.64798272   1.27094175]
 [  0.5060743   21.64798272  11.06117677]
 [  4.51983105  21.64798272   1.74626435]
 [  7.50948795  21.64798272   3.46497629]
 [ 11.10199877  21.64798272   1.78047269]
 [-10.15444935  22.64798272   1.47486166]
 [ -6.26274479  22.64798272   4.73707852]
 [ -3.45440904  22.64798272   1.72516012]
 [  0.52759064  22.64798272  12.58470433]
 [  4.22258017  22.64798272   2.63827535]
 [  7.03480033  22.64798272   3.506412  ]
 [ 10.63560314  22.64798272   3.56076386]
 [ -5.95693623  23.64798272   2.97403863]
 [ -3.66261423  23.64798272   2.31667236]
 [  0.52051366  23.64798272  12.5526344 ]
 [  4.21083787  23.64798272   1.95794387]
 [  6.82438636  23.64798272   4.77995659]
 [ 10.18138299  23.64798272   5.21836205]
 [ -9.94629932  24.64798272   0.4074823 ]
 [ -5.74101948  24.64798272   2.60992238]
 [  0.52987226  24.64798272  10.68846987]
 [  6.29981921  24.64798272   3.56204471]
 [  9.96431168  24.64798272   2.85079129]
 [ -9.64229717  25.64798272   0.4503241 ]
 [ -5.579063    25.64798272   0.64475469]
 [  0.52053534  25.64798272  10.05046667]
 [  5.79167815  25.64798272   0.92797027]
 [ 10.05116919  25.64798272   2.52194933]
 [ -8.55286247  26.64798272   0.94447148]
 [  0.45065604  26.64798272  10.97432823]
 [  5.50068393  26.64798272   2.39645232]
 [ 10.08992273  26.64798272   2.77716257]
 [-16.62381217  27.64798272   0.2021621 ]
 [ -9.62146213  27.64798272   0.62245778]
 [ -7.66905507  27.64798272   2.84466396]
 [  0.38656111  27.64798272  10.74369366]
 [  5.76925402  27.64798272   1.13362978]
 [  9.83525197  27.64798272   1.18241147]
 [-15.64874512  28.64798272   0.18279302]
 [ -7.52932494  28.64798272   2.94012191]
 [  0.32171219  28.64798272  10.73770466]
 [  9.4062684   28.64798272   1.41714298]
 [-12.71287717  29.64798272   0.70268073]
 [ -7.59473877  29.64798272   2.16183026]
 [  0.20748772  29.64798272  12.97312987]
 [  3.92952496  29.64798272   1.54987681]
 [  9.05148017  29.64798272   2.40563748]
 [ 14.96021523  29.64798272   0.55258241]
 [-12.14428813  30.64798272   0.36365363]
 [ -7.12360666  30.64798272   2.54312163]
 [  0.40594038  30.64798272  12.64839117]
 [  4.59465757  30.64798272   1.23496581]
 [  8.54333134  30.64798272   2.18912857]
 [-10.6296531   31.64798272   1.4839259 ]
 [ -7.09532763  31.64798272   2.0113838 ]
 [  0.37037733  31.64798272  12.2071139 ]
 [  3.01253349  31.64798272   3.01591777]
 [  4.64523695  31.64798272   3.50267541]
 [  8.39369696  31.64798272   2.53195817]
 [ -7.07947026  32.64798272   1.01324147]
 [  0.39269437  32.64798272   9.67368625]
 [  8.58669997  32.64798272   1.00475646]
 [ 12.02329114  32.64798272   0.50782399]
 [-10.13060786  33.64798272   0.31475653]
 [ -7.30360407  33.64798272   0.35065243]
 [  0.49556923  33.64798272   9.66608818]
 [ -5.37822311  34.64798272   0.38727401]
 [  0.4958055   34.64798272   7.5415026 ]
 [  6.07719006  34.64798272   0.63012453]
 [ -4.64579055  35.64798272   0.39990249]
 [  0.46323666  35.64798272   4.60449213]
 [  4.72819312  35.64798272   0.98050594]
 [ -4.62418372  36.64798272   0.64160709]
 [  0.48866236  36.64798272   4.29331656]
 [  5.06493722  36.64798272   0.59888608]
 [  0.49730481  37.64798272   1.32828464]
 [ -1.31849217  38.64798272   0.70780886]
 [  1.70966455  38.64798272   0.88052135]
 [  0.06305774  39.64798272   0.47366487]
 [  2.13639356  39.64798272   0.67971461]
 [ -0.84726354  40.64798272   0.63787522]
 [  0.55723562  40.64798272   0.62855097]
 [  2.22359779  40.64798272   0.33884894]
 [  0.77309816  41.64798272   0.4605534 ]
 [  0.56144565  42.64798272   0.43678788]]

Cuando se trazan, mi pointcloud normalmente se ve como este, con un color que representa la "idoneidad" de un punto (cuánto de su barrio es simétrica, de acuerdo a un determinado criterio):

EDIT: los puntos de color se obtiene de la siguiente manera:

1. Para cada nivel horizontal (puntos con la misma coordenada Y), casi un segmento horizontal se toma de la parte posterior de la superficie del modelo. Para cada punto en este sector, una tangente dibujada y el sector se divide en lado derecho y lado izquierdo. Entonces, la diferencia entre las distancias de cada bilateral par de puntos para que la línea se mide y se suman. Entonces, para cada punto tenemos un peso que indica la cantidad del lado derecho y el lado izquierdo de la corte en ese punto son iguales. Las imágenes siguientes ilustran:

Gráfico que muestra una porción de la superficie de la espalda a nivel lumbar, donde cada punto recibe un valor por su asimetría. Los mínimos locales seleccionadas (puntos rojos) con un peso (triángulos verdes) proporcional a cuánto es "mejor" que la de sus alrededores (qué tan profundo su asimetría valor es, a nivel local).

Me pregunto cual de ajuste/método de regresión debo usar, con el fin de obtener:

• Local de la diferenciabilidad (podría ser aproximada);
• Capacidad para interpolar para resolución arbitraria;
• No se rebasen o impredecible oscilaciones;
• Si es posible, la eliminación de valores atípicos.

Sospecho que algo se podría hacer con menos plazas, funciones de base radial, RANSAC, etc. pero también sospecho yo no tengo el conocimiento suficiente para seleccionar un método apropiado, así que por eso estoy pidiendo a los expertos.

EDIT: la realización de un promedio ponderado de interpolación spline, obtengo resultados similares a los de abajo, pero es demasiado sensible a los "valores atípicos", y sujeto a las oscilaciones obviamente lejos de lo que debería ser el "derecho" de la línea central:

Cualquier ayuda es bienvenida, gracias por la lectura"

Answer 1

(suponiendo que usted desea utilizar python)

La forma más fácil (dado lo que está disponible actualmente) sería el uso de un polinomio y métodos robustos de statsmodels.

algo como

endog = y  # observed points, one dimensional array
#polynomial array as explanatory variable:
#assuming x contains the vertical points
x = x / float(x.max() - x.min()) * 2 - 1 #optional rescaling 
exog = np.vander(x, 5) 

import statsmodels.api as sm
#robust estimation
res = sm.RLM(endog, exog).fit()
print res.summary()
poly = np.poly1d(res.params)
#now we have a polynomial, and we can use differentiation, ... with it

Me acaba de escribir esto, y no lo hizo. También no sé cual es el orden del polinomio sería útil en este caso.

Me gustaría utilizar, chebvander de numpy.polinomio, pero no recuerdo los detalles, especialmente el manejo del cambio en el dominio.

RLM.ajuste tiene varias opciones diferentes, si los "outliers" o la gran dispersión de los puntos no son manejados correctamente con el valor predeterminado. http://statsmodels.sourceforge.net/stable/rlm.html

Cuantil de regresión sería una buena alternativa a RLM. Alguien contribuido de una secuencia de comandos para statsmodels, pero no está incluido aún. cuantil de regresión podría ser utilizado para la estimación de una "línea media" que puede funcionar bien en este caso.

Los splines que se acoplan en lugar de polinomios podría ser mejor, pero no está disponible en python sin un poco de trabajo, y creo que no va a hacer mucha diferencia, dado el número relativamente pequeño de la vertical de puntos.

actualización:

utilizando los datos de las columnas 0 y 1, he probado RLM. No he podido conseguir las opciones obvias para RLM para el trabajo, pero todavía hay algunos más infrecuente de las opciones que no he probado.

Lo he intentado en su lugar fue para poner directamente todo el peso en el centro de los puntos. Esto es esencialmente lo menos recortan cuadrados, con una inicial robusta (M-) estimación de los puntos centrales. La primera etapa se ajusta un polinomio de orden bajo. A continuación, seleccione los puntos con un pequeño residual. He seleccionado el umbral manualmente para obtener los puntos de derecho. Entonces me ajuste a un polinomio de orden superior con mínimos cuadrados ponderados y poner el peso solo en el centro de los puntos.

Esto parece funcionar para la estimación de la línea por el centro de los puntos, pero aún no sé cómo se iba a realizar para obtener más general de problemas.

Excepto para la importación de numpy y matplotlib.pyplot y para la carga de los datos, mi script para producir la gráfica es:

#use data in column 1
#rescaling didn't matter much
d = 0.1
x_rs = a[:,1]  #a is original data
x_rs = (x_rs - x_rs.min())
x_rs = x_rs / x_rs.max() #* (1-d) *2 - 1 + d
print x_rs.min(), x_rs.max()

degree = 3  #preliminary regression
exog = x_rs[:,None] ** np.arange(degree+1)
degree2 = 6  #trimmed regression
exog2 = x_rs[:,None] ** np.arange(degree2+1)

endog = a[:,0]

import statsmodels.api as sm
#fit low order polynomial with robust estimator
res = sm.RLM(endog, exog).fit()

#fit on center points, threshold 1.31 chosen by inspection of residuals
resw = sm.WLS(endog, exog2, weights=(np.abs(res.resid) < 1.31)).fit()

plt.plot(a[:,1], a[:,0], 'o')
plt.plot(a[:,1], res.fittedvalues,  '-', lw=2, label='RLM')
idx = np.nonzero(np.abs(res.resid) < 1.31)[0]
plt.plot(a[idx,1], a[idx,0], 'ro', alpha=0.5)
plt.plot(a[:,1], resw.fittedvalues,  '-', lw=2, label='WLS center')
plt.legend()

plt.show()

actualización 2

Uno de los problemas con las opciones para reforzar la estimación de cuánto control que tenemos sobre la varianza de la estimación. RLM actualizaciones de la varianza y en este caso es downweights la periferia de las observaciones no es suficiente.

Si tenemos una estimación previa de la varianza, se puede forzar a RLM para usarlo y no para su actualización.

La varianza de la estimación por el tapizado de los mínimos cuadrados, WLS anteriormente, es sólo acerca de 0.014. Si utilizamos esta información,

scale = resw.scale
resw2 = sm.RLM(endog, exog2).fit(init=dict(scale=scale), update_scale=False)

a continuación, el conjunto de la curva es muy similar a la de WLS arriba la curva, aunque las estimaciones de los parámetros para el polinomio son diferentes.

(Advertencia: Esto no funciona, sin embargo en statsmodels, init opción sólo está disponible en una oficina donde yo trabajaba en la solidez de la estimación y mejoras a RLM.)