摄影测量之内定向程序
1. 内定向概念
内定向的目的是将相片纠正到相片坐标，通常方法是相片的周边有一系列的框标点，通常有4个或8个，它们的相片坐标是事先经过严格校正过的，利用这些点构成一个仿射变换的模型（或多项式），把象素纠正到相片坐标系。通过这一步基本上消除了相片因扫描、压平等因素导致的变形。
2. 内定向原理
参考课件
3. 内定向程序设计
1、首先编写一个处理矩阵的模块文件：MatrixOperations.h和MatrixOperations.cpp
MatrixOperations.h
#pragma once
//transpose数组将from进行转置，转置结果存储在to里
//from为row行、col列
void transpose(double* from, double* to, int row, int col);
//矩阵a与矩阵b相乘，相乘的结果放在矩阵c里
//矩阵a为m行n列，矩阵b为n行、l列
void multiply(double* a, double* b, double* c, int m, int n, int l);
//矩阵a与矩阵b相加，结果存放在矩阵a中
void mat_add(double* a, double* b, int m, int n);
void make2array(double**& a, int n);
void deletarray(double**& a, int n);
//矩阵求逆
//ip为要求逆的矩阵，rp是返回的逆矩阵，矩阵维数为n
//行列式为0，返回0；否则返回值为1
int converse(double* ip, double* rp, int n);
MatrixOperations.cpp
#include "MatrixOperations.h"
//transpose数组将from进行转置，转置结果存储在to里
//from为row行、col列
void transpose(double* from, double* to, int row, int col)
{
int i, j;
for (i = 0; i < row; i++)
{
for (j = 0; j < col; j++)
{
to[j * row + i] = from[i * col + j];
}
}
}
//矩阵a与矩阵b相乘，相乘的结果放在矩阵c里
//矩阵a为m行n列，矩阵b为n行、l列
void multiply(double* a, double* b, double* c, int m, int n, int l)
{
int i, j, k;
for (i = 0; i < m * l; i++)
{
c[i] = 0;
}
for (i = 0; i < m; i++)
{
for (j = 0; j < l; j++)
{
for (k = 0; k < n; k++)
{
c[i * l + j] += a[i * n + k] * b[k * l + j];
}
}
}
}
//矩阵a与矩阵b相加，结果存放在矩阵a中
void mat_add(double* a, double* b, int m, int n)
{
int i, j;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
a[i * n + j] += b[i * n + j];
}
}
}
void make2array(double**& a, int n)
{
int i, j;
a = new double* [n];
for (i = 0; i < n; i++)
{
a[i] = new double[2 * n];
}
for (i = 0; i < n; i++)
{
for (j = 0; j < 2 * n; j++)
{
a[i][j] = 0;
}
}
}
void deletarray(double**& a, int n)
{
int i;
for (i = 0; i < n; i++)
{
delete[]a[i];
}
delete[]a;
}
//矩阵求逆
//ip为要求逆的矩阵，rp是返回的逆矩阵，矩阵维数为n
//行列式为0，返回0；否则返回值为1
int converse(double* ip, double* rp, int n)
{
double** mat; //做行变换的矩阵
int i, j, r;
double k, temp;
//初始化mat为0,大小为n*2n
make2array(mat, n);
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
mat[i][j] = ip[i * n + j];
}
mat[i][n + i] = 1; //初始化右侧的单位矩阵
}
//做行变换化为上三角阵
for (i = 0; i < n; i++)
{
if (mat[i][i] == 0) //第i行i列为0
{
for (j = i + 1; j < n; j++)
{
if (mat[j][i] != 0) //找一个非0行
{
for (r = i; r < 2 * n; r++)
{
temp = mat[j][r];
mat[j][r] = mat[i][r];
mat[i][r] = temp;
}
break; //跳出j循环
}
}
}
if (mat[i][i] == 0) return 0; //行列式为0则返回
for (j = i + 1; j < n; j++)
{
if (mat[j][i] != 0) //i行i列下方的j行i列元素不为0
{
k = -mat[j][i] / mat[i][i]; //做行变换
for (r = i; r < 2 * n; r++)
mat[j][r] = mat[j][r] + k * mat[i][r];
}
}
}
//化成单位矩阵
for (i = n - 1; i >= 0; i--)
{
k = mat[i][i];
for (r = i; r < 2 * n; r++)
mat[i][r] = mat[i][r] / k;
for (j = i - 1; j >= 0; j--)
{
k = -mat[j][i];
for (r = i; r < 2 * n; r++)
mat[j][r] = mat[j][r] + k * mat[i][r];
}
}
//将结果输出
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
rp[i * n + j] = mat[i][j + n];
//mat矩阵释放
deletarray(mat, n);
return 1;
}
2、然后编写主程序main.cpp，实现对内定向参数的求解
main.cpp
#include<stdio.h>
#include <stdlib.h>
#include"MatrixOperations.h"
#define M1 8
#pragma warning(disable : 4996)
int main()
{
int i=0,j,k=0,q=0;
double A[2*M1][6]={0},A1[2][6]={0},A2[2*M1][1],A3[M1][2],A4[2][1],B[6][2*M1],C[6][6],D[6][6],E[6][2*M1],L[2*M1][1],X[6][1];
//其中A1,A4在计算单点的像平面坐标时使用
double h1,h2,width=0.045;
double test[M1][2];
double hls[M1][2];
FILE *fin;
fin = fopen("matlab.txt","r"); // 打开文件matlab.txt，按读的方式打开
for (i=0;i<M1;i++)
fscanf(fin,"%lf %lf", &hls[i][0], &hls[i][1]); // 循环读
fclose(fin); //关闭文件
double v[100],L0[100]; //开一个足够大的数组。
double x0,y0,z0;
i=0;
FILE *fp; //文件指针
fp = fopen("rc30.cmr", "r"); //以文本方式打开文件rc30.cmr。
if(fp == NULL) //打开文件出错。
return -1;
while(fscanf(fp, "%lf", &v[i]) != EOF) //读取数据到数组，直到文件结尾(返回EOF)
i++;
fclose(fp); //关闭文件
x0=v[0];
y0=v[1];
z0=v[2];
for(j=3;j<i;j++){
if(j%3==0) continue;
L0[q++]=v[j];}
for(i=0;i<M1;i++)
for(j=0;j<2;j++){
test[i][j]=hls[i][j]*width;} //test
for(q=0,i=0;i<2*M1;i++)
L[k++][0]=L0[q++]; //test-L 一维矩阵的转置
/***test-L-系数矩阵A***/
for(i=0;i<2*M1;i=i+2)
{
A[i][0]=1;
A[i][1]=test[i/2][0];
A[i][2]=test[i/2][1];
}
for(i=1;i<2*M1;i=i+2)
{
A[i][3]=1;
A[i][4]=test[(i-1)/2][0];
A[i][5]=test[(i-1)/2][1];
}
/*矩阵运算,求定向转换参数X*/
transpose(*A,*B,2*M1,6);
multiply(*B,*A,*C,6,2*M1,6);
converse(*C,*D,6);
multiply(*D,*B,*E,6,6,2*M1);
multiply(*E,*L,*X,6,2*M1,1);
multiply(*A,*X,*A2,2*M1,6,1);
for(i=0,j=0;i<2*M1;i=i+2,j++)
A3[j][0]=A2[i][0]-x0;
for(i=1,j=0;i<2*M1;i=i+2,j++)
A3[j][1]=A2[i][0]-y0; //考虑主点偏移量
printf("内定向转换参数X:\n");
for(i=0;i<6;i++){
printf("\n");
for(j=0;j<1;j++)
printf("%12.6lf\t",X[i][j]);} //内定向转换参数X
printf("\n\n");
printf("各个框标的像平面坐标A3:(用于对比验证)\n");
for(i=0;i<M1;i++){
printf("\n");
for(j=0;j<2;j++)
printf("%12.6lf\t",A3[i][j]);} //主点偏移量后像平面坐标A3
printf("\n\n");
ccc:
printf("请输入一个像点的行列值\n\n");
scanf("%lf%lf",&h1,&h2);
test[0][0]=h1*width;
test[0][1]=h2*width; //test2
/***test-L-系数矩阵A1***/
{
A1[0][0]=1;
A1[0][1]=test[0][0];
A1[0][2]=test[0][1];
A1[1][3]=1;
A1[1][4]=test[0][0];
A1[1][5]=test[0][1];
}
multiply(*A1,*X,*A4,16,6,1);
A4[0][0]=A4[0][0]-x0;
A4[1][0]=A4[1][0]-y0;
printf("\n该点的像空间坐标\n\n");
for(i=0;i<2;i++)
printf("%12.6lf\t",A4[i][0]); //点的像平面坐标
printf("%12.6lf",-z0);
printf("\n\n");
goto ccc;
return 0;
}
3、需要读入的的数据是相机文件rc30.cmr和量测的框标点的行列号文件matlab.txt
rc30.cmr
-0.004
-0.008
152.72
1
106.003
-106.005
2
-106.003
-106.003
3
-106.001
106.003
4
106.001
106.001
5
0.002
-109.998
6
-110.001
0.002
7
-0.004
110.004
8
109.997
0.000
matlab.txt
4862.5 4932.5
155.5 4889.5
197.5 181.5
4905.5 223.5
2508 5000
87.5 2535
2552 114
4972 2578.5
4. 结果输出
1926
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
