沿加权点云的非代数曲线拟合（如果可能，使用 python）

机器算法验证 scipy 插值 回归 曲线拟合 最小二乘
2022-04-17 02:57:16

我有一个从人体背面的对称分析中提取的加权二维点列表。我应该找到代表描述椎骨位置（实际上是椎骨的脊柱过程）的最可能路径的“中线”。

编辑：一个有代表性的数据集如下：

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

绘制时，我的点云通常看起来像这样，颜色代表一个点的“适合度”（根据某个标准，它的邻域有多少是对称的）：

在此处输入图像描述

编辑：彩色点的获得如下：

  1. 对于每个水平级别（具有相同 Y 坐标的点），从背面模型中获取几乎水平的切片。对于该切片中的每个点，绘制一条切线并将该切片分为右侧和左侧。然后，测量每对双边点到该线的距离之差并求和。然后，对于每个点，我们都有一个权重，它告诉切片的右侧和左侧在该点处相等的程度。下面的图片说明了这一点：

图表显示了腰椎水平的后表面切片，其中每个点都接收到一个不对称值。选择局部最小值（红色点），权重（绿色三角形）与它比附近“更好”的程度（局部不对称值的深度）成正比。

在此处输入图像描述

我想知道我应该使用哪种拟合/回归方法，以获得：

  • 局部可微性（可以近似）；
  • 能够内插到任意分辨率；
  • 没有过冲或不可预测的振荡；
  • 如果可能，去除异常值。

我怀疑可能会用最小二乘法、径向基函数、RANSAC 等来做一些事情，但我也怀疑我没有足够的知识来选择合适的方法，所以这就是我问专家的原因。

编辑：执行加权样条插值，我得到类似于下面的结果，但它对“异常值”过于敏感，并且明显远离应该是“正确”中心线的振荡：

在此处输入图像描述

非常欢迎任何帮助，感谢阅读”

2个回答

（假设你想使用 python）

最简单的（考虑到目前可用的）是使用 statsmodels 中的多项式和稳健的方法。

就像是

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

我只是输入了这个，并没有运行它。此外，我不知道在这种情况下多项式的什么阶会有用。

我会使用chebvandernumpy.polynomial，但我不记得细节，尤其是处理域的变化。

RLM.fit 有几个不同的选项，如果默认情况下没有正确处理“异常值”或点的大散布。 http://statsmodels.sourceforge.net/stable/rlm.html

分位数回归将是 RLM 的一个很好的替代方案。有人为 statsmodels 贡献了一个脚本，但它还没有被包含在内。分位数回归可用于估计在这种情况下可能运作良好的“中线”。

拟合样条曲线而不是多项式可能会更好，但是如果没有一些工作，它在 python 中是不可用的，而且我认为考虑到垂直点的数量相对较少，它不会有太大的不同。

更新：

使用第 0 列和第 1 列的数据，我尝试了 RLM。我没有设法让 RLM 工作的明显选项，但还有一些更不常见的选项我没有尝试过。

相反，我尝试将所有重量直接放在中心点上。这本质上是最小修剪平方，具有对中心点的初始稳健 (M-) 估计。第一阶段适合低阶多项式。然后我们选择残差小的点。我手动选择了阈值以获得正确的分数。然后我用加权最小二乘拟合一个高阶多项式，并只在中心点上加权。

2 阶段拟合，仅在中心点上使用 WLS

这似乎适用于估计中心点的线，但我还不知道这将如何解决更一般的问题。

除了导入 numpy 和 matplotlib.pyplot 以及加载数据，我生成图表的脚本是：

#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()

更新 2

稳健估计选项的一个问题是我们对方差估计有多少控制。RLM 更新方差，在这种情况下，它对离群观测值的权重不足。

如果我们有方差的先验估计，那么我们可以强制 RLM 使用它而不是更新它。

上面的 WLS 修剪最小二乘法估计的方差仅为 0.014 左右。如果我们使用这些信息，

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

那么拟合曲线看起来与上面的 WLS 曲线非常相似，尽管多项式的参数估计值不同。

（警告：这在 中还不起作用statsmodels，该init选项仅在我从事稳健估计和改进 RLM 的分支中可用。）

这是一篇旧文章，但这里有一篇值得考虑的论文：

“通过非线性最小化对无组织数据点进行多维曲线拟合”，Lian Fang 和 David C Gossard

http://www.cs.jhu.edu/~misha/Fall05/Papers/fang95.pdf

对于无序点，这篇论文很有趣。它描述了一种找到彼此自然相邻点的顺序的方法：“基于布朗运动的自然距离的点排序”，Philsu Kim1 和 Hyoungseok Kim2

http://downloads.hindawi.com/journals/mpe/2010/450460.pdf

其它你可能感兴趣的问题
上一篇R中的随机块设计ANOVA 下一篇剂量反应和致死剂量 50 分析