魚眼鏡頭標定及畸變校正

龍嘯wyh發表於2020-10-03

魚眼攝像頭畸變校正的方法:

  1. 棋盤矯正法
  2. 經緯度矯正法。

相機為什麼會出現畸變?
當前相機的畸變主要分為徑向畸變和切向畸變兩種。
徑向畸變產生的原因:相機的光學鏡頭厚度不均勻,離鏡頭越遠場景的光線就越彎曲從而產生徑向畸變。
切向畸變產生的原因:鏡頭與影像感測器不完全平行造成的。
在這裡插入圖片描述
在這裡插入圖片描述

相機引數有哪些?
相機內參:主要包括相機矩陣(包括焦距,光學中心,這些都是相機本身屬性)和畸變係數(畸變數學模型的5個引數 D = {K1,K2,K3,P1,P2})。
相機外參:通過旋轉和平移將實際場景3D對映到相機的2D座標過程中的旋轉和平移就是外參。(他描述的是世界座標轉化成相機座標的過程)

相機的標定流程
機標定流程就是4個座標系在轉換過程中求出計算機的外參和內參的過程。四個座標系分別是:世界座標系(真實的實際場景),相機座標系(攝像頭鏡頭中心),影像座標系(影像感測器成像中心,圖片中心,影布中心),畫素座標系(影像左上角為原點)。如圖三所示,O1是影像座標系,O0 是畫素座標系,兩者之間的區別只是原點發生了變化。
世界座標系 ->相機座標系 求解外參(旋轉和平移)
相機座標系 ->影像座標系 求解內參(攝像頭矩陣,畸變係數)
影像座標系 ->畫素座標系 求解畫素轉化矩陣(可簡單理解為原點從圖片中心到左上角,單位釐米變行列)
在這裡插入圖片描述

一、基於棋盤的相機畸變校正方法
列印棋盤並採集魚眼攝像頭下的棋盤圖片:
1.棋盤獲取:連結: https://pan.baidu.com/s/14qB3kQ_MbWORay1i0GFm1A 提取碼: ksqw 複製這段內容後開啟百度網盤手機App,操作更方便哦
2.採集圖片,採集圖片如下圖所示,可以多采集一些,標註的會更加準確。
在這裡插入圖片描述

3.使用採集的圖片求出相機的內參和矯正係數(DIM, K, D),然後使用得到的(DIM, K, D)再進行測試,程式碼如下。

import cv2
import numpy as np
import glob

def get_K_and_D(checkerboard, imgsPath):
 
    CHECKERBOARD = checkerboard
    subpix_criteria = (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 0.1)
    calibration_flags = cv2.fisheye.CALIB_RECOMPUTE_EXTRINSIC+cv2.fisheye.CALIB_CHECK_COND+cv2.fisheye.CALIB_FIX_SKEW
    objp = np.zeros((1, CHECKERBOARD[0]*CHECKERBOARD[1], 3), np.float32)
    objp[0,:,:2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2)
    _img_shape = None
    objpoints = []
    imgpoints = []
    images = glob.glob(imgsPath + '/*.png')
    for fname in images:
        img = cv2.imread(fname)
        if _img_shape == None:
            _img_shape = img.shape[:2]
        else:
            assert _img_shape == img.shape[:2], "All images must share the same size."
 
        gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
        ret, corners = cv2.findChessboardCorners(gray, CHECKERBOARD,cv2.CALIB_CB_ADAPTIVE_THRESH+cv2.CALIB_CB_FAST_CHECK+cv2.CALIB_CB_NORMALIZE_IMAGE)
        if ret == True:
            objpoints.append(objp)
            cv2.cornerSubPix(gray,corners,(3,3),(-1,-1),subpix_criteria)
            imgpoints.append(corners)
    N_OK = len(objpoints)
    K = np.zeros((3, 3))
    D = np.zeros((4, 1))
    rvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
    tvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
    rms, _, _, _, _ = cv2.fisheye.calibrate(
        objpoints,
        imgpoints,
        gray.shape[::-1],
        K,
        D,
        rvecs,
        tvecs,
        calibration_flags,
        (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 1e-6)
    )
    DIM = _img_shape[::-1]
    print("Found " + str(N_OK) + " valid images for calibration")
    print("DIM=" + str(_img_shape[::-1]))
    print("K=np.array(" + str(K.tolist()) + ")")
    print("D=np.array(" + str(D.tolist()) + ")")
    return DIM, K, D
 
 
def undistort(img_path,K,D,DIM,scale=0.6,imshow=False):
    img = cv2.imread(img_path)
    dim1 = img.shape[:2][::-1]  #dim1 is the dimension of input image to un-distort
    assert dim1[0]/dim1[1] == DIM[0]/DIM[1], "Image to undistort needs to have same aspect ratio as the ones used in calibration"
    if dim1[0]!=DIM[0]:
        img = cv2.resize(img,DIM,interpolation=cv2.INTER_AREA)
    Knew = K.copy()
    if scale:#change fov
        Knew[(0,1), (0,1)] = scale * Knew[(0,1), (0,1)]
    map1, map2 = cv2.fisheye.initUndistortRectifyMap(K, D, np.eye(3), Knew, DIM, cv2.CV_16SC2)
    undistorted_img = cv2.remap(img, map1, map2, interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT)
    if imshow:
        cv2.imshow("undistorted", undistorted_img)
    return undistorted_img
 
if __name__ == '__main__':
 
   # 開始使用圖片來獲取內參和畸變係數
    DIM, K, D = get_K_and_D((6,9), '')
 
   # 得到內參和畸變係數畸變矯正進行測試
    '''
    DIM=(2560, 1920)
    K=np.array([[652.8609862494474, 0.0, 1262.1021584894233], [0.0, 653.1909758659955, 928.0871455436396], [0.0, 0.0, 1.0]])
    D=np.array([[-0.024092199861108887], [0.002745976275100771], [0.002545415522352827], [-0.0014366825722748522]])
    img = undistort('../imgs/pig.jpg',K,D,DIM)
    cv2.imwrite('../imgs/pig_checkerboard.jpg', img)
    '''

二、基於經緯度的矯正方法
1.演算法原理:經緯度矯正法, 可以把魚眼圖想象成半個地球, 然後將地球展開成地圖,經緯度矯正法主要是利用幾何原理, 對影像進行展開矯正。(此演算法操作簡單不需要使用棋盤來進行資料採集)
2.程式碼實現

import cv2
import numpy as np
import time
 
# 魚眼有效區域擷取
def cut(img):
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    (_, thresh) = cv2.threshold(img_gray, 20, 255, cv2.THRESH_BINARY)
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    cnts = sorted(contours, key=cv2.contourArea, reverse=True)[0]
    x,y,w,h = cv2.boundingRect(cnts)
    r = max(w/ 2, h/ 2)
    # 提取有效區域
    img_valid = img[y:y+h, x:x+w]
    return img_valid, int(r)
 
# 魚眼矯正
def undistort(src,r):
    # r: 半徑, R: 直徑
    R = 2*r
    # Pi: 圓周率
    Pi = np.pi
    # 儲存對映結果
    dst = np.zeros((R, R, 3))
    src_h, src_w, _ = src.shape
 
    # 圓心
    x0, y0 = src_w//2, src_h//2
 
    # 陣列, 迴圈每個點
    range_arr = np.array([range(R)])
 
    theta = Pi - (Pi/R)*(range_arr.T)
    temp_theta = np.tan(theta)**2
 
    phi = Pi - (Pi/R)*range_arr
    temp_phi = np.tan(phi)**2
 
    tempu = r/(temp_phi + 1 + temp_phi/temp_theta)**0.5
    tempv = r/(temp_theta + 1 + temp_theta/temp_phi)**0.5
 
    # 用於修正正負號
    flag = np.array([-1] * r + [1] * r)
 
    # 加0.5是為了四捨五入求最近點
    u = x0 + tempu * flag + 0.5
    v = y0 + tempv * np.array([flag]).T + 0.5
 
    # 防止陣列溢位
    u[u<0]=0
    u[u>(src_w-1)] = src_w-1
    v[v<0]=0
    v[v>(src_h-1)] = src_h-1
 
    # 插值
    dst[:, :, :] = src[v.astype(int),u.astype(int)]
    return dst
 
if __name__ == "__main__":
    t = time.perf_counter()
    frame = cv2.imread('../imgs/pig.jpg')
    cut_img,R = cut(frame)
    t = time.perf_counter()
    result_img = undistort(cut_img,R)
    cv2.imwrite('../imgs/pig_vector_nearest.jpg',result_img)
    print(time.perf_counter()-t)

影像演算法中會經常用到攝像機的畸變校正,有必要總結分析OpenCV中畸變校正方法,其中包括普通針孔相機模型和魚眼相機模型fisheye兩種畸變校正方法。
普通相機模型畸變校正函式針對OpenCV中的cv::initUndistortRectifyMap(),魚眼相機模型畸變校正函式對應OpenCV中的cv::fisheye::initUndistortRectifyMap()。兩種方法算出對映Mapx和Mapy後,統一用cv::Remap()函式進行插值得到校正後的影像。
1. FishEye模型的畸變校正。
方便起見,直接貼出OpenCV原始碼,我在裡面加了註釋說明。建議參考OpenCV官方文件看畸變模型原理會更清楚:https://docs.opencv.org/3.0-beta/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html#fisheye-initundistortrectifymap

簡要流程就是:
1.求內參矩陣的逆,由於攝像機座標系的三維點到二維影像平面,需要乘以旋轉矩陣R和內參矩陣K。那麼反向投影回去則是二維影像座標乘以K*R的逆矩陣。
2.將目標影像中的每一個畫素點座標(j,i),乘以1中求出的逆矩陣iR,轉換到攝像機座標系(_x,_y,_w),並歸一化得到z=1平面下的三維座標(x,y,1);
3.求出平面模型下畫素點對應魚眼半球模型下的極座標(r, theta)。
4.利用魚眼畸變模型求出擁有畸變時畫素點對應的theta_d。
在這裡插入圖片描述

5.利用求出的theta_d值將三維座標點重投影到二維影像平面得到(u,v),(u,v)即為目標影像對應的畸變影像中畫素點座標
6.使用cv::Remap()函式,根據mapx,mapy取出對應座標位置的畫素值賦值給目標影像,一般採用雙線性插值法,得到畸變校正後的目標影像。

#include <opencv2\opencv.hpp>

void cv::fisheye::initUndistortRectifyMap( InputArray K, InputArray D, InputArray R, InputArray P,
    const cv::Size& size, int m1type, OutputArray map1, OutputArray map2 )
{
    CV_Assert( m1type == CV_16SC2 || m1type == CV_32F || m1type <=0 );
    map1.create( size, m1type <= 0 ? CV_16SC2 : m1type );
    map2.create( size, map1.type() == CV_16SC2 ? CV_16UC1 : CV_32F );

    CV_Assert((K.depth() == CV_32F || K.depth() == CV_64F) && (D.depth() == CV_32F || D.depth() == CV_64F));
    CV_Assert((P.empty() || P.depth() == CV_32F || P.depth() == CV_64F) && (R.empty() || R.depth() == CV_32F || R.depth() == CV_64F));
    CV_Assert(K.size() == Size(3, 3) && (D.empty() || D.total() == 4));
    CV_Assert(R.empty() || R.size() == Size(3, 3) || R.total() * R.channels() == 3);
    CV_Assert(P.empty() || P.size() == Size(3, 3) || P.size() == Size(4, 3));

    //從內參矩陣K中取出歸一化焦距fx,fy; cx,cy
    cv::Vec2d f, c;
    if (K.depth() == CV_32F)
    {
        Matx33f camMat = K.getMat();
        f = Vec2f(camMat(0, 0), camMat(1, 1));
        c = Vec2f(camMat(0, 2), camMat(1, 2));
    }
    else
    {
        Matx33d camMat = K.getMat();
        f = Vec2d(camMat(0, 0), camMat(1, 1));
        c = Vec2d(camMat(0, 2), camMat(1, 2));
    }
    //從畸變係數矩陣D中取出畸變係數k1,k2,k3,k4
    Vec4d k = Vec4d::all(0);
    if (!D.empty())
        k = D.depth() == CV_32F ? (Vec4d)*D.getMat().ptr<Vec4f>(): *D.getMat().ptr<Vec4d>();

    //旋轉矩陣RR轉換資料型別為CV_64F,如果不需要旋轉,則RR為單位陣
    cv::Matx33d RR  = cv::Matx33d::eye();
    if (!R.empty() && R.total() * R.channels() == 3)
    {
        cv::Vec3d rvec;
        R.getMat().convertTo(rvec, CV_64F);
        RR = Affine3d(rvec).rotation();
    }
    else if (!R.empty() && R.size() == Size(3, 3))
        R.getMat().convertTo(RR, CV_64F);
    
    //新的內參矩陣PP轉換資料型別為CV_64F
    cv::Matx33d PP = cv::Matx33d::eye();
    if (!P.empty())
        P.getMat().colRange(0, 3).convertTo(PP, CV_64F);

    //關鍵一步:新的內參矩陣*旋轉矩陣,然後利用SVD分解求出逆矩陣iR,後面用到
    cv::Matx33d iR = (PP * RR).inv(cv::DECOMP_SVD);

    //反向對映,遍歷目標影像所有畫素位置,找到畸變影像中對應位置座標(u,v),並分別儲存座標(u,v)到mapx和mapy中
    for( int i = 0; i < size.height; ++i)
    {
        float* m1f = map1.getMat().ptr<float>(i);
        float* m2f = map2.getMat().ptr<float>(i);
        short*  m1 = (short*)m1f;
        ushort* m2 = (ushort*)m2f;

        //二維影像平面座標系->攝像機座標系
        double _x = i*iR(0, 1) + iR(0, 2),
               _y = i*iR(1, 1) + iR(1, 2),
               _w = i*iR(2, 1) + iR(2, 2);

        for( int j = 0; j < size.width; ++j)
        {
            //歸一化攝像機座標系,相當於假定在Z=1平面上
            double x = _x/_w, y = _y/_w;

            //求魚眼半球體截面半徑r
            double r = sqrt(x*x + y*y);
            //求魚眼半球面上一點與光心的連線和光軸的夾角Theta
            double theta = atan(r);
            //畸變模型求出theta_d,相當於有畸變的角度值
            double theta2 = theta*theta, theta4 = theta2*theta2, theta6 = theta4*theta2, theta8 = theta4*theta4;
            double theta_d = theta * (1 + k[0]*theta2 + k[1]*theta4 + k[2]*theta6 + k[3]*theta8);
            //利用有畸變的Theta值,將攝像機座標系下的歸一化三維座標,重投影到二維影像平面,得到(j,i)對應畸變影像中的(u,v)
            double scale = (r == 0) ? 1.0 : theta_d / r;
            double u = f[0]*x*scale + c[0];
            double v = f[1]*y*scale + c[1];

            //儲存(u,v)座標到mapx,mapy
            if( m1type == CV_16SC2 )
            {
                int iu = cv::saturate_cast<int>(u*cv::INTER_TAB_SIZE);
                int iv = cv::saturate_cast<int>(v*cv::INTER_TAB_SIZE);
                m1[j*2+0] = (short)(iu >> cv::INTER_BITS);
                m1[j*2+1] = (short)(iv >> cv::INTER_BITS);
                m2[j] = (ushort)((iv & (cv::INTER_TAB_SIZE-1))*cv::INTER_TAB_SIZE + (iu & (cv::INTER_TAB_SIZE-1)));
            }
            else if( m1type == CV_32FC1 )
            {
                m1f[j] = (float)u;
                m2f[j] = (float)v;
            }

            //這三條語句是上面 ”//二維影像平面座標系->攝像機座標系“的一部分,是矩陣iR的第一列,這樣寫能夠簡化計算
            _x += iR(0, 0);
            _y += iR(1, 0);
            _w += iR(2, 0);
        }
    }
}

opencv-python 實現魚眼矯正 棋盤矯正法
魚眼矯正有很多的方法, 比較常用的有:
1、棋盤標定法
2、經緯度法
opencv自帶魚眼矯正演算法, 也就是第一種, 棋盤矯正法。

第一步:製作棋盤格
用A4紙列印一張棋盤格, 固定到硬紙板上, 然後用魚眼鏡頭對著拍攝。 保留拍到的圖片, 如下圖所示:
在這裡插入圖片描述

可以從不同角度拍攝, 多儲存一些。

第二步: 計算內參和矯正係數
棋盤標定法, 必須要先計算出魚眼的內參和矯正係數, 可直接呼叫以下函式計算

import cv2
assert cv2.__version__[0] == '3'
import numpy as np
import os
import glob

def get_K_and_D(checkerboard, imgsPath):
    CHECKERBOARD = checkerboard
    subpix_criteria = (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 0.1)
    calibration_flags = cv2.fisheye.CALIB_RECOMPUTE_EXTRINSIC+cv2.fisheye.CALIB_CHECK_COND+cv2.fisheye.CALIB_FIX_SKEW
    objp = np.zeros((1, CHECKERBOARD[0]*CHECKERBOARD[1], 3), np.float32)
    objp[0,:,:2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2)
    _img_shape = None
    objpoints = [] 
    imgpoints = [] 
    images = glob.glob(imgsPath + '/*.png')
    for fname in images:
        img = cv2.imread(fname)
        if _img_shape == None:
            _img_shape = img.shape[:2]
        else:
            assert _img_shape == img.shape[:2], "All images must share the same size."
        
        gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
        ret, corners = cv2.findChessboardCorners(gray, CHECKERBOARD,cv2.CALIB_CB_ADAPTIVE_THRESH+cv2.CALIB_CB_FAST_CHECK+cv2.CALIB_CB_NORMALIZE_IMAGE)
        if ret == True:
            objpoints.append(objp)
            cv2.cornerSubPix(gray,corners,(3,3),(-1,-1),subpix_criteria)
            imgpoints.append(corners)
    N_OK = len(objpoints)
    K = np.zeros((3, 3))
    D = np.zeros((4, 1))
    rvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
    tvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
    rms, _, _, _, _ = \
    cv2.fisheye.calibrate(
                                objpoints,
                                imgpoints,
                                gray.shape[::-1],
                                K,
                                D,
                                rvecs,
                                tvecs,
                                calibration_flags,
                                (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 1e-6)
                                )
    DIM = _img_shape[::-1]
    print("Found " + str(N_OK) + " valid images for calibration")
    print("DIM=" + str(_img_shape[::-1]))
    print("K=np.array(" + str(K.tolist()) + ")")
    print("D=np.array(" + str(D.tolist()) + ")")
    
    return DIM, K, D

# 計算內參和矯正係數
'''
# checkerboard: 棋盤格的格點數目
# imgsPath: 存放魚眼圖片的路徑
'''
get_K_and_D((6,9), 'fisheyeImage')

第三步: 根據計算的K和D, 矯正魚眼圖
儲存第二步計算的DIM, K, D, 利用一下程式, 矯正魚眼圖

def undistort(img_path):
    img = cv2.imread(img_path)
    img = cv2.resize(img, DIM)
    map1, map2 = cv2.fisheye.initUndistortRectifyMap(K, D, np.eye(3), K, DIM,cv2.CV_16SC2)
    undistorted_img = cv2.remap(img, map1, map2, interpolation=cv2.INTER_LINEAR,borderMode=cv2.BORDER_CONSTANT)    
    cv2.imwrite('unfisheyeImage.png', undistorted_img)

DIM, K, D是固定不變的,因此,map1和map2也是不變的, 當你有大量的資料需要矯正時, 應當避免map1和map2的重複計算, 只需要計算remap即可。

矯正效果
在這裡插入圖片描述

後續
矯正結束, 對於上面的矯正結果, 應該適用於很多場景了, 但是矯正之後, 很多有效區域被截掉了, 能否矯正出更大的有效面積呢? 答案當然是有的, 具體的細節, 去github上看專案裡面的pdf吧, 這裡不再介紹了。
https://github.com/HLearning/fisheye

相關文章