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