时间戳

首页 > 解决方案 > 正确实现 scipy.optimize 以求解平面方程的所有参数

python - 正确实现 scipy.optimize 以求解平面方程的所有参数

问题描述

一个三维平面方程定义为： 

a * x + b * y + c * z + d = 0

其中 x、y 和 z 是数据点的坐标，a、b、c 和 d 定义平面参数。我有一堆 x,y,z 数据点，我知道它们大致形成了一个平面。现在我想将这些数据点拟合到一个平面上并找到 a、b、c、d 参数。到目前为止我有这个功能：

import scipy.optimize as op
def _plane(ws,x,y,z):
    cost = np.sum(ws[0]*x + ws[1]*y + ws[2]*z + ws[3])
    print cost, np.sum(ws)
    return cost
out = op.minimize(_plane,ws,args=(x,y,z),method='SLSQP',options={'maxiter':1000, 'disp':1})

ws = plane([np.mean(x),np.mean(y),np.mean(z),0.001],x,y,z)

输出：

130.78467 -0.3011288588643074
130.78467 -0.3011288588643074
130.78467 -0.3011288439631462
130.78467 -0.3011288439631462
130.78468 -0.3011288439631462
130.78467 -0.3011288439631462
-119765.375 -1024.3011288588643
-119765.375 -1024.3011288588643
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -119765.375
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
out:      fun: -119765.375
     jac: array([0., 0., 0., 0.])
 message: 'Optimization terminated successfully.'
    nfev: 12
     nit: 2
    njev: 2
  status: 0
 success: True
       x: array([-4.07030255e-01, -8.12448636e-02, -1.02381385e+03,  1.00000000e-03])

正如你所看到cost的，它总是在减少，d无论我最初给出什么，它的值都不会改变。这个实现正确吗？

更新：

使用cost = abs(np.sum(ws[0]*x + ws[1]*y + ws[2]*z + ws[3]))可解决优化问题，设置method='cg'可提供最佳结果。

标签： pythonscipy

解决方案

对于将平面调整为一组点的具体问题，可以使用更简单的方法：

>>> import numpy as np
>>> points = np.array([(3,0,0), (0,4,0), (3,0,1),
...                    (0,4,1), (6,-4,0), (6,-4,1)])

计算质心并从每个点中减去它。

>>> centroid = np.average(points, axis=0)
>>> centered = points - centroid

计算中心点数组的单值分解。第三个奇异值对应的向量垂直于平面，范数为1，所以可以直接用它作为平面的a,b,c系数。

>>> _, _, vh = np.linalg.svd(centered)
>>> a, b, c = normal = vh[2]

质心必须位于平面上，因此可以从法线向量和质心计算出最后一个系数 (d)。

>>> d = -normal @ centroid
>>> plane_coefficients = a, b, c, d
>>> print(plane_coefficients)
(0.8, 0.6, 0.0, -2.4000000000000004)

推荐阅读