# -*- coding: utf-8 -*- # for interpreting unicode literals from numpy import * # for the array handling code from visual.filedialog import get_file # I mostly use visual python for its nice file menus import os.path # used for the filename routines import cv # get the opencv package ##import sys # Camera Calibration & Undistortion routines # Suzanne Amador Kane 11-15-11 http://www.haverford.edu/physics-astro/Amador/index.php # SOURCES FOR THE CODE & LEARNING ABOUT ALL THIS: # The code itself was adapted from the original C code from the OpenCV bible: http://opencv.willowgarage.com/wiki/OpenCVBooks # Learning OpenCV: Computer Vision with the OpenCV Library by Gary Bradski and Adrian Kaehler, Published by O'Reilly Media, October 3, 2008. # and translations from C to python from these sites: # https://github.com/abidrahmank/OpenCV-Python/blob/master/Other_Examples/camera_calibration.py # http://www.neuroforge.co.uk/index.php/camera-calibration # plus very helpful postings at: # Martin Peris's blog: http://blog.martinperis.com/2011/01/opencv-stereo-camera-calibration.html (the Download button leads to sample code & images) # Vezhnevets Vladimir site: He wrote the original checkerboard finding code! # http://www.graphicon.ru/oldgr/en/research/calibration/opencv.html # has helpful commentary, examples, & samples images (I used these to check out my code was working and to compare my values with his; they did # not seem rich enough to get truly undistorted images, though. Maybe because this version of OpenCV doesn't allow skew?) # Interesting commentary about doing this in Matlab from Bouguet: http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/example.html # (though I'm not able to use his posted example images--they appear to be too low res to find corners more than 1/4 of the time!) # # Build instructions: I am running all this on a Windows 7 machine using the following python build: # # python 2.7.1 http://www.python.org/getit/releases/2.7/ # visual python 2.7a0 http://vpython.org/index.html # opencv 2.2 http://opencv.willowgarage.com/wiki/PythonInterface # OpenCV documentation for python: http://opencv.willowgarage.com/documentation/python/index.html # numpy-1.5.1-py2.7 http://sourceforge.net/projects/numpy/files/NumPy/1.5.1/ # # # This program uses the regular python bindings for OpenCV (i.e., not ctypes_opencv or pyopencv. (I had problems building and finding the DLL's for ctypes opencv; # pyopencv has renamed and restructured all function calls, without providing documentation about correspondences!) # Purpose: First read in a series of checkerboard images taken with your camera, then find the corners of each # Then use these to find the camera distortion parameters, as well as the camera's intrinsic and extrinsic matrices # Finally, use to undistort images and print out the relevant matrices and other values # Inputs: image stacks in any readable image format, numbered consecutively with four digit numbers at the end of each file corename (before the extension) # for example: image0000.png is an acceptable filename # See the original code from Learning OpenCV pg 398 and the original python code above for how to stream video from a live camera feed # Outputs: .xml files Intrinsics & Distortion hold the relevant matrices found; these matrix values are also printed out to the python shell window # Units: if you properly input the size of your checkerboard square sizes in real world units (e.g. square_size=2.54 for 2.54 cm squares) # then the extrinsic matrix values will be in real world units (here, cm) # Focal lengths and image centers (cx,cy) are in pixels, as noted. # The units of the distortion coefficients (called k1, k2, etc. in the opencv book) are cryptic since they apply to # when you are using homogeneous coordinates for x and y. So, they are not in pixels and hence work for any size # image taken using the same lens! That's great, but how to think about their units. See this website for clues about # how to think about their native units: http://imagej.588099.n2.nabble.com/RE-LENS-DISTORTION-CORRECTION-td631005.html # (see second posting) # Using it: Look throughout for the triple asterisks: *** for modifications you'll have to make for your own use. My defaults work well for the example # images found via the Download button at Martin Peris' website: http://blog.martinperis.com/2011/01/opencv-stereo-camera-calibration.html # Getting good images: # You can test out your code using the sample images at: http://www.graphicon.ru/oldgr/en/research/calibration/opencv.html # Remember--you'll have to search for *** & change things for this to work! # Taking your own images: # > Use a white and black checkerboard image like the one I've provided at this website. # > It really needs to be black and white! Be sure to leave a wide white border around all four edges. # > Secure it so it's truly flat! I use masking tape, paper that's flat (not puckered from humidity) and a truly flat surface. # Be sure to double-check flatness & straightness using high-quality metal rulers/machinists flats/carpenters squares, etc. # > Use an assymetrical board with one edge odd and one even in number of corners. This breaks the symmetry and improves performance. # > Most people seem to use from 9 to 5 corners per edge. I found these values work well. More corners are problematic and harder to find. # (Make sure you input the number of corners below--not the number of square!) # > Take at least 10 or somewhat more images (Too many can overdetermine the finicky searches, while too few do not sufficiently determine them--read the book.) # > Try for an uncluttered background, especially without posters, windows, other corner-forming objects. A very flat painted wall or floortiles are good. # > Don't use really distorted perspectives (i.e.,very edge-on views) since that's hard for the checkerboard corner finder(search on FILTER_QUADS flag) # > Make sure your board isn't too small on the images--see suggestions below # > Use high enough resolution images--if it can't find the gradient, you are doomed. Actually look at your images first! # > Try for uniform lighting, with no strong reflected highlights on the checkerboard (though the adaptive thresholding flag can help here--search on ADAPTIVE_THRESH flag) # > The images need to be "balanced"--that is, evenly distributed through all parts of the image plane, with no bias toward any one region # > Suggested distribution of where the checkerboard ought to be: (from: http://stackoverflow.com/questions/3197523/camera-and-projector-calibration) # * Take one that has target taking up entire image # * four with target taking up each quadrant # * then another four with target taking up ~1/9 the area in each of the four corners. # > Note: bias to corners is often useful since lens distortion is usually bad at the corners and minimal in the middle. # # Troubleshooting: # > I recommend first trying to get this working with Martin Peris' sample images first (see http://blog.martinperis.com/2011/01/opencv-stereo-camera-calibration.html) # > Every time an image is posted, the program hangs, waiting for you to close it before proceeding. So, just click the image window closed and off you go. # > Make sure you've changed the values indicated by *** for your setup # > Look for the commented out print statements. I recommend at least uncommenting these and viewing the original images # > are original images OK? # > Is it finding the checkerboard corners OK? # > If not, see comments about images below and/or try different flag combinations # Note I've added code that indicates whether it found the correct # of corners, whether it found corners at all, and a file-by-file review # telling you which ones worked and which ones didn't work. That should help pinpoint what's going wrong. # # Everyone agrees that it's a bit of an art getting this all to work properly. Use the sample image sets to make sure you have the code working, # then work on getting good images for your camera. board_w = 9 # *** number of corners wide: I used the checkerboard from Martin Peris' website (see above) which has 9 X 6 corners board_h = 6 # *** number of corners high: set for your chessboard's dimensions--number of CORNERS NOT SQUARES!!! board_n = board_w*board_h # total number of corners board_sz = (board_w,board_h) # size of board square_size = 0.254/8. # *** size of the checkerboard squares in the real world--measure yours! (in my case, in units of meters) # *** this was what I got when printing out the checkerboard from Martin Peris' website (see above) on 11" x 17" paper win_subpix = 8 # *** the window size (2 x this value to an edge) used for the subpixel edge finding routine # *** look at your images to see over what region your light-to-dark transitions take place fish_eye_lens = 0 # *** set to 1 if using a fisheye lens; otherwise use FALSE (0) tangent_undistort = 0 # *** set to 1 if using tangential undistort--only for very distorted lenses--set to 0 most of the time! first_order_in_k = 0 # *** only use the first order correction in k (makes sense only for very high end lens, not wide angle!) fit_only_f_center = 0 # *** set = 1 to fit only the focal length and center fit_only_k1_k2 = 0 # *** if set = 1 what it sounds like--fix f and cx,cy, fit only k1, k2 f_x = 1080 # *** initial guesses at the x and y focal lengths (currently, I fix their aspect ratio at 1--see later flag f_y = 1080 # *** num_usable_images=0 # number of images where we could find the corners successes=0 # number of usable images with the correct number of corners num_ext= 2 # *** each image filename must end in four digits and be consecutively numbered (you will have to rewrite some code # *** the present code works if your file numbers are zero-padded at the left like this: 00, 01, 02, etc. # *** to deactivate this zeropadding (1,2,3,4,..., 10, 11, ...) then, remove all instances of .zfill(num_ext) from the code ######end of initialization ###################################################################### print "Get the first image file name" fd = get_file() print "Get the last image file name" fd_last = get_file() if fd and fd_last: # we got the filenames OK image_file = os.path.abspath(fd.name) # get full path name of input file # get number of digits at end; strip these off the corename! image_file_corename=os.path.splitext(image_file) # get only the core (no extension) of input file for use in naming output files image_file_last = os.path.abspath(fd_last.name) # get full path name of input file # get number of digits at end; strip these off the corename! image_file_last_corename=os.path.splitext(image_file_last) # get only the core (no extension) of input file for use in naming output files start_num = image_file_corename[0][-num_ext:] # starting number for the files to analyze end_num = image_file_last_corename[0][-num_ext:] # ending number for the files to analyze frame_spacing = int(input("Filenumber increment between images used (skip this many files in between images): ")) # we will skip this many frames in between each read else: print "could not find this image file" exit() n_boards= len(arange(int(start_num),int(end_num)+1,frame_spacing)) # no of images of checkerboards used to calibrate : # creation of memory storage arrays image_points=cv.CreateMat(n_boards*board_n,2,cv.CV_32FC1) object_points=cv.CreateMat(n_boards*board_n,3,cv.CV_32FC1) point_counts=cv.CreateMat(n_boards,1,cv.CV_32SC1) intrinsic_matrix=cv.CreateMat(3,3,cv.CV_32FC1) distortion_coefficient=cv.CreateMat(5,1,cv.CV_32FC1) for i in arange(int(start_num),int(end_num)+1,frame_spacing): # read in the different chessboard views, find the checkerboard corners # add in a part to locate the files to use *** ## image = cv.LoadImage("checkerboard"+str(i).zfill(2)+".jpg",cv.CV_LOAD_IMAGE_GRAYSCALE) # don't delete these three lines--useful if you ever need to debug ## image_corners = cv.LoadImage("checkerboard"+str(i).zfill(2)+".jpg",cv.CV_LOAD_IMAGE_COLOR) image = cv.LoadImage(image_file_corename[0][:-num_ext]+str(i).zfill(num_ext)+image_file_corename[1],cv.CV_LOAD_IMAGE_GRAYSCALE) image_corners = cv.LoadImage(image_file_corename[0][:-num_ext]+str(i).zfill(num_ext)+image_file_corename[1],cv.CV_LOAD_IMAGE_COLOR) (image_width,image_height) = cv.GetSize(image) # figure out how many pixels high and wide the images is--use for image center later ## print "# pixels wide by high",cv.GetSize(image) # don't delete these three lines--useful if you ever need to debug ## cv.NamedWindow("chessboard") # so you can uncomment out and look at the images ## cv.ShowImage("chessboard",image) ## cv.WaitKey(0) # this makes HighGui happy! ## (found,corners)=cv.FindChessboardCorners(image,board_sz,0) # no flags set -- works OK for uniform lighting, not very distorted squares (found,corners)=cv.FindChessboardCorners(image,board_sz,cv.CV_CALIB_CB_ADAPTIVE_THRESH | cv.CV_CALIB_CB_FILTER_QUADS) # *** maybe you want to set some flags? read da book! corners=cv.FindCornerSubPix(image,corners,(win_subpix,win_subpix),(-1,-1),(cv.CV_TERMCRIT_EPS+cv.CV_TERMCRIT_ITER,30,0.1)) # *** set win_subpix large enough for high def images # (11,11) (original code) is a bit large for my images. Each entry is how wide (on either side) a window you # utilize to deterimine the transition from dark to light for each edge. 8 seemed OK to me. Too big # and you are into the next square! if found: # it found chessboard corners OK! corner_count=len(corners) print "found ",corner_count," chessboard corners in file ",i num_usable_images=num_usable_images+1 ## cv.DrawChessboardCorners( image_corners, board_sz, corners, found ) # don't delete these three lines--useful if you ever need to debug ## cv.NamedWindow("chessboard corners") # don't delete these three lines--useful if you ever need to debug ## cv.ShowImage("chessboard corners", image_corners) # so you can uncomment out and look at the images ## cv.WaitKey() # this makes HighGui happy! if len(corners)==board_n: # if it found all the corners, add to matrix step=successes*board_n # this is the starting point in the image_points & object_points arrays to insert the new set of board_n values k = step for j in range(board_n): cv.Set2D(image_points,k,0,corners[j][0]) # this arrange holds the corner's coordinates in the images cv.Set2D(image_points,k,1,corners[j][1]) cv.Set2D(object_points,k,0,square_size*float(j)/float(board_w)) # x = row number (actual physical coordinates of the checkerboard corners) cv.Set2D(object_points,k,1,square_size*float(j)%float(board_w)) # y = column number cv.Set2D(object_points,k,2,0.0) # z = 0.0 for the real world coordinates in the checkerboard frame # if I multiply object_points by a scaling factor, this can give me real world units! k=k+1 cv.Set2D(point_counts,successes,0,board_n) successes=successes+1 ## time.sleep(2) print "-------------------------------------------------" print "\n" else: # why didnt' it find all the corners? take a look! cv.DrawChessboardCorners( image_corners, board_sz, corners, found ) cv.NamedWindow("chessboard corners") cv.ShowImage("chessboard corners", image_corners) else: print "did not find corners in file ",i cv.NamedWindow("chessboard") # what's wrong with this image? take a look! cv.ShowImage("chessboard",image) cv.DestroyWindow("chessboard corners") cv.DestroyWindow("chessboard") print "checkerboard corners finding finished" if successes == 0: print "did not find any correct checkerboard corner counts--exiting" exit() ########################################################################################################### # now assigning new matrices according to view_count object_points2=cv.CreateMat(successes*board_n,3,cv.CV_32FC1) # this second set of arrays is just big enough to hold the results for valid corners image_points2=cv.CreateMat(successes*board_n,2,cv.CV_32FC1) # (the first set of arrays above was big enough for all possible image corner data) point_counts2=cv.CreateMat(successes,1,cv.CV_32SC1) #transfer points to matrices for all the various images for i in range(successes*board_n): cv.Set2D(image_points2,i,0,cv.Get2D(image_points,i,0)) # the Learning OpenCV book says use CV_MAT_ELEM but this uses cv.Get2D--works the same cv.Set2D(image_points2,i,1,cv.Get2D(image_points,i,1)) cv.Set2D(object_points2,i,0,cv.Get2D(object_points,i,0)) cv.Set2D(object_points2,i,1,cv.Get2D(object_points,i,1)) cv.Set2D(object_points2,i,2,cv.Get2D(object_points,i,2)) for i in range(successes): cv.Set2D(point_counts2,i,0,cv.Get2D(point_counts,i,0)) cv.Set2D(intrinsic_matrix,0,0,f_x) # Initialize the intrinsic matrix so the center is centered & the focal lengths have a value of 1.0 <-- change ??? *** cv.Set2D(intrinsic_matrix,0,1,0) # Initialize the intrinsic matrix so the center is centered & the focal lengths have a value of 1.0 <-- change ??? *** cv.Set2D(intrinsic_matrix,0,2,int(float(image_width)/2.)) cv.Set2D(intrinsic_matrix,1,0,0) cv.Set2D(intrinsic_matrix,1,1,f_y) cv.Set2D(intrinsic_matrix,1,2,int(float(image_height)/2.)) cv.Set2D(intrinsic_matrix,2,0,0) cv.Set2D(intrinsic_matrix,2,1,0) cv.Set2D(intrinsic_matrix,2,2,1) rcv = cv.CreateMat(successes, 3, cv.CV_64FC1) # Initialize the rotation vectors and translation vectors that define the extrinsic matrix tcv = cv.CreateMat(successes, 3, cv.CV_64FC1) # in the python code I found, this was originally first dimension n_boards, but that's wrong--we only want it to be successes print "checking camera calibration............." if fit_only_k1_k2: print "fixing focal length and image center, fitting only k1 and k2" cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_PRINCIPAL_POINT|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_ZERO_TANGENT_DIST|cv.CV_CALIB_FIX_K3) else: # *** since I only use high quality cameras, I always fix the focal lengths fx = fy by setting flag cv.CV_CALIB_FIX_ASPECT_RATIO cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_FIX_PRINCIPAL_POINT|cv.CV_CALIB_ZERO_TANGENT_DIST|cv.CV_CALIB_FIX_K1|cv.CV_CALIB_FIX_K2|cv.CV_CALIB_FIX_K3) # first solve only for the focal length cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_ZERO_TANGENT_DIST|cv.CV_CALIB_FIX_K1|cv.CV_CALIB_FIX_K2|cv.CV_CALIB_FIX_K3) # next, solve for both the focal length and the image center pixels if fish_eye_lens and tangent_undistort: # poor quality camera--use all possible corrections (k1,k2,k3 and tangential terms != 0) cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO) elif fish_eye_lens: # allows k1, k2 and k3 != 0 cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_ZERO_TANGENT_DIST) elif first_order_in_k: # no tangential terms, only k1 != 0, k2,k3 = 0 cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_ZERO_TANGENT_DIST|cv.CV_CALIB_FIX_K2|cv.CV_CALIB_FIX_K3) elif fit_only_k1_k2: print "not fitting to the distortion model" else: # currently, by default (at the book's suggestion) k1, k2 != 0, k3, tangential terms = 0 cv.CalibrateCamera2(object_points2,image_points2,point_counts2,cv.GetSize(image),intrinsic_matrix,distortion_coefficient,rcv,tcv,cv.CV_CALIB_USE_INTRINSIC_GUESS|cv.CV_CALIB_FIX_ASPECT_RATIO|cv.CV_CALIB_ZERO_TANGENT_DIST|cv.CV_CALIB_FIX_K3) # *** Choose your flags wisely--read the book! This part is important: there are many, many possible parameters. Creeping up (maybe fixing f if you know it, fixing # *** the centers, finding k1, k2, or maybe zeroing out the k's and finding f and the center are a good strategy. print "camera calibration went OK " # storing results in xml files cv.Save("Intrinsics.xml",intrinsic_matrix) cv.Save("Distortion.xml",distortion_coefficient) # Loading from xml files intrinsic = cv.Load("Intrinsics.xml") distortion = cv.Load("Distortion.xml") print "saved & reloaded all distortion parameters" mapx = cv.CreateImage( cv.GetSize(image), cv.IPL_DEPTH_32F, 1 ) # initialize the undistortion maps used to correct for radial and tangential camera distortions mapy = cv.CreateImage( cv.GetSize(image), cv.IPL_DEPTH_32F, 1 ) cv.InitUndistortMap(intrinsic,distortion,mapx,mapy) # build the undistortion maps print "distortion mapping completed" ########################################################################################################### print "summary so far: # usable images (corners found) = ",num_usable_images," # images with correct # corners = ",successes print "focal length (fx,fy) in pixels = ",intrinsic_matrix[0,0],intrinsic_matrix[1,1] print "center of images (x,y) in pixels = ",intrinsic_matrix[0,2],intrinsic_matrix[1,2] print "distortion parameters:" distort_params = ['k_1', 'k_2', 'p_1', 'p_2', 'k_3'] for j in arange(5): print distort_params[j],' ',distortion_coefficient[j,0] print "image size (pixels) = ",(image_width,image_height)," dead-center would be (pixels) = ",(float(image_width)/2.,float(image_height)/2.) ########################################################################################################### print "Now undistorting the original images." print "The three windows show the original, the attempt at undistoring it, and their pixel-by-pixel difference" # *** if you are confident your distortion code is doing its job, you can comment out the rest. cv.NamedWindow( "Original calibraction image" ) cv.NamedWindow( "Undistorted image" ) for i in arange(int(start_num),int(end_num)+1,frame_spacing): # loop over all input images and apply undistortion print "working on image ",image_file_corename[0][:-num_ext]+str(i).zfill(num_ext)+image_file_corename[1] image=cv.LoadImage(image_file_corename[0][:-num_ext]+str(i).zfill(num_ext)+image_file_corename[1],cv.CV_LOAD_IMAGE_GRAYSCALE) orig_image = cv.CloneImage(image) cv.ShowImage( "Original calibraction image", orig_image) # show original image (before undistortion correction) cv.Remap( orig_image, image, mapx, mapy ) cv.ShowImage("Undistorted image", image) # show remapped undistorted version image_diff = cv.CloneImage(image) # now show the image difference--is it really undistorting? cv.AbsDiff( orig_image, image, image_diff ); cv.NamedWindow( "original minus undistort") cv.ShowImage( "original minus undistort", image_diff ) c = cv.WaitKey(0) # this makes HighGui happy!