照相機標定——棋盤格

照相機標記步驟：

1.制作棋盤格（每個格子的大小可測量），最好是打印出來，貼在平面上（實在不行就用我剛才演示的黑白棋盤格方法）

2.根據棋盤格，采集10-20張圖片，提取角點

3.解算出內外參數，內參截圖放在博客中，外部參數最好能可視化

一、棋盤格選定:

規格:10cm×10cm（5×5）  每個方格2cm×2cm

手機型號:vivoy67  相機分辨率:1300×1600 

二、棋盤格的Harris角點檢測

     首先Harris 角點檢測是基於圖像像素灰度值變化梯度的，灰度值圖像的角點附近，是其像素灰度值變化非常大的區域，其梯度也非常大。換句話說，在非角點位置鄰域里，各點的像素值變化不大，甚至幾乎相等，其梯度相對也比較小。如果鄰域內點的灰度值與中心點Image (i,j) 的灰度值之差 的絕對值 在一個閾值t 范圍 內，那就認為這個點與中心點是相似的。與此同時，屬於該 Image(i,j) 點的相似點計數器nlike(i,j) 也隨之加1。在 Image (i,j) 點的n 鄰域 全部被遍歷一邊之后，就能得到在這個鄰域范圍內與中心點相似的 點個數的統計值nlike(i,j) 。根據nlike(i,j) 的大小，就可以判斷這個中心點是否可能為角點。由於選擇3*3的檢測窗口，所以，對於中心像素點 , 在下面的討論中只考慮其8 鄰域內像素點的相似度。算該范圍的像素點與中心像素點的灰度值之差的絕對值 (記為 Δ ) , 如果該值小於等於設定的閾值 ( 記為 t) , 則認為該像素點與目標像素點相似。

nlike(i,j)=sum(R(i+x,j+y))(-1≤x≤1,-1≤y≤1,且 x≠0 , y≠0) ,其中 : R(i+x, j+y)=1 , 當Δ(i+x , j+y)≤t

        R(i+x,j+y)=0 , 當Δ(i+x , j+y)>t

 從定義中可以看出 : 0≤nlike (i,j)≤8。現在討論nlike( i , j) 值的含義 。

   (1) nlike (i , j) = 8 ,表示當前中心像素點的 8鄰域范圍內都是與之相似的像素點 , 所以該像素點鄰域范圍內的梯度不會很大 , 因此角點檢測時 , 應該排除此類像素點，不將其作為角點的候選點 。

   (2) nlike (i , j) = 0 ,表示當前中心像素點的 8鄰域范圍內沒有與之相似的像素點 , 所以該像素點為孤立像素點或者是噪聲點 , 因此角點檢測時 , 也應該排除此類像素點 。

   (3) nlike (i , j) = 7 ,可以歸結為以下的兩者情況 ,其他情形都可以通過旋轉來得到 ( 圖中黑色區域僅表示與中心像素相似 , 而兩個黑色區域像素可能是相似的 , 也可能不相似 ) 

   (4) nlike (i , j) = 1 , 中心像素點也不可能為角點 。　 

   (5) 2≤ nlike ( i , j) ≤ 6 , 情況比較復雜 , 無法確認像素點准確的性質。我采取的方法是先將其列入候選角點之列，對其進行計算角點響應CRF等后續操作。

   然后計算各候選角點的CBF，進行非極大抑制，再次由於Harris算子不具有尺度不變性，為確保角點檢測的准確度和完整性，我們采用Mikolajczyk和Schmid提出的Harris - Laplace檢測方法[8]，將Harris算子與高斯尺度空間相結合，從而克服Harris 算子只能在單一尺度下檢測角點的缺點

代碼:

#include <iostream>
#include <opencv2/opencv.hpp>
int main(int argc, char *argv[])
{
    int blockSize = 2;
    int kSize = 3;
    double k = 0.04;
    double thread_value = 85;  // R的閾值
    
    cv::Mat  orignal_image = cv::imread (argv[1]);     
    cv::Mat gray_image , after_harris_image;
    cv::Mat norm_image; //歸一化后的圖
    cv::Mat scaled_image; //線性變換后的八位無符號整形的圖
    cv::cvtColor (orignal_image,gray_image,CV_BGR2GRAY);     // 灰度變換
    cv::cornerHarris(gray_image, after_harris_image, blockSize, kSize, k, cv::BORDER_DEFAULT); //harris corner detect
    cv::normalize(after_harris_image, norm_image, 0, 255, cv::NORM_MINMAX, CV_32FC1, cv::Mat()); //歸一化處理
    convertScaleAbs(norm_image, scaled_image);
    for (int i = 0; i < norm_image.rows; i++)
    {
        for (int j = 0; j < norm_image.cols; j++)
        {
            if ((int)norm_image.at<float>(i, j) > thread_value)
            {
                cv::circle(orignal_image, cv::Point(j, i), 10, cv::Scalar(10, 10, 255), 1, 8, 0);  //tag the corner
            }
        }
    }
    cv::imshow("after_harris_corner", orignal_image);
    cv::imwrite("after_harris_corner.png",orignal_image);
    cv::waitKey(0);

    return 0;
}

檢測結果:

 

小結:

    基於圖像灰度的角點檢測和圖像邊緣提取的角點檢測，其中后一種依賴於圖像邊緣提取的好壞，實用性不是很強。其中Harris的定位性能和魯棒性較好，加之Harris算法原理簡單，不受攝像機姿態和光照的影響，現已成為了使用最為廣泛的角點提取算法之一。因此，針對Harris算法現在也已提出了許多改進算法：亞像素級的角點檢測算法與尺度空間相結合的算法，自適應角點提取算法根據棋盤格的幾何特性提出的檢測算法等，但大部分都是以增加計算量或者低效率為代價，同時，用 Harris 算法進行檢測，有三點不足：（1 ）該算法不具有尺度不變性；（2 ）該算法提取的角點是像素級的；（3 ）該算法檢測時間不是很令人滿意。

三、求解內參數矩陣:

單應性：在計算機視覺中被定義為一個平面到另一個平面的投影映射。首先確定，圖像平面與標定物棋盤格平面的單應性。設三維世界坐標的點為二維相機平面像素坐標為，所以標定用的棋盤格平面到圖像平面的單應性關系為:(其中，K為相機的內參矩陣，R為外部參數矩陣（旋轉矩陣），T為平移向量。令，設棋盤格位於Z=0的平面，定義旋轉矩陣R的第i列為 ri, 則有：

於是空間到圖像的映射可改為：H=λK[r1 r2 t]

其中H 是描述Homographic矩陣,可通過最小二乘，從角點世界坐標到圖像坐標的關系求解

根據步驟1中的式子，令H為H =[h1 h2 h3]，則[h1 h2 h3]=λK[r1 r2 t],再根據正交和歸一化的約束可以得到等式:

即每個單應性矩陣能提供兩個方程，而內參數矩陣包含5個參數，要求解，至少需要3個單應性矩陣。為了得到三個不同的單應性矩陣，我們使用至少三幅棋盤格平面的圖片進行標定。通過改變相機與標定板之間的相對位置來得到三個不同的圖片。為了方便計算，我們定義:

 

B中的未知量可表示為6D向量b,設H中的第列為h;根據b的定義，可以推導出公式

最后推到出通過上式，我們可知當觀測平面n≥3時，即至少3幅棋盤格圖像，可以得到b的唯一解，求得相機內參數矩陣K。

 3、計算外參數矩陣:

外部參數可通過Homography求解，由 H = [h1 h2 h3] = λA[r1 r2 t]，可推出

代碼:

# coding=utf-8

import numpy as np

import cv2

import glob

# 終止標准

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

#w = 6

#h = 6

#准備對象點，如（0,0,0），（1,0,0），（2,0,0）......，（6,5,0）

#objp = np.zeros((w*h,3), np.float32)

objp = np.zeros((4*4,3), np.float32)

objp[:,:2] = np.mgrid[0:4,0:4].T.reshape(-1,2)

# 用於存儲所有圖像中的對象點和圖像點的數組。

objpoints = [] # 在現實世界空間的3d點

imgpoints = [] # 圖像平面中的2d點。

#glob是個文件名管理工具

images = glob.glob('C:/B1/*.jpg')

print('...loading')

for fname in images:

    #對每張圖片，識別出角點，記錄世界物體坐標和圖像坐標

    print('processing img:{fname}')

    img = cv2.imread(fname)

    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)#轉灰度

    print('grayed')

    #尋找角點，存入corners，ret是找到角點的flag

    #ret, corners = cv2.findChessboardCorners(gray, (w, h),None)

    ret, corners = cv2.findChessboardCorners(gray, (4, 4),None)

    # 如果找到，添加對象點，圖像點（精煉后）

    if ret == True:

        print('chessboard detected')

        objpoints.append(objp)

        #執行亞像素級角點檢測

        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)

        imgpoints.append(corners2)

        # 繪制並顯示角點

        #img = cv2.drawChessboardCorners(img, (w,h), corners2,ret)

        img = cv2.drawChessboardCorners(img, (4,4), corners2,ret)

        cv2.namedWindow('img',0)

        cv2.resizeWindow('img', 500, 500)

        cv2.imshow('img',img)

        cv2.waitKey(500)

        cv2.destroyAllWindows()

'''

傳入所有圖片各自角點的三維、二維坐標，相機標定。

每張圖片都有自己的旋轉和平移矩陣，但是相機內參和畸變系數只有一組。

mtx，相機內參；dist，畸變系數；revcs，旋轉矩陣；tvecs，平移矩陣。

'''

img2 = cv2.imread("C:/B1/13.jpg")

print("type objpoints:{objpoints[0].shape}")

print("type imgpoints:{imgpoints[0].shape}")

ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

h,  w = img2.shape[:2]

'''

優化相機內參（camera matrix），這一步可選。

參數1表示保留所有像素點，同時可能引入黑色像素，

設為0表示盡可能裁剪不想要的像素，這是個scale，0-1都可以取。

'''

newcameramtx, roi=cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))

#糾正畸變

dst = cv2.undistort(img2, mtx, dist, None, newcameramtx)

# 裁剪圖像，輸出糾正畸變以后的圖片

x,y,w,h = roi

dst = dst[y:y+h, x:x+w]

cv2.imwrite('calibresult.png',dst)

#打印我們要求的兩個矩陣參數

print ("newcameramtx外參:\n",newcameramtx)

print ("dist畸變值:\n",dist)

print ("newcameramtx旋轉（向量）外參:\n",rvecs)

print ("dist平移（向量）外參:\n",tvecs)

#計算誤差

tot_error = 0

for i in range(len(objpoints)):

    imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)

    error = cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2)

    tot_error += error

print ("total error: ", tot_error/len(objpoints))

運行結果:

實驗總結:

　  實驗時要保證視場大小與標定板大小相當，或者稍微大一些，以保證在標定板平移時仍然在視場內。先把光圈調到最大進行對焦，對焦時要考慮到全場的清晰程度，達到一種全場的均衡。對焦完成后，把光圈減小，增加曝光，同時補光增加亮度，這樣做是為了增加拍攝時的景深。關於標定板亮度，看灰度直方圖，平均亮度在150可以，在180~200之間為宜，不能出現過曝。  其次是相機水平和對稱問題，根據實驗結果分析。要減小誤差需要讓雙攝像機盡量與水平面平行，且兩相機的y坐標盡量相同，左右攝像機關於中軸線對稱。