Compute camera ray (direction vector) for a pixel

[geo_enth3]

Hi,

So I am trying to compute 3D world coordinates for certain pixels within my aligned photo. However, I do not want to intersect the camera ray with the 3D model (as the pickPoint() function does) but I want to intersect it with a custom plane (which I have computed from the sparse point cloud). The reason is that the 3D model is not available or sometimes has holes.

For intersecting the camera ray of a certain pixel (X/Y) with the plane I need to compute the direction vector (cameraCenter to Pixel). As I understand the upper right 3X3 matrix of

Code: [Select]

rotation = np.asarray(segmentCam.transform).reshape(4,4)[0:3,0:3]    contains the rotation matrix that allows the transformation of the pixel coordinates in the 3D world coordinates. On this I then apply the chunk transformation matrix to derive the real world coordinates:

Code: [Select]

M = np.asarray(chunk.transform.matrix).reshape(4,4)[0:3,0:3]                 # extracting the transformation from internal Metashape to world coordinates
rotation_M = np.matmul(M,rotation)   
and this Matrix I finally apply on the internal camera coordinates which I picked (and corrected using the camera calibration parameters):

Code: [Select]

b = np.matmul(np.transpose(rotation_M), b_loc)       
However, the results are not correct and I strongly suspect that I am misunderstanding the two transformation matrices. Could somebody help me with estimating the direction vector (cameraCenter to Pixel)?

Thanks!

Here is the complete code for a better understanding:

Code: [Select]

segmentCam = chunk.cameras[0]
R, C = getProjectionPlane(chunk)                                             # derive plane from sparse point cloud (works fine)
                                                                             # R....3x3 matrix containting the plane defining direction vectors
                                                                             # C...barymetric center of the plane (EPSG:31256)

C_camera = chunk.crs.project(chunk.transform.matrix.mulp(segmentCam.center)) # 3D coordinates of Camera center [m] (EPSG: 31256)
row = np.asarray([3000,1000])                                                # image coordinates [px]
b_temp1 = np.append(row,0)                                                   # homogeneous image coordinates [px]
b_temp2 = np.array([segmentCam.calibration.cx, segmentCam.calibration.cy,    # calibration parameters of the camera
                    segmentCam.calibration.f])
b_loc = b_temp1-b_temp2                                                      # image coordinates in the "local" camera coordinate system

rotation = np.asarray(segmentCam.transform).reshape(4,4)[0:3,0:3]            # extractin the camera transformation parameters (3X3 matrix)
M = np.asarray(chunk.transform.matrix).reshape(4,4)[0:3,0:3]                 # extracting the transformation from internal Metashape to world coordinates
rotation_M = np.matmul(M,rotation)                                           # combining the matrices by multiplying
b = np.matmul(np.transpose(rotation_M), b_loc)                               # applying the transformations to the picked camera coordinates
b_norm = b / np.linalg.norm(b)                                               # normalizing (not really necessary)

C_plane = Metashape.Vector(C)                                                # barycenter of the plane
a1 = np.asarray(R).reshape(3,3)[0,0:3]                                       # first plane direction vector
a2 = np.asarray(R).reshape(3,3)[1,0:3]                                       # second plane direction vector



A = np.transpose(np.array([a1,a2,-b]))                                       # solve equation (intersect plane with image ray)
param = np.matmul(np.linalg.inv(A),(C_camera-C_plane))

S1 = a1*param[0] + a2*param[1] + C_plane                                     # apply parameters to the equation and retrieve the
                                                                             # intersection of the image ray with the plane -> S1==S2
S2 = b*param[2] + C_camera

[geo_enth3]

I solved the issue myself.

Just in case somebody faces the same issue here is the code which solved the problem explained above:

Code: [Select]

pt_2d  = Metashape.Vector([3091,2034])
img_3d = camera.unproject(pt_2d)
ray = chunk.crs.project(chunk.transform.matrix.mulp(img_3d)) - chunk.crs.project(chunk.transform.matrix.mulp(camera.center))
b = ray/ray.norm()

C_plane = Metashape.Vector(C)
a1 = np.asarray(R).reshape(3,3)[0,0:3]
a2 = np.asarray(R).reshape(3,3)[1,0:3]



A = np.transpose(np.array([a1,a2,-b]))
param = np.matmul(np.linalg.inv(A),(C_camera-C_plane))

S1 = a1*param[0] + a2*param[1] + C_plane

S2 = b*param[2] + C_camera