professorfabioandrade / georef

6 stars 1 forks source link

Georeferencing using DJI Drones #4

Open peakyquest opened 6 months ago

peakyquest commented 6 months ago

I am performing the same procedure on both the DJI Mavic 2 Pro and the Mavic Mini, but I am unable to achieve consistent positional accuracy of the object. I have calculated the K matrix with the help of Agisoft Lens.

import cv2, math
import numpy as np
import utm
class GeoReference:
    def __init__(self, image_path,
                 gimbal_roll, gimbal_pitch, gimbal_yaw,
                 uas_roll, uas_pitch, uas_yaw,
                 latitude, longitude, altitude):
        # Input the image path for geo-reference
        self.image_path = image_path
        # values of gimbal's roll pitch and yaw with w.r.t uas
        self.gimbal_roll = math.radians(gimbal_roll)
        self.gimbal_pitch = math.radians(gimbal_pitch)
        self.gimbal_yaw = math.radians(gimbal_yaw)
        # values of the uas with w.r.t ned frame
        self.uas_roll = math.radians(uas_roll)
        self.uas_pitch = math.radians(uas_pitch)
        self.uas_yaw = math.radians(uas_yaw)
        # latitude longitude and altitude values of the uas
        self.latitude = latitude
        self.longitude = longitude
        self.UTM_z = altitude
        # Intrinsic parameters camera matrix, Calculated by the resizing or the Mavic 2 Pro Images
        self.K = np.array([
            [1664.59, 0, 935.755],
            [0, 1326.18, 487.819],
            [0, 0, 1]
        ])
        # Rotation and Translation Matrix from Camera frame to Gimbal frame
        self.R_C_to_G = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
        self.T_C_to_G = np.array([[0], [0], [0]])
        # Rotation and Translation Matrix from Gimbal frame to UAS frame
        self.R_G_to_UAS = math.nan
        self.T_G_to_UAS = np.array([[-0.002], [0.023], [0.002]])
        # Rotation and Translation Matrix from UAS frame to NED frame
        self.R_UAS_to_NED = math.nan
        self.T_UAS_to_NED = np.array([[0], [0], [0]])
        # Rotation and Translation Matrix from NED frame to ENU frame
        self.R_NED_to_ENU = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
        self.T_NED_to_ENU = math.nan

    def open_image_and_capture_coordinates(self):
        # Read the image
        image = cv2.imread(self.image_path)
        # Check if the image was successfully loaded
        if image is not None:
            # Display the image
            cv2.imshow('Image', image)
            # Set the callback function for mouse events
            cv2.setMouseCallback('Image', self.on_mouse_click)
            # Wait for a key event
            cv2.waitKey(0)
            cv2.destroyAllWindows()
        else:
            print(f"Error: Unable to load the image at '{self.image_path}'. Make sure the file path is correct.")

    def on_mouse_click(self, event, x, y, flags, param):
        if event == cv2.EVENT_LBUTTONDOWN:
            # Get the pixel values i-e x and y
            print(f"Clicked at pixel coordinates: ({x}, {y})")
            P_C_prime = self.image_to_camera_frame(self.K, x, y)
            # Camera frame to gimbal frame
            P_G_prime = self.camera_to_gimbal_frame(P_C_prime, self.R_C_to_G, self.T_C_to_G)
            # Gimbal frame to UAS frame
            self.R_G_to_UAS = self.rotation_matrix_gimbal_to_uas(self.gimbal_yaw, self.gimbal_pitch, self.gimbal_roll)
            P_UAS_prime = self.gimbal_to_uas_frame(P_G_prime, self.R_G_to_UAS, self.T_G_to_UAS)
            # UAS frame to NED frame
            self.R_UAS_to_NED = self.rotation_matrix_uas_to_ned(self.uas_yaw, self.uas_pitch, self.uas_roll)
            P_NED_prime = self.uas_to_ned_frame(P_UAS_prime, self.R_UAS_to_NED, self.T_UAS_to_NED)
            # NED frame to ENU frame
            UTM_x, UTM_y, letter, zone = self.latlon_to_utm(self.latitude, self.longitude)
            self.T_NED_to_ENU = np.array([[UTM_x], [UTM_y], [self.UTM_z]])
            P_ENU_prime = self.ned_to_enu_frame(P_NED_prime, self.R_NED_to_ENU, self.T_NED_to_ENU)
            z_ENU = 0
            T = self.T_NED_to_ENU + self.R_NED_to_ENU @ self.T_UAS_to_NED + self.R_NED_to_ENU @ self.R_UAS_to_NED @ self.T_G_to_UAS + self.R_NED_to_ENU @ self.R_UAS_to_NED @ self.R_G_to_UAS @ self.T_C_to_G
            z_T = T[2]
            z_ENU_prime = P_ENU_prime[2]
            z_C = self.calculate_z_C(z_ENU, z_T, z_ENU_prime)

            # Calculating P_ENU
            P_ENU = self.calculate_P_ENU(z_C, P_ENU_prime, T)
            P_ENU = list(map(str, P_ENU))
            # Use list comprehension to apply the conversion function to each element
            float_list = [self.convert_to_float(s) for s in P_ENU]
            target_latitude, target_longitude = self.utm_to_latlon(float_list[0], float_list[1], letter, zone)
            print(target_latitude, target_longitude)

    def convert_to_float(self, s):
        # Strip the square brackets and convert to float
        return float(s.strip('[]'))

    def image_to_camera_frame(self, K, u, v):
        i = np.array([[u], [v], [1]])
        P_C_prime = np.linalg.inv(K) @ i
        return P_C_prime

    def camera_to_gimbal_frame(self, P_C_prime, R_C_to_G, T_C_to_G):
        P_G_prime = R_C_to_G @ P_C_prime + T_C_to_G
        return P_G_prime

    def gimbal_to_uas_frame(self, P_G_prime, R_G_to_UAS, T_G_to_UAS):
        P_UAS_prime = R_G_to_UAS @ P_G_prime + T_G_to_UAS
        return P_UAS_prime

    def uas_to_ned_frame(self, P_UAS_prime, R_UAS_to_NED, T_UAS_to_NED):
        P_NED_prime = R_UAS_to_NED @ P_UAS_prime + T_UAS_to_NED
        return P_NED_prime

    def ned_to_enu_frame(self, P_NED_prime, R_NED_to_ENU, T_NED_to_ENU):
        P_ENU_prime = R_NED_to_ENU @ P_NED_prime + T_NED_to_ENU
        return P_ENU_prime

    def calculate_z_C(self, z_ENU, z_T, z_ENU_prime):
        return (z_ENU - z_T) / (z_ENU_prime - z_T)

    def calculate_P_ENU(self, z_C, P_ENU_prime, T):
        P_ENU = z_C * P_ENU_prime - z_C * T + T
        return P_ENU.round(4)

    def rotation_matrix_gimbal_to_uas(self, yaw, pitch, roll):
        R_G_to_UAS = np.array([
            [math.cos(yaw) * math.cos(pitch),
             -math.sin(yaw) * math.cos(roll) + math.cos(yaw) * math.sin(pitch) * math.sin(roll),
             math.sin(yaw) * math.sin(roll) + math.cos(yaw) * math.cos(roll) * math.sin(pitch)],
            [math.sin(yaw) * math.cos(pitch),
             math.cos(yaw) * math.cos(roll) + math.sin(roll) * math.sin(pitch) * math.sin(yaw),
             -math.cos(yaw) * math.sin(roll) + math.sin(pitch) * math.sin(yaw) * math.cos(roll)],
            [-math.sin(pitch),
             math.cos(pitch) * math.sin(roll),
             math.cos(pitch) * math.cos(roll)]
        ])
        return R_G_to_UAS

    def rotation_matrix_uas_to_ned(self, yaw, pitch, roll):
        R_UAS_to_NED = np.array([
            [math.cos(yaw) * math.cos(pitch),
             -math.sin(yaw) * math.cos(roll) + math.cos(yaw) * math.sin(pitch) * math.sin(roll),
             math.sin(yaw) * math.sin(roll) + math.cos(yaw) * math.cos(roll) * math.sin(pitch)],
            [math.sin(yaw) * math.cos(pitch),
             math.cos(yaw) * math.cos(roll) + math.sin(roll) * math.sin(pitch) * math.sin(yaw),
             -math.cos(yaw) * math.sin(roll) + math.sin(pitch) * math.sin(yaw) * math.cos(roll)],
            [-math.sin(pitch),
             math.cos(pitch) * math.sin(roll),
             math.cos(pitch) * math.cos(roll)]
        ])
        return R_UAS_to_NED

    def latlon_to_utm(self, latitude, longitude):
        # Convert latitude and longitude to UTM coordinates
        utm_coords = utm.from_latlon(latitude, longitude)
        return utm_coords

    def utm_to_latlon(self, utm_easting, utm_northing, utm_zone, utm_letter):
        # Convert UTM coordinates to latitude and longitude
        # zone number is 43, zone letter is 'S'

        lat, lon = utm.to_latlon(easting=utm_easting,
                                 northing=utm_northing,
                                 zone_number=utm_zone,
                                 zone_letter=utm_letter
                                 )
        return lat, lon

def main():
    geo = GeoReference(image_path='images/robo-flow/DJI_0753.MP4.jpg',
                       longitude=-117.1697468,latitude=33.03167344,altitude=78.7000025,
                       gimbal_yaw= 318.8, gimbal_roll=0, gimbal_pitch=-90,
                       uas_pitch=-7.9, uas_roll= -4.9, uas_yaw= 334.5)
    geo.open_image_and_capture_coordinates()

if __name__ == '__main__':
    main()

I am experiencing a 30 to 40 meter error with some angle when attempting to determine the UTM coordinates of an object by clicking on the image to obtain the (u, v) values. The process involves converting these coordinates into UTM and then into latitude and longitude. This error persists across both DJI Mavic 2 Pro and Mavic Mini drones.

Am I doing something wrong ?

professorfabioandrade commented 6 months ago

Hi,

First thing you should do is to check if the gimbal in the DJI drones already give the attitude (roll, pitch and yaw) with respect to the world frame. Therefore you do not need the rotation from gimbal to UAS neither the UAS attitude.

You can investigate that here: https://www.phantomhelp.com/LogViewer/Upload/

Regards, Fabio

On Thu, May 23, 2024, 9:46 AM Haroon Rasheed @.***> wrote:

I am performing the same procedure on both the DJI Mavic 2 Pro and the Mavic Mini, but I am unable to achieve consistent positional accuracy of the object. I have calculated the K matrix with the help of Agisoft Lens.

import cv2, math import numpy as np import utm class GeoReference: def init(self, image_path, gimbal_roll, gimbal_pitch, gimbal_yaw, uas_roll, uas_pitch, uas_yaw, latitude, longitude, altitude):

Input the image path for geo-reference

    self.image_path = image_path
    # values of gimbal's roll pitch and yaw with w.r.t uas
    self.gimbal_roll = math.radians(gimbal_roll)
    self.gimbal_pitch = math.radians(gimbal_pitch)
    self.gimbal_yaw = math.radians(gimbal_yaw)
    # values of the uas with w.r.t ned frame
    self.uas_roll = math.radians(uas_roll)
    self.uas_pitch = math.radians(uas_pitch)
    self.uas_yaw = math.radians(uas_yaw)
    # latitude longitude and altitude values of the uas
    self.latitude = latitude
    self.longitude = longitude
    self.UTM_z = altitude
    # Intrinsic parameters camera matrix, Calculated by the resizing or the Mavic 2 Pro Images
    self.K = np.array([
        [1664.59, 0, 935.755],
        [0, 1326.18, 487.819],
        [0, 0, 1]
    ])
    # Rotation and Translation Matrix from Camera frame to Gimbal frame
    self.R_C_to_G = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
    self.T_C_to_G = np.array([[0], [0], [0]])
    # Rotation and Translation Matrix from Gimbal frame to UAS frame
    self.R_G_to_UAS = math.nan
    self.T_G_to_UAS = np.array([[-0.002], [0.023], [0.002]])
    # Rotation and Translation Matrix from UAS frame to NED frame
    self.R_UAS_to_NED = math.nan
    self.T_UAS_to_NED = np.array([[0], [0], [0]])
    # Rotation and Translation Matrix from NED frame to ENU frame
    self.R_NED_to_ENU = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
    self.T_NED_to_ENU = math.nan

def open_image_and_capture_coordinates(self):
    # Read the image
    image = cv2.imread(self.image_path)
    # Check if the image was successfully loaded
    if image is not None:
        # Display the image
        cv2.imshow('Image', image)
        # Set the callback function for mouse events
        cv2.setMouseCallback('Image', self.on_mouse_click)
        # Wait for a key event
        cv2.waitKey(0)
        cv2.destroyAllWindows()
    else:
        print(f"Error: Unable to load the image at '{self.image_path}'. Make sure the file path is correct.")

def on_mouse_click(self, event, x, y, flags, param):
    if event == cv2.EVENT_LBUTTONDOWN:
        # Get the pixel values i-e x and y
        print(f"Clicked at pixel coordinates: ({x}, {y})")
        P_C_prime = self.image_to_camera_frame(self.K, x, y)
        # Camera frame to gimbal frame
        P_G_prime = self.camera_to_gimbal_frame(P_C_prime, self.R_C_to_G, self.T_C_to_G)
        # Gimbal frame to UAS frame
        self.R_G_to_UAS = self.rotation_matrix_gimbal_to_uas(self.gimbal_yaw, self.gimbal_pitch, self.gimbal_roll)
        P_UAS_prime = self.gimbal_to_uas_frame(P_G_prime, self.R_G_to_UAS, self.T_G_to_UAS)
        # UAS frame to NED frame
        self.R_UAS_to_NED = self.rotation_matrix_uas_to_ned(self.uas_yaw, self.uas_pitch, self.uas_roll)
        P_NED_prime = self.uas_to_ned_frame(P_UAS_prime, self.R_UAS_to_NED, self.T_UAS_to_NED)
        # NED frame to ENU frame
        UTM_x, UTM_y, letter, zone = self.latlon_to_utm(self.latitude, self.longitude)
        self.T_NED_to_ENU = np.array([[UTM_x], [UTM_y], [self.UTM_z]])
        P_ENU_prime = self.ned_to_enu_frame(P_NED_prime, self.R_NED_to_ENU, self.T_NED_to_ENU)
        z_ENU = 0
        T = self.T_NED_to_ENU + self.R_NED_to_ENU @ self.T_UAS_to_NED + self.R_NED_to_ENU @ self.R_UAS_to_NED @ self.T_G_to_UAS + self.R_NED_to_ENU @ self.R_UAS_to_NED @ self.R_G_to_UAS @ self.T_C_to_G
        z_T = T[2]
        z_ENU_prime = P_ENU_prime[2]
        z_C = self.calculate_z_C(z_ENU, z_T, z_ENU_prime)

        # Calculating P_ENU
        P_ENU = self.calculate_P_ENU(z_C, P_ENU_prime, T)
        P_ENU = list(map(str, P_ENU))
        # Use list comprehension to apply the conversion function to each element
        float_list = [self.convert_to_float(s) for s in P_ENU]
        target_latitude, target_longitude = self.utm_to_latlon(float_list[0], float_list[1], letter, zone)
        print(target_latitude, target_longitude)

def convert_to_float(self, s):
    # Strip the square brackets and convert to float
    return float(s.strip('[]'))

def image_to_camera_frame(self, K, u, v):
    i = np.array([[u], [v], [1]])
    P_C_prime = np.linalg.inv(K) @ i
    return P_C_prime

def camera_to_gimbal_frame(self, P_C_prime, R_C_to_G, T_C_to_G):
    P_G_prime = R_C_to_G @ P_C_prime + T_C_to_G
    return P_G_prime

def gimbal_to_uas_frame(self, P_G_prime, R_G_to_UAS, T_G_to_UAS):
    P_UAS_prime = R_G_to_UAS @ P_G_prime + T_G_to_UAS
    return P_UAS_prime

def uas_to_ned_frame(self, P_UAS_prime, R_UAS_to_NED, T_UAS_to_NED):
    P_NED_prime = R_UAS_to_NED @ P_UAS_prime + T_UAS_to_NED
    return P_NED_prime

def ned_to_enu_frame(self, P_NED_prime, R_NED_to_ENU, T_NED_to_ENU):
    P_ENU_prime = R_NED_to_ENU @ P_NED_prime + T_NED_to_ENU
    return P_ENU_prime

def calculate_z_C(self, z_ENU, z_T, z_ENU_prime):
    return (z_ENU - z_T) / (z_ENU_prime - z_T)

def calculate_P_ENU(self, z_C, P_ENU_prime, T):
    P_ENU = z_C * P_ENU_prime - z_C * T + T
    return P_ENU.round(4)

def rotation_matrix_gimbal_to_uas(self, yaw, pitch, roll):
    R_G_to_UAS = np.array([
        [math.cos(yaw) * math.cos(pitch),
         -math.sin(yaw) * math.cos(roll) + math.cos(yaw) * math.sin(pitch) * math.sin(roll),
         math.sin(yaw) * math.sin(roll) + math.cos(yaw) * math.cos(roll) * math.sin(pitch)],
        [math.sin(yaw) * math.cos(pitch),
         math.cos(yaw) * math.cos(roll) + math.sin(roll) * math.sin(pitch) * math.sin(yaw),
         -math.cos(yaw) * math.sin(roll) + math.sin(pitch) * math.sin(yaw) * math.cos(roll)],
        [-math.sin(pitch),
         math.cos(pitch) * math.sin(roll),
         math.cos(pitch) * math.cos(roll)]
    ])
    return R_G_to_UAS

def rotation_matrix_uas_to_ned(self, yaw, pitch, roll):
    R_UAS_to_NED = np.array([
        [math.cos(yaw) * math.cos(pitch),
         -math.sin(yaw) * math.cos(roll) + math.cos(yaw) * math.sin(pitch) * math.sin(roll),
         math.sin(yaw) * math.sin(roll) + math.cos(yaw) * math.cos(roll) * math.sin(pitch)],
        [math.sin(yaw) * math.cos(pitch),
         math.cos(yaw) * math.cos(roll) + math.sin(roll) * math.sin(pitch) * math.sin(yaw),
         -math.cos(yaw) * math.sin(roll) + math.sin(pitch) * math.sin(yaw) * math.cos(roll)],
        [-math.sin(pitch),
         math.cos(pitch) * math.sin(roll),
         math.cos(pitch) * math.cos(roll)]
    ])
    return R_UAS_to_NED

def latlon_to_utm(self, latitude, longitude):
    # Convert latitude and longitude to UTM coordinates
    utm_coords = utm.from_latlon(latitude, longitude)
    return utm_coords

def utm_to_latlon(self, utm_easting, utm_northing, utm_zone, utm_letter):
    # Convert UTM coordinates to latitude and longitude
    # zone number is 43, zone letter is 'S'

    lat, lon = utm.to_latlon(easting=utm_easting,
                             northing=utm_northing,
                             zone_number=utm_zone,
                             zone_letter=utm_letter
                             )
    return lat, lon

def main(): geo = GeoReference(image_path='images/robo-flow/DJI_0753.MP4.jpg', longitude=-117.1697468,latitude=33.03167344,altitude=78.7000025, gimbal_yaw= 318.8, gimbal_roll=0, gimbal_pitch=-90, uas_pitch=-7.9, uas_roll= -4.9, uas_yaw= 334.5) geo.open_image_and_capture_coordinates()

if name == 'main': main()

I am experiencing a 30 to 40 meter error with some angle when attempting to determine the UTM coordinates of an object by clicking on the image to obtain the (u, v) values. The process involves converting these coordinates into UTM and then into latitude and longitude. This error persists across both DJI Mavic 2 Pro and Mavic Mini drones.

Am I doing something wrong ?

— Reply to this email directly, view it on GitHub https://github.com/professorfabioandrade/georef/issues/4, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACEWBLAZB6KVAAHNG7LNPRTZDWNEHAVCNFSM6AAAAABIFCQB4CVHI2DSMVQWIX3LMV43ASLTON2WKOZSGMYTEMJZHEZDQNA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

peakyquest commented 6 months ago

image

I am clicking on the cross as shown in the picture, to get the values of that specific point, and I have gather the information, of the UAS from Air Data logs. The results that I get is:

Clicked at pixel coordinates: (924, 463)
33.03161185118164 -117.16964745153744

When I plot the values I get the following results, and these results are not constant. image

peakyquest commented 6 months ago

the gimbal in the DJI drones give the attitude (roll, pitch and yaw) with respect to the body of the drone (UAS)

peakyquest commented 6 months ago

I was trying to implement georeferencing using this repo, and it works very accurate when the gimbal is looking straight downward, but in the repo there is no calculation how to get the location when gimbal have some angle

https://github.com/roboflow/dji-aerial-georeferencing

professorfabioandrade commented 6 months ago

Could you send me an image and the metadata so I can try here?

On Thu, May 23, 2024, 10:01 AM Haroon Rasheed @.***> wrote:

I was trying to implement georeferencing using this repo, and it works very accurate when the gimbal is looking straight downward, but in the repo there is no calculation how to get the location when gimbal have some angle

https://github.com/roboflow/dji-aerial-georeferencing

— Reply to this email directly, view it on GitHub https://github.com/professorfabioandrade/georef/issues/4#issuecomment-2126476863, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACEWBLD75GNDPHQB6E73TVDZDWO43AVCNFSM6AAAAABIFCQB4CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMRWGQ3TMOBWGM . You are receiving this because you commented.Message ID: @.***>

peakyquest commented 6 months ago

Yes sure, but what I am not sure about few parameters, which I am using in my code, for example "intrinsic parameters" there is no function in DJI SDK to get the value of K matrix, so I am using Agisoft to calculate the intrinsic parameters. (Not sure the values that I am getting are correct or not)

Secondly, I am not sure about the coordinate frame of DJI (ENU or NED). In the documentation of the DJI SDK it is mentioned that they have used two aircraft coordinate system which is GROUND and BODY. The body coordinate system is relative to the aircraft itself. Three perpendicular axes are defined such that the origin is the center of mass, and the X axis is directed through the front of the aircraft and the Y axis through the right of the aircraft. Using the coordinate right hand rule, the Z axis is then through the bottom of the aircraft. The Ground coordinate system, is actually NED frame, and in the documentation the origin for a NED coordinate system is usually a point in the world frame (like take-off location). How to calculate the rotational matrix and translation matrix ?

peakyquest commented 6 months ago

DJI_0107

image

Clicked at pixel coordinates: (996, 631) u and v GPS Position: 33.6627975833333, 73.039458833333 Altitude = 43 Flight Roll, Pitch, Yaw (4.8, -5.8, -122.3) Gimbal Roll Pitch Yaw (0, -90, -130) T_C_to_G = [0, 0, 0] T_G_to_UAS = [1.2, 0, 0] TUAS NED = [0, 0, 0] zENU = 0

peakyquest commented 5 months ago

any updates ?

professorfabioandrade commented 5 months ago

Hi,

I am in a work trip so I didn't have time to check but it is in my to-do list.

Regards, Fabio

On Thu, May 30, 2024, 9:39 AM Haroon Rasheed @.***> wrote:

any updates ?

— Reply to this email directly, view it on GitHub https://github.com/professorfabioandrade/georef/issues/4#issuecomment-2138875605, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACEWBLA5RRQ5TOZDAKZAAFTZE3JRXAVCNFSM6AAAAABIFCQB4CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMZYHA3TKNRQGU . You are receiving this because you commented.Message ID: @.***>

peakyquest commented 5 months ago

Thankyou for looking into this.

peakyquest commented 5 months ago

image

If I assume that the origin of the NED frame is the point where the drone takes off, the translation between the NED frame and the UAS frame will be the x and y components of the distance from the home position to its current position. Am I correct ?

peakyquest commented 5 months ago

in your paper you have mentioned an example with ROS, do you have any repository or ros package related to it ?

professorfabioandrade commented 5 months ago

Hi, I started to process this.

Why is your camera matrix self.K = np.array([

        [1664.59, 0, 935.755],
        [0, 1326.18, 487.819],
        [0, 0, 1]
    ])

instead of self.K = np.array([ [4950.24, 7.95966, 2633.58], [0, 4948.64, 1704.28], [0, 0, 1] ])

?

You mentioned something about resizing. But the image that you sent me is 5472 x 3078. The same as the metadata file says.

Em seg., 3 de jun. de 2024 às 09:52, Haroon Rasheed < @.***> escreveu:

in your paper you have mentioned an example with ROS, do you have any repository or ros package related to it ?

— Reply to this email directly, view it on GitHub https://github.com/professorfabioandrade/georef/issues/4#issuecomment-2144512145, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACEWBLDSEBIE3KUIC3HWPP3ZFQOFJAVCNFSM6AAAAABIFCQB4CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNBUGUYTEMJUGU . You are receiving this because you commented.Message ID: @.***>

--

Fabio Andrade, PhD Associate Professor at University of South-Eastern Norway Research Scientist at NORCE Norwegian Research Center professorfabioandrade.com https://www.professorfabioandrade.com

professorfabioandrade commented 5 months ago

Hi,

I get the expected results here when I only consider the gimbal angles and ignore the flight angles. This is correct as according to the documentation of DJI the gimbal angles are already with respect to the NED frame ( https://developer.dji.com/mobile-sdk/documentation/introduction/flightController_concepts.html). This is also confirmed here: https://forum.dji.com/forum.php?mod=viewthread&tid=305582. Therefore, we do not need to use the flight angles. I also included the undistortion (using ChatGPT) as you provided k1,k2,k3,k4,p1 and p2. Please confirm if these equations from ChatGPT are correct or just remove the undistortion function.

import numpy as np import utm import math

Drone GPS position and altitude

latitude = 33.6627975833333 longitude = 73.039458833333 drone_altitude = 571.47 # drone altitude in meters (from metadata) relative_altitude = 43.5 # difference between drone and ground altitude in meters (from metadata) ground_altitude = drone_altitude - relative_altitude # ground altitude in meters (assuming flat ground - I checked that's the cade in google earth)

Constants from metadata

fx = 4950.24 fy = 4948.64 cx = 2633.58 cy = 1704.28 skew = 7.95966

K = np.array([ [fx, skew, cx], [0, fy, cy], [0, 0, 1] ])

Distortion coefficients

k1 = 0.0232456 k2 = -0.00062289 k3 = -0.0848578 k4 = 0.141931 p1 = -0.00614786 p2 = -0.00273619

Gimbal angles (convert to radians)

gimbal_yaw = np.deg2rad(-130.2) gimbal_pitch = np.deg2rad(-90) gimbal_roll = np.deg2rad(0)

Flight angles (convert to radians)

uas_yaw = np.deg2rad(0) uas_pitch = np.deg2rad(0) uas_roll = np.deg2rad(0)

Prompt user for pixel coordinates

u = float(input("Enter u coordinate: ")) v = float(input("Enter v coordinate: "))

def undistort_point(u, v, k1, k2, k3, k4, p1, p2, fx, fy, cx, cy):

Normalize the coordinates

x_norm = (u - cx) / fx
y_norm = (v - cy) / fy

# Compute the radius squared
r2 = x_norm**2 + y_norm**2
r4 = r2 * r2
r6 = r4 * r2
r8 = r6 * r2

# Compute radial distortion
radial_distortion = 1 + k1 * r2 + k2 * r4 + k3 * r6 + k4 * r8

# Compute tangential distortion
x_tangential = 2 * p1 * x_norm * y_norm + p2 * (r2 + 2 * x_norm**2)
y_tangential = p1 * (r2 + 2 * y_norm**2) + 2 * p2 * x_norm * y_norm

# Apply the distortion to the normalized coordinates
x_undistorted = x_norm * radial_distortion + x_tangential
y_undistorted = y_norm * radial_distortion + y_tangential

# Convert back to pixel coordinates
u_undistorted = fx * x_undistorted + cx
v_undistorted = fy * y_undistorted + cy

return u_undistorted, v_undistorted

def image_to_camera_frame(K, u, v): i = np.array([[u], [v], [1]]) P_C_prime = np.linalg.inv(K) @ i return P_C_prime

def camera_to_gimbal_frame(P_C_prime, R_C_to_G, T_C_to_G): P_G_prime = R_C_to_G @ P_C_prime + T_C_to_G return P_G_prime

def gimbal_to_uas_frame(P_G_prime, R_G_to_UAS, T_G_to_UAS): P_UAS_prime = R_G_to_UAS @ P_G_prime + T_G_to_UAS return P_UAS_prime

def uas_to_ned_frame(P_UAS_prime, R_UAS_to_NED, T_UAS_to_NED): P_NED_prime = R_UAS_to_NED @ P_UAS_prime + T_UAS_to_NED return P_NED_prime

def ned_to_enu_frame(P_NED_prime, R_NED_to_ENU, T_NED_to_ENU): P_ENU_prime = R_NED_to_ENU @ P_NED_prime + T_NED_to_ENU return P_ENU_prime

def calculate_z_C(z_ENU, z_T, z_ENU_prime): return (z_ENU - z_T) / (z_ENU_prime - z_T)

def calculate_P_ENU(z_C, P_ENU_prime, T): P_ENU = z_C P_ENU_prime - z_C T + T return P_ENU.round(4)

def rotation_matrix_gimbal_to_uas(yaw, pitch, roll): R_G_to_UAS = np.array([ [math.cos(yaw) math.cos(pitch), -math.sin(yaw) math.cos(roll) + math.cos(yaw) math.sin(pitch) math.sin(roll), math.sin(yaw) math.sin(roll) + math.cos(yaw) math.cos(roll)

def rotation_matrix_uas_to_ned(yaw, pitch, roll): R_UAS_to_NED = np.array([ [math.cos(yaw) math.cos(pitch), -math.sin(yaw) math.cos(roll) + math.cos(yaw) math.sin(pitch) math.sin(roll), math.sin(yaw) math.sin(roll) + math.cos(yaw) math.cos(roll)

def latlon_to_utm(latitude, longitude):

Convert latitude and longitude to UTM coordinates

utm_coords = utm.from_latlon(latitude, longitude)
return utm_coords

def utm_to_latlon(utm_easting, utm_northing, utm_zone, utm_letter):

Convert UTM coordinates to latitude and longitude

# zone number is 43, zone letter is 'S'

lat, lon = utm.to_latlon(easting=utm_easting,
                            northing=utm_northing,
                            zone_number=utm_zone,
                            zone_letter=utm_letter
                            )
return lat, lon

R_C_to_G = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) R_G_to_UAS = rotation_matrix_gimbal_to_uas(gimbal_yaw, gimbal_pitch, gimbal_roll) R_UAS_to_NED = rotation_matrix_uas_to_ned(uas_yaw, uas_pitch, uas_roll) R_NED_to_ENU = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]]) T_C_to_G = np.array([[0], [0], [0]]) T_G_to_UAS = np.array([[0], [0], [0]]) T_UAS_to_NED = np.array([[0], [0], [0]])

UTM_x, UTM_y, letter, zone = latlon_to_utm(latitude, longitude) UTM_z = drone_altitude T_NED_to_ENU = np.array([[UTM_x], [UTM_y], [UTM_z]])

u_undistorted, v_undistorted = undistort_point(u, v, k1, k2, k3, k4, p1, p2, fx, fy, cx, cy) P_C_prime = image_to_camera_frame(K, u_undistorted, v_undistorted) P_G_prime = camera_to_gimbal_frame(P_C_prime, R_C_to_G, T_C_to_G) P_UAS_prime = gimbal_to_uas_frame(P_G_prime, R_G_to_UAS, T_G_to_UAS) P_NED_prime = uas_to_ned_frame(P_UAS_prime, R_UAS_to_NED, T_UAS_to_NED) P_ENU_prime = ned_to_enu_frame(P_NED_prime, R_NED_to_ENU, T_NED_to_ENU)

z_ENU = ground_altitude T = T_NED_to_ENU + R_NED_to_ENU @ T_UAS_to_NED + R_NED_to_ENU @ R_UAS_to_NED @ T_G_to_UAS + R_NED_to_ENU @ R_UAS_to_NED @ R_G_to_UAS @ T_C_to_G z_T = T[2] z_ENU_prime = P_ENU_prime[2] z_C = calculate_z_C(z_ENU, z_T, z_ENU_prime)

Calculating P_ENU

P_ENU = calculate_P_ENU(z_C, P_ENU_prime, T) lat_ground, lon_ground = utm_to_latlon(P_ENU[0], P_ENU[1], letter, zone)

print("Ground coordinates (lat, lon, alt):", lat_ground[0], lon_ground[0], ground_altitude)

Em sex., 7 de jun. de 2024 às 13:54, Fabio Andrade @.***> escreveu:

Hi, I started to process this.

Why is your camera matrix self.K = np.array([

        [1664.59, 0, 935.755],
        [0, 1326.18, 487.819],
        [0, 0, 1]
    ])

instead of self.K = np.array([ [4950.24, 7.95966, 2633.58], [0, 4948.64, 1704.28], [0, 0, 1] ])

?

You mentioned something about resizing. But the image that you sent me is 5472 x 3078. The same as the metadata file says.

Em seg., 3 de jun. de 2024 às 09:52, Haroon Rasheed < @.***> escreveu:

in your paper you have mentioned an example with ROS, do you have any repository or ros package related to it ?

— Reply to this email directly, view it on GitHub https://github.com/professorfabioandrade/georef/issues/4#issuecomment-2144512145, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACEWBLDSEBIE3KUIC3HWPP3ZFQOFJAVCNFSM6AAAAABIFCQB4CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNBUGUYTEMJUGU . You are receiving this because you commented.Message ID: @.***>

--

Fabio Andrade, PhD Associate Professor at University of South-Eastern Norway Research Scientist at NORCE Norwegian Research Center professorfabioandrade.com https://www.professorfabioandrade.com

--

Fabio Andrade, PhD Associate Professor at University of South-Eastern Norway Research Scientist at NORCE Norwegian Research Center professorfabioandrade.com https://www.professorfabioandrade.com

professorfabioandrade commented 5 months ago

ChatGPT forgot to add the skew factor in the undistortion function.

Here is the corrected one:

def undistort_point(u, v, k1, k2, k3, k4, p1, p2, fx, fy, cx, cy, skew):
    # Normalize the coordinates
    x_norm = (u - cx) / fx
    y_norm = (v - cy - skew * x_norm) / fy

    # Compute the radius squared
    r2 = x_norm**2 + y_norm**2
    r4 = r2 * r2
    r6 = r4 * r2
    r8 = r6 * r2

    # Compute radial distortion
    radial_distortion = 1 + k1 * r2 + k2 * r4 + k3 * r6 + k4 * r8

    # Compute tangential distortion
    x_tangential = 2 * p1 * x_norm * y_norm + p2 * (r2 + 2 * x_norm**2)
    y_tangential = p1 * (r2 + 2 * y_norm**2) + 2 * p2 * x_norm * y_norm

    # Apply the distortion to the normalized coordinates
    x_undistorted = x_norm * radial_distortion + x_tangential
    y_undistorted = y_norm * radial_distortion + y_tangential

    # Convert back to pixel coordinates
    u_undistorted = fx * x_undistorted + cx
    v_undistorted = fy * y_undistorted + cy + skew * x_undistorted

    return u_undistorted, v_undistorted
peakyquest commented 5 months ago

Thanks for helping me out. I will update my code according to the modification you have done, and check whether the error is reduced or not.