Inversing brown-conrady equations (for camera distortion)

[elliebowler]

Hi,

I am working on an application where i need to project 2d image pixel coordinates into the 3D camera coordinate system. I understand that this can be done with the metashape .unproject method, for example:

Code: [Select]

chunk = Metashape.app.document.chunk
camera = chunk.cameras[0] # test using first camera
cam_calib = camera.sensor.calibration # camera 0 calibration values

test_pixel = Metashape.Vector([100, 100]) # random point in image pixel coordinates
X, Y, Z = cam_calib.unproject(test_pixel) # result using unproject function
However I would like to understand the equations behind this reverse projection. In particular i have radial and tangental distortion parameters which i need to include. I understand and can replicate the forward projection from camera crs to pixel, as below, using the equations for frame cameras at the end of the documentation https://www.agisoft.com/pdf/photoscan-pro_1_4_en.pdf. Using this i get the point (100, 100) as expected.

Code: [Select]

test_pixel = Metashape.Vector([100, 100]) # random point in image pixel coordinates
X, Y, Z = cam_calib.unproject(test_pixel) # result using unproject function

# Test projection from camera crs to test pixel...
print('Inputing ' + str(X) + ',' + str(Y) + ',' + str(Z))
print('Expecting ' + str(test_pixel))

x = X/Z
y = Y/Z

### distortion corrections
r = math.sqrt(x**2 + y**2)

x_rad = x*(1 + cam_calib.k1*r**2 + cam_calib.k2*r**4 + cam_calib.k3*r**6 + cam_calib.k4*r**8)
x_tang = (cam_calib.p1*(r**2 + 2*x**2) + 2*cam_calib.p2*x*y)*(1 + cam_calib.p3*r**2 + cam_calib.p4*r**4)
y_rad = y*(1 + cam_calib.k1*r**2 + cam_calib.k2*r**4 + cam_calib.k3*r**6 + cam_calib.k4*r**8)
y_tang = (cam_calib.p1*(r**2 + 2*y**2) + 2*cam_calib.p2*x*y)*(1 + cam_calib.p3*r**2 + cam_calib.p4*r**4)

x_dash = x_rad + x_tang
y_dash = y_rad + y_tang
###

# Equivalent to multiplying x_dash, y_dash by intrinsic matrix K
u = cam_calib.width*0.5 + cam_calib.cx + x_dash*cam_calib.f + x_dash*cam_calib.b1 + y_dash*cam_calib.b2
v = cam_calib.height *0.5 + cam_calib.cy + y_dash*cam_calib.f
# i.e K = [[f + b1, b2, cx + width*0.5], [0, fy, cy + height*0.5], [0,0,1]]
# so [u,v] = K.mulp(x_dash, y_dash)

print('Got: ' + str(u) + ',' + str(v))
I also understand that i can invert the intrinsic matrix K in the final stage, and have also separately saved depth maps so i can later use the z values. Really it is the section where you apply the brown-conrady distortion equations which i'm unsure how to perform in reverse. Would it be possible to give any insight into how this is performed in the .unproject method? I have parameters f, cx, cy, k1, k2, k3, p1 and p2 for this particular camera model.

Thanks in advance anyone who might be able to advise!

[Paulo]

Hi EllieBowler,

going from x', y' to x,  y you can use a first approximation as:

Code: [Select]

x₀ ≈ x'  - x' * (k1 * r'² + k2 * r'⁴ + k3 * r'⁶ + k4 * r'⁸) - p1 * (r'² + 2 * x'²) - 2 * p2 * x' * y' = x' + Δx₀
y₀ ≈ y'  - y' * (k1 * r'² + k2 * r'⁴ + k3 * r'⁶ + k4 * r'⁸) - p2 * (r'² + 2 * y'²) - 2 * p1 * x' * y' = y' + Δy₀  where r' = sqrt( x'² + y'²)

and then iterate x, y as

Code: [Select]

x₁ = x'  - x₀ * (k1 * r₀² + k2 * r₀⁴ + k3 * r₀⁶ + k4 * r₀⁸) - p1 * (r₀² + 2 * x₀²) - 2 * p2 * x₀ * y₀ = x' + Δx₁
y₁ = y'  - y₀ * (k1 * r₀² + k2 * r₀⁴ + k3 * r₀⁶ + k4 * r₀⁸) - p2 * (r₀² + 2 * y₀²) - 2 * p1 * x₀ * y₀ = y' + Δy₁  where r₀ = sqrt( x₀² + y₀²)
and after a few iterations it will converge to final x, y as in :

x'  y'  Δx₀  Δy₀ = -0.416170108   0.187124851   -0.001347323   -0.000846066
         
x₀  y₀  Δx₁  Δy₁ = -0.417517431   0.186278785   -0.001365358   -0.00084913
         
x₁  y₁  Δx₂  Δx₂ = -0.417535466   0.186275721   -0.001365583   -0.000849169
         
x₂  y₂  Δx₃  Δy₃ = -0.417535691   0.186275682   -0.001365586   -0.00084917
   
x₃  y₃  Δx₄  Δy₄ = -0.417535694   0.186275681   -0.001365586   -0.00084917

where k1, k2, k3, k4, p1, p2  = -0.0435999,   0.0398564,   -0.0155917,   0,   -0.00109523,   0.00746545

Maybe this can help,