undo rotation of measured vector components after affine transformation of coordinates

Question

Problem

I have a vector field with three dimensions. Each direction of the vector I have to sample seperately and therefore the three grids are slightly unaligned within the measured field (the sample moves, not the measurement grid).

When I realign my measurements my vector components are no longer orthogonal to each other since each has a different transformation. Because the realignment also has a translation, I think each sample point has a slightly different rotation. Basically I want to rotate my sample points but keep my vector directions and make sure they are orthogonal.

Question

How do I correctly 'unrotate' my vectors?

Example

Example MatLab code (note that in my real data I only want to perform the correction in MatLab, I dont perform the transformation in MatLab, just creating fake data to give an example):

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')

The three realigned coordinates, each component of the vector field has been translated and rotated slightly differently.

The original field, each component of one vector was not really measured at the same location, which I correct for by transforming the coordinates and then interpolating my measurements.

The transformed field, each component of one vector is no longer really along the axis it represents, and they are no longer orthogonal to each other.

Trying to show the problem in 2d with paint. Each arrow shows a measured component, each cross shows a coordinate system. The two blue arrows are the two components that I measure, the two gray arrows are my result after realigning my measurements, the two orange arrows are what I need after somehow 'unrotating' and combining them.

Answer 1

I calculated the displacement for each sampling point D = X - XX and then calculated the rotation of that field with curl(). That showed a constant rotation for all locations and all directions. So probably I was wrong in assuming that I needed a different correction for each sampling point.

Then its probably possible to invert/transpose the rotation matrix and apply it to the measurement (instead of the coordinates of the measurement). Then add the three resulting components of the three measurements together would then give the corrected data.

I added that code to the loop and plotted the displacement and the 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

Displacement of vectors (color as a function of how close to the edge).

Rotational part of displacement.

Which is equal everywhere:

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

ans =

0.1617
0.1617
0.1617
0.1617

So I think I can just apply R' to each measurement and then add them together as the corrected data.

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

Then just reshape it to an array again.

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

Now I think the three data sets are aligned and the vector components are orthogonal and aligned with the new coordinate system. Just not sure how to check that. It seems easy to mix up the different components of the measurements and the dimensions of the rotations.

Answer 2

This isn't a very complete answer, and it probably isn't an answer that merits votes; but, at the time the answer is given, there aren't any others, so I had better post what I know.

If your problem were an academic problem, then it would be a graduate-level academic vector-algebra problem. Someone may already have prefabricated a neat formula for it; but, lacking such a formula, if I were in your place, before rotating, I might try to reset all measurements on a regular orthonormal coordinate grid. (I recently instructed a junior-level electromagnetics course in which I had the students do something vaguely similar, but simpler, to discretize Laplace's equation.) Before resetting, one would need to expand each orthonormal component of the continuous function in a multidimensional Taylor series, then fit the nearby measured values to the undetermined coefficients of the series.... In other words, the analysis before you might be painful.

Splines work on a related idea.

Your comment mentioned an affine analysis. Unfortunately, it is not obvious to me how that would help.

Fortunately, once you had done the analysis and coding, the associated calculations should be pretty quick for the computer to complete.

This answer doesn't answer anything, of course. Given 16 hours or so, I could probably do the analysis; then given another 40 hours or so, I could probably come up with some code—by which time I might discover where in the literature someone had already solved the problem in a neater, simpler way. Good luck.