Bug on pcl::PCA when processing two dimentional dataset

[Danfoa]

Your Environment

• Operating System and version: Ubuntu Linux 16.4

Context

When using pcl::PCA class on a 2-dimensional dataset of points, the third principal component does not always follow the right-hand-rule. I found this issue while processing a plane-projected point cloud, and noticed that some of the returned Affine3d homogenous matrices, when converted to Quaternion, returned an orientation identity. i.e. when extracting the orientation from the Eigenvectors returned, even do the rotation matrix was not the Identity, the quaternion was initialized as Identity.

Expected Behavior

Principal components should, besides of being orthogonal to each other, follow the right-hand-rule, since most geometric libraries are assuming this as a fact

Current Behavior

When pcl::PCA is being used on a two-dimensional dataset (i.g. Z coordinate for all points is the same) the third principal component does not always follow the right-hand rule, triggering unexpected behavior when using the eigenvectors as rotation matrix, in Eigen, and in ROS framework.

Code to Reproduce

The way I am using PCA is fairly normal

      pcl::PCA<PointT> pca;
      pca.setInputCloud( projected_cloud );
      Eigen::Matrix3f eigen_vectors = pca.getEigenVectors();
  
      percentage_of_variance = pca.getEigenValues()/pca.getEigenValues().sum() * 100;
      grasping_pose.translation() = pca.getMean().head(3).cast<double>();
      grasping_pose.linear() = eigen_vectors.cast<double>();

And I noticed the bug when visualizing all the homogenous transformations in space and noticed that a bunch of them had orientation equals to Identity. I am sure about this issue because if I force the rotation matrix to follow the right hand rule Z = X (x) Y --- ((x) = cross product):

grasping_pose.linear().col(2) = grasping_pose.linear().col(0).cross(grasping_pose.linear().col(1));

Then the issues dissapear.

Possible Solution

Look above.

[taketwo]

Do you have a proposal how to fix this issue inside PCL? In general case we don't know whether the dataset is 2D or not. Can we detect this by examining the smallest eigenvalue? Can you post some numbers (values/vectors) that you get so we can decide on the threshold?

[SergioRAgostinho]

In general case we don't know whether the dataset is 2D or not. Can we detect this by examining the smallest eigenvalue?

That's exactly the way to do it.

[Danfoa]

Sorry for the late answer:

This, for example, is one result that illustrates the error:

Eigen Vectors:
 0.243846  -0.969814         0
-0.969814  -0.243846         0
        0         -0         1
Eigenvalues: [0.0209367   0.0019931     0]

Considering the three eigenvectors as a rotation matrix, this matrix does not follow the right-hand rule, and thus triggers errors when using its rotation information.

It this case, we are talking of a dataset projected to the XY plane, hence it can be easily perceived that the new X-axis (or first component) will be on the 4th quadrant and the new Y axis (second component) will be on the 3th, and in order to follow the right-hand rule the new Z axis (third component) should be pointing down (be negative).

Probably this is an edge case in the algorithm, but it can be easily processed by detecting the 3rd eigenvalue is 0, and forcing the Z-axis or the "third" axis to align properly.

if (eigenvalues[2] == 0)  
    eigenvectors.col(2) = eigenvectors.col(0).cross(eigenvectors.col(1));

PS: An eigen value will always be 0 if the data has no variance in the correspondant eigenvector axis/dimention. This can be a way to detect the edge case, not only in the 3D case.

[taketwo]

OK, last question is, should we compare with exactly 0, or some small epsilon?

[Danfoa]

@taketwo Quite frankly, I did not look at the implementation you have of PCA, so I am not sure if the error can occur also on 3D sparse data, but, if nobody has experienced this issue before (assuming a lot of people are using pcl::PCA on 3D data) I am inclined to think that this issue occurs only when the input data has "fewer" dimensions than the ones expected, i.e. since PCL only accepts 3D datasets, the input that triggers the error is a two-dimensional one (one dimensional data does not make sense here). Having this in mind you should compare directly to 0 since such and eigenvalue occurs only when the dataset does not have variance in one dimension (fewer dimensions than expected).

Note that you should not modify the mean since there might be the case where the dimension with no variance is a constant not equal to 0 for all points in the dataset.

[taketwo]

OK, then we can use the snippet you posted before. A pull request is welcome, ideally with a unit test.

[Danfoa]

@taketwo I will have a look at the code and post a PR, by the end of the week.

[SergioRAgostinho]

OK, last question is, should we compare with exactly 0, or some small epsilon?

Despite @Danfoa 's remark, the check should be done with a small epsilon, data type dependent. Please see the info on numpy's linalg.matrix_rank text. It's a similar concept.

since PCL only accepts 3D datasets, the input that triggers the error is a two-dimensional one (one dimensional data does not make sense here)

The 3D assumption is prevalent because of our 3d oriented point types, but any group of points which only reside in a 2D- linear subspace are subject to the same problem you experienced. A table surface scan is an example of that, for instance.

[taketwo]

Thanks for the link, 👍 for data type dependent threshold (to be determined empirically).

[Danfoa]

Would not then make more sense to always enforce the correct "sign" of the third component/axis?.
Considering the third component is always unique in orientation (but not in sign), maybe is more efficient to just enforce this instead of detecting the missing dimension?.

To detect the missing dimension we should:

1. Check is the third eigenvalue is below the data type dependent threshold.
2. Enforce the correct vector orientation by a cross product operation.

If we consider that PCL input data can only be 2 or 3-dimensional, we can skip step 1. and do always step 2.?

PS: I am not sure whether there is a more efficient way to check if the eigenvectors follow the right-hand rule.

[SergioRAgostinho]

Would not then make more sense to always enforce the correct "sign" of the third component/axis?.

Agreed.

Assuming your eigen vectors are v1, v2, v3 normalized and in decreasing order of eigen value magnitude, you can construct

v3 = cross(v1, v2)

and this is valid for every 3D symmetric matrix (and implicitly PCA).

On a side note, I do find it curious that you use this to describe a rotation since there's at least 4 rotation matrices you can generate from these eigen vectors, due to the sign ambiguity in v1 and v2.

[Danfoa]

On a side note, I do find it curious that you use this to describe a rotation since there are at least 4 rotation matrices you can generate from these eigenvectors, due to the sign ambiguity in v1 and v2.

You are right, all components/eigenvectors are strictly ambiguous in sign, but indeed the eigenvectors can be interpreted (and I believe quite often they are) as a base for a given vector space (PCA is mostly used for vector space dimensions reduction).

Now, since PCL implies a 3D space, for me personally it makes sense that the eigenvectors as a base of a 3-dimensional space follow the right-hand rule, and thus can be used directly as a rotation matrix. Indeed, there are 4 possible bases/results, all of them collinear, and only two of them follow the right-hand-rule.

Having this in mind, the result of PCA does not need strictly to follow the right-hand rule. So, now the question should be, shall we enforce this behavior? shall we only document it? or just let it be?

Advantages of enforcing it:

• Results can be directly used to form homogenous transformations, which are, I believe, the most common use case of pcl::PCA.

Disadvantages:

• Leaving it like it is, user might have the same problem that I faced, and spend a lot of time and effort in knowing where does the problem come from :P.

@taketwo, @SergioRAgostinho Any thoughts?

[taketwo]

This behavior seems reasonable to me. Enforcing it costs nearly nothing at runtime and it can save someone a headache.

[SergioRAgostinho]

This behavior seems reasonable to me. Enforcing it costs nearly nothing at runtime and it can save someone a headache.

Same opinion from my side. Sorry for the headache it already caused you.

[Danfoa]

What do you guys think of providing an additional function to obtain a Homogenous Matrix holding the mean as translation vector and eigenvectors as a rotation matrix?.

In my opinion, if the user is aware he is obtaining a homogenous transformation, with a rotation matrix inside, then the implication of considering the eigenvectors as a rotation matrix becomes somewhat clearer.

PS: I will do the PR soon

[SergioRAgostinho]

What do you guys think of providing an additional function to obtain a Homogenous Matrix holding the mean as translation vector and eigenvectors as a rotation matrix?.

I'm inclined to reject it without having a strong opinion on the matter. I just don't think it adds much to the interface.

[Danfoa]

@SergioRAgostinho @taketwo, sorry for the delay in the PR.

With the PR #2946 I believe this issue is solved, I added the line of code that we discussed to always enforce the third axis/component proper orientation to follow the right-hand rule.