Open Ryan-Lima opened 5 years ago
Here you will find my images, and what I did to try to rectify them in: what_ive_done.py https://github.com/Ryan-Lima/sandbar_image_rectification
Your code sample is a bit too long and hard to follow. But what I can guess is that your homography and/or camera parameters are not being defined properly.
I wrote a script in python to calibrate the camera and they will generate the coefficients in the form Flamingo needs. Have a look here:
https://github.com/caiostringari/pywavelearn/blob/master/scripts/calibrate_camera.py
then you can load the arrays using this function (I will translate this to JSON I/O at some point this year):
def camera_parser(fname):
"""
Read the camera matrix and distortion coefficients generated from
"../scripts/calibrate.py"
----------
Args:
fname [Mandatory (str)]: Path to the camera calibration file.
----------
Returns:
K [ (np.ndarray)]: Camera matrix (OpenCV conventions)
DC [ (np.ndarray)]: Vector with distortion coefficients
[k_1, k_2, p_1, p_2(, k_3(, k_4, k_5, k_6))]
of 4, 5, or 8 elements.
"""
# open the file
f = open(fname, "r")
# read all lines
lines = f.readlines()
# build the camera matrix
K = np.zeros([3, 3])
for l, line in enumerate(lines):
if "Camera Matrix:" in line:
break
K[0, 0] = float(lines[l+2].strip("\n").split(",")[0])
K[0, 1] = float(lines[l+2].strip("\n").split(",")[1])
K[0, 2] = float(lines[l+2].strip("\n").split(",")[2])
K[1, 0] = float(lines[l+3].strip("\n").split(",")[0])
K[1, 1] = float(lines[l+3].strip("\n").split(",")[1])
K[1, 2] = float(lines[l+3].strip("\n").split(",")[2])
K[2, 0] = float(lines[l+4].strip("\n").split(",")[0])
K[2, 1] = float(lines[l+4].strip("\n").split(",")[1])
K[2, 2] = float(lines[l+4].strip("\n").split(",")[2])
# read the distortion coeficients
DC = []
for l, line in enumerate(lines):
if "Distortion Coeficients:" in line:
break
DC.append(float(lines[l+2].split(":")[1])) # K1
DC.append(float(lines[l+3].split(":")[1])) # K2
DC.append(float(lines[l+4].split(":")[1])) # P1
DC.append(float(lines[l+5].split(":")[1])) # P2
DC.append(float(lines[l+6].split(":")[1])) # K3
return K, np.array(DC)
Then, try this (Replace USERINPUT with your own data):
# read frame
img = os.path.isfile(USERINPUT)
if img:
Ir = skimage.io.imread(USERINPUT)
h, w = Ir.shape[:2]
else:
raise ValueError("Could not find input frame.")
# camera matrix
K,DC = ipwl.camera_parser(output from calibrate_camera.py)
# read XYZ coords
dfxyz = pd.read_csv(USERINPUT) # this is a csv file with x, y, z columns of data
XYZ = dfxyz[["x","y","z"]] .values
gcp_x = XYZ[0,0]
gcp_y = XYZ[0,1]
# read UV coords
dfuv = pd.read_csv(USERINPUT) # this is a csv file with u, v columns of data
UV = dfuv[["u","v"]].values
# horizon
hor = float(USERINPUT)
# rotation Angle
theta = float(USERINPUT)
# translations
xt = float(USERINPUT)
yt = float(USERINPUT)
# projection height
z = float(USERINPUT)
# undistort image
Kn,roi = cv2.getOptimalNewCameraMatrix(K,DC,(w,h),1,(w,h))
Ir = cv2.undistort(Ir, K, DC, None, Kn)
# homography
H = ipwl.find_homography(UV, XYZ, K, z=z, distortion=0)
# rectify coordinates
X,Y = ipwl.rectify_image(Ir, H)
# find the horizon limits
horizon = ipwl.find_horizon_offset(X, Y, max_distance=hor)
# rotate and translate
Xr, Yr = ipwl.rotate_translate(X, Y, rotation=theta, translation=[xt,yt])
# final arrays
if hor == -999:
Xc = Xr
Yc = Yr
Ic = Ir
else:
Xc = Xr[horizon:,:]
Yc = Yr[horizon:,:]
Ic = Ir[horizon:,:,:]
When I follow your instructions, I do get slightly different values in the K matrix. However when I tried to use plot_rectified(X,Y,imgs) I still get the same value error in the title of this post. I wonder if it has something to do with the horizona, rotation, and translations. in your example above you have the lines
hor = float(USERINPUT)
# rotation Angle
theta = float(USERINPUT)
# translations
xt = float(USERINPUT)
yt = float(USERINPUT)
# projection height
z = float(USERINPUT)
I am not sure what to put in for these values. My images are in a canyon, so there is no horizon line. the elevation is somewhere near 860m I used z = np.mean(XYZ[:,2]). I am not sure what to input for rotation or translation though....
I am able to get X and Y values using ipwl.rectify_image(Ir,H) however, I want to plot the image to see what it looks like and see if the rectification worked...., and that is where I am stuck. I hoped the ipwl.plot_rectified(X,Y,imgs) would work but it does not. In that it returns
ValueError: Zero-size array to reduction operation minimum which has no identity
Try something like this:
# rotation angle
theta = 0.0
# translations
xt = 0.0
yt = 0.0
# projection height
z = 0.0 # you can vary this value
# undistort image
Kn,roi = cv2.getOptimalNewCameraMatrix(K,DC,(w,h),1,(w,h))
Ir = cv2.undistort(Ir, K, DC, None, Kn)
# homography
H = find_homography(UV, XYZ, K, z=z, distortion=0)
# rectify coordinates
X,Y = rectify_image(Ir, H)
Now, don't use plot_rectified directly.
Do:
fig, ax = plt.subplots()
ax.pcolormesh(X, Y, Ir.T, cmap="Greys_r") # not sure if you need Ir or Ir.T here.
plt.show()
Have a look and see if you get any sensible results.
For some reason the K and DC values were wrong, and caused cv2.undistort to fail. However, using the K and DC values I extracted from the matlab camera calibration software, I maged to get it to work.
Also, the image has to be grayscale to use ax.pcolormesh()
Thanks!
I'll post the full code that worked for me soon, so others can follow along.
You can also add a RGB overlay to your plot if you want (as originally in plot_rectified):
# plot as before
im = ax.pcolormesh(X, Y, np.mean(Ir.T, axis=2), cmap="Greys_r")
# add colors
rgba = _construct_rgba_vector(Ir, n_alpha=n_alpha) # from flamingo.rectification import _construct_construct_rgba_vector
im.set_array(None) # remove the array
im.set_edgecolor('none')
im.set_facecolor(rgba)
For some reason the K and DC values were wrong, and caused cv2.undistort to fail. However, using the K and DC values I extracted from the matlab camera calibration software, I maged to get it to work. Also, the image has to be grayscale to use ax.pcolormesh()
Thanks!
I'll post the full code that worked for me soon, so others can follow along.
Could you post the opencv error? I may work to fix it in the future if I find some time...
Here is what I did. And the real problem (I think) is that the camera calibration is not work, so when I used cv2.undistort() it did not undistort the image properly
first I generated the camera calibration on the chessboard images
python calibrate.py -c "Canon_Rebel_30mile" -i H:\Hypsomtery_paper\Rectification\chessboard_images -fmt "JPG" -p 4 5 -o "Calibration_parameters.txt"
I provided 14 images but it only found the chess board in 4 of those 14. Also, Even though my chess board was 6x6 with 5x5 internal corners, it only managed to find corners if I provided 4x5 as the pattern
Mean Error 18.673 pixels (seems like a large error rate)
Camera : Canon_rebel_RC0307Rf
Camera Matrix:
8975.4345,0.0000,1973.8829
0.0000,3525.1319,1230.9077
0.0000,0.0000,1.0000
Pixel Centers:
U0 : 1973.8829088961977 U0 : 1230.907733571387
Distortion Coeficients:
k1 : -186.98709935770518 k2 : 9817.00603175839 p1 : 3.4189540645420666 p2 : -0.6669167341407937 k3 : 62.097644731053244
then I extracted those parameters and tried to undistort
import numpy as np
import cv2
import os
def camera_parser(fname):
"""
Read the camera matrix and distortion coefficients generated from
"../scripts/calibrate.py"
----------
Args:
fname [Mandatory (str)]: Path to the camera calibration file.
----------
Returns:
K [ (np.ndarray)]: Camera matrix (OpenCV conventions)
DC [ (np.ndarray)]: Vector with distortion coefficients
[k_1, k_2, p_1, p_2(, k_3(, k_4, k_5, k_6))]
of 4, 5, or 8 elements.
"""
# open the file
f = open(fname, "r")
# read all lines
lines = f.readlines()
# build the camera matrix
K = np.zeros([3, 3])
# had to adjust these because the provided code was actually grabbing the \n lines in between and not getting the camera matrix values
K[0, 0] = float(lines[6].strip("\n").split(",")[0])
K[0, 1] = float(lines[6].strip("\n").split(",")[1])
K[0, 2] = float(lines[6].strip("\n").split(",")[2])
K[1, 0] = float(lines[8].strip("\n").split(",")[0])
K[1, 1] = float(lines[8].strip("\n").split(",")[1])
K[1, 2] = float(lines[8].strip("\n").split(",")[2])
K[2, 0] = float(lines[10].strip("\n").split(",")[0])
K[2, 1] = float(lines[10].strip("\n").split(",")[1])
K[2, 2] = float(lines[10].strip("\n").split(",")[2])
# read the distortion coeficients
DC = []
DC.append(float(lines[20].split(":")[1])) # K1
DC.append(float(lines[21].split(":")[1])) # K2
DC.append(float(lines[22].split(":")[1])) # P1
DC.append(float(lines[23].split(":")[1])) # P2
DC.append(float(lines[24].split(":")[1])) # K3
return K, np.array(DC)
K, P = camera_parser("Calibration_parameters.txt")
K array([[8.9754345e+03, 0.0000000e+00, 1.9738829e+03], [0.0000000e+00, 3.5251319e+03, 1.2309077e+03], [0.0000000e+00, 0.0000000e+00, 1.0000000e+00]]) P array([-1.86987099e+02, 9.81700603e+03, 3.41895406e+00, -6.66916734e-01, 6.20976447e+01])
K = np.array([[3707.28817166849,2.73440717127128,1853.53725535792],[0,3718.83680802725,1196.97037378192],[0,0,1]]) P = np.array([-0.154149563752979,0.147525526652382,0.00161089070708720,0.00200212660140590]) note that K3 = 0 in matlab calibration
GCPpoints_rw = [(219564.143,611772.081,856.506),(219545.827,611789.288,856.596),(219517.774,611728.37,857.519),
(219513.184,611704.139,865.029),(219540.082,611757.827,859.451),(219524.911,611831.131,876.81),(219541.506,611806.246,858.974),
(219570.485,611856.435,857.012),(219590.036,611862.6,857.551)]
## reshape GCP real world points
rwpts = np.array(GCPpoints_rw).reshape(9,3)
## UV location of GCP panels in image RC0307Rf_20171002_1129.JPG
GCP_imgpoints = [(1104,1295),(1543,1090.5),(642.25,997.75),(322.25,1017),(968.25,1060.5),(2350.75,471),(1883.75,974.25),
(3036.25,1158.25),(3310.75,1356)]
## reshape UV image points
impts = np.array(GCP_imgpoints).reshape(9,2)
## I downloaded the rectification.py, plot.py, and __init__.py from Openearth/flamingo
from rectification import *
from plot import *
## find homography
H = find_homography(impts,rwpts, K, z = 0, distortion = 0)
img = "RC0307Rf_20091012_0942_regis.JPG"
## seperate the subject image into filename, extention, and relative path
filename, ext = os.path.splitext(img)
ir = cv2.imread(img)
## convert to grayscale
ir_gray = cv2.cvtColor(ir,cv2.COLOR_RGB2GRAY)
h,w = ir.shape[:2]
## undistort image
Kn, roi = cv2.getOptimalNewCameraMatrix(K,P,(w,h),1,(w,h))
ir_undst = cv2.undistort(ir,K,P,None,Kn)
ir_gray_undst = cv2.undistort(ir_gray,K,P,None,Kn)
plt.imshow(ir_undst)
plt.show()
K = np.array([[3707.28817166849,2.73440717127128,1853.53725535792],[0,3718.83680802725,1196.97037378192],[0,0,1]])
P = np.array([-0.154149563752979,0.147525526652382,0.00161089070708720,0.00200212660140590])
## note that K3 = 0 in matlab calibration
## find homography
H = find_homography(impts,rwpts, K, z = 0, distortion = 0)
img = "RC0307Rf_20091012_0942_regis.JPG"
## seperate the subject image into filename, extention, and relative path
filename, ext = os.path.splitext(img)
ir = cv2.imread(img)
## convert to grayscale
ir_gray = cv2.cvtColor(ir,cv2.COLOR_RGB2GRAY)
h,w = ir.shape[:2]
## undistort image
Kn, roi = cv2.getOptimalNewCameraMatrix(K,P,(w,h),1,(w,h))
ir_undst = cv2.undistort(ir,K,P,None,Kn)
ir_gray_undst = cv2.undistort(ir_gray,K,P,None,Kn)
plt.imshow(ir_undst)
plt.show()
[https://drive.google.com/drive/folders/1JrSDR1MTAQB0d6z5VVwc8_74tWnTCT5u?usp=sharing]
Here is the calibrate.py that I used
#------------------------------------------------------------------------
#------------------------------------------------------------------------
#
#
# SCRIPT : calibrate.py
# POURPOSE : Find Intrinsic and Extrinsic camera parameters.
# AUTHOR : Caio Eadi Stringari
# EMAIL : Caio.EadiStringari@uon.edu.au
#
# V1.0 : 29/07/2016 [Caio Stringari]
#
#
# The algorithm is based on [Zhang2000] and [BouguetMCT]
#
# Zhang. A Flexible New Technique for Camera Calibration. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(11):1330-1334, 2000.
# Y.Bouguet. MATLAB calibration tool. http://www.vision.caltech.edu/bouguetj/calib_doc/
#
#------------------------------------------------------------------------
#------------------------------------------------------------------------
# system
import os
# arguments
import argparse
# files
from glob import glob
# numpy
import numpy as np
# OpenCV
import cv2
# I/O
from pandas import DataFrame
# Matplotlib
import matplotlib.pyplot as plt
def main():
# camera name
camera = args.camera[0]
# input image path
path = os.path.abspath(args.images[0])
# pattern
n1 = int(args.pattern[0])
n2 = int(args.pattern[1])
images = sorted(glob(path+"/*."+args.format[0]))
# prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
objp = np.zeros((n1*n2,3), np.float32)
objp[:,:2] = np.mgrid[0:n2,0:n1].T.reshape(-1,2)
# termination criteria
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 300, 0.0001)
# arrays to store object points and image points from all the images.
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.
nrets = 0 # counter
# loop over the images trying to identify the chessboard
for k,fname in enumerate(images):
print ("Loading image {} of {}".format(k+1,len(images)))
# read the image as 8 bit
img = cv2.imread(fname)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
# find the chess board corners
ret, corners = cv2.findChessboardCorners(gray, (n1,n2),None)
# if found, add object points, image points (after refining them)
if ret == True:
print (" Chessboard pattern found in image {} of {}".format(k+1,len(images)))
# refine corners
corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
# append
objpoints.append(objp)
imgpoints.append(corners2)
# store the last processed images
I = gray
nrets+=1
# Draw and display the corners
img = cv2.drawChessboardCorners(img, (n1,n2), corners2,ret)
cv2.imshow('img',img)
cv2.waitKey(500)
if nrets > 10: break
cv2.destroyAllWindows()
# calculating the calibration
h, w = I.shape[:2]
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)
# get the intrinsic camera variables [See Holman et al 1997] for terminology
fx_px = mtx[0,0] # Focal lenght in pixels [x]
fy_px = mtx[1,1] # Focal lenght in pixels [y]
u0 = mtx[0,2]
v0 = mtx[1,2]
k1 = dist[0][0]
k2 = dist[0][1]
p1 = dist[0][2]
p2 = dist[0][3]
k3 = dist[0][4]
ud = np.arange(1,w+1,1)
vd = np.arange(1,h+1,1)
# calculating the error
error = []
for i in range(len(objpoints)):
imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
error.append(cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2))
error = np.array(error)
merror = round(error.mean(),3)
print ("Mean error is {} pixels".format(merror))
# write the parameters to file
f = open(args.output[0],"w")
f.write("### Camera Calibration Results ###\n")
f.write("\n")
f.write("Camera : {}\n".format(camera))
f.write("\n")
f.write("Camera Matrix:\n")
f.write("\n")
df = DataFrame(mtx)
df.to_csv(f,mode="a",header=False,index=False,float_format="%.4f")
f.write("\n")
f.write("Pixel Centers:\n")
f.write("\n")
f.write("U0 : {}\n".format(u0))
f.write("U0 : {}\n".format(v0))
f.write("\n")
f.write("Distortion Coeficients:\n")
f.write("\n")
f.write("k1 : {}\n".format(k1))
f.write("k2 : {}\n".format(k2))
f.write("p1 : {}\n".format(p1))
f.write("p2 : {}\n".format(p2))
f.write("k3 : {}\n".format(k3))
if __name__ == '__main__':
parser = argparse.ArgumentParser()
# Camera name
parser.add_argument('--camera','-c',
nargs = 1,
action = 'store',
default = ["My Camera"],
dest = 'camera',
required = False,
help = "Camera name.",)
# Input images location
parser.add_argument('--images','-i',
nargs = 1,
action = 'store',
dest = 'images',
required = True,
help = "Input images folder.",)
# Input images location
parser.add_argument('--format','-fmt',
nargs = 1,
action = 'store',
dest = 'format',
required = True,
help = "Input images format. Usually is \'jpg\', \'JPG\', or \'png\'.",)
# Define the pattern to search
parser.add_argument('--pattern','-p',
nargs = 2,
action = 'store',
dest = 'pattern',
required = True,
help = "Number of rows and colums in the chessboard. Ususally 1-ncols and 1-nrows.",)
# Output file
parser.add_argument('--output','-o',
nargs = 1,
action = 'store',
dest = 'output',
required = True,
help = "Output file name.",)
# handle inputs
args = parser.parse_args()
# main call
main()
Yes, your camera calibration with my script does not look correct. I usually calibrate the camera in my lab and not in the field. Also, I use a 8x8 calibration board recommended by OpenCV: https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_calib3d/py_calibration/py_calibration.html. I would recommend to use a bit more images if possible. I usually use no less than 30.
Yes, I realize the camera calibration is not Ideal. These are stationary cameras that have been in fixed locations since 1990, though the cameras have been changed out over time. They were not calibrated before they were deployed. The Cameras are perched on cliffs so getting in front of them is tricky, if not impossible in some cases. They are fitted with zoom lenses and we adjust the zoom based on the field of view and feasible camera location so its difficult to know what focal length we will use in the field beforehand in the lab to calibrate. At this particular site (in the bottom of the Grand Canyon) my field tech took the images of the checkerboard pattern, obviously not enough, and some are partially out of the frame. I did not originally intend to use cv2 for the camera calibration, and Matlab calibration did not seem to require a 9x9. In any case, this is what I have to work with. Any suggestions on where I went wrong with your script?
I am relatively new to python. I am trying to utilize your rectification.py and plot.py to Geo-rectify oblique images in the Grand Canyon and reproject them in a birdseye view projected in latitude and longitude.
I used matlabs camera calibration software to extract the K and P matricies, then I provided ground control points in both UV and XYZ to generate the homography H.
However, when I try to plot the lists of X,Y values and the list of images using fig, ax = plot.plot_rectified(Xlist,Ylist,imlist)
I receive the following error
ValueError: Zero-size array to reduction operation minimum which has no identity
Here is my code `` import numpy as np import matplotlib.pyplot as plt from rectification import* import plot
camera intrinsic matrix, manually exported from matlab camera calibration tool
K = np.array([[3707.28817166849,2.73440717127128,1853.53725535792],[0,3718.83680802725,1196.97037378192],[0,0,1]])
in matlab this matrix is transposed so I will try it with both...
Kt = np.matrix.transpose(K)
tangential distortion parameters
P = np.array([-0.154149563752979,0.147525526652382,0.00161089070708720,0.00200212660140590])
locations of the GCP panels in real world coordinates
GCPpoints_rw = [(219564.143,611772.081,856.506),(219545.827,611789.288,856.596),(219517.774,611728.37,857.519), (219513.184,611704.139,865.029),(219540.082,611757.827,859.451),(219524.911,611831.131,876.81),(219541.506,611806.246,858.974), (219570.485,611856.435,857.012),(219590.036,611862.6,857.551)]
rwpts = np.array(GCPpoints_rw).reshape(9,3)
location of GCP panels in image RC0307Rf_20171002_1129.JPG
GCP_imgpoints = [(1104,1295),(1543,1090.5),(642.25,997.75),(322.25,1017),(968.25,1060.5),(2350.75,471),(1883.75,974.25), (3036.25,1158.25),(3310.75,1356)]
impts = np.array(GCP_imgpoints).reshape(9,2)
find homography using camera intrinsic matrix or solve for extrinsics
H = find_homography(impts,rwpts,K,P) HKt = find_homography(impts,rwpts,Kt,P)
import registered image
img1 = cv2.imread("RC0307Rf_20091012_0942_regis.jpg") img2 = cv2.imread("RC0307Rf_20130304_1048_regis.jpg") img3 = cv2.imread("RC0307Rf_20140527_1335_regis.JPG") img4 = cv2.imread("RC0307Rf_20150509_1446_regis.JPG")
X1, Y1 = rectify_image(img1,H) X2, Y2 = rectify_image(img2,H) X3, Y3 = rectify_image(img3,H) X4, Y4 = rectify_image(img4,H)
Xlist = [X1,X2,X3,X4] Ylist = [Y1,Y2,Y3,Y4] imlist =[img1,img2,img3,img4]
fig, ax = plot_rectified(Xlist,Ylist,imlist) # does not work.....
ValueError Traceback (most recent call last)