撤销坐标仿射变换后测量矢量分量的旋转
undo rotation of measured vector components after affine transformation of coordinates
问题
我有一个三维向量场。矢量的每个方向我都必须单独采样，因此三个网格在测量场内略微未对齐（sample 移动，而不是测量网格）。
当我重新调整我的测量值时，我的矢量分量不再相互正交，因为每个分量都有不同的变换。因为重新对齐也有平移，我 认为 每个样本点的旋转略有不同。基本上我想旋转我的样本点但保持我的矢量方向并确保它们是正交的。
问题
如何正确 'unrotate' 我的载体？
例子
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')
三个重新对齐的坐标，矢量场的每个分量都进行了稍微不同的平移和旋转。
原始场，一个向量的每个分量并不是在同一位置真正测量的，我通过转换坐标然后对我的测量结果进行插值来校正。
变换后的场，一个向量的每个分量不再真正沿着它代表的轴，它们不再相互正交。
试图用绘画在二维中显示问题。每个箭头表示一个被测组件，每个十字表示一个坐标系。两个蓝色箭头是我测量的两个组件，两个灰色箭头是我重新调整测量后的结果，两个橙色箭头是我以某种方式 'unrotating' 并将它们组合后所需要的。
这不是一个非常完整的答案，也可能不是一个值得投票的答案；但是，在给出答案时，没有其他答案，所以我最好 post 我所知道的。
如果您的问题是学术问题，那么它将是研究生水平的学术矢量代数问题。有人可能已经为它预制了一个简洁的公式；但是，如果没有这样的公式，如果我处在你的位置，在旋转之前，我可能会尝试在规则的正交坐标网格上重置所有测量值。 （我最近指导了一门初级电磁学课程，在课程中我让学生做一些类似但更简单的事情来离散化拉普拉斯方程。）在重置之前，需要展开每个 正交 多维泰勒级数中连续函数的分量，然后将附近的测量值拟合到该级数的待定系数......换句话说，你面前的分析可能会很痛苦。
Splines 致力于相关的想法。
您的评论提到了仿射分析。不幸的是，我不清楚这会有什么帮助。
幸运的是，一旦你完成了分析和编码，相关的计算应该很快就能让计算机完成。
当然，这个答案没有回答任何问题。给定 16 个小时左右，我可能可以进行分析；然后再给我 40 个小时左右的时间，我可能会想出一些代码——到那时我可能会发现文献中的哪些地方已经有人以更简洁、更简单的方式解决了问题。祝你好运。
我计算了每个采样点的位移D = X - XX，然后用curl()计算了那个场的旋转。这表明所有位置和所有方向都在不断旋转。所以我假设我需要对每个采样点进行不同的校正可能是错误的。
然后可能 invert/transpose 旋转矩阵并将其应用于测量（而不是测量的坐标）。然后将三个测量的三个结果分量加在一起，然后给出校正后的数据。
我将该代码添加到循环中并绘制了位移和卷曲：
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
矢量位移（ 作为与边缘距离的函数）。
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
现在我觉得三个数据集是对齐的，向量分量是正交的，对齐新的坐标系。只是不确定如何检查。混淆测量的不同分量和旋转的维度似乎很容易。
问题
我有一个三维向量场。矢量的每个方向我都必须单独采样，因此三个网格在测量场内略微未对齐（sample 移动，而不是测量网格）。
当我重新调整我的测量值时，我的矢量分量不再相互正交，因为每个分量都有不同的变换。因为重新对齐也有平移，我 认为 每个样本点的旋转略有不同。基本上我想旋转我的样本点但保持我的矢量方向并确保它们是正交的。
问题
如何正确 'unrotate' 我的载体？
例子
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')
三个重新对齐的坐标，矢量场的每个分量都进行了稍微不同的平移和旋转。
原始场，一个向量的每个分量并不是在同一位置真正测量的，我通过转换坐标然后对我的测量结果进行插值来校正。
变换后的场，一个向量的每个分量不再真正沿着它代表的轴，它们不再相互正交。
试图用绘画在二维中显示问题。每个箭头表示一个被测组件，每个十字表示一个坐标系。两个蓝色箭头是我测量的两个组件，两个灰色箭头是我重新调整测量后的结果，两个橙色箭头是我以某种方式 'unrotating' 并将它们组合后所需要的。
这不是一个非常完整的答案，也可能不是一个值得投票的答案；但是，在给出答案时，没有其他答案，所以我最好 post 我所知道的。
如果您的问题是学术问题，那么它将是研究生水平的学术矢量代数问题。有人可能已经为它预制了一个简洁的公式；但是，如果没有这样的公式，如果我处在你的位置，在旋转之前，我可能会尝试在规则的正交坐标网格上重置所有测量值。 （我最近指导了一门初级电磁学课程，在课程中我让学生做一些类似但更简单的事情来离散化拉普拉斯方程。）在重置之前，需要展开每个 正交 多维泰勒级数中连续函数的分量，然后将附近的测量值拟合到该级数的待定系数......换句话说，你面前的分析可能会很痛苦。
Splines 致力于相关的想法。
您的评论提到了仿射分析。不幸的是，我不清楚这会有什么帮助。
幸运的是，一旦你完成了分析和编码，相关的计算应该很快就能让计算机完成。
当然，这个答案没有回答任何问题。给定 16 个小时左右，我可能可以进行分析；然后再给我 40 个小时左右的时间，我可能会想出一些代码——到那时我可能会发现文献中的哪些地方已经有人以更简洁、更简单的方式解决了问题。祝你好运。
我计算了每个采样点的位移D = X - XX，然后用curl()计算了那个场的旋转。这表明所有位置和所有方向都在不断旋转。所以我假设我需要对每个采样点进行不同的校正可能是错误的。
然后可能 invert/transpose 旋转矩阵并将其应用于测量（而不是测量的坐标）。然后将三个测量的三个结果分量加在一起，然后给出校正后的数据。
我将该代码添加到循环中并绘制了位移和卷曲：
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
矢量位移（
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
现在我觉得三个数据集是对齐的，向量分量是正交的，对齐新的坐标系。只是不确定如何检查。混淆测量的不同分量和旋转的维度似乎很容易。