[Registration] Rotation epsilon enhancement : (sin(angle))^2 instead of cos(angle)

Describe the bug

Currently, transformation rotation epsilon (maximum allowable rotation difference between two consecutive transformations) must
be set as cos(angle) value.

(in registration.h)

  /** \brief Set the transformation rotation epsilon (maximum allowable rotation
   * difference between two consecutive transformations) in order for an optimization to
   * be considered as having converged to the final solution. \param[in] epsilon the
   * transformation rotation epsilon in order for an optimization to be considered as
   * having converged to the final solution (epsilon is the cos(angle) in a axis-angle
   * representation).
   */
  inline void
  setTransformationRotationEpsilon(double epsilon)
  {
    transformation_rotation_epsilon_ = epsilon;
  }

But in common, cos(angle) is not good for using as threshold value around angle is near to 0.
I think sin(angle) is better for it.

Context

When angle is close to 0, cos(angle) is close to 1.0, means 0.9999999....
It waists significand bits of floating point value.
In addition, the slope of the Cosine function approaches 0 if angle is close to 0.
It means we'll need many significand bits to distinguish cos(angle) value if angle is very close to 0.

https://en.wikipedia.org/wiki/Sine_and_cosine

Instead of cos(angle) we can easily calculate (sin(angle))^2 value from rotation matrix.
I think it should be much better than cos(angle).

I modified NormalDistributionsTransform::computeTransformation() like as follows.

https://github.com/PointCloudLibrary/pcl/blob/master/registration/include/pcl/registration/impl/ndt.hpp

    delta *= delta_norm;

    // Convert delta into matrix form
    convertTransform(delta, transformation_);

    transform += delta;

    // Update Visualizer (untested)
    if (update_visualizer_)
      update_visualizer_(output, pcl::Indices(), *target_, pcl::Indices());

    const double cos_angle =
        0.5 * (transformation_.template block<3, 3>(0, 0).trace() - 1);
    const double translation_sqr =
        transformation_.template block<3, 1>(0, 3).squaredNorm();

    const double sin_angle_sqr = ((transformation_(2, 1) - transformation_(1, 2)) * (transformation_(2, 1) - transformation_(1, 2))
                                + (transformation_(0, 2) - transformation_(2, 0)) * (transformation_(0, 2) - transformation_(2, 0))
                                + (transformation_(1, 0) - transformation_(0, 1)) * (transformation_(1, 0) - transformation_(0, 1))) / 4.0;

    // FOR DEBUG
    printf ("%g %g %g %g\n", score, translation_sqr, cos_angle, sin_angle_sqr);

    nr_iterations_++;

sin_angle_sqr was added and output log to compare other values. (translation_sqr, cos_angle)
This calculation formula can be easily derived from Rodrigues' Rotation Formula. (https://mathworld.wolfram.com/RodriguesRotationFormula.html)

I tried NDT matching using following example, I posted before.
https://github.com/PointCloudLibrary/pcl/issues/5056#issuecomment-988864151
I changed setTransformationEpsilon to 0 to check more iteration loop.

87862.4 0.988927 0.994469 0.0110306
126489 0.0891851 0.997194 0.00560402
128543 0.0202943 0.999676 0.000648835
129053 0.00153224 0.999974 5.24327e-05
129071 6.64301e-06 1 3.91958e-07
129071 5.71179e-14 1 2.97089e-15
129071 5.65053e-14 1 2.93572e-15
129071 5.58576e-14 1 2.90353e-15
129071 5.51443e-14 1 2.86564e-15
129071 5.44038e-14 1 2.82411e-15
129071 5.37615e-14 1 2.79186e-15
129071 5.29458e-14 1 2.74898e-15
129071 5.22264e-14 1 2.70875e-15
129071 8.35703e-38 1 4.33727e-39
129071 8.35703e-38 1 4.33727e-39
129071 8.35703e-38 1 4.33727e-39
129071 8.35703e-38 1 4.33727e-39
...

cos(angle) became 1.0 in early stage, but (sin(angle))^2 can show small value.
So it will be useful if the user would like to use more small rotation epsilon value, I supposed.

Your Environment (please complete the following information):

• OS: Windows10
• Compiler: Compiler: Visual Studio 2019 Build Tools
• PCL Version: HEAD

Start with registration.h to understand the transformation rotation epsilon contract, then read NormalDistributionsTransform::computeTransformation() in registration/include/pcl/registration/impl/ndt.hpp. Compare the existing cosine-based convergence calculation with the reported sin_angle_sqr values from the linked NDT example. Done means the rotation epsilon behavior is consistently defined for small angles and validated against the example's convergence behavior.