Transformation Matrix (4 x 4) in lidR - Issue during export

[Sylva3D]

I'm trying to use a transformation matrix like the one implemented in Cloud Compare.
When I'm applying the matrix, the X,Y and Z coordinates look fine. Yet when I'm exporting the laz point cloud using the WriteLas function, I observe a significant offset in the exported point cloud (i.e., the coordinates are different than the one I had in R - see the two pictures in attachment).

I first thought that it was a shift caused by the coordinate system but this is not the source of the issue.
Maybe it's related to the large number od decimal


value in the transformation matrix?
Potentially a header issue?

Any idea how to resolve this issue? Thanks a lot!

Here is the code with a small subset of the las file available on Gdrive:
https://drive.google.com/drive/folders/1gtf8jQjw82jWm7MjnDD97KNbL9KbH9fP?usp=sharing

library(lidR)
las <- readTLSLAS("input")

# Function to apply a 4x4 transformation matrix to the point cloud
apply_transformation <- function(las, transformation_matrix) {
  # Extract point coordinates (X, Y, Z) and add a column of 1s for affine transformation
  coords <- cbind(las@data$X, las@data$Y, las@data$Z, rep(1, nrow(las@data)))
  
  # Apply the transformation matrix (matrix multiplication)
  transformed_coords <- coords %*% t(transformation_matrix)
  
  # Update the LAS object with transformed coordinates
  las@data$X <- transformed_coords[, 1]
  las@data$Y <- transformed_coords[, 2]
  las@data$Z <- transformed_coords[, 3]
  
  return(las)
}

# Define the first transformation matrix 
transformation_matrix_1 <- matrix(c(
  0.967502420747, -0.252091995625, -0.019715262882, 305599.051770884835,
  0.252486279972, 0.967372280401, 0.021013080159, 5093252.171847742982,
  0.013774769502, -0.025308039304, 0.999584793238, 203.270126236990,
  0.000000000000, 0.000000000000, 0.000000000000, 1.000000000000
), nrow = 4, byrow = TRUE)

# Apply the first transformation matrix to the point cloud
las <- apply_transformation(las, transformation_matrix_1)

# Optionally, set the projection to EPSG:2959 (NAD83(CSRS) / MTM zone 9)
#projection(las) <- "EPSG:2959"

# Save the transformed LAS point cloud
writeLAS(las, "output")

[Jean-Romain]

Try not to use las@data. The operator $ is here to protect the object from bad operation. Use las$X. Using the $ operator you would have had an error:

Erreur : Trying to store values ranging in [-682709.88887492, -682706.253151878] but storable range is [-214748.3647, 214748.3647]

Before to modify the point cloud, modify the header such as your LAS object can store your coordinates.

# Function to apply a 4x4 transformation matrix to the point cloud
apply_transformation <- function(las, transformation_matrix) 
{
  # Extract point coordinates (X, Y, Z) and add a column of 1s for affine transformation
  coords <- cbind(las@data$X, las@data$Y, las@data$Z, rep(1, nrow(las@data)))
  
  # Apply the transformation matrix (matrix multiplication)
  transformed_coords <- coords %*% t(transformation_matrix)
  
  # Update the LAS object with transformed coordinates
  Xt <- transformed_coords[, 1]
  Yt <- transformed_coords[, 2]
  Zt <- transformed_coords[, 3]
  
  xoffset = floor(min(Xt))
  yoffset = floor(min(Yt))
  zoffset = floor(min(Zt))
  
  las@header$`X offset` = xoffset
  las@header$`Y offset` = yoffset
  las@header$`Z offset` = zoffset
  
  # Update the LAS object with transformed coordinates
  las$X <- transformed_coords[, 1]
  las$Y <- transformed_coords[, 2]
  las$Z <- transformed_coords[, 3]
  
  return(las)
}