采用的方法為最小二乘法：

首先我們要構(gòu)建以下方程：

我們討論角點(diǎn)的情況：

q是我們要求的角點(diǎn)

p0和p1為q周?chē)狞c(diǎn)

（q-pi）為一個(gè)向量

Gi為pi處的梯度

所以滿足一下公式

Gi*(q-pi)=0

有以下兩種情況：

（1）p0處的梯度為0,雖然（q-pi）不為0

（2）p1處（q-pi）和p1處的梯度垂直，因此乘積為0.

Gi*(q-pi)=0

我們寫(xiě)成最小二乘的形式：

Gi*q = Gi*pi

根據(jù)最小二乘解：

同理可得：

?

代碼:

// 最大迭代次數(shù)為100次，誤差精度為eps*eps，也就是0.1*0.1。const int MAX_ITERS = 100;int win_w = win.width * 2 + 1, win_h = win.height * 2 + 1;int i, j, k;int max_iters = (criteria.type & CV_TERMCRIT_ITER) ? MIN(MAX(criteria.maxCount, 1), MAX_ITERS) : MAX_ITERS;double eps = (criteria.type & CV_TERMCRIT_EPS) ? MAX(criteria.epsilon, 0.) : 0;eps *= eps; // use square of error in comparsion operations /* 然后是高斯權(quán)重的計(jì)算，如下所示，窗口中心附近權(quán)重高，越往窗口邊界權(quán)重越小。如果設(shè)置的有“零區(qū)域”，則權(quán)重值設(shè)置為0。計(jì)算出的權(quán)重分布如下圖： */Mat maskm(win_h, win_w, CV_32F), subpix_buf(win_h+2, win_w+2, CV_32F);float* mask = maskm.ptr<float>();for( i = 0; i < win_h; i++ ){float y = (float)(i - win.height)/win.height;float vy = std::exp(-y*y);for( j = 0; j < win_w; j++ ){float x = (float)(j - win.width)/win.width;mask[i * win_w + j] = (float)(vy*std::exp(-x*x));}}// make zero_zoneif( zeroZone.width >= 0 && zeroZone.height >= 0 &&zeroZone.width * 2 + 1 < win_w && zeroZone.height * 2 + 1 < win_h ){for( i = win.height - zeroZone.height; i <= win.height + zeroZone.height; i++ ){for( j = win.width - zeroZone.width; j <= win.width + zeroZone.width; j++ ){mask[i * win_w + j] = 0;}}}/* ① 代碼中CI2為本次迭代獲取的亞像素角點(diǎn)位置，CI為上次迭代獲取的亞像素角點(diǎn)位置，CT是初始的整數(shù)角點(diǎn)位置。② 每次迭代結(jié)束計(jì)算CI與CI2之間的歐式距離err，如果兩者之間的歐式距離err小于設(shè)定的閾值，或者迭代次數(shù)達(dá)到設(shè)定的閾值，則停止迭代。③停止迭代后，需要再次判斷最終的亞像素角點(diǎn)位置和初始整數(shù)角點(diǎn)之間的差異，如果差值大于設(shè)定窗口尺寸的一半，則說(shuō)明最小二乘計(jì)算中收斂性不好，丟棄計(jì)算得到的亞像素角點(diǎn)，仍然使用初始的整數(shù)角點(diǎn)。 */// do optimization loop for all the pointsfor( int pt_i = 0; pt_i < count; pt_i++ ){Point2f cT = corners[pt_i], cI = cT;int iter = 0;double err = 0;do{Point2f cI2;double a = 0, b = 0, c = 0, bb1 = 0, bb2 = 0;getRectSubPix(src, Size(win_w+2, win_h+2), cI, subpix_buf, subpix_buf.type());const float* subpix = &subpix_buf.at<float>(1,1);// process gradientfor( i = 0, k = 0; i < win_h; i++, subpix += win_w + 2 ){double py = i - win.height;for( j = 0; j < win_w; j++, k++ ){double m = mask[k];double tgx = subpix[j+1] - subpix[j-1];double tgy = subpix[j+win_w+2] - subpix[j-win_w-2];double gxx = tgx * tgx * m;double gxy = tgx * tgy * m;double gyy = tgy * tgy * m;double px = j - win.width;a += gxx;b += gxy;c += gyy;bb1 += gxx * px + gxy * py;bb2 += gxy * px + gyy * py;}}double det=a*c-b*b;if( fabs( det ) <= DBL_EPSILON*DBL_EPSILON )break;// 2x2 matrix inversiondouble scale=1.0/det;cI2.x = (float)(cI.x + c*scale*bb1 - b*scale*bb2);cI2.y = (float)(cI.y - b*scale*bb1 + a*scale*bb2);err = (cI2.x - cI.x) * (cI2.x - cI.x) + (cI2.y - cI.y) * (cI2.y - cI.y);cI = cI2;if( cI.x < 0 || cI.x >= src.cols || cI.y < 0 || cI.y >= src.rows )break;}while( ++iter < max_iters && err > eps );// if new point is too far from initial, it means poor convergence.// leave initial point as the resultif( fabs( cI.x - cT.x ) > win.width || fabs( cI.y - cT.y ) > win.height )cI = cT;corners[pt_i] = cI;}

?

?

總結(jié)

以上是生活随笔為你收集整理的OpenCV中cornerSubPixel（）亚像素求精原理的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。