鱼眼图以“镜头球面ROI”方式还原

Posted 2020-07-16

大致实现思路：

【】引入原始的鱼眼图片，并显示；

【】创建另外几个Mat的大小；

【】确定镜头的参数（焦距，最大半径，最小半径，r = f*theta。1rad = CV_PI/180）；

【】原始img的坐标参考点左上角，变化到img中心处，原点O（ox，oy）；

【】创建一个二维动态数组（   vector<Point2f> points ）,以原点O为参考，创建一个直角坐标系的网格；

【】创建两个数组xp，yp，分别存 points的x，y坐标

【】进行直角坐标到笛卡尔坐标的转换，转换的结果存到对应的xp，yp中;  [ cartToPolar(xp, yp, rou, phi) ]

【】求原图的每个扇形ROI区域夹角边，对应在，球面上的ROI的夹角边，度数；

【】求视点的角度；（这个角度 = 每个扇形ROI两条夹角边之和的一半）

【】选择一个原鱼眼图的ROI，变换成展开的正常图片，并显示；

【】设计算法把原鱼眼图的极坐标映射到球坐标：二维->三维；

【】遍历原鱼眼图的二维数组points里的每一个值，求解它映射到球面后的theta（线面角）、phi（光线与三维坐标里x轴的夹角值），rou；

【】用三个Mat分别存取三维坐标里的值；

【】再用一个3层的Mat合成数据；

【】使用旋转矩阵，实现球面ROI的位置变化；

【】显示原图ROI常规矫正后的图片、显示映射到球面的图片、显示球面图片旋转过后矫正的图片

【源码：】注：以下代码中，变量名与上面举例时列出的没关系

#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/core/types_c.h>

using namespace cv;
using namespace std;

Mat fish, fishROI, nor, nor2; //鱼眼图（原图），ROI区域，正常图（效果图）
Mat mx, my; //remap参数，变换矩阵
Mat Rotation; //旋转后的坐标矩阵
Mat rou, phi; //rou 和 phi

IplImage *i1, *i2;

int ox, oy; //原点xy坐标
const int num_View = 6;

double f = 1, DR = 275, dr = 20, thetaMax = 80, thetaMin = 10, r = 1;//275
double viewPoint = CV_PI / 4; //虚拟视点度数
double view_phi[num_View + 1]; //分割区域的度数
double points[num_View + 1]; //虚拟视点的投影方向

void setO(Mat& m)
{
ox = m.cols / 2;
oy = m.rows / 2;
}

int main()
{
fish = imread("1.jpg", 1);
fishROI.create(fish.size(), fish.type());
//nor.create(fish.rows*2,fish.cols*2,fish.type());
nor.create(fish.size(), fish.type());
nor2.create(fish.size(), fish.type());
mx.create(fish.size(), CV_32FC1);
my.create(fish.size(), CV_32FC1);
setO(fish);
f = DR / (thetaMax / 180 * CV_PI);
thetaMin = dr / f;
//cout<<f;

imshow("【原鱼眼图】",fish);
waitKey(1000);
/*--------------------------------------------------------------------------------------------------------------------*/
vector<Point2f> point;//QQQQQQQQQQQQQQQQQQQ
//i对应x j对应y
for (int j = 0; j<fish.rows; j++)
{
for (int i = 0; i<fish.cols; i++)
{
point.push_back(Point2f(i + 1 - ox, oy - j + 1));
}
}

//极坐标与直角坐标互相转换的xy数组
Mat xp(point.size(), 1, CV_32F, &point[0].x, 2 * sizeof(float));
Mat yp(point.size(), 1, CV_32F, &point[0].y, 2 * sizeof(float));

//极坐标rou与phi
cartToPolar(xp, yp, rou, phi);
//ps:角度转弧度 π/180×角度,弧度变角度 180/π×弧度
//ps:依行遍历
/*--------------------------------------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/
//每个区域的度数
for (int i = 0; i<num_View; i++)
{
view_phi[i] = i * 2 * CV_PI / num_View;//弧度制
}
view_phi[num_View] = view_phi[0] + 2 * CV_PI;

//视点的位置角度
for (int i = 0; i<num_View; i++)
{
points[i] = ((view_phi[i + 1] + view_phi[i]) / 2);//角度制
}
/*--------------------------------------------------------------------------------------------------------------------*/
for (int tmpnum = 0; tmpnum<num_View; tmpnum++)
{
//ROI
Mat xxp(fish.size(), CV_32FC1);
Mat yyp(fish.size(), CV_32FC1);
for (int j = 0; j < fish.rows; j++)
{
for (int i = 0; i < fish.cols; i++)
{
if (phi.at<float>(i + j * fish.cols) >= view_phi[tmpnum] && phi.at<float>(i + j * fish.cols) <= view_phi[tmpnum + 1])
{
xxp.at<float>(j, i) = static_cast<float>(rou.at<float>(i + j * fish.cols) * cos(phi.at<float>(i + j * fish.cols)) + ox);
yyp.at<float>(j, i) = static_cast<float>(rou.at<float>(i + j * fish.cols) * -sin(phi.at<float>(i + j * fish.cols)) + oy);
}
else
{
xxp.at<float>(j, i) = static_cast<float>(fish.cols);
yyp.at<float>(j, i) = static_cast<float>(fish.rows);
}
}
}
remap(fish, fishROI, xxp, yyp, CV_INTER_LINEAR);
imshow("原图里的ROI", fishROI);
/*--------------------------------------------------------------------------------------------------------------------*/
//极坐标映射到球坐标：二维->三维 ***
Mat xx(point.size(), 1, CV_32F);
Mat yy(point.size(), 1, CV_32F);
Mat zz(point.size(), 1, CV_32F);
double *tphi = new double[point.size()];
double *ttheta = new double[point.size()];
for (int i = 0; i<point.size(); i++)
{
ttheta[i] = (CV_PI / 2 - rou.at<float>(i) / f);
tphi[i] = (phi.at<float>(i)-CV_PI);
}
for (int j = 0; j<fish.rows; j++)
{
for (int i = 0; i<fish.cols; i++)
{
if ((CV_PI / 2 - ttheta[i + j*fish.cols])*f <= DR)/**********************DR dr********************************/
{
xx.at<float>(i + j*fish.cols) = static_cast<float>(r*cos(ttheta[i + j*fish.cols])*-cos(tphi[i + j*fish.cols]));
yy.at<float>(i + j*fish.cols) = static_cast<float>(r*cos(ttheta[i + j*fish.cols])*sin(tphi[i + j*fish.cols]));
zz.at<float>(i + j*fish.cols) = static_cast<float>(r*sin(ttheta[i + j*fish.cols]));
}
else
{
xx.at<float>(i + j*fish.cols) = static_cast<float>(0);
yy.at<float>(i + j*fish.cols) = static_cast<float>(0);
zz.at<float>(i + j*fish.cols) = static_cast<float>(0);
}
}
}
/*--------------------------------------------------------------------------------------------------------------------*/
//[x,y,z]
Mat xyz(point.size(), 3, CV_32F);
for (int i = 0; i<point.size(); i++)
{
xyz.at<float>(i * 3) = xx.at<float>(i);
xyz.at<float>(i * 3 + 1) = yy.at<float>(i);
xyz.at<float>(i * 3 + 2) = zz.at<float>(i);
}
//旋转矩阵
double mtheta = -(view_phi[2] - view_phi[tmpnum + 1]), mphi = CV_PI / 4;
double r1[3][3] = { cos(mtheta), -sin(mtheta), 0, sin(mtheta), cos(mtheta), 0, 0, 0, 1 };//绕Z轴
double r2[3][3] = { 1, 0, 0, 0, cos(mphi), sin(mphi), 0, -sin(mphi), cos(mphi) };//绕X轴
double **db_xyz = new double *[point.size()];
double **db_xyz_2 = new double *[point.size()];
double **db_xyz_3 = new double *[point.size()];
for (int i = 0; i<point.size(); i++)
{
db_xyz[i] = new double[3];
db_xyz_2[i] = new double[3];
db_xyz_3[i] = new double[3];
db_xyz[i][0] = xyz.at<float>(i * 3);
db_xyz[i][1] = xyz.at<float>(i * 3 + 1);
db_xyz[i][2] = xyz.at<float>(i * 3 + 2);
}
for (int i = 0; i<point.size(); i++)
{
db_xyz_2[i][0] = db_xyz[i][0] * r2[0][0] + db_xyz[i][1] * r2[1][0] + db_xyz[i][2] * r2[2][0];
db_xyz_2[i][1] = db_xyz[i][0] * r2[0][1] + db_xyz[i][1] * r2[1][1] + db_xyz[i][2] * r2[2][1];
db_xyz_2[i][2] = db_xyz[i][0] * r2[0][2] + db_xyz[i][1] * r2[1][2] + db_xyz[i][2] * r2[2][2];
}
for (int i = 0; i<point.size(); i++)
{
db_xyz_3[i][0] = db_xyz_2[i][0] * r1[0][0] + db_xyz_2[i][1] * r1[1][0] + db_xyz_2[i][2] * r1[2][0];
db_xyz_3[i][1] = db_xyz_2[i][0] * r1[0][1] + db_xyz_2[i][1] * r1[1][1] + db_xyz_2[i][2] * r1[2][1];
db_xyz_3[i][2] = db_xyz_2[i][0] * r1[0][2] + db_xyz_2[i][1] * r1[1][2] + db_xyz_2[i][2] * r1[2][2];
}
for (int i = 0; i<point.size(); i++)
{
xx.at<float>(i) = static_cast<float>(db_xyz_3[i][0]);
yy.at<float>(i) = static_cast<float>(db_xyz_3[i][1]);
zz.at<float>(i) = static_cast<float>(db_xyz_3[i][2]);
}
/*--------------------------------------------------------------------------------------------------------------------*/
//球坐标映射到平面坐标：三维->二维
Mat nrou(point.size(), 1, CV_32F);
Mat nphi(point.size(), 1, CV_32F);
double PY = 0;

for (int j = 0; j<fish.rows; j++)
{
for (int i = 0; i<fish.cols; i++)
{
double ntheta = (CV_PI / 2 - atan((zz.at<float>(i + j*fish.cols) / (sqrt(xx.at<float>(i + j*fish.cols)*xx.at<float>(i + j*fish.cols) + yy.at<float>(i + j*fish.cols)*yy.at<float>(i + j*fish.cols))))));
nrou.at<float>(i + j*fish.cols) = static_cast<float>((ntheta)*f);
nphi.at<float>(i + j*fish.cols) = static_cast<float>(-(atan2(yy.at<float>(i + j*fish.cols), xx.at<float>(i + j*fish.cols)) - CV_PI));
mx.at<float>(j, i) = static_cast<float>(ox - nrou.at<float>(i + j*fish.cols)*cos(nphi.at<float>(i + j*fish.cols)));
my.at<float>(j, i) = static_cast<float>(oy - nrou.at<float>(i + j*fish.cols)*-sin(nphi.at<float>(i + j*fish.cols)));

if (PY<nrou.at<float>(i + j*fish.cols)*-sin(nphi.at<float>(i + j*fish.cols)))
{
PY = nrou.at<float>(i + j*fish.cols)*-sin(nphi.at<float>(i + j*fish.cols));
}
// double range=sqrt((i-ox)*(i-ox)+(j-oy)*(j-oy));
// mx.at<float>(j,i) = static_cast<float>((ox+((i-ox)/((range/f)/atan(range/f)))));
// my.at<float>(j,i) = static_cast<float>((oy+((j-oy)/((range/f)/atan(range/f)))));
}
}
cout << PY << endl;
remap(fishROI, nor, mx, my, CV_INTER_LINEAR);
imshow("【常规矫正过的图】", nor);

Mat ans;
ans.create(nor.size(), nor.type());

for (int j = 0; j<ans.rows; j++)
{
for (int i = 0; i<ans.cols; i++)
{
double range = sqrt((i - ox)*(i - ox) + (j - oy)*(j - oy));
mx.at<float>(j, i) = static_cast<float>((ox + ((i - ox) / ((range / f) / atan(range / f)))));
my.at<float>(j, i) = static_cast<float>((oy + ((j - oy) / ((range / f) / atan(range / f)))));
}
}
remap(nor, ans, mx, my, CV_INTER_LINEAR);
imshow("【映射后的球面图】", ans);
/*--------------------------------------------------------------------------------------------------------------------*/
Point2f center = Point2f(ans.cols / 2, ans.rows / 2);
double angle = 60 * (tmpnum - 1);
double scale = 1;
int nThresholdEdge = 1;
Mat rotateMat = getRotationMatrix2D(center, angle, scale);
Mat rotateImg;
warpAffine(ans, rotateImg, rotateMat, ans.size());
imshow("【球面旋转过后矫正结果图】", rotateImg);
/*--------------------------------------------------------------------------------------------------------------------*/
delete[](tphi);
delete[](ttheta);
delete[](db_xyz);
delete[](db_xyz_2);
delete[](db_xyz_3);
cvWaitKey();
}

//cvWaitKey();
return 0;
}