Calculating world coordinates of a pixel from a camera picture

Question

Sorry for this question. I know there are many similar questions but I really know nothing about math that involves this case.

I have a picture taken from a camera. I know some information about this camera and this picture:

Camera:

* FOV       : 61 deg.
* Altitude  : 91 meters
* Long.     : -43.17687898427574
* Lat.      : -22.89925324277859
* Azimuth   : 109 deg.
* Pitch     : 91 deg. ( 1 degree below the horizon )

These numbers came from Cesium. I have no access to the camera specifications so they can be with approximate values. However, Cesium was extremely accurate in calculating the actual pixel coordinates of the image, so I think these numbers are pretty close to reality.

Picture:

* Height    : 485px
* Width     : 496px

Now, I need to know the real world geographic coordinates of a specific pixel in this image.

I've made a test with Cesium, and I was able to check some results to serve as test parameters:

* Pixel at 220 x 322 (w/h)
* Real world coordinates of that pixel: 
    Lat: -22.89974712930635
    Lon: -43.167470162955375

Cesium is a very good tool to do this job, but I need some backend resource like Java or Python and can't find any ( maybe due to my math limitations I can't determine what works and what doesn't ).

I found a material I think is very close to what I need, but can't implement it because it is very theoretical: https://medium.com/swlh/ray-tracing-from-scratch-in-python-41670e6a96f9

Here is the picture from the test case. The red arrow is where the pixel is (approximately, of course. I'm not so precise with the mouse). Have in mind that this picture is a screenshot and may not have same dimensions that the real one I have.:

... and this is the screenshot of my Cesium map to the same area ( but in 3D view ):

... and 2D view ( the red mark is the point used in the test case. The Blue mark is where the camera is):

Resuming, I need to take the geographic coordinates of a pixel from a picture using Python or Java ( or any agnostic math pseudo code for laymen that I could implement ).

EDIT:

Well... I have some help from ChatGPT and now I have this code in Python (sorry for PT_BR comments):

import numpy as np

# FOV       : 61 deg.
# Altitude  : 91 meters
# Long.     : -43.17687898427574
# Lat.      : -22.89925324277859
# Azimuth   : 109 deg.
# Pitch     : 91 deg. ( 1 degree below the horizon )

# Width     : 859.358px
# Height    : 484.983px

# Pixel at 850.1561889648438 x 475.18054962158203 (w/h)

camera_fov = 61            # FOV
image_width = 859.358      # largura da imagem em pixels
image_height = 484.983     # altura da imagem em pixels
camera_height = 90.0       # altura da câmera em metros
pitch_angle = 91.0         # ângulo de inclinação da câmera em graus
yaw_angle = 109.0          # ângulo de guinada da câmera em graus
camera_position = np.array([-43.17687898427574, -22.89925324277859, camera_height])  # posição da câmera no sistema de coordenadas do mundo

# coordenadas 2D do pixel na imagem ( u=w  v=h)
u = 850.1561889648438     # w
v = 475.18054962158203    # h

# cálculo das coordenadas 3D do pixel no espaço do mundo
aspect_ratio = image_width / image_height
fov = np.radians( camera_fov ) # campo de visão em radianos
y_fov = 2 * np.arctan(np.tan(fov / 2) * np.cos(np.radians(pitch_angle)))  # campo de visão vertical em radianos
y_angle = np.radians((v / image_height - 0.5) * y_fov * 180 / np.pi + pitch_angle)  # ângulo vertical em radianos
x_fov = 2 * np.arctan(np.tan(fov / 2) * np.cos(y_angle))  # campo de visão horizontal em radianos
x_angle = np.radians((u / image_width - 0.5) * x_fov * 180 / np.pi)  # ângulo horizontal em radianos
direction = np.array([-np.sin(x_angle), np.cos(x_angle) * np.sin(y_angle), np.cos(x_angle) * np.cos(y_angle)])
rotation_matrix = np.array([
    [np.cos(np.radians(yaw_angle)), -np.sin(np.radians(yaw_angle)), 0],
    [np.sin(np.radians(yaw_angle)), np.cos(np.radians(yaw_angle)), 0],
    [0, 0, 1]
])

direction = rotation_matrix @ direction
world_coords = camera_position + direction * camera_height / np.tan(y_angle)

# Real world coordinates of that pixel: 
#    Lat: -22.901614465334827
#    Lon: -43.17465101364515


print(world_coords[1],"/",world_coords[0] )

I've found the numbers very close to what Cesium gave me, then I've made some tests. First I found the pixel near the horizon. All pixels in height less than it will be sky. Then I take all four pixels from horizon to the image bottom:

The result is ( pixel pos ) LAT , LON : ( coordinates calculated by the code )

p1 : (   1,264)  -22.40093040092381 , -41.77409577136955
p2 : ( 858,254)  -22.42828595321902 , -41.76467637515944
p3 : (   1,481)  -22.68142494320101 , -42.55301747856779 
p4 : ( 858,481)  -22.68681453552334 , -42.55116168224547

Now I need to check these coordinates the code gave against those the Cesium calculated to see if the results are the same.

Don't, but was close.

I'm curious about what coordinates the code are giving me, so I take a WKT Polygon (GeoJSON) and use the site geojson.io to draw it.

The gray area is the polygon representing the "viewport" ( to the horizon ) that I made from the camera picture ( p1 to p4 ) and calculated by the code.

The red arrow is where the camera is.

The Blue triangle is what I think the camera picture are actualy seeing.

As you can see I'm very close. I also suspect that input values ​​such as pitch and yaw (azimuth) orientation may be different than what Cesium uses. Any help to check the math in the code will be wellcome!!