Display units for translations in registration.txt

[meyer-usda]

Would it be possible to display the units for the translations in registration.txt? We are trying to apply a transformation found for a pair of DSMs onto an orthomosaic derived from the AOI DSM (and pointcloud). It seems like the units could be pixels or perhaps cm/mm but so far haven't found a conversion that fits the observed shift.

When viewing DSM and DSM_registered there appears to be a 10m shift in x and ~0m in Y. Understanding how this shift maps from the provided translations would be great, thanks!

Transformation matrix:
[[ 9.99599931e-01 -1.46808327e-04 -5.44040176e-03 1.14434870e+03]
[ 1.51526200e-04 9.99614360e-01 8.66455776e-04 2.57916565e+03]
[ 5.44027240e-03 -8.67267615e-04 9.99599567e-01 4.11771693e+03]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
Transformation Parameters:
Omega = -0.050 degrees
Phi = -0.312 degrees
Kappa = 0.009 degrees
X Translation = 1144.349
Y Translation = 2579.166
Z Translation = 4117.717
Scale = 0.999615
Number of pairs = 339092
RMSEs:
X = +/-0.535,
Y = +/-0.576,
Z = +/-0.576,
3D = +/-0.975
Propagated Error:
Assumed global 3D error = +/-3m
3D_RMSE = +/-0.975
Total Error is computed as √(global_3d_error²+3D_RMSE²)
Total Error = +/-3.154384032226434

[hobu]

The units should be whatever the units of the coordinate system of the data are provided with. CODEM shouldn't do any unit translation. Could the DSM vs DTM differences be some kind of aliasing effect?

[meyer-usda]

The CRS is WGS84/UTM so should be in meters.

X Translation = 1144.349
Y Translation = 2579.166
Z Translation = 4117.717

This is an example transect from our dataset (a forested area). You can see that the alignment is only off by a few meters at most. This seems inconsistent with the large (~1km) translation values provided. Do you have any insight on what I might be missing?

(Note: that CODEM produces a registered cloud/dsm that is well aligned, so those values work under the hood. I'm trying to understand how to then apply that shift to a related orthomosaic).

Thanks!

[hobu]

Without some data and example commands that demonstrate the process you used to produce the output that isn't behaving, anything we might have to comment about it is just wild ass guesses.

[meyer-usda]

Makes sense and thanks for the prompt reply. Here is more information (I can also send the actual files directly if needed).

I ran codem using the default settings
Here is the log and registration .txt from the run including the settings used:
logCopy.txt
registration.txt

Here is the foundation DSM

Here is the unregistered DSM overlaid (green box highlights an example of the mismatch measured as ~-11m in x and ~2m in y)

Here is the registered DSM after CODEM runs showing successful alignment.

The AOI DSM came from a photogrammetry workflow. I'd like to apply the same transformation to the orthomosaic produced with that DSM (This ignores projection/distortion issues, but we want to test this result first). What I'm not sure I understand is: How do the registration.txt values lead to this shift? It seems like those values are much larger than the observed shift.

Let me know if this information helps.

[sergeikharchenko]

Hello, Meyer!
I am trying to solve the same issue. Here is the code that basically applies the copied matrix from registration.txt to a 3-channel photoplan raster. As far as I can understand, there are no gross errors. But also in my case the photoplane lays down a bit differently than the DEM after codem.

import rasterio
import numpy as np
from rasterio.transform import Affine

# 0. Source orthophotoplan (check the path)
with rasterio.open(r'E:\Data\Apsny\Statsionar\DEM_coregistration\photoplan2025_hq.tif') as src:
    data = src.read() #1
    transform = src.transform
    crs = src.crs
    profile = src.profile

# 1. Copypaste here your matrix, check the commas
transform_matrix = np.array([
    [ 9.98519026e-01, -4.45550829e-03, -4.04216143e-03,  2.22110198e+04], 
    [ 4.44359565e-03,  9.98522912e-01, -2.94701400e-03,  4.32315233e+03], 
    [ 4.05525349e-03,  2.92897247e-03,  9.98524618e-01, -1.63850218e+04], 
    [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.00000000e+00]
])

# 2. Trasformation for georeference (X, Y)
a, b, c, d, e, f = transform[0], transform[1], transform[2], transform[3], transform[4], transform[5]

# New transformation matrix
new_transform = Affine(
    transform_matrix[0, 0] * a + transform_matrix[0, 1] * d,
    transform_matrix[0, 0] * b + transform_matrix[0, 1] * e,
    transform_matrix[0, 0] * c + transform_matrix[0, 1] * f + transform_matrix[0, 3],
    transform_matrix[1, 0] * a + transform_matrix[1, 1] * d,
    transform_matrix[1, 0] * b + transform_matrix[1, 1] * e,
    transform_matrix[1, 0] * c + transform_matrix[1, 1] * f + transform_matrix[1, 3]
)

# 3. Update metadata
profile.update({
    'transform': new_transform,
    'height': data.shape[1],
    'width': data.shape[2],
    'dtype': data.dtype
})

# 4. Saving georeferences orthophotoplan
with rasterio.open('transformed_opp.tif', 'w', **profile) as dst:
    dst.write(data,)

Difference between corrected DEM (using CODEM) and non-corrected orthophotoplan = 3 m. DEM is semi-transparent over the photoplan.

Difference between corrected DEM (using CODEM) and corrected orthophotoplan (using the script above and trasformation matrix from CODEM) = 0.5 m

I am not an expert in linear algebra etc., but probably the transformation matrix obtained from three-dimensional data (DEM), especially its Z component, distorts X and Y. Probably it is somehow possible to collapse a 4*4 matrix with height into a 3*3 matrix without height by not simply extracting the necessary cells.

[chambbj]

This can be misleading and a source for confusion for sure. But remember that the translation and the rotation are both components of the final composed transformation, and they not should be considered independently. CODEM does not currently have a translation only registration option.

[meyer-usda]

Makes sense. Are there plans to add a translation only option?

[sergeikharchenko]

Hello, Meyer. I've solved the problem (accouting chambbj's commentary), but on R. It'll take some time to rewrite it on Python.
I gonna upload R-script as soon as I'll put it in order (mostly, nessesary comments to it).

The general principle - add to your orthomosaic (3 or 4 bands) another band: elevation from DEM/DSM, then apply matrices production to XYZ of each pixel and save related RGB of it. Then you'll get irregular data (because of transformation), you need project it on GRID again.

Look at the fallen trees for example.

[meyer-usda]

Nice! Good work.

[sergeikharchenko]

If you can use R - try it:
https://github.com/sergeikharchenko/AATMCO

[disarticulate]

it'd be helpful if the registration.txt produced was in yaml, along with the noted units, for easy to ingest.

But a feature request state above would be to use the registration on other files. an api like:

run registration forward

codem -forward -registration=registration.txt unregistered.(laz|tif...)

run registration reverse

codem -reverse -registration=registration.txt registered.(laz|tif...)

This would generally match similar tools for converting between coordinate systems. I can open up a new issue and possibly get a initial pull request going. The registration.txt should start by being parseable.