圖像去畸變的思路
對于目標(biāo)圖像(無畸變圖像)上的每個像素點，轉(zhuǎn)換到normalize平面，再進行畸變變換，進行投影，得到這個像素點畸變后的位置，然后將這個位置的源圖像（畸變圖像）的像素值作為目標(biāo)圖像該點的像素值。
通常我們得到的原圖是畸變后的圖像(x_distort,y_distort)，要計算畸變之前的真實圖像(x,y)，不是用逆運算，而是計算真實圖像畸變后會投影在哪，對應(yīng)過去。先把原圖像設(shè)置為一個空的圖像，把一個個像素畸變投影過去，找到和畸變后圖像像素點的對應(yīng)關(guān)系
魚眼相機模型
世界坐標(biāo)系內(nèi)的點轉(zhuǎn)換到相機坐標(biāo)系
X_c = RX_w + t
轉(zhuǎn)到歸一化平面
?
對歸一化平面內(nèi)的點，進行畸變變換
?
畸變模型
徑向畸變：桶形或者魚眼
切向畸變：鏡片安裝和成像平面不平行引起的。
undistort
?注意，OpenCV畸變系數(shù)矩陣的默認(rèn)順序 [k1, k2, p1, p2, k3]
#include <opencv2/opencv.hpp>
#include <chrono>
using namespace std;
int main(int argc, char **argv)
{
// 內(nèi)參
double fx = 458.654, fy = 457.296, cx = 367.215, cy = 248.375;
/**內(nèi)參矩陣K
* fx 0 cx
* 0 fy cy
* 0 0 1
*/
// 畸變參數(shù)
double k1 = -0.28340811, k2 = 0.07395907, p1 = 0.00019359, p2 = 1.76187114e-05;
cv::Mat image = cv::imread(argv[1], 0); // 圖像是灰度圖，CV_8UC1
int rows = image.rows, cols = image.cols;
cv::Mat image_undistort = cv::Mat(rows, cols, CV_8UC1); // 方法1去畸變以后的圖
cv::Mat image_undistort2 = cv::Mat(rows, cols, CV_8UC1); // 方法2 OpenCV去畸變以后的圖
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
//! 方法1. 自己寫計算去畸變后圖像的內(nèi)容
for (int v = 0; v < rows; v++)
{
for (int u = 0; u < cols; u++)
{
double x = (u - cx) / fx, y = (v - cy) / fy; //要求解的真實圖，歸一化平面上的坐標(biāo)
double r = sqrt(x * x + y * y);
double x_distorted = x * (1 + k1 * r * r + k2 * r * r * r * r) + 2 * p1 * x * y + p2 * (r * r + 2 * x * x); //畸變后歸一化坐標(biāo)
double y_distorted = y * (1 + k1 * r * r + k2 * r * r * r * r) + p1 * (r * r + 2 * y * y) + 2 * p2 * x * y;
double u_distorted = fx * x_distorted + cx; //畸變后像素坐標(biāo)，即原圖
double v_distorted = fy * y_distorted + cy;
// 投影賦值
if (u_distorted >= 0 && v_distorted >= 0 && u_distorted < cols && v_distorted < rows) //真實圖畸變后仍然在圖上的
{
image_undistort.at<uchar>(v, u) = image.at<uchar>((int)v_distorted, (int)u_distorted);
}
else
{
image_undistort.at<uchar>(v, u) = 0; //這里最好用插值法
}
}
}
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "time = " << time_used.count() << endl;
//! 方法2. OpenCV自帶的undistort函數(shù)，更快速
cv::Mat K = cv::Mat::eye(3, 3, CV_32FC1); //內(nèi)參矩陣
K.at<float>(0, 0) = fx;
K.at<float>(1, 1) = fy;
K.at<float>(0, 2) = cx;
K.at<float>(1, 2) = cy;
cv::Mat distort_coeffs = cv::Mat::zeros(1, 5, CV_32FC1); //畸變系數(shù)矩陣 順序是[k1, k2, p1, p2, k3]
distort_coeffs.at<float>(0, 0) = k1;
distort_coeffs.at<float>(0, 1) = k2;
distort_coeffs.at<float>(0, 2) = p1;
distort_coeffs.at<float>(0, 3) = p2;
cout << "K = " << endl
<< K << endl;
cout << "distort_coeffs = " << endl
<< distort_coeffs << endl;
t1 = chrono::steady_clock::now();
cv::undistort(image, image_undistort2, K, distort_coeffs); //去畸變
t2 = chrono::steady_clock::now();
time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "time = " << time_used.count() << endl;
// 展示去畸變后圖像
cv::imshow("distorted", image);
cv::imshow("undistorted", image_undistort);
cv::imshow("image_undistort2", image_undistort2);
cv::waitKey(0);
return 0;
}
undistortPoints
?
if (_distCoeffs)
{
// compensate tilt distortion
cv::Vec3d vecUntilt = invMatTilt * cv::Vec3d(x, y, 1);
double invProj = vecUntilt(2) ? 1. / vecUntilt(2) : 1;
x0 = x = invProj * vecUntilt(0);
y0 = y = invProj * vecUntilt(1);
double error = 1.7976931348623158e+308;
// compensate distortion iteratively
//注1： 開始循環(huán)迭代去除畸變
for (int j = 0;; j++)
{
//注2：一個是最大迭代次數(shù)條件
if ((criteria.type & cv::TermCriteria::COUNT) && j >= criteria.maxCount)
break;
//注3：一個是迭代最大誤差條件
if ((criteria.type & cv::TermCriteria::EPS) && error < criteria.epsilon)
break;
//注4：在只有k1,k2,k3,p1,p2共5個畸變參數(shù)時，對應(yīng)的k數(shù)組中只有k[0]~k[4]有值，其他都為0
double r2 = x*x + y*y;
double icdist = (1 + ((k[7] * r2 + k[6])*r2 + k[5])*r2) / (1 + ((k[4] * r2 + k[1])*r2 + k[0])*r2);
double deltaX = 2 * k[2] * x*y + k[3] * (r2 + 2 * x*x) + k[8] * r2 + k[9] * r2*r2;
double deltaY = k[2] * (r2 + 2 * y*y) + 2 * k[3] * x*y + k[10] * r2 + k[11] * r2*r2;
//注5：形如式（8）的迭代
x = (x0 - deltaX)*icdist;
y = (y0 - deltaY)*icdist;
if (criteria.type & cv::TermCriteria::EPS)
{
double r4, r6, a1, a2, a3, cdist, icdist2;
double xd, yd, xd0, yd0;
cv::Vec3d vecTilt;
//注6：將第k次計算的無畸變點代入畸變模型，得到于當(dāng)前有畸變的偏差
r2 = x*x + y*y;
r4 = r2*r2;
r6 = r4*r2;
a1 = 2 * x*y;
a2 = r2 + 2 * x*x;
a3 = r2 + 2 * y*y;
cdist = 1 + k[0] * r2 + k[1] * r4 + k[4] * r6;
icdist2 = 1. / (1 + k[5] * r2 + k[6] * r4 + k[7] * r6);
xd0 = x*cdist*icdist2 + k[2] * a1 + k[3] * a2 + k[8] * r2 + k[9] * r4;
yd0 = y*cdist*icdist2 + k[2] * a3 + k[3] * a1 + k[10] * r2 + k[11] * r4;
vecTilt = matTilt*cv::Vec3d(xd0, yd0, 1);
invProj = vecTilt(2) ? 1. / vecTilt(2) : 1;
xd = invProj * vecTilt(0);
yd = invProj * vecTilt(1);
double x_proj = xd*fx + cx;
double y_proj = yd*fy + cy;
error = sqrt(pow(x_proj - u, 2) + pow(y_proj - v, 2));
}
}
}
void cvUndistortPointsInternal( const CvMat* _src, CvMat* _dst, const CvMat* _cameraMatrix,
const CvMat* _distCoeffs,
const CvMat* matR, const CvMat* matP, cv::TermCriteria criteria)
{
// 判斷迭代條件是否有效
CV_Assert(criteria.isValid());
// 定義中間變量--A相機內(nèi)參數(shù)組，和matA共享內(nèi)存；RR-矯正變換數(shù)組，和_RR共享內(nèi)存
// k-畸變系數(shù)數(shù)組
double A[3][3], RR[3][3], k[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
CvMat matA=cvMat(3, 3, CV_64F, A), _Dk;
CvMat _RR=cvMat(3, 3, CV_64F, RR);
cv::Matx33d invMatTilt = cv::Matx33d::eye();
cv::Matx33d matTilt = cv::Matx33d::eye();
// 檢查輸入變量是否有效
CV_Assert( CV_IS_MAT(_src) && CV_IS_MAT(_dst) &&
(_src->rows == 1 || _src->cols == 1) &&
(_dst->rows == 1 || _dst->cols == 1) &&
_src->cols + _src->rows - 1 == _dst->rows + _dst->cols - 1 &&
(CV_MAT_TYPE(_src->type) == CV_32FC2 || CV_MAT_TYPE(_src->type) == CV_64FC2) &&
(CV_MAT_TYPE(_dst->type) == CV_32FC2 || CV_MAT_TYPE(_dst->type) == CV_64FC2));
CV_Assert( CV_IS_MAT(_cameraMatrix) &&
_cameraMatrix->rows == 3 && _cameraMatrix->cols == 3 );
cvConvert( _cameraMatrix, &matA );// _cameraMatrix <--> matA / A
// 判斷輸入的畸變系數(shù)是否有效
if( _distCoeffs )
{
CV_Assert( CV_IS_MAT(_distCoeffs) &&
(_distCoeffs->rows == 1 || _distCoeffs->cols == 1) &&
(_distCoeffs->rows*_distCoeffs->cols == 4 ||
_distCoeffs->rows*_distCoeffs->cols == 5 ||
_distCoeffs->rows*_distCoeffs->cols == 8 ||
_distCoeffs->rows*_distCoeffs->cols == 12 ||
_distCoeffs->rows*_distCoeffs->cols == 14));
_Dk = cvMat( _distCoeffs->rows, _distCoeffs->cols,
CV_MAKETYPE(CV_64F,CV_MAT_CN(_distCoeffs->type)), k);// _Dk和數(shù)組k共享內(nèi)存指針
cvConvert( _distCoeffs, &_Dk );
if (k[12] != 0 || k[13] != 0)
{
cv::detail::computeTiltProjectionMatrix<double>(k[12], k[13], NULL, NULL, NULL, &invMatTilt);
cv::detail::computeTiltProjectionMatrix<double>(k[12], k[13], &matTilt, NULL, NULL);
}
}
if( matR )
{
CV_Assert( CV_IS_MAT(matR) && matR->rows == 3 && matR->cols == 3 );
cvConvert( matR, &_RR );// matR和_RR共享內(nèi)存指針
}
else
cvSetIdentity(&_RR);
if( matP )
{
double PP[3][3];
CvMat _P3x3, _PP=cvMat(3, 3, CV_64F, PP);
CV_Assert( CV_IS_MAT(matP) && matP->rows == 3 && (matP->cols == 3 || matP->cols == 4));
cvConvert( cvGetCols(matP, &_P3x3, 0, 3), &_PP );// _PP和數(shù)組PP共享內(nèi)存指針
cvMatMul( &_PP, &_RR, &_RR );// _RR=_PP*_RR 放在一起計算比較高效
}
const CvPoint2D32f* srcf = (const CvPoint2D32f*)_src->data.ptr;
const CvPoint2D64f* srcd = (const CvPoint2D64f*)_src->data.ptr;
CvPoint2D32f* dstf = (CvPoint2D32f*)_dst->data.ptr;
CvPoint2D64f* dstd = (CvPoint2D64f*)_dst->data.ptr;
int stype = CV_MAT_TYPE(_src->type);
int dtype = CV_MAT_TYPE(_dst->type);
int sstep = _src->rows == 1 ? 1 : _src->step/CV_ELEM_SIZE(stype);
int dstep = _dst->rows == 1 ? 1 : _dst->step/CV_ELEM_SIZE(dtype);
double fx = A[0][0];
double fy = A[1][1];
double ifx = 1./fx;
double ify = 1./fy;
double cx = A[0][2];
double cy = A[1][2];
int n = _src->rows + _src->cols - 1;
// 開始對所有點開始遍歷
for( int i = 0; i < n; i++ )
{
double x, y, x0 = 0, y0 = 0, u, v;
if( stype == CV_32FC2 )
{
x = srcf[i*sstep].x;
y = srcf[i*sstep].y;
}
else
{
x = srcd[i*sstep].x;
y = srcd[i*sstep].y;
}
u = x; v = y;
x = (x - cx)*ifx;//轉(zhuǎn)換到歸一化圖像坐標(biāo)系（含有畸變）
y = (y - cy)*ify;
//進行畸變矯正
if( _distCoeffs ) {
// compensate tilt distortion--該部分系數(shù)用來彌補沙氏鏡頭畸變？？
// 如果不懂也沒管，因為普通鏡頭中沒有這些畸變系數(shù)
cv::Vec3d vecUntilt = invMatTilt * cv::Vec3d(x, y, 1);
double invProj = vecUntilt(2) ? 1./vecUntilt(2) : 1;
x0 = x = invProj * vecUntilt(0);
y0 = y = invProj * vecUntilt(1);
double error = std::numeric_limits<double>::max();// error設(shè)定為系統(tǒng)最大值
// compensate distortion iteratively
// 迭代去除鏡頭畸變
// 迭代公式 x′= (x?2p1 xy?p2 (r^2 + 2x^2))∕( 1 + k1*r^2 + k2*r^4 + k3*r^6)
// y′= (y?2p2 xy?p1 (r^2 + 2y^2))∕( 1 + k1*r^2 + k2*r^4 + k3*r^6)
for( int j = 0; ; j++ )
{
if ((criteria.type & cv::TermCriteria::COUNT) && j >= criteria.maxCount)// 迭代最大次數(shù)為5次
break;
if ((criteria.type & cv::TermCriteria::EPS) && error < criteria.epsilon)// 迭代誤差閾值為0.01
break;
double r2 = x*x + y*y;
double icdist = (1 + ((k[7]*r2 + k[6])*r2 + k[5])*r2)/(1 + ((k[4]*r2 + k[1])*r2 + k[0])*r2);
double deltaX = 2*k[2]*x*y + k[3]*(r2 + 2*x*x)+ k[8]*r2+k[9]*r2*r2;
double deltaY = k[2]*(r2 + 2*y*y) + 2*k[3]*x*y+ k[10]*r2+k[11]*r2*r2;
x = (x0 - deltaX)*icdist;
y = (y0 - deltaY)*icdist;
// 對當(dāng)前迭代的坐標(biāo)加畸變，計算誤差error用于判斷迭代條件
if(criteria.type & cv::TermCriteria::EPS)
{
double r4, r6, a1, a2, a3, cdist, icdist2;
double xd, yd, xd0, yd0;
cv::Vec3d vecTilt;
r2 = x*x + y*y;
r4 = r2*r2;
r6 = r4*r2;
a1 = 2*x*y;
a2 = r2 + 2*x*x;
a3 = r2 + 2*y*y;
cdist = 1 + k[0]*r2 + k[1]*r4 + k[4]*r6;
icdist2 = 1./(1 + k[5]*r2 + k[6]*r4 + k[7]*r6);
xd0 = x*cdist*icdist2 + k[2]*a1 + k[3]*a2 + k[8]*r2+k[9]*r4;
yd0 = y*cdist*icdist2 + k[2]*a3 + k[3]*a1 + k[10]*r2+k[11]*r4;
vecTilt = matTilt*cv::Vec3d(xd0, yd0, 1);
invProj = vecTilt(2) ? 1./vecTilt(2) : 1;
xd = invProj * vecTilt(0);
yd = invProj * vecTilt(1);
double x_proj = xd*fx + cx;
double y_proj = yd*fy + cy;
error = sqrt( pow(x_proj - u, 2) + pow(y_proj - v, 2) );
}
}
}
// 將坐標(biāo)從歸一化圖像坐標(biāo)系轉(zhuǎn)換到成像平面坐標(biāo)系
double xx = RR[0][0]*x + RR[0][1]*y + RR[0][2];
double yy = RR[1][0]*x + RR[1][1]*y + RR[1][2];
double ww = 1./(RR[2][0]*x + RR[2][1]*y + RR[2][2]);
x = xx*ww;
y = yy*ww;
if( dtype == CV_32FC2 )
{
dstf[i*dstep].x = (float)x;
dstf[i*dstep].y = (float)y;
}
else
{
dstd[i*dstep].x = x;
dstd[i*dstep].y = y;
}
}
}
畸變矯正
通過calibrate Camera()得到的內(nèi)參和畸變系數(shù)
圖像的畸變矯正需要相機的內(nèi)參和畸變系數(shù), 在opencv中, 有以下兩個函數(shù)可以實現(xiàn):
- initUndistortRectifyMap() + remap()函數(shù)
- undistort()函數(shù)
推薦使用第一種方法：initUndistortRectifyMap只用運行一次，remap讀取一次圖像運行一次。
void initUndistortRectifyMap( InputArray cameraMatrix, InputArray distCoeffs,
InputArray R, InputArray newCameraMatrix, Size size, int m1type, OutputArray map1, OutputArray map2 );
參數(shù)說明：
cameraMatrix——輸入的攝像頭內(nèi)參數(shù)矩陣（3X3矩陣）
distCoeffs——輸入的攝像頭畸變系數(shù)矩陣（5X1矩陣）
R——輸入的第一和第二攝像頭坐標(biāo)系之間的旋轉(zhuǎn)矩陣
newCameraMatrix——輸入的校正后的3X3攝像機矩陣
size——攝像頭采集的無失真圖像尺寸
m1type——map1的數(shù)據(jù)類型，可以是CV_32FC1或CV_16SC2
map1——輸出的X坐標(biāo)重映射參數(shù)
map2——輸出的Y坐標(biāo)重映射參數(shù)
?
void remap(InputArray src, OutputArray dst, InputArray map1, InputArray map2, int interpolation, int borderMode=BORDER_CONSTANT, const Scalar& borderValue=Scalar())
參數(shù)說明：
src——輸入圖像，即原圖像，需要單通道8位或者浮點類型的圖像
dst(c++)——輸出圖像，即目標(biāo)圖像，需和原圖形一樣的尺寸和類型
map1——它有兩種可能表示的對象：（1）表示點（x,y）的第一個映射；（2）表示CV_16SC2，CV_32FC1等
map2——有兩種可能表示的對象：（1）若map1表示點（x,y）時，這個參數(shù)不代表任何值；（2）表示 CV_16UC1，CV_32FC1類型的Y值
intermap2polation——插值方式，有四中插值方式：
（1）INTER_NEAREST——最近鄰插值
（2）INTER_LINEAR——雙線性插值（默認(rèn)）
（3）INTER_CUBIC——雙三樣條插值（默認(rèn)）
（4）INTER_LANCZOS4——lanczos插值（默認(rèn)）intborderMode——邊界模式，默認(rèn)BORDER_CONSTANT
borderValue——邊界顏色，默認(rèn)Scalar()黑色
?
projectPoints
先將世界坐標(biāo)系內(nèi)的點轉(zhuǎn)換到相機坐標(biāo)系、圖像坐標(biāo)系、像素坐標(biāo)系。
然后再利用undistortPoint對這個點畸變矯正。
文章來源：http://www.zghlxwxcb.cn/news/detail-469507.html
文章來源地址http://www.zghlxwxcb.cn/news/detail-469507.html
到了這里，關(guān)于opencv圖像去畸變的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
