Tuesday, May 24, 2016

Use Rodrigues' Rotation to Convert 3D data to 2D

In this blog, we introduce how Rodtrigues' rotation formula can be use to convert 3D data to 2D. Let us assume that we have a bunch of 3D points, and we want to project them to x/y axis plane for easier analysis. The Matlab code below implements this function.

clear all; close all;

% collect 3D points 
xyz = [-0.485 0.130 1.489
    0.008 0.134 1.464
    -0.008 0.425 1.458
    -0.465 0.414 1.470
    -0.489 0.136 1.483
    -0.465 0.412 1.470
    0.006 0.137 1.464];

% Model them as a plane in 3D using least square method
[h,l]=size(xyz);
coef = inv(xyz'*xyz)*(xyz'*ones(h,1));
coef_unit = (coef/sqrt(coef'*coef)).';

% Rodrigues' rotation formula
in0 = coef_unit;
in1 = [0 0 1];
% Reference
in0Norm = in0 ./ sqrt(in0*in0.');
in1Norm = in1 ./ sqrt(in1*in1.');
axisRaw = cross(in0Norm, in1Norm);
axis = axisRaw./sqrt(axisRaw*axisRaw.');
angle = acos(dot(in0Norm, in1Norm));
M = [axis(1)^2+(1-axis(1)^2)*cos(angle) axis(1)*axis(2)*(1-cos(angle))-axis(3)*sin(angle) axis(1)*axis(3)*(1-cos(angle))+axis(2)*sin(angle);
axis(1)*axis(2)*(1-cos(angle))+axis(3)*sin(angle) axis(2)^2+(1-axis(2)^2)*cos(angle) axis(2)*axis(3)*(1-cos(angle))-axis(1)*sin(angle);
axis(1)*axis(3)*(1-cos(angle))-axis(2)*sin(angle) axis(2)*axis(3)*(1-cos(angle))+axis(1)*sin(angle) axis(3)^2+(1-axis(3)^2)*cos(angle)];
    
% Project 3D points to 2D plane with z=0 
xyz_ortho = xyz - xyz*coef_unit.'*coef_unit;   
xyz_prj = M*xyz_ortho.';

Let us walk through the code step-by-step.

Step 1: A plane in 3D space can be characterized as \( a x + b y + c z = 1 \). Assuming that matrix \(X\) has all 3D points: \[ X = \begin{bmatrix}
    x_{0}       & y_{0} & z_{0}  \\
    x_{1}       & y_{1} & z_{1} \\
    ...              & ...        &...\\
    x_{N}       & y_{N} & z_{N}
\end{bmatrix}\]
Then least square estimator for [a b c] is \( coef = (X'X)^{-1}(X'h)\) where \(h = [1 1 ... 1]'\) and \(X'\) is the transpose of \(X\). Normalized version of coef is called coef_unit.

Step 2: Use Rodrigues' rotation formula to find the matrix \(M\) which can rotate coef_unit to [0 0 1].

Step 3: First we project all 3D points to a plane characterized by [x y z]*coef_unit = 0. This is achieved by:
xyz_ortho = xyz - xyz*coef_unit.'*coef_unit;   
The last step is to use rotation matrix \(M\) to reduce the dimension of xyz_ortho to 2D. After the rotation, all z-axis elements in xyz_prj are 0's.




No comments:

Post a Comment