[Registration] gauss_c2 calculation formula in NormalDistributionsTransform

Describe the bug

Thanks a lot for great library.
I found a suspicious code in ndt.hpp to calculate "gauss_c2" value.
But I'm not sure that what is the correct answer, so I would like to hear professional member's opinion...

Context

Following code is in ndt.hpp to calculate NDT related parameters.

  // Initializes the gaussian fitting parameters (eq. 6.8) [Magnusson 2009]
  const double gauss_c1 = 10.0 * (1 - outlier_ratio_);
  const double gauss_c2 = outlier_ratio_ / pow(resolution_, 3); // *** HERE ***
  const double gauss_d3 = -std::log(gauss_c2);
  gauss_d1_ = -std::log(gauss_c1 + gauss_c2) - gauss_d3;
  gauss_d2_ =
      -2 * std::log((-std::log(gauss_c1 * std::exp(-0.5) + gauss_c2) - gauss_d3) /
                    gauss_d1_);

But I wondered why resolution^3 is used to calculate gauss_c2.
For example, if resolution value is small, gauss_c2 can become very big value, and
NDT tends to work incorrectly.
Several paper pointed out that the user need to decide resolution carefully, but I supposed that
this code may be one of the cause of difficulty to decide best resolution parameter.
I checked Dr. M.Magnusson's paper, but I think this calculation formula is not written clearly, I supposed.

I found a following code, seems to be for D2D (Distribution-to-Distribution) NDT paper, Dr. M.Magnusson is one of the author.
https://gitsvn-nt.oru.se/software/ndt_core_public/-/blob/master/ndt_registration/src/ndt_matcher_p2d.cpp

On this code, resolution is not used to calculate gauss_c2.

    support_size = 4; //???
    ...
    lfc2 = outlier_ratio/pow(support_size,3);

On this code, "support_size" parameter is used instead of resolution, and support_size = 4 fixed value.

I tentatively checked using NDT tutorial (https://pointclouds.org/documentation/tutorials/normal_distributions_transform.html)
using following parameter. (Note that resolution is small (0.5m = 50cm))

  // Setting scale dependent NDT parameters
  // Setting minimum transformation difference for termination condition.
  ndt.setTransformationEpsilon (0.000001);
  // Setting maximum step size for More-Thuente line search.
  ndt.setStepSize (1.0);
  //Setting Resolution of NDT grid structure (VoxelGridCovariance).
  ndt.setResolution (0.5);
  // Setting max number of registration iterations.
  ndt.setMaximumIterations (35);

Current PCL code. (gauss_c2 = outlier_ratio_ / pow(resolution_, 3);)
NDT matching failed.

Using modified code. (gauss_c2 = outlier_ratio_ / pow(4.0, 3);)
NDT matching succeeded.

I'm not sure that above code (gauss_c2 = outlier_ratio_ / pow(4.0, 3);) is best and correct.
But I think it's better than current PCL code. (gauss_c2 = outlier_ratio_ / pow(resolution_, 3);)

What do you think?
Again, I would like to hear professional member's opinion to improve this point...

Your Environment (please complete the following information):

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

Start in ndt.hpp and reproduce the tutorial configuration with resolution 0.5, then compare its gauss_c2 calculation with the cited ndt_matcher_p2d.cpp implementation and Magnusson's paper. Done means establishing the correct formula and validating its effect on the reported NDT matching behavior.