Non-constant transformation of photogrammetric point clouds in R?

Question

I have a process which attempts to register two forested point clouds using CloudCompare's Iterative Closest Point (ICP) algorithm. Often times, warping (especially in the vertical direction) occurs in the generation of UAV derived photogrammetric clouds.

Using a moving window, my process attempts to match the sample clouds found in each moving window (top of figure). ICP returns a matrix of transformation values, of which the estimated X,Y and Z transformation values are then stored and associated with the center coords of each window (bottom of figure - the coloured points).

With the resulting point dataset from the moving window ICP process, we can fit models to predict the X, Y and Z shift values as a function of 2D space. I have experimented with low order polynomial functions for example:

Z shift = b0(UTMx) + b1(UTMx)^2 + b2(UTMx)^3 + b3(UTMy) + b4(UTMy)^2 + intercept

where b0-b4 are model coefficients and UTMx/y refers to the 2D coordinates. 2 other models are fit for X and Y shift. So far the models prove to be accurate (rsquared ~0.9) and applying the shifts to the sample subsections of the original cloud show very close alignment visually.

With these models we can go back to the original Photogrammetric Cloud and apply X,Y and Z shifts to correct for warping (at least so it matches the LiDAR cloud which we know has better georeferencing procedure than UAV). These corrections ideally would be applied to each point of the original Photogrammetric Cloud... I am looking for an efficient way to apply the shifts (which vary according to the models) in a computationally efficient way because in some cases we deal with up to a billion points. My first instinct is to using some sort of tiling procedure, running these in parallel... I prefer using R and the lidR package however I'm open to other suggestions/comments/questions regarding any part of this process!

Answer 1

I'm answering for a lidR solution based on what I understand.

If I understand well, you have 3 rasters or at least 3 lattices of points. The first one gives you the X translation to apply on a given window. The second the Y translation and so on.

Also you have 3 models one for each coordinates shift that give you, for a given XY the shift to apply to X, Y or Z. This allows you to apply a continuous transformation instead of a discretized transformation.

I would start by a function that apply the shift:

lasshift = function(las, ...)
{
  las$Z <- las$Z + b0(las$X) + b1(las$X)^2 + b2(las$X)^3 + b3(las$Y) + b4(las$Y)^2
  # the same for X and Y
  return(las)
}

This function will apply the computed shift to each point. Because you did not give more information about your model I can't write more specific code.

Now, for a "chunked" version of this function you should look at this vignette. Basically, it looks like the code draft below (which won't work as it is), but you must first read the vignette and the doc of catalog_apply to finish this draft. Also, for speed I recommend you to read this vignette.

lasshift = function(las, ...)
{
  UseMethod("lasshift", las)
}

lasshift.LAS = function(las, ...)
{
  # The function above
}

lasshift.LAScluster = function(las, ...)
{
  las <- readLAS(las)                        
  if (is.empty(las)) return(NULL)           

  las <- lasshift(las, ...)
  return(las)
}

lasshift.LAScatalog = function(las, ...)
{
   # Force some options
   opt_select(las) <-  "*" 
   opt_chunck_buffer(las) <- 0

   options <- list(need_output_file = TRUE)

   output <- catalog_apply(las, lasshift, ..., .options = options)
   output <- unlist(output)
   output <- catalog(output)
   return(output)
}

ctg = catlog("path")
opt_chunk_size(ctg) <- 100
opt_output_file(ctg) <- "/templated/path/"
new_ctg = lasshift(ctg, xshift, yshift, zshift)