The Result of Certification Problem

[banbanhank]

Hi, when I tried the performance of certification, the η which the certification estimated was not good. I think the η value is between 0 and 1, but the result is 1233.39. I would like to know if I run the code in wrong way. Is the theta value I gave right?
The program of cetification is shown below:

  teaser::DRSCertifier::Params params2; //default
  params2.noise_bound = 0.05;
  params2.cbar2 = 1;
  params2.max_iterations = 100;
  Eigen::Matrix<bool, 1, Eigen::Dynamic> theta;
  theta = solver.getRotationInliersMask();
  teaser::DRSCertifier cert(params2);
  teaser::CertificationResult result;
  result = cert.certify(solution.rotation, src, tgt, theta);
  cout << "==========================certifiable===================================" << endl;
  cout << "The certification result is " << endl;
  cout << "The optimal result is ";
  if (result.is_optimal)
    cout << "success" << endl;
  else 
    cout << "fail" << endl;
  cout << "The best_suboptimality is " << result.best_suboptimality << endl;

---

The steps how I find the transformation

I implemented scan matching with the lidar sensor and the dataset is from the nuScence. The source and target point clouds are identical. The target point cloud is obtained by applying a certain transformation on the source point cloud. And then I added normal distributed noise on the target. The correspondence is given as the indices are unchanged. The total number of correspondence is 1991 and the number of outlier is 1000. The performance is good.
The experiment is shown below:

  teaser::RobustRegistrationSolver::Params params;
  params.noise_bound = 0.05;
  params.cbar2 = 1;
  params.estimate_scaling = false;
  params.rotation_max_iterations = 100;
  params.rotation_gnc_factor = 1.4;
  params.rotation_estimation_algorithm =
  teaser::RobustRegistrationSolver::ROTATION_ESTIMATION_ALGORITHM::GNC_TLS;
  params.rotation_cost_threshold = 0.005;
  teaser::RobustRegistrationSolver solver(params);
  solver.solve(src, tgt);
  auto solution = solver.getSolution();

(The red point cloud and the green point cloud are the source and target, respectively.)

(The matching of source and target is perfect)

[hankyang94]

Hi, in order to use the certification module, some preprocessing is required to pass the correct problem formulation to the solver.

Essentially, the certification module tries to certify the “rotation” sub-problem (aka the robust Wahba problem). Therefore, instead of passing the entire point clouds, one needs to pass the translation-invariant measurements (TIMs) for certification:

R_est = solution.rotation
A_TIMs = solver.getMaxCliqueSrcTIMs()
B_TIMs = solver.getMaxCliqueDstTIMs()
theta_TIMs_raw = solver.getRotationInliersMask()
theta_TIMs = getBinaryTheta(theta_TIMs_raw) 

Where the getBinaryTheta function simply converts the “1-0” mask to “+1” and “-1” binary values:

def getBinaryTheta(theta):
    N = theta.shape[0]
    output = -1.0 * np.ones(N)
    output[theta] = 1.0 
    return output

Now the certifier should work:

certificate = certifier.certify(R_est,A_TIMs,B_TIMs,theta_TIMs)

[banbanhank]

@hankyang94, thanks a lot for your reply. I had corrected my input. I replaced the source cloud and target cloud with the TIMs of source and target like below:

  teaser::DRSCertifier::Params params2; //default
  params2.noise_bound = 0.05;
  params2.cbar2 = 1;
  params2.max_iterations = 100;
  Eigen::Matrix<bool, 1, Eigen::Dynamic> theta;
  Eigen::Matrix<double, 3, Eigen::Dynamic> src_TIM, tgt_TIM;
  src_TIM = solver.getMaxCliqueSrcTIMs();
  tgt_TIM = solver.getMaxCliqueDstTIMs();
  theta = solver.getRotationInliersMask();

  teaser::DRSCertifier cert(params2);
  teaser::CertificationResult result;
  result = cert.certify(solution.rotation, src_TIM, tgt_TIM, theta);
  cout << "==========================certifiable===================================" << endl;
  cout << "The certification result is " << endl;
  cout << "The optimal result is ";
  if (result.is_optimal)
    cout << "success" << endl;
  else 
    cout << "fail" << endl;
  cout << "The best_suboptimality is " << result.best_suboptimality << endl;

And the result is shown below:

The best suboptimality is 72.7375 and it is much smaller than before. However, I think the value should be between 0 and 1. I would like to know if there are still other wrongs with my program. By the way, I set max_iterations = 100. It seemed that the suboptimality might converge more smaller if I set max_iterations bigger.

[jingnanshi]

@banbanhank Try manually convert theta to an Eigen matrix with double type. There might be a bug in the interface that takes in the bool matrix.

[banbanhank]

@jingnanshi, I have tried to manually convert theta to an Eigen matrix with double type. However, the problem is still the same.
Here is my program shown below:

  teaser::DRSCertifier::Params params2; //default
  params2.noise_bound = 0.05;
  params2.cbar2 = 1;
  params2.max_iterations = 200;
  Eigen::Matrix<bool, 1, Eigen::Dynamic> theta;
  Eigen::Matrix<double, 3, Eigen::Dynamic> src_TIM, tgt_TIM;
  src_TIM = solver.getMaxCliqueSrcTIMs();
  tgt_TIM = solver.getMaxCliqueDstTIMs();
  theta = solver.getRotationInliersMask();
  //convert theta to an Eigen matrix with double type
  Eigen::Matrix<double, 1, Eigen::Dynamic> theta_double(1, theta.cols());
  for (size_t i = 0; i < theta.cols(); ++i) {
    if (theta(i)) {
      theta_double(i) = 1;
    } else {
      theta_double(i) = -1;
    }
  }
  teaser::DRSCertifier cert(params2);
  teaser::CertificationResult result;
  result = cert.certify(solution.rotation, src_TIM, tgt_TIM, theta_double);
  cout << "==========================certifiable result===================================" << endl;
  cout << "The certification result is " << endl;
  cout << "The optimal result is ";
  if (result.is_optimal)
    cout << "success" << endl;
  else 
    cout << "fail" << endl;
  cout << "The best_suboptimality is " << result.best_suboptimality << endl;

Do you mind to try my program and my data if you have time, please? I would like to know how to let the certification module run successfully. Thanks. teaser_demo.zip

[hankyang94]

@banbanhank I have tried your data, and I believe the reason why DRS certifier failed on your problem is because your input point cloud is not normalized.

Below I attached a visualization of the translation-invariant measurements (TIMs) from your data, notice that because the input point cloud is not normalized, the range of the coordinates is about 200, which does not match the assumption in the Wahba problem https://arxiv.org/abs/1905.12536, where we assume the TIMs are typically unit-norm vectors. Moreover, because the TIMs are not normalized, it looks like they all lie on a single straight line.

So there are two ways of fixing the issue:

• One way is that one can manually normalize the TIMs to have unit norm, and then use the GNC solver to solve for the rotation. I also attached a visualization of the TIMs after normalization, and you can see that the TIMs now lie on a unit-sphere. I have tried GNC+DRS certifier on this set of normalized TIMs and it works perfectly. However, currently the TEASER++ implementation does not support manually normalizing the TIMs and calling the GNC-TLS solver independently. Therefore, I tried it using our Matlab prototype code. (@jingnanshi Maybe consider adding a flag in TEASER++ to normalize the TIMs before passing to GNC-TLS?)

• However, there is an easier fix to this issue -- we can directly normalize the input point cloud to fit inside a unit cube, that is, we scale the input point cloud so that the range of the coordinates is precisely 1.0. I have tried this method too, using TEASER++'s python binding, and I have attached the code for your reference. This method also works perfectly and the DRS certifier is able to return a suboptimality gap below 0.1%. (A sample result with 90% outliers, and 200 TIMs is shown below.)

Give the code a try and let us know if you have further questions. BTW, this is a very good issue report -- we did not realize this potential problem before. Thanks for sharing!

example_lidar_scan_cert.zip

TEASER Registration@0.9 outlier ratio: R_err=2.8920013776999873[deg], t_err=0.002105632500156549.
Number of TIMs: 200.
Starting linear inverse map calculation.
Linear proj at i=0
Linear proj at i=10
Linear proj at i=20
Linear proj at i=30
Linear proj at i=40
Linear proj at i=50
Linear proj at i=60
Linear proj at i=70
Linear proj at i=80
Linear proj at i=90
Linear proj at i=100
Linear proj at i=110
Linear proj at i=120
Linear proj at i=130
Linear proj at i=140
Linear proj at i=150
Linear proj at i=160
Linear proj at i=170
Linear proj at i=180
Linear proj at i=190
Finalizing A_inv ...
A_inv finalized.
Obtained linear inverse map.
Linear projection time: 0.303704
Obtained Q_cost matrix.
Obtained D_omega matrix.
Starting Douglas-Rachford Splitting.
Iteration: 0
PSD time: 0.367166
Dual Projection time: 0.0429138
Sub Optimality Gap time: 0.350968
Current sub-optimality gap: 2.31014
PSD time: 0.322108
Dual Projection time: 0.0398719
Sub Optimality Gap time: 0.340229
Current sub-optimality gap: 1.71097
PSD time: 0.314531
Dual Projection time: 0.0392528
Sub Optimality Gap time: 0.334258
Current sub-optimality gap: 1.53477
PSD time: 0.306462
Dual Projection time: 0.0348892
Sub Optimality Gap time: 0.332215
Current sub-optimality gap: 0.936499
PSD time: 0.295047
Dual Projection time: 0.0360362
Sub Optimality Gap time: 0.33864
Current sub-optimality gap: 0.769344
PSD time: 0.310949
Dual Projection time: 0.0369563
Sub Optimality Gap time: 0.344102
Current sub-optimality gap: 0.478863
PSD time: 0.302425
Dual Projection time: 0.0380352
Sub Optimality Gap time: 0.345331
Current sub-optimality gap: 0.0973418
PSD time: 0.325942
Dual Projection time: 0.0411443
Sub Optimality Gap time: 0.356996
Current sub-optimality gap: 0.0533233
PSD time: 0.335501
Dual Projection time: 0.03823
Sub Optimality Gap time: 0.33796
Current sub-optimality gap: 0.0364319
PSD time: 0.317635
Dual Projection time: 0.0388267
Sub Optimality Gap time: 0.334922
Current sub-optimality gap: 0.032993
Iteration: 10
PSD time: 0.316812
Dual Projection time: 0.0359963
Sub Optimality Gap time: 0.348936
Current sub-optimality gap: 0.0320276
PSD time: 0.312475
Dual Projection time: 0.0365663
Sub Optimality Gap time: 0.344132
Current sub-optimality gap: 0.0280966
PSD time: 0.307729
Dual Projection time: 0.0396079
Sub Optimality Gap time: 0.347025
Current sub-optimality gap: 0.0196238
PSD time: 0.321662
Dual Projection time: 0.0389538
Sub Optimality Gap time: 0.350755
Current sub-optimality gap: 0.00951329
PSD time: 0.317571
Dual Projection time: 0.03746
Sub Optimality Gap time: 0.35406
Current sub-optimality gap: 0.00570769
PSD time: 0.331516
Dual Projection time: 0.0355681
Sub Optimality Gap time: 0.354861
Current sub-optimality gap: 0.00572525
PSD time: 0.320576
Dual Projection time: 0.0360841
Sub Optimality Gap time: 0.360485
Current sub-optimality gap: 0.0068241
PSD time: 0.31892
Dual Projection time: 0.0353581
Sub Optimality Gap time: 0.348707
Current sub-optimality gap: 0.00953713
PSD time: 0.313584
Dual Projection time: 0.0363902
Sub Optimality Gap time: 0.3497
Current sub-optimality gap: 0.00192235
PSD time: 0.317518
Dual Projection time: 0.0407545
Sub Optimality Gap time: 0.358829
Current sub-optimality gap: 1.89059e-15
Suboptimality condition reached in 20 iterations. Stopping DRS.
TEASER Certifier: best subopt: 1.8905874965858548e-13%, iterations: 20.

[banbanhank]

@hankyang94, thank you very much. I have a question. When I'm trying to do registration on lidar scans from two different scanning time, the certification process is confusing. I realized that I have to normalize the point cloud before calculating the η value. However, the upper and lower bound of two scans are different. I wonder which value should I use on normalizing the input of teaser++, or I should just normalize them separately.

[hankyang94]

I think a good idea is to use the maximum of the two upper bounds. Basically you want to apply the same scaling/normalization to both scans so that both of them fit inside a unit cube.