Re: [eigen] Geometry module - Quaternion fitting alternative

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Review:

* in the initial comment, remove "Eigen itself is part of the KDE project."

+  const Scalar one_over_n = 1 / static_cast<Scalar>(n);

here, the day we'll have optimized complex numbers support, it would
be useful in the complex case (if you want to support it) that
one_over_n be a RealScalar. This will it make it faster to multiply by
it.

+  // computation of mean
+  VectorType src_mean = VectorType::Zero(m);
+  VectorType dst_mean = VectorType::Zero(m);
+  for (int i=0; i<n; ++i)
+  {
+    src_mean += src.col(i).cast<Scalar>();
+    dst_mean += dst.col(i).cast<Scalar>();
+  }

OK, this is a rowwise().sum() but you didn't want to #include Array
just for that, right? makes sense.

There's one thing i'm not sure about. Why do you cast<> all the time?
It seems to me that it's better to just enforce with a static
assertion that the types are the same. Indeed, casting is very costly
and kills vectorization, so it can be a source of inefficiency; so if
the user does something that requires casting, better force him to do
it himself explicitly.

+* \todo Does MatrixBase::cast<T>() become a NOP when T == MatrixBase::Scalar?
+*       In case it does not incur additional cost, mixed types seem to be fine.

yes it does since a recent commit by Gael, but still, better require
the user to cast explicitly.

+  // demeaning of src and dst points
+  MatrixType src_demean(m,n);
+  MatrixType dst_demean(m,n);
+  for (int i=0; i<n; ++i)
+  {
+    src_demean.col(i) = src.col(i).cast<Scalar>() - src_mean;
+    dst_demean.col(i) = dst.col(i).cast<Scalar>() - dst_mean;
+  }

yet more casting.again i understand that you prefer to avoid including Array.

+  // Eq. (36)-(37)
+  Scalar src_var = 0;
+  Scalar dst_var = 0;
+  for (int i=0; i<m; ++i)
+  {
+    src_var += src_demean.row(i).dot(src_demean.row(i));
+    dst_var += dst_demean.row(i).dot(dst_demean.row(i));
+  }

instead of x.dot(x), better use x.squaredNorm().

+  int rank = 0; for (int i=0; i<m; ++i) if (d(i)>0.0) ++rank;

if a singular value d(i) is 1e-15, you probably don't want to consider
it nonzero and increase the rank.
is what matters here is if it's "much smaller than" the biggest
singular value. With such a definition you would ensure that
multiplying a nonzero matrix by a nonzero scalar doesn't change its
rank.

By the way, this should definitely be done in the SVD api.

+    if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) {

oh so right away you assume that Scalar is real, right? otherwise the
numbers here are complex.

This is OK but must be stated in the dox.

*** unit-test: great.

*** general remark ***

inside Eigen, in inner loops, we tend to avoid operator() and use
coeff() and coeffRef() instead. This ensures that enabling asserts
doesn't have a too big effect on performance. here it's mostly ok, the
critical part is in SVD and matrix product outside of that file
anyway. also when everything is known at compile time (like Vector3f
v; v(1)) the assert goes away at compilation.

Feel free to push after these small remarks are adressed! Your code is great.

Benoit


2009/5/22 Hauke Heibel <hauke.heibel@xxxxxxxxxxxxxx>:
> On Thu, Apr 30, 2009 at 2:12 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> Maybe let's first migrate to git and then see how to integrate this.
>
> A test implementation has been added.
>
> http://bitbucket.org/hauke/umeyama-integration
>
> I am happy to discuss the implementation details.
>
> - Hauke
>
>
>