找到一组 2D 轮廓的平均表示

Question

我在实验室工作，使用小鼠进行行为研究。我想展示的部分内容是已知的刺激如何以可预测的方式引起耳朵的运动。我有一个脚本可以在前例中分割耳朵。 21 个视频帧并生成轮廓。

当轮廓被叠加并编码时[红色 = 刺激前；蓝色=刺激后]，人们可以欣赏到耳朵形状/方向的整体差异。

我想以更友好的方式展示这些视觉数据。我的想法是找到红色和蓝色轮廓的“平均”表示，然后比较两者，而不是使用 21 个轮廓进行这种混乱的业务。如果从这个“平均”轮廓中，我们还可以显示红色或蓝色轮廓集中固有方差的标准偏差，那就更好了。

为了进一步阐明，我想将这一挑战与可视化一维时间序列数据的任务联系起来。例如，本文的图 4A 将时间序列数据可视化，如下所示。许多神经荧光记录被平均为一条绿线，方差显示为灰色背景。我想知道在 2D 空间中我的轮廓是否可能出现类似的情况？

主要的挑战是，与时间序列数据相比，一个轮廓的坐标与另一个轮廓的坐标之间不存在一一对应的关系。对于一维数据，我们只需对某个 t 的所有值进行平均即可找到该点的均值和方差，然后在整个域上应用该操作。没有索引 t 关联轮廓之间的坐标对。此外，由于尺寸/形状的差异，轮廓甚至不会具有相同数量的坐标。

但是如果 2D 轮廓中的坐标有一个粗略索引 t 该怎么办？假设我们计算蓝色轮廓的质心 (CoM)？从 CoM 中，也许我们可以迭代封装轮廓的非常大的圆（r=1000，d_theta=1）的小扇区。该扇区的作用类似于参考t，因此该扇区内的所有坐标将被平均为单个坐标。如果我们迭代这个圆的所有扇区，那么我假设我们会得到“平均”二维轮廓的相当好的表示。此外，我们可以获得每个扇区的方差度量，并将其绘制为从新“平均”轮廓的每个新“平均”点切向延伸的灰色条。问题是这些条可能不会像一维时间序列情况那样看起来非常漂亮和平滑。下面是我设想的输出的粗略可视化。

如果您有任何想法，请告诉我！如果有人想玩/测试东西，我在下面留下示例代码来生成大致对齐的圆形轮廓。

import numpy as np
import matplotlib.pyplot as plt

#generate some ovals
num_ovals = 4
points = []

disp = np.random.rand(num_ovals, 2) * 0.5 # random displacements

for _ in range(num_ovals):
    center = np.random.rand(2) * 0.05
    d_x = 1 + np.random.rand() * 0.5  # Dilation noise x
    d_y = 1 + np.random.rand() * 0.5  # Dilation noise y 
    angle = np.random.rand() * 0.2  # Angle offset
    
    # generate with random dilation
    t = np.linspace(0, 2 * np.pi, 100)
    x = center[0] + d_x * 2 * np.cos(t)
    y = center[1] + d_y * np.sin(t)
    
    # apply random displacement
    x += disp[_][0]
    y += disp[_][1]
    
    # apply random rotations
    rotation_matrix = np.array([[np.cos(angle), -np.sin(angle)],
                                [np.sin(angle), np.cos(angle)]])
    points.append(np.dot(np.vstack((x, y)).T - center, rotation_matrix) + center)
    
comX, comY = np.mean(np.vstack(points), axis=0).T #center of mass

#plot figure
plt.figure(figsize=(8, 8))

for p in points:
    plt.plot(p[:, 0], p[:, 1], marker='.', markersize = 2)

plt.plot(comX, comY, marker='x', markersize=10, color='red', label='Center of Mass')
    
plt.legend()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
 

也许以与时间序列数据平均相同的方式找到平均轮廓可能是不可行的。在这种情况下，我很快就会提出我的另一个想法/好奇心。让我们忘记连接构成轮廓的所有坐标对的线，只需将轮廓视为形成 2D 点云的坐标对的集合。是否有可能找到最接近该二维点云的单个轮廓？我认为这会更容易实现，即使我没有看到任何方法以简单的方式合理地推断方差。

如果您已经阅读完毕，感谢您的阅读！

Answer 1

很快，一种缺乏准确性的解决方案是计算密度等值线，如下所示：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# stack all points
all_points = np.vstack(points)

# Get density estimation
kde = multivariate_normal(all_points.mean(axis=0), np.cov(all_points.T))
density_level = 0.03  # Set arbitrary level

x, y = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xy = np.column_stack([x.flatten(), y.flatten()])

density_values = kde.pdf(xy).reshape(x.shape)


# Plot the point cloud & average contour
plt.scatter(all_points[:, 0], all_points[:, 1], s=5, alpha=0.5)
plt.contour(x, y, density_values, levels=[density_level], colors='red')

plt.gca().set_aspect('equal', adjustable='box')
plt.show()

产生

但是，没有明确的指标表明哪种密度水平最接近点云。此外，对于我正在处理的稍微复杂的耳朵形状，轮廓图并不是最准确的轮廓。尽管对于椭圆形轮廓示例来说这是完全可以接受的。