如何计算方向的轴?
以前,我根据解剖结构计算了方向轴,例如脚趾中的脚趾。
但是我发现,当我无法很好地区分脚趾或“脚跟”(蓝色方块)离开时,这不起作用。 所以我决定寻找更好的选择,我决定尝试计算惯性轴。
这个页面给出了如何计算它的一个很好的解释,但是我很难理解从质心(或者在我的情况下是压力)到一个角度的步骤。
解释归结为: 它使用压力中心和值p,其中我不知道它是什么。
我有权访问计算人体脚的这个轴的Matlab代码,并尽我所能将它转换为Python:
x = 0.508 # sensor size in the x-direction
y = 0.762 # sensor size in the y-direction
Ptot = 0 # total pressure
Px = 0 # first order moment(x)
Py = 0 # first order moment(y)
Pxx = 0 # second order moment (y)
Pyy = 0 # second order moment (x)
Pxy = 0 # second order moment (xy)
for row in range(rows): # y-direction
for col in range(cols): # x-direction
if data[row,col] > 0.0: # If not zero
temp = 1
else:
temp = 0
Ptot = Ptot + temp # Add 1 for every sensor that is nonzero
Px = Px + (x * col + x / 2) * temp
Py = Py + (y * row + y / 2) * temp
Pxx = Pxx + (x * y * y * y / 12 + x * y * (row * y + y / 2) * (row * y + y / 2) ) * temp
Pyy = Pyy + (y * x * x * x / 12 + x * y * (col * x + x / 2) * (col * x + x / 2) ) * temp
Pxy = Pxy + (x * y * (row * y + y / 2) * (row * x + x / 2)) * temp
CoPY = Py / Ptot
CoPX = Px / Ptot
CoP = [CoPX, CoPY]
Ixx = Pxx - Ptot * self.x * self.y * CoPY * CoPY
Iyy = Pyy - Ptot * self.x * self.y * CoPX * CoPX
Ixy = Pxy - Ptot * self.x * self.y * CoPY * CoPX
angle = (math.atan(2 * Ixy / (Iyy - Ixx))) / 2
Ixp = Ixx * math.cos(angle) * math.cos(angle) + Iyy * math.sin(angle) * math.sin(angle) - 2 * Ixy * math.sin(angle) * math.cos(angle)
Iyp = Iyy * math.cos(angle) * math.cos(angle) + Ixx * math.sin(angle) * math.sin(angle) + 2 * Ixy * math.sin(angle) * math.cos(angle)
RotationMatrix = [[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]]
所以就我所知,使用来自RotationMatrix的sin(角度)和cos(角度)来确定轴。 但我真的不明白如何使用这些值来绘制通过爪子的轴并围绕它旋转。
任何想法我做错了什么和/或我应该怎么做来解决它?
如果有人觉得有必要进行实验,这里有一个包含所有切片阵列的文件,其中包含每个爪子的压力数据。 clarfiy:walk_sliced_data是包含['ser_3','ser_2','sel_1','sel_2','ser_1','sel_3']的字典,它们是测量的名称。 每个度量包含另一个字典,[0,1,2,3,4,5,6,7,8,9,10](例如'sel_1'),它们表示已提取的影响。
那么,这里有一个实现与上面的代码做同样的事情(并以相关角度旋转图像)。
然而,就你的爪子而言,我不确定它会像人脚那样工作。
首先,对于一只狗的爪子来说,这样定义的“长”轴沿着爪的宽度而不是爪的长度。 只要一致,这并不重要,因为我们可以简单地旋转计算的角度而不是90° - 计算的角度。
但是,狗爪接近圆形这一事实给我们带来了更多问题。
基本上,这可能不会像狗一样对人类有用。 由图像的第二个中心矩形成的图像的协方差矩阵推导出的“长”轴的旋转(这是您的代码在上面做了什么(我认为))不太可能是对方向的精确测量爪子。
换句话说,狗的爪子接近圆形,并且它们的大部分重量都显示在脚趾上,所以在这种计算中,“后背”脚趾的重量要小于字体的重量。 因此,我们得到的轴不会一直与“后”脚趾与前脚趾的位置有关系。 (希望这是有道理的......我是一个可怕的作家......这就是为什么我回答这个问题而不是写在我应该工作的论文上......)
无论如何,足够漫不经心......这是一个例子:
import cPickle
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
def main():
measurements = cPickle.load(open('walk_sliced_data', 'r'))
plot(measurements['ser_1'].values())
plt.show()
def raw_moment(data, iord, jord):
nrows, ncols = data.shape
y, x = np.mgrid[:nrows, :ncols]
data = data * x**iord * y**jord
return data.sum()
def intertial_axis(data):
data_sum = data.sum()
m10 = raw_moment(data, 1, 0)
m01 = raw_moment(data, 0, 1)
x_bar = m10 / data_sum
y_bar = m01 / data_sum
u11 = (raw_moment(data, 1, 1) - x_bar * m01) / data_sum
u20 = (raw_moment(data, 2, 0) - x_bar * m10) / data_sum
u02 = (raw_moment(data, 0, 2) - y_bar * m01) / data_sum
angle = 0.5 * np.arctan(2 * u11 / (u20 - u02))
return x_bar, y_bar, angle
def plot(impacts):
def plot_subplot(pawprint, ax):
x_bar, y_bar, angle = intertial_axis(pawprint)
ax.imshow(pawprint)
plot_bars(x_bar, y_bar, angle, ax)
return angle
fig1 = plt.figure()
fig2 = plt.figure()
for i, impact in enumerate(impacts[:9]):
ax1 = fig1.add_subplot(3,3,i+1)
ax2 = fig2.add_subplot(3,3,i+1)
pawprint = impact.sum(axis=2)
angle = plot_subplot(pawprint, ax1)
pawprint = ndimage.rotate(pawprint, np.degrees(angle))
plot_subplot(pawprint, ax2)
fig1.suptitle('Original')
fig2.suptitle('Rotated')
def plot_bars(x_bar, y_bar, angle, ax):
def plot_bar(r, x_bar, y_bar, angle, ax, pattern):
dx = r * np.cos(angle)
dy = r * np.sin(angle)
ax.plot([x_bar - dx, x_bar, x_bar + dx],
[y_bar - dy, y_bar, y_bar + dy], pattern)
plot_bar(1, x_bar, y_bar, angle + np.radians(90), ax, 'wo-')
plot_bar(3, x_bar, y_bar, angle, ax, 'ro-')
ax.axis('image')
if __name__ == '__main__':
main()
在这些图中,中心点是图像的质心,红线定义“长”轴,而白线定义“短”轴。
原始(未旋转)爪子:
旋转爪子:
有一件事要注意这里......我只是围绕它的中心旋转图像。 (另外, scipy.ndimage.rotate
也适用于ND阵列,与2D相同,您可以轻松旋转原始3D“pawprint-over-time”数组。
如果你想旋转一点(比如质心),并将该点移动到新图像上的新位置,可以通过一些技巧在scipy的ndimage
模块中很容易地完成。 如果你愿意,我可以举个例子。 这个例子有点长,不过......