如何在MATLAB中找到矩阵中波的任意方向[关闭]

Posted 2023-03-24

【中文标题】如何在MATLAB中找到矩阵中波的任意方向[关闭]

【英文标题】：How to find the arbitrary direction of waves in a matrix in MATLAB [closed]

【发布时间】：2020-06-19 01:53:07

【问题描述】：

我有一个二维数据，M，类似于这段代码的结果：

xv = 1:.2:10;
yv = 1:.2:5;
[X, Y] = meshgrid(xv, yv);
M = sin(X - 4 * Y) + 0.5 * rand(size(X));
figure;
imagesc(xv, yv, M)

M 显示平面的周期性变化。我想找到这波的方向。

我的方法是使用circshift，每列数据线性增加移位，并调整移位直到波对齐。但是对于最佳对齐，我没有一个好的数学定义。

是否有任何替代方法来定义和查找波浪的方向？

【问题讨论】：

这将受益于editing 中的minimal reproducible example，您的图像显示了您的最终目标应用程序，但是例如 5*5 矩阵的示例会更容易构建.

@Wolfie 我添加了一个小例子。

您的预期输出是什么？假设您平均沿 45 度或 60 度？否则，很难验证任何提议的解决方案是否按预期工作。您的示例仍然很大，无法在这种情况下进行任何手动计算以进行验证...您提出了一种方法-这行不通吗？

我想我的观点是，当角度不是 45 度的倍数时，如何通过矩阵定义角度？当它是 90 的倍数时，这是微不足道的，对于（沿对角线）之间的 45 度是可以理解的，但是对于您要求的任意角度，您是否希望在多个计算中包含元素，或者是否应该没有重叠？当您因为跨越多个而没有列/行计数时，您如何定义应该计算多少平均值？如果没有预期的输出示例，这个简介就太宽泛了...

查看氡变换。它会计算您在许多不同角度的投影，然后您可以简单地选择对比度最高的角度。

【参考方案1】：

我认为最简单的方法是使用imrotate，然后使用正常平均值。在这里，我编写了一个小脚本，它通过蛮力迭代角度并尝试通过最大化该平均向量的 tv 范数来找到正确的：

xv = 1:.2:10;
yv = 1:.2:5;
[X, Y] = meshgrid(xv, yv);
M = sin(X - 4 * Y) + 0.5 * rand(size(X));
imagesc(xv, yv, M)

offset = -min(M(:))+1;
angles = 0:1:180; % use any step size you want
tv = zeros(size(angles));
for k = 1:numel(angles)
    % use offset to ensure we can set all artificial zeros to NaN
     imr = imrotate(M +  offset, angles(k));
     imr(imr==0) = NaN;
     imr = imr - offset;
     % measure tv norm of the mean vector
     mean_vector = nanmean(imr, 1);
     tv(k) = mean(abs(diff(mean_vector(~isnan(mean_vector)))));

     subplot(2,2,1);
     imagesc(xv, yv, imr)
     subplot(2,2,2);
     plot(mean_vector);
     subplot(2,2,3);
     plot(angles(1:k), tv(1:k),'-r')
     ylim([0, 1])
     drawnow
end

% find angle with maximal tv norm of means
[~, k] = max(tv);
angle = angles(k);
IMR = imrotate(M, angle);
imagesc(xv, yv, IMR);
subplot(2,2,4);
plot(nanmean(IMR));

【讨论】：

顺便说一句：这实际上基本上是 CrisLuengo 建议的蛮力方法，我们使用 tv norm 测量对比度。

【参考方案2】：

卷积有用吗？对于这种特殊情况，

size = 10;
a = -4; % y=-4x, direction
x = repmat([-size:1:size],2*size+1,1);
y = x';
d = 1-abs(a*x-y)/sqrt(a^2+1);
d(d<0.293)=0;
d2 = sum(d,1);
d(:,d2==0) = [];
d = d/sum(d2);
LM = conv2(M,d,'valid');
LM = [0,0,0,LM,0,0,0];

给你 46 表示方向 y=-4x，对应于 y=3 上的 46 个网格点。零分接边框无效，因此设置为 0。

但这对于角度接近 45 度的人来说并不好。

我也可以想象使用插值，如果每行上一些样本点的坐标保存到矩阵xq和yq的一列中，那么这将起作用

V = interp2(X,Y,M,xq,yq);
LM = mean(V,1);

希望这会有所帮助。