晶须的匹配轨迹

Question

我正在进行胡须追踪实验。我有老鼠拍打物体的高速视频 (500fps)。在每个这样的视频中，我都跟踪了老鼠鼻子和胡须的形状。由于跟踪存在噪声，因此每帧中的胡须数量可能不同（请参阅附图中的 2 个连续帧，请注意黄色假阳性胡须出现在左侧帧中，而不是右侧帧中）。

参见示例 1：

作为跟踪的最终结果，对于每一帧，我得到了不同数量的可变长度向量；每个向量对应一个晶须。在这一点上，我想匹配胡须在帧之间。我曾尝试使用 Matlab 的示例 align 来执行此操作，但它只能正常工作。它的结果附在下面（在附图中显示了超过 227 帧的所有晶须的基点）。

参见示例 2：

我想运行一些算法来正确聚类胡须，这样每个胡须都被识别为自己，并在许多帧的过程中与其他胡须分开。换句话说，我希望将第二张图像中的每条略微正弦的轨迹都识别为一个轨迹。无论我使用哪种排序算法，都应该考虑到胡须可能会在连续帧之间消失并重新出现。不幸的是，我完全没有想法......

有什么帮助吗？

再次请记住，对于附图 2 中的每个点，我都有很多数据点，因为这只是胡须基点的图，而实际上我有整个胡须长度的数据。

Answer 1

这就是我处理问题的方式。假设不同大小的数据向量在称为dataVectors 的cell 类型中，并且知道晶须的数量（nSignals），我会尝试将数据扩展到从原始数据派生的第二维，然后执行二维上的 k 均值。

所以，首先我会获得向量的最大大小，以便将数据转换为矩阵并执行NaN-padding。

maxSize = -Inf;
for k = 1:nSignals
    if length(dataVectors{k}.data) > maxSize
        maxSize = length(dataVectors{k}.data);
    end
end

现在，我将通过将数据提升到 2 次方（或 3 次方，您的选择）来制作 2D 数据。这只是一个非常简单的转换。但是您也可以在这里使用kernel methods 并将每个向量与其余向量进行投影；但是，我认为这没有必要，如果您的数据真的很大，它可能效率低下。目前，将数据提高到二的幂应该可以解决问题。结果存储在第二维中。

projDegree = 2;
projData = zeros(nSignals, maxSize, 2).*NaN;
for k = 1:nSignals
    vecSize = length(dataVectors{k}.data);
    projData(k, 1:vecSize, 1) = dataVectors{k}.data;
    projData(k, 1:vecSize, 2) = dataVectors{k}.data.*projDegree;
end
projData = reshape(projData, [], 2);

在这里，projData 将在行 1 和列 1 中具有第一个晶须的原始数据（或我在这里称之为信号），而列 2 将具有新维度。假设您总共有8 晶须，那么projData 将拥有1、9、17 等行中第一个晶须的数据。 2、10、18 等行中的第二个晶须的数据。如果您想以自己的方式返回原始数据，这一点很重要。此外，您可以尝试使用不同的projDegrees，但我怀疑它会产生很大的不同。

现在我们对二维数据执行k-means；但是，我们提供了初始点，而不是让它用 k-means++ 确定它们。正如我在此提出的，初始点是每个晶须的每个向量的第一个数据点。以这种方式，k-means 将从那里离开并相应地移动到集群均值。我们将结果保存在idxK。

idxK = kmeans(projData,nSignals, 'Start', projData(1:nSignals, :));

你有它。变量idxK 会告诉你哪个数据点属于哪个集群。

以下是我提出的解决方案的工作示例。第一部分只是尝试生成看起来像您的数据的数据，您可以跳过它。

rng(9, 'twister')
nSignals = 8;   % number of whiskers
n = 1000;       % number of data points
allData = zeros(nSignals, n);   % all the data will be stored here

% this loop will just generate some data that looks like yours
for k = 1:nSignals
    x = sort(rand(1,n));
    nPeriods = round(rand*9)+1;     % the sin can have between 1-10 periods
    nShiftAmount = round(randn*30);     % shift between ~ -100 to +100
    y = sin(x*2*pi*nPeriods) + (randn(1,n).*0.5);
    y = y + nShiftAmount;
    allData(k, :) = y;
end
nanIdx = round(rand(1, round(n*0.05)*nSignals).*((n*nSignals)-1))+1;   
allData(nanIdx) = NaN;      % about 5% of the data is now missing

figure(1);
for k = 1:nSignals
    nanIdx = ~isnan(allData(k, :));
    dataVectors{k}.data = allData(k, nanIdx);
    plot(dataVectors{k}.data, 'kx'), hold on;
end

% determine the max size 
maxSize = -Inf;
for k = 1:nSignals
    if length(dataVectors{k}.data) > maxSize
        maxSize = length(dataVectors{k}.data);
    end
end

% making the data now into two dimensions and NaN pad
projDegree = 2;
projData = zeros(nSignals, maxSize, 2).*NaN;
for k = 1:nSignals
    vecSize = length(dataVectors{k}.data);
    projData(k, 1:vecSize, 1) = dataVectors{k}.data;
    projData(k, 1:vecSize, 2) = dataVectors{k}.data.*projDegree;
end
projData = reshape(projData, [], 2);
figure(2); plot(projData(:,1), projData(:,2), 'kx');

% run k-means using the first points of all measure as the initial points
idxK = kmeans(projData,nSignals, 'Start', projData(1:nSignals, :));
figure(3);
liColors = [{'yx'},{'mx'},{'cx'},{'bx'},{'kx'},{'gx'},{'rx'},{'gd'}];
for k = 1:nSignals
    plot(projData(idxK==k,1), projData(idxK==k,2), liColors{k}), hold on;
end

% plot results on original data
figure(4);
for k = 1:nSignals
    plot(projData(idxK==k,1), liColors{k}), hold on;
end

如果这有帮助，请告诉我。