坐标仿射变换后撤消测量矢量分量的旋转

Question

问题

我有一个具有三个维度的向量场。我必须单独采样矢量的每个方向，因此三个网格在测量场内略微未对齐（sample 移动，而不是测量网格）。

当我重新调整我的测量值时，我的矢量分量不再相互正交，因为每个分量都有不同的变换。因为重新对齐也有平移，我认为每个采样点的旋转略有不同。基本上我想旋转我的样本点，但保持我的矢量方向并确保它们是正交的。

问题

如何正确“取消旋转”我的向量？

示例

示例 MatLab 代码（注意，我的真实数据中我只想在 MatLab 中进行校正，我不在 MatLab 中进行转换，只是创建假数据来举例）：

n1=10; n2=10; n3=10; %10x10x10 samples for each direction / vector component measurement
sig=5; %smoothness of measured vector field
vec = struct; A = struct; %define structs
%define measurement grid
[vec(1).X,vec(2).X,vec(3).X] = ndgrid(1:n1,1:n2,1:n3); figure; 
for itrans = 1:3
    vec(itrans).x = imgaussfilt3(rand(n1,n2,n3), sig);%make random smooth vector field
    t = rand(1,3); %random translations
    r = rand(1,3)*2*pi/50; %random pitch roll yaw
    %seperate rotation matrices
    R1 = [1,         0,         0;...
          0, cos(r(1)), sin(r(1));...
          0,-sin(r(1)), cos(r(1))];
    R2 = [cos(r(2)), 0,-sin(r(2));...
                  0, 1,         0;...
          sin(r(2)), 0, cos(r(2))];
    R3 = [cos(r(3)), sin(r(3)), 0;...
         -sin(r(3)), cos(r(3)), 0;...
                  0,         0, 1];
    %make affine component matrices
    T = eye(4); T(1,4) = t(1); T(2,4) = t(2); T(3,4) = t(3);%make translation matrix
    R = R1*R2*R3; R(4,4) = 1;%combine rotations
    S = eye(4); %no skew or scaling
    %compose affine transformation
    A(itrans).mat = T*R*S;

    %apply transformation to coordinates and plot alignment
    X = [vec(1).X(:)'; vec(2).X(:)'; vec(3).X(:)';ones(1,numel(vec(3).X))]; XX = A(itrans).mat*X;
    subplot(1,3,itrans); scatter3(vec(1).X(:),vec(2).X(:),vec(3).X(:)); hold on; title('displacement measurement 3'); scatter3(XX(1,:)',XX(2,:)',XX(3,:)', 'r'); legend('original', 'displaced')
end
%apply tranformation to data
for itrans = 1:3
    vec(itrans).xtrans = interp3(vec(itrans).x ,XX(1,:)',XX(2,:)',XX(3,:)','cubic',0);
end
%plot new vectors
figure; subplot(1,2,1); quiver3(vec(1).X,vec(2).X,vec(3).X,vec(1).x,vec(2).x,vec(3).x); title('original field')
subplot(1,2,2); quiver3(vec(1).X,vec(2).X,vec(3).X,vec(1).xtrans,vec(2).xtrans,vec(3).xtrans); legend('displaced field')

三个重新对齐的坐标，矢量场的每个分量都被平移和旋转，略有不同。

原始场，一个向量的每个分量都没有真正在同一位置测量，我通过变换坐标然后对我的测量值进行插值来进行校正。

变换后的场，一个向量的每个分量不再真正沿着它所代表的轴，它们不再相互正交。

试图用油漆在 2d 中显示问题。每个箭头表示一个测量分量，每个十字表示一个坐标系。两个蓝色箭头是我测量的两个组件，两个灰色箭头是我重新调整测量后的结果，两个橙色箭头是我在以某种方式“不旋转”并将它们组合后需要的结果。

Answer 1

这不是一个非常完整的答案，也可能不是一个值得投票的答案；但是，在给出答案的时候，还没有其他人，所以我最好把我知道的发布出来。

如果您的问题是学术问题，那么它将是研究生级别的学术向量代数问题。有人可能已经为它预制了一个简洁的公式；但是，缺少这样的公式，如果我在你的位置，在旋转之前，我可能会尝试在常规正交坐标网格上重置所有测量值。 （我最近指导了一门初级电磁学课程，在该课程中我让学生做一些类似但更简单的事情来离散拉普拉斯方程。）在重置之前，需要扩展多维泰勒级数中的连续函数，然后将附近的测量值拟合到该级数的未确定系数....换句话说，您之前的分析可能会很痛苦。

样条线致力于一个相关的想法。

您的评论提到了仿射分析。不幸的是，这对我有什么帮助并不明显。

幸运的是，一旦你完成了分析和编码，相关的计算应该很快就可以让计算机完成。

这个答案当然没有回答任何东西。给定 16 个小时左右的时间，我可能可以进行分析；然后再给我 40 小时左右的时间，我可能会想出一些代码——到那时我可能会发现在文献中已经有人以更简洁、更简单的方式解决了问题。祝你好运。

Answer 2

我计算了每个采样点D = X - XX 的位移，然后用curl() 计算了该场的旋转。这表明所有位置和所有方向都在不断旋转。所以我假设我需要对每个采样点进行不同的校正可能是错误的。

然后它可能会反转/转置旋转矩阵并将其应用于测量（而不是测量的坐标）。然后将三个测量结果的三个分量相加，即可得到校正后的数据。

我将该代码添加到循环中并绘制了位移和卷曲：

n1=10; n2=10; n3=10; sig=5; vec = struct; A = struct; D = struct; C = struct; [vec(1).X,vec(2).X,vec(3).X] = ndgrid(1:n1,1:n2,1:n3); fig1 = figure; fig2 = figure; fig3 = figure;
for itrans = 1:3
    t = rand(1,3); r = rand(1,3)*2*pi/50; vec(itrans).x = imgaussfilt3(rand(n1,n2,n3), sig);
    R1 = [1,     0,     0; 0, cos(r(1)), sin(r(1)); 0,-sin(r(1)), cos(r(1))];
    R2 = [cos(r(2)), 0,-sin(r(2)); 0, 1,         0; sin(r(2)), 0, cos(r(2))];
    R3 = [cos(r(3)), sin(r(3)), 0; -sin(r(3)), cos(r(3)), 0; 0,        0, 1];
    %make affine component matrices
    T = eye(4); T(1,4) = t(1); T(2,4) = t(2); T(3,4) = t(3); R = R1*R2*R3; R(4,4) = 1;S = eye(4);
    A(itrans).mat = T*R*S;

    %apply transformation to coordinates and plot alignment
    X = [vec(1).X(:)'; vec(2).X(:)'; vec(3).X(:)';ones(1,numel(vec(3).X))]; XX = A(itrans).mat*X;
    dcenter = sqrt(sum((X(1:3,:)-5).^2)); dcenter = dcenter./max(dcenter(:));
    figure(fig1); subplot(1,3,itrans); scatter3(vec(1).X(:),vec(2).X(:),vec(3).X(:)); hold on; title('displacement measurement 3'); scatter3(XX(1,:)',XX(2,:)',XX(3,:)', 'r'); legend('original', 'displaced')

    %calculate and plot displacement due to transformation
    D(itrans).d = X-XX;
    figure(fig2); subplot(1,3,itrans); q = quiver3(vec(1).X(:)',vec(2).X(:)',vec(3).X(:)',D(itrans).d(1,:),D(itrans).d(2,:),D(itrans).d(3,:)); title(['displacement of each voxel measurement ' num2str(itrans)])
    currentColormap = colormap(gca); [~, ~, ind] = histcounts(dcenter, size(currentColormap, 1)); cmap = uint8(ind2rgb(ind(:), currentColormap) * 255); cmap(:,:,4) = 255; cmap = permute(repmat(cmap, [1 3 1]), [2 1 3]); set(q.Head, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:3,:,:), [], 4).'); set(q.Tail, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:2,:,:), [], 4).');

    %calculate and plot rotational part of displacement field
    [C(itrans).x1,C(itrans).x2,C(itrans).x3,C(itrans).av] = curl(reshape(vec(1).X(:)', n1, n2, n3),reshape(vec(2).X(:)', n1, n2, n3),reshape(vec(3).X(:)', n1, n2, n3),reshape(D(itrans).d(1,:), n1, n2, n3),reshape(D(itrans).d(2,:), n1, n2, n3),reshape(D(itrans).d(3,:), n1, n2, n3));
    figure(fig3); subplot(1,3,itrans); q = quiver3(vec(1).X(:)',vec(2).X(:)',vec(3).X(:)',C(itrans).x1(:)',C(itrans).x2(:)',C(itrans).x3(:)'); title(['curl of each voxel measurement ' num2str(itrans)])
    currentColormap = colormap(gca); [~, ~, ind] = histcounts(dcenter, size(currentColormap, 1)); cmap = uint8(ind2rgb(ind(:), currentColormap) * 255); cmap(:,:,4) = 255; cmap = permute(repmat(cmap, [1 3 1]), [2 1 3]); set(q.Head, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:3,:,:), [], 4).'); set(q.Tail, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:2,:,:), [], 4).');
    %rotation is constant, so not changed by translation or position of sample?
end

向量的位移（color 是与边缘接近程度的函数）。

Rotational部分位移。

处处相等：

unique(C(1).x1(:))

ans =

0.1617
0.1617
0.1617
0.1617

所以我想我可以将R' 应用于每次测量，然后将它们加在一起作为校正数据。

rotmeas = zeros(3, numel(vec(1).x));
for itrans = 1:3
    %apply transformation to coordinates and interpolate
    vec(itrans).xtrans = interp3(vec(itrans).x(:) ,XX(1,:)',XX(2,:)',XX(3,:)','cubic',0);
    %get inverse rotation
    Rinv = A(itrans).mat(1:3,1:3)'
    curmeas = vec(itrans).xtrans;
    compmeas = zeros(3, numel(vec(1).x));
    compmeas(itrans,:) = curmeas;
    %apply inverse rotation to spread old component over new axes
    rotmeas = rotmeas + Rinv*compmeas;
end

然后将其重新整形为数组。

for itrans = 1:3
    vec(itrans).xcor = reshape(rotmeas(itrans,:), n1,n2,n3);
end

现在我认为三个数据集是对齐的，向量分量是正交的，并且与新的坐标系对齐。只是不知道如何检查。混淆测量的不同组成部分和旋转的维度似乎很容易。