How to find the best translation available to match two sets of related points that do not exactly match?

Question

I'm working with lidar data for forestry analysis.

I'm trying to match some tree measurements that were taken on the ground. Circle areas were defined, and the center coordinates of those circles were measured using a GPS that was only accurate to 2-3 meters and gives a systematic error to all trees' position.

From the center of the circle, every tree's position in a radius of 15m was recorded, with azimuth and distance.

The lidR package in R automatically detects the top of the trees in an area, and returns their XYZ position. Some trees that are below the canopy are not detected, and some points that are detected as tree tops are not actual trees.

Now, I'm trying to match the trees XY_obs from the observed area on the ground with the detected trees XY_lid from the lidar point cloud, accounting for the systematic error caused by the poor accuracy of the GPS that measured the center coordinates XY_plac of the observed area, and the random error caused by tree tops not being exactly at the same XY coordinates as their base.

The question is : does an algorithm that finds the best 2D translation to match two related point sets exist ? How would you proceed ?

Here's a reproducible script with a brute force algorithm that I implemented :

library(sf)

XY_lid <- structure(c(948163.83, 948164.21, 948177.15, 948164.68, 948177.23, 
                      948170.3, 948161.59, 948175.5, 948170.57, 948170.26, 948168.4, 
                      948157.87, 948164.95, 948166.13, 948179.01, 948159.62, 948156.52, 
                      948171.06, 948181.58, 948169.07, 948153.4, 948158.94, 948163.87, 
                      948155.16, 948175.63, 948166.31, 6635593.7, 6635603.3, 6635613.56, 
                      6635622.33, 6635603.46, 6635604.54, 6635619.97, 6635596.58, 6635621.05, 
                      6635609.86, 6635601.03, 6635596.77, 6635605.17, 6635609.83, 6635599.31, 
                      6635616.89, 6635615.83, 6635593.86, 6635607.01, 6635613.42, 6635614.45, 
                      6635601.49, 6635597.17, 6635601.63, 6635605.68, 6635619.43), dim = c(26L, 
                                                                                           2L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", 
                                                                                                                  "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", 
                                                                                                                  "20", "21", "22", "23", "24", "25", "26"), c("X", "Y")))
XY_obs <- structure(c(948165.392030324, 948175.061077685, 948174.673817546, 
                      948170.855442077, 948174.342852745, 948173.474049551, 948174.777193994, 
                      948174.213948658, 948182.594583396, 948168.682365432, 948169.914549858, 
                      948166.336999996, 948166.338676136, 948163.678106303, 948163.090351684, 
                      948163.024286944, 948161.365684326, 948158.606039234, 948159.763234718, 
                      948161.599215024, 948161.07742998, 948170.002990815, 6635611.40733171, 
                      6635615.65017971, 6635612.4194865, 6635608.44601593, 6635608.84904698, 
                      6635608.5873978, 6635608.6633901, 6635605.56821887, 6635586.96396967, 
                      6635603.82751194, 6635600.86259593, 6635587.01052903, 6635600.71554885, 
                      6635599.22824923, 6635602.07525215, 6635604.68862755, 6635603.91066638, 
                      6635603.7185643, 6635606.30449629, 6635608.14719949, 6635608.66388758, 
                      6635602.05110294), dim = c(22L, 2L), dimnames = list(c("1", "2", 
                                                                             "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", 
                                                                             "15", "16", "17", "18", "19", "20", "21", "22"), c("X", "Y")))

performance <- function(a, b) {
  
  N = dim(a)[1]
  M = dim(b)[1]
  
  perf = vector()
  dist = matrix(nrow = N, ncol = M) #matrice des distances
  for (i in 1:N) {
    for (j in 1:M) {
      dist[i,j] = sqrt((a[i,1]-b[j,1])^2 + (a[i,2]-b[j,2])^2)
    }
  }
  
  best_dist = numeric(length = min(N,M))
  k = 1
  while(length(dim(dist)) == 2) {
    ind = which.min(dist)
    best_dist[k] = dist[ind]
    # Delete points that matched
    ind = arrayInd(ind, dim(dist))
    # Delete corresponding row and column. Different case when there's just 2 rows or column left
    dist = dist[-ind[1],]
    if (!is.null(dim(dist))) {
      dist = dist[,-ind[2]]
    }
    else {
      dist = dist[-ind[2]]
    }
    k = k+1
  }
  
  return(sum(best_dist))
}

offsets <- matrix(ncol=3)
colnames(offsets) <- c("X", "Y", "Perf")

# Brute force algorithm

XY_obs_temp = XY_obs
for (x in seq(-10, 10, by = 0.2)) {
  for (y in seq(-10, 10, by = 0.2)) {
    # Set offset to observed points
    XY_obs_temp[,"X"] = XY_obs[,"X"]+x
    XY_obs_temp[,"Y"] = XY_obs[,"Y"]+y
    
    # Measure sum of distances to closest points
    perf = performance(XY_obs_temp, XY_lid)

    # Save performance for this xy offset
    offsets <- rbind(offsets, c(x,y,perf))
  }
}
  
  # Find the indice of xy with the best performance
  ind = which.min(offsets[,"Perf"])
  
  # Return best offset with performance
  offsets[ind,]

Thank you !