数字摄影测量实习报告(DOC) - 图文(5)

x = pic[i].x; y = pic[i].y; X = 0; Y = 0; Z = 0;

for (size_t j = 0; j < obj.size(); j++) {

if (obj[j].index == pic[i].index) {

X = obj[j].X; Y = obj[j].Y; Z = obj[j].Z; break; } }

if (abs(x) < 1) { char *arg1 = \像点\; char *arg2 = \没有找到对应的控制点\; ui->textBrowser->append(QString(\).arg(QString::fromLocal8Bit(arg1)).arg(pic[i].index).arg(QString::fromLocal8Bit(arg2))); continue; }

//求观测值的改正

r_r = (x - x0)*(x - x0) + (y - y0)*(y - y0);

dx = (x - x0)*(k1*r_r + k2*r_r*r_r) + p1*(r_r + 2 * (x - x0)*(x - x0)) - 2 * p2*(x - x0)*(y - y0);

dy = -(y - y0)*(k1*r_r + k2*r_r*r_r) + p2*(r_r + 2 * (y - y0)*(y - y0)) - 2 * p1*(x - x0)*(y - y0);

//经过畸变改正后的像点坐标观测值 x = (x - x0) - dx; y = -(y - y0) - dy; //求矩阵R各元素的值

a1 = cos(phi)*cos(kappa) - sin(phi)*sin(omega)*sin(kappa); a2 = -cos(phi)*sin(kappa) - sin(phi)*sin(omega)*cos(kappa); a3 = -sin(phi)*cos(omega); b1 = cos(omega)*sin(kappa); b2 = cos(omega)*cos(kappa); b3 = -sin(omega);

c1 = sin(phi)*cos(kappa) + cos(phi)*sin(omega)*sin(kappa); c2 = -sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa); c3 = cos(phi)*cos(omega); //求X_，Y_，Z_

X_ = a1*(X - Xs) + b1*(Y - Ys) + c1*(Z - Zs); Y_ = a2*(X - Xs) + b2*(Y - Ys) + c2*(Z - Zs); Z_ = a3*(X - Xs) + b3*(Y - Ys) + c3*(Z - Zs); //求x，y的近似值

double x_ = -fx*X_ / Z_; double y_ = -fy*Y_ / Z_; //求矩阵L各元素的值 lx = x - x_; ly = y - y_;

//求矩阵A各元素的值

a11 = 1 / Z_*(a1*fx + a3*x); a12 = 1 / Z_*(b1*fx + b3*x); a13 = 1 / Z_*(c1*fx + c3*x); a21 = 1 / Z_*(a2*fy + a3*y); a22 = 1 / Z_*(b2*fy + b3*y); a23 = 1 / Z_*(c2*fy + c3*y);

a14 = y*sin(omega) - (x / fx*(x*cos(kappa) - y*sin(kappa)) + fx*cos(kappa))*cos(omega);

a15 = -fx*sin(kappa) - x / fx*(x*sin(kappa) + y*cos(kappa)); a16 = y;

a24 = -x*sin(omega) - (y / fy*(x*cos(kappa) - y*sin(kappa)) - fy*sin(kappa))*cos(omega);

a25 = -fy*cos(kappa) - y / fy*(x*sin(kappa) + y*cos(kappa)); a26 = -x; //合成矩阵

double a[12] = { a11, a12, a13, a14, a15, a16, a21, a22, a23, a24, a25,

19

a26 };

double l[2] = { lx, ly };

cv::Mat Ai(2, 6, CV_64FC1, a); cv::Mat Li(2, 1, CV_64FC1, l); A.push_back(Ai); L.push_back(Li); }

//计算

dX = (A.t() * A).inv() * A.t() * L; Xs+=dX.at(0,0); Ys+=dX.at(1,0); Zs+=dX.at(2,0); phi+=dX.at(3,0); omega+=dX.at(4,0); kappa+=dX.at(5,0); ui->textBrowser->append(QString(\).arg(index).arg(dX.at(0)).arg(dX.at(1)).arg(dX.at(2)).arg(dX.at(3)).arg(dX.at(4)).arg(dX.at(5))); index++; cv::Mat M0 = (A*dX - L).t()*(A*dX - L); m0 = sqrt(M0.at(0) / (2 * pic.size() - 6)); A.release(); L.release();

if (index > 100) break; }

while(fabs(dX.at(0))>0.001||fabs(dX.at(3))>e); ui->textBrowser->append(QString(\)); ui->textBrowser->append(QString(\ Ys:%2 Zs:%3 phi:%4 omega:%5 kappa:%6\).arg(Xs).arg(Ys).arg(Zs).arg(phi).arg(omega).arg(kappa)); char *log = \单位权中误差：\; ui->textBrowser->append(QString(\).arg(QString::fromLocal8Bit(log)).arg(m0)); }

void HFJHDialog::on_pushButton_clicked() { QString qpath = QFileDialog::getOpenFileName(this); string path = qpath.toLocal8Bit(); if(readxy(path) == -1) { char *log =\读取像点坐标错误XXXXXX\; ui->textBrowser->append(QString::fromLocal8Bit(log)); return; } char *log = \读取像点坐标结束OOOOOO\; ui->textBrowser->append(QString::fromLocal8Bit(log)); ui->lineEdit->setText(qpath); }

void HFJHDialog::on_pushButton_2_clicked() { QString qpath = QFileDialog::getOpenFileName(this); string path = qpath.toLocal8Bit(); if(readXYZ(path) == -1) { char *log = \读取控制点坐标错误XXXXXX\; ui->textBrowser->append(QString::fromLocal8Bit(log)); return; } char *log = \读取控制点坐标结束OOOOOO\; ui->textBrowser->append(QString::fromLocal8Bit(log)); ui->lineEdit_2->setText(qpath); }

void HFJHDialog::on_pushButton_4_clicked() { HFJH(); }

void HFJHDialog::on_pushButton_5_clicked() { HFJHDialog::close();

20

}

图1 控件后方交会

2、特征匹配

#include \ #include #include #include using namespace std; using namespace cv; //构造函数

Mmatch::Mmatch(int length, int flag) { Mlength = length; Mflag = flag; match(); }

//构建子窗体

void Mmatch::findsubwin(Mat &src, Mat &dst, int i, int j, int length) { dst = Mat::zeros(length, length, CV_8UC1); int xl = i - length / 2; int yl = j - length / 2; int xr = i + length / 2; int yr = j + length / 2; int xdl = 0; int ydl = 0; int xdr = length - 1; int ydr = length - 1; if (xl<0) { xdl = 0 - xl; xl = 0; } if (yl < 0) { ydl = 0 - yl; yl = 0; } if (xr >= src.cols) {

21

xdr = (dst.cols -1)- (xr - (src.cols-1)); xr = src.cols - 1; } if (yr >= src.rows) { ydr = (dst.rows-1) - (yr - (src.rows-1)); yr = src.rows - 1; } Mat sROI = src(Range(yl+1, yr+1), Range(xl+1, xr+1)); sROI.copyTo(dst(Range(ydl+1, ydr+1), Range(xdl+1, xdr+1))); }

//匹配算法主体函数 int Mmatch::match() { string path1 = \我的课程\\\\数字摄影测量影像数据\\\\匹配\\\\1BAK.bmp\; string path2 = \我的课程\\\\数字摄影测量影像数据\\\\匹配\\\\2BAK.bmp\; string path = \我的课程\\\\数字摄影测量影像数据\\\\匹配\\\\1BAK_points.txt\; ifstream fil; vector pointssrc; fil.open(path); if (!fil) { cout << \ << endl; return -1; } int index; point pt; while (!fil.eof()) { fil >> index >> pt.x >> pt.y; pointssrc.push_back(pt); } Mat src1 = imread(path1,CV_LOAD_IMAGE_GRAYSCALE); Mat src2 = imread(path2,CV_LOAD_IMAGE_GRAYSCALE); if (src1.empty() || src2.empty()) { cout << \ << endl; return -1; } vector pointsdst; Mat subwinsrc, subwindst; for (int i = 0; i

22

} imshow(\, dst1); imshow(\, dst2); waitKey(); return 0; }

//寻找同名点函数

void Mmatch::findsamename(Mat &src, Mat &subwin, point &dst) { Mat subwind; int length = subwin.cols; int helf_length = length / 2; int flag = Mflag; float relate = 0; for (int i = 0; i < src.cols; ++i) { for (int j = 0; j < src.rows; ++j) { findsubwin(src, subwind, i, j, length); if (i == helf_length && j == helf_length) { relate = compare(subwin, subwind, flag); dst.x = i; dst.y = j; } float tmp = compare(subwin, subwind, flag); if (tmp > relate) { relate = tmp; dst.x = i; dst.y = j; } } } }

//特征向量比较函数

float Mmatch::compare(Mat &subwin, Mat &subwind,int flag) { float sum = 0.0; float avg1 = 0.0, avg2 = 0.0; float sdv1 = 0.0, sdv2 = 0.0; int nr = subwin.rows; int nc = subwin.cols; if(subwin.isContinuous()&&subwind.isContinuous()) { nc = nc*nr; nr = 1; } switch (flag) { case 1: //相关函数法 for (int i = 0; i < nr; ++i) { uchar *subwin_ptr = subwin.ptr(i); uchar *subwind_ptr= subwind.ptr(i); for (int j = 0; j < nc; j++) { sum += (float)(*(subwin_ptr+j)* (*(subwind_ptr+j))); } } return sum; case 2://协方差函数法 for (int i = 0; i < nr; ++i) { uchar *subwin_ptr = subwin.ptr(i); uchar *subwind_ptr= subwind.ptr(i); for (int j = 0; j < nc; j++) { avg1 += *(subwin_ptr+j);

23