Applying colinearity function to dem with no flat plane

Happy new year,
Trying to obtain source undistorted images x,y coordinates relative to dem tif file,and im pretty close thanks of course to @pierotofy and his orthorectify script.I realized that colinearinity function apllied by piero used a flat plane at dem min vallue and i have big errors on different elevation pixels.
Pierotify suggests in his script that is possible to use dtm valid elevation by ‘‘scanning the whole dem’’ ,but i cant find the way to apply it.The sample code:

shot_image = udata.load_undistorted_image(shot.id)

r = shot.pose.get_rotation_matrix()
print(r)
Xs, Ys, Zs = shot.pose.get_origin()
print(Xs, Ys, Zs)
a1 = r[0][0]
b1 = r[0][1]
c1 = r[0][2]
a2 = r[1][0]
b2 = r[1][1]
c2 = r[1][2]
a3 = r[2][0]
b3 = r[2][1]
c3 = r[2][2]

print(“Camera pose: (%f, %f, %f)” % (Xs, Ys, Zs))

img_h, img_w, num_bands = shot_image.shape
print(“Image dimensions: %sx%s pixels” % (img_w, img_h))
f = shot.camera.focal * max(img_h, img_w)
print(‘num of bands {}’.format(num_bands))
has_nodata = dem_raster.profile.get(‘nodata’) is not None
#! Compute bounding box of image coverage
#! assuming a flat plane at Z = min Z
#! (Otherwise we have to scan the entire DEM)
#! The Xa,Ya equations are just derived from the colinearity equations
#! solving for Xa and Ya instead of x,y
def dem_coordinates(cpx, cpy):
“”"
:param cpx principal point X (image coordinates)
:param cpy principal point Y (image coordinates)
“”"
Za = dem_min_value
m = (a3b1cpy - a1b3cpy - (a3b2 - a2b3)cpx - (a2b1 - a1b2)f)
Xa = dem_offset_x + (m
Xs + (b3
c1cpy - b1c3cpy - (b3c2 - b2c3)cpx - (b2c1 - b1c2)f)Za - (b3c1cpy - b1c3cpy - (b3c2 - b2c3)cpx - (b2c1 - b1c2)f)Zs)/m
Ya = dem_offset_y + (m
Ys - (a3
c1
cpy - a1c3cpy - (a3c2 - a2c3)cpx - (a2c1 - a1c2)f)Za + (a3c1cpy - a1c3cpy - (a3c2 - a2c3)cpx - (a2c1 - a1c2)*f)*Zs)/m

    y, x = dem_raster.index(Xa, Ya)
    return (x, y)

dem_ul = dem_coordinates(-(img_w - 1) / 2.0, -(img_h - 1) / 2.0)
dem_ur = dem_coordinates((img_w - 1) / 2.0, -(img_h - 1) / 2.0)
dem_lr = dem_coordinates((img_w - 1) / 2.0, (img_h - 1) / 2.0)
dem_ll = dem_coordinates(-(img_w - 1) / 2.0, (img_h - 1) / 2.0)

The code is part of orthorectify contrib ODM/orthorectify.py at master · OpenDroneMap/ODM · GitHub

Best regards

1 Like

Thx to this paper https://www.asprs.org/wp-content/uploads/pers/2005journal/jul/2005_jul_863-871.pdf IP method can be used to define correct Za (except extreme cases).Still working on it having issues on corner of images that are outside the dem.Any advise would help.

This topic was automatically closed 30 days after the last reply. New replies are no longer allowed.