我想使用python找到位於圓柱形表面上的點的最佳擬合axis
。圓柱形表面上的點的最佳擬合軸
好像是scipy.linalg.svd
是查找的函數。
所以要測試,我決定從這個線程How to generate regular points on cylindrical surface產生一些點,功能makeCylinder
,並估計軸。
這是代碼:
def rotMatrixAxisAngle(axis, theta, theta2deg=False):
# Load
from math import radians, cos, sin
from numpy import array
# Convert to radians
if theta2deg:
theta = radians(theta)
#
a = cos(theta/2.0)
b, c, d = -array(axis)*sin(theta/2.0)
# Rotation matrix
R = array([ [a*a+b*b-c*c-d*d, 2.0*(b*c-a*d), 2.0*(b*d+a*c)],
[2.0*(b*c+a*d), a*a+c*c-b*b-d*d, 2.0*(c*d-a*b)],
[2.0*(b*d-a*c), 2.0*(c*d+a*b), a*a+d*d-b*b-c*c] ])
return R
def makeCylinder(radius, length, nlength, alpha, nalpha, center, orientation):
# Load
from numpy import array, allclose, linspace, tile, vstack
from numpy import pi, cos, sin, arccos, cross
from numpy.linalg import norm
# Create the length array
I = linspace(0, length, nlength)
# Create alpha array avoid duplication of endpoints
if int(alpha) == 360:
A = linspace(0, alpha, num=nalpha, endpoint=False)/180.0*pi
else:
A = linspace(0, alpha, num=nalpha)/180.0*pi
# Calculate X and Y
X = radius * cos(A)
Y = radius * sin(A)
# Tile/repeat indices so all unique pairs are present
pz = tile(I, nalpha)
px = X.repeat(nlength)
py = Y.repeat(nlength)
# Points
points = vstack((pz, px, py)).T
## Shift to center
points += array(center) - points.mean(axis=0)
# Orient tube to new vector
ovec = orientation/norm(orientation)
cylvec = array([1,0,0])
if allclose(cylvec, ovec):
return points
# Get orthogonal axis and rotation
oaxis = cross(ovec, cylvec)
rot = arccos(ovec.dot(cylvec))
R = rotMatrixAxisAngle(oaxis, rot)
return points.dot(R)
from numpy.linalg import norm
from numpy.random import rand
from scipy.linalg import svd
for i in xrange(100):
orientation = rand(3)
orientation[0] = 0
orientation /= norm(orientation)
# Generate sample points
points = makeCylinder(radius = 3.0,
length = 20.0, nlength = 20,
alpha = 360, nalpha = 30,
center = [0,0,0],
orientation = orientation)
# Least Square
uu, dd, vv = svd(points - points.mean(axis=0))
asse = vv[0]
assert abs(abs(orientation.dot(asse)) - 1) <= 1e-4, orientation.dot(asse)
正如你可以看到,我生成多個氣缸的軸線是隨機的(rand(3)
)。
有趣的是,svd
返回一個絕對完美的座標軸,如果orientation
的第一個分量爲零(orientation[0] = 0
)。
如果我評論這一行,估計的軸是關閉的。
更新1: 即使使用leastsq在圓筒式返回相同的行爲:
def bestLSQ1(points):
from numpy import array, sqrt
from scipy.optimize import leastsq
# Expand
points = array(points)
x = points[:,0]
y = points[:,1]
z = points[:,2]
# Calculate the distance of each points from the center (xc, yc, zc)
# http://geometry.puzzles.narkive.com/2HaVJ3XF/geometry-equation-of-an-arbitrary-orientated-cylinder
def calc_R(xc, yc, zc, u1, u2, u3):
return sqrt((x-xc)**2 + (y-yc)**2 + (z-zc)**2 - ((x-xc)*u1 + (y-yc)*u2 + (z-zc)*u3)**2)
# Calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc, zc)
def dist(c):
Ri = calc_R(*c)
return Ri - Ri.mean()
# Axes - Minimize residu
xM, yM, zM = points.mean(axis=0)
# Calculate the center
center, ier = leastsq(dist, (xM, yM, zM, 0, 0, 1))
xc, yc, zc, u1, u2, u3 = center
asse = u1, u2, u3
return asse
我們沒有'angoli'模塊,因此我們無法驗證'rotMatrixAxisAngle(oaxis,腐爛) '按預期工作。 –
添加rotMatrixAxisAngle(oaxis,rot)。 – yellowhat