数字摄影测量实习报告
数字摄影测量实习报告书
学 号：20111000684
班级序号： 113112-05
姓 名： 舒 超
指导老师： 宋 妍
成 绩：
中国地质大学（武汉）
信息工程学院遥感科学技术系
20##年6月
目录
实习一：Moravec算子点特征提取... 3
1.1 实习目的：... 3
1.2 实习原理：... 3
1.3 实习步骤以及代码分析：... 4
1.4 结果分析：... 8
实习二：边缘提取算法... 10
2.1 实习目的：... 10
2.2 实习原理：... 10
2.3 实习步骤以及代码：... 10
2.4 结果分析：... 12
实习总结... 13
实习一：Moravec算子点特征提取
1.1 实习目的：
用程序设计语言（VisualC++或者C语言）编写一个完整的提取点特征的程序，通过对提供的图像数据进行特征点提取，输出提取出的点特征坐标。本实验的目的在于让学生深入理解Moravec算子原理。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
1.2 实习原理：
Moravec角点检测算法是最早的角点检测算法之一。该算法将角点定义为具有低“自相关性”的点。算法会检测图像的每一个像素，将像素周边的一个邻域作为一个patch，并检测这个patch和周围其他patch的相关性。这种相关性通过两个patch间的平方差之和(SSD)来衡量，SSD值越小则相似性越高。
如果像素位于平滑图像区域内，周围的patch都会非常相似。如果像素在边缘上，则周围的patch在与边缘正交的方向上会有很大差异，在与边缘平行的方向上则较为相似。而如果像素是各个方向上都有变化的特征点，则周围所有的patch都不会很相似。
Moravec会计算每个像素patch和周围patch的SSD最小值作为强度值，取局部强度最大的点作为特征点。
1.3 实习步骤以及代码分析：
步骤流程图如下：
程序实现以及相关关键代码：
voidCMy2010302590183cylView::OnMoravec()//读取图像以及相关算法
{
//TODO:Addyourcommandhandlercodehere
CMmoravecDlgdlg;
dlg.DoModal();
CMy2010302590183cylDoc*pDoc=GetDocument();
LPSTRm_pDIB=(LPSTR)::GlobalLock((HGLOBAL)pDoc->hdib);//得到句柄内存起始地址存放位图数据hdib句柄变量存放BMP位图
::GlobalUnlock((HGLOBAL)pDoc->hdib);
LPBITMAPINFOm_pBMP;//指向BITMAPINFO结构的指针
m_pBMP=(LPBITMAPINFO)::GlobalLock(pDoc->hdib);
//获取指向BITMAPINFO结构的指针
::GlobalUnlock((HGLOBAL)pDoc->hdib);
intBitCount=m_pBMP->bmiHeader.biBitCount;
DWORDWidth=::DIBWidth(m_pDIB);//获取位图宽
DWORDHeight=::DIBHeight(m_pDIB);//获取位图高
LPBYTElpData=(LPBYTE)::FindDIBBits(m_pDIB);//定义字符指针变量,原位图指针
intWidthBytes=WIDTHBYTES(Width*BitCount);//获取字节
DWORDpixelCount=WidthBytes*Height;
intck1=dlg.c1;
intck2=dlg.c2;
doubleyz=dlg.m_yuzhi;
DWORDr,c;
INTh;
double*xx=newdouble[Width*Height];
intk;
k=INT(ck1/2);
for(r=ck1/2;r<Height-ck1/2;r++)
for(c=ck1/2;c<Width-ck1/2;c++)
{
doublemin,v[4]={0.0};
for(h=0;h<=ck1-1;h++) {
v[0]+=pow((double)(*((BYTE*)(lpData+r*WidthBytes+(c-k+h)))-*((BYTE*)(lpData+(r)*WidthBytes+(c-k+1+h)))),2);//0°方向
v[1]+=pow((double)(*((BYTE*)(lpData+(r-k+h)*WidthBytes+(c+k-h)))-*((BYTE*)(lpData+(r-k+h+1)*WidthBytes+(c+k-h-1)))),2);
//45°方向
v[2]+=pow((double)(*((BYTE*)(lpData+(r-k+h)*WidthBytes+(c)))-*((BYTE*)(lpData+(r-k+1+h)*WidthBytes+(c)))),2);
//90°方向
v[3]+=pow((double)(*((BYTE*)(lpData+(r-k+h)*WidthBytes+(c-k+h)))-*((BYTE*)(lpData+(r-k+1+h)*WidthBytes+(c-k+h+1)))),2);
//135°方向
}
min=min(min(min(v[0],v[1]),v[2]),v[3]);
//求出v1,v2,v3,v4中的最小值
if(min>yz)
xx[r*Width+c]=min;
}
bool*bMatrix=newbool[Width*Height];
memset(bMatrix,0,Width*Height*sizeof(bool));
DWORDx,y;
doublemax2;
boolb=false;
inttempX(0),tempY(0);
for(x=ck2/2;x<Height-ck2/2;x+=ck2)
{
for(y=ck2/2;y<Width-ck2/2;y+=ck2)
{
max2=0;
for(DWORDm=(x-ck2/2);m<(x+ck2/2);m++)
{
for(DWORDn=(y-ck2/2);n<(y+ck2/2);n++)
if(xx[m*Width+n]>max2)
{
max2=xx[m*Width+n];
tempY=m;
tempX=n;
b=true;
}
}
if(b)
{ bMatrix[tempY*Width+tempX]=1; }
}
}
intsum=0;//特征点总数
for(DWORDi=0;i<Height;i++)
for(DWORDj=0;j<Width;j++)
{
if(bMatrix[i*Width+j])
{
*((BYTE*)(lpData+i*WidthBytes+j))=0;
*((BYTE*)(lpData+i*WidthBytes+j+1))=0;
*((BYTE*)(lpData+i*WidthBytes+j-1))=0;
*((BYTE*)(lpData+(i+1)*WidthBytes+j))=0;
*((BYTE*)(lpData+(i-1)*WidthBytes+j))=0;
*((BYTE*)(lpData+i*WidthBytes+j+2))=0;
*((BYTE*)(lpData+i*WidthBytes+j-2))=0;
*((BYTE*)(lpData+(i+2)*WidthBytes+j))=0;
*((BYTE*)(lpData+(i-2)*WidthBytes+j))=0;
sum++;
}
}
if(sum<4000)
{
CStringstrInfo;
strInfo.Format("特征点数%d\n",sum);
MessageBox(strInfo,"提示",MB_OK);
}
else
{
CStringstrInfo;
strInfo.Format("特征点数较多，请设置合理参数");
MessageBox(strInfo,"提示",MB_OK);
}
Invalidate();
}
1.4 结果分析：
按照提示，对老师所给数据进行分析，当窗口大小设置为5*5,，阈值设置为5000的时候，对右核线影像进行分析，得到特征点43个，同时图像分析，得出如下结果：
调整阈值和窗口大小，程序能够正常运行，且经过测试，结果精确度有较好的保证。
实习二：边缘提取算法
2.1 实习目的：
熟悉Matlab环境下的编程，熟悉边缘提取算法。
2.2 实习原理：
Sobel算子实现思路如下：对输入图像分别使用水平和垂直模板做卷积计算，对得到的两个处理结果求平方和，该平方和与阈值的平方比较。只有当某点的两种卷积的平方大于阈值的平方，且水平占优（水平模板卷积结果大于垂直模板卷积的结果，且该点的卷积平方大于其左右两点的卷积平方和）或者垂直占优（垂直模板卷积的结果大于水平模板卷积的结果，且该点的卷积平方和大于其上下两点的卷积平方和）时，该点的输出结果为255，否则为0。输出的结果为二值图像。第一行和最后一行本来就是图像边界，不包括可用信息，因此相应的输出为0，按照这个思路课题编写了相应的Sobel算子实现程序
2.3 实习步骤以及代码：
2.4 结果分析：
原图像 sobel边缘提取
实习总结
本次实习过程中，根据自身实际情况，我选择使用vc环境下的编程完成实习，而没有采用Matlab环境下的编程。在实习过程中，我熟悉了sobel算法以及Moravec算子，在程序调试的过程中，我认识到任何算法都有其局限性，比如说本次实习过程中，sobel算子的边缘提取就将许多的噪音提取了出来，导致边缘特征提取的不准确性。本次实习让我认识到了编程能力的重要性，学会编写基本的代码来实现基本的算法，能让我们摆脱软件已有算法的束缚，更多的按照需要来实现一些步骤。
总体来说，本次实习还是很成功的，让我认识到，在以后的学习生活中,我认识到,应该把理论和实践结合起来,多锻炼自己的动手能力,好好把握住每一次实习的机会.
第二篇：摄影测量实习报告翟智佳
摄影测量实习报告
——单张影像空间后方交会程序设计
实习时间 2013.5.20-2013.5.24
学生班级 学生姓名
学生学号所在院系 矿业工程学院
指导教师
一、实习目的
1．深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
2．利用Visual C++或者Matlab（或其他熟悉的计算机语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并进行评定精度。
3．通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的动手能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。
二、实习环境
1．硬件环境：Window操作系统
2．软件环境：VC++或Matlab或其他计算机语言
三、实习内容
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素并进行精度评定。
四、实习原理
1．共线方程
x??fa1(X?Xs)?b1(Y?Ys)?c1(Z?Zs)
a3(X?Xs)?b3(Y?Ys)?c3(Z?Zs)
a(X?Xs)?b2(Y?Ys)?c2(Z?Zs)y??f2
a3(X?Xs)?b3(Y?Ys)?c3(Z?Zs)
2．精度评定
mi?m0?ii 其中m0?? VV2n?6
五、实习数据
1．模拟像片一对：左片号23 右片号24
2．像片比例尺： 1/30000
3．航摄机主距：f=150mm
4．每张像片有4个控制点
5．点位略图
6．各片像点坐标及其地面坐标
六、验证数据
1．已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。 2．已知4对点的影像坐标和地面坐标
3．解算结果
?0.997710.067530.00399??? YS?27476.46m R??-0.067530.99772-0.00221??ZS?7572.69m?-0.004120.001840.9999??0XS?39795.45m
七、程序过程框图：
八.实现程序
#include "iostream.h"
#include"stdio.h"
#include "stdlib.h"
#include<math.h>
#define N 4
void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2)//矩阵相乘
{
int i,j,k;
for(i=0;i<i_1;i++)
for(j=0;j<j_2;j++)
{
result[i*j_2+j]=0.0;
for(k=0;k<j_12;k++)
result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2]; }
return;
}
int invers_matrix(double *m1,int n)//矩阵求逆
{
int *is,*js;
int i,j,k,l,u,v;
double temp,max_v;
is=(int *)malloc(n*sizeof(int));
js=(int *)malloc(n*sizeof(int));
if(is==NULL||js==NULL){
printf("out of memory!\n");
return(0);
}
for(k=0;k<n;k++){
max_v=0.0;
for(i=k;i<n;i++)
for(j=k;j<n;j++){
temp=fabs(m1[i*n+j]);
if(temp>max_v){
max_v=temp; is[k]=i; js[k]=j;
}
}
if(max_v==0.0){ free(is); free(js); printf("invers is not availble!\n"); return(0); } if(is[k]!=k) for(j=0;j<n;j++){ u=k*n+j; v=is[k]*n+j; temp=m1[u]; m1[u]=m1[v]; m1[v]=temp; } if(js[k]!=k) for(i=0;i<n;i++){ u=i*n+k; v=i*n+js[k]; temp=m1[u]; m1[u]=m1[v]; m1[v]=temp; } l=k*n+k; m1[l]=1.0/m1[l]; for(j=0;j<n;j++) if(j!=k){ u=k*n+j; m1[u]*=m1[l]; } for(i=0;i<n;i++) if(i!=k) for(j=0;j<n;j++) if(j!=k){ u=i*n+j; m1[u]-=m1[i*n+k]*m1[k*n+j]; } for(i=0;i<n;i++) if(i!=k){ u=i*n+k; m1[u]*=-m1[l]; } } for(k=n-1;k>=0;k--){ if(js[k]!=k) for(j=0;j<n;j++){ u=k*n+j; v=js[k]*n+j; temp=m1[u]; m1[u]=m1[v]; m1[v]=temp; } if(is[k]!=k) for(i=0;i<n;i++){ u=i*n+k; v=i*n+is[k];
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
}
free(is); free(js);
return(1);
}
void transpose(double *m1,double *m2,int m,int n) //矩阵转置 {
int i,j;
for(i=0;i<m;i++)
for(j=0;j<n;j++)
m2[j*m+i]=m1[i*n+j];
return;
}
void main()
{
double Xs,Ys,Zs,q,w,k;
double a[3],b[3],c[3];
double x0,y0,f;
double x[N],y[N];
double X[N],Y[N],Z[N];
double x1[N],y1[N];
double m;
double L[2*N];
double XX[6];
double A[2*N][6];
double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1]; int i,n=0;
double sum=0,m0;
/*---------------输入点地面坐标---------------------*/
for(i=0;i<N;i++)
{
printf("请输入第%d个点的地面坐标：",i+1);
scanf("%lf%lf%lf",&X[i],&Y[i],&Z[i]);
}
/*---------------输入点像片坐标---------------------*/
for(i=0;i<N;i++)
{
printf("请输入第%d个点的像片坐标：",i+1);
scanf("%lf%lf",&x[i],&y[i]); } cout<<endl; /*-----------------设定外方位元素初始值--------------*/ x0=0;y0=0;f=150;m=30000; Xs=0;Ys=0;Zs=f*m/1000; q=0;w=0;k=0; XX[3]=1; /*------------------迭代计算--------------------------*/ while((XX[3]>0.00001 || XX[4]>0.00001 || XX[5]>0.00001)&&n<100) { /*----------------旋转矩阵R-----------------------*/ a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a[2]=-sin(q)*cos(w); b[0]=cos(w)*sin(k); b[1]=cos(w)*cos(k); b[2]=-sin(w); c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c[2]=cos(q)*cos(w); /*-----------------像点坐标计算值------------------*/ for(i=0;i<N;i++) { X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs); Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs); Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs); x1[i]=x0-f*X0[i]/Z0[i]; y1[i]=y0-f*Y0[i]/Z0[i]; } /*-------------误差方程中各偏导数的值--------------*/ for(i=0;i<N;i++) { A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i]; A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i]; A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i]; A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin(k))/f+f*cos(k))
*cos(w);
A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f; A[2*i][5]=y[i]-y0;
L[2*i]=x[i]-x1[i];
A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][1]=((b[1]*f+b[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][2]=((c[1]*f+c[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][3]=-(x[i]-x0)*sin(w)-((y[i]-y0)*((x[i]-x0)*cos(k)-(y[i]-y0)*sin(k))/f-f*sin(k))
*cos(w);
A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;
A[1+2*i][5]=-x[i]+x0;
L[1+2*i]=y[i]-y1[i];
}
/*-------------------解法方程--------------------*/
transpose(&A[0][0],&At[0][0],2*N,6);
mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);
invers_matrix(&result1[0][0],6);
mult(&At[0][0],L,&result2[0][0],6,2*N,1);
mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);
Xs+=XX[0];
Ys+=XX[1];
Zs+=XX[2];
q+=XX[3];
w+=XX[4];
k+=XX[5];
n++;
}
/*----------------旋转矩阵R-----------------------*/
a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);
a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);
a[2]=-sin(q)*cos(w);
b[0]=cos(w)*sin(k);
b[1]=cos(w)*cos(k);
b[2]=-sin(w);
c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);
c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);
c[2]=cos(q)*cos(w);
cout<<"迭代次数为:"<<n<<endl;
printf("\n像片外方位元素的解\n");
} cout<<"航向顷角q:"<<q<<endl; cout<<"旁向倾角w:"<<w<<endl; cout<<"像片旋角k:"<<k<<endl; printf("Xs=%lf\t\tYs=%lf\t\tZs=%lf\n",Xs,Ys,Zs); cout<<endl; printf("旋转矩阵R:\n"); for(i=0;i<3;i++) printf("%lf\t",a[i]); printf("\n"); for(i=0;i<3;i++) printf("%lf\t",b[i]); printf("\n"); for(i=0;i<3;i++) printf("%lf\t",c[i]); printf("\n"); /*-------------------计算单位权中误差---------------*/ for(i=0;i<2*N;i++) sum+=L[i]*L[i]; m0=sqrt(sum/(2*N-6)); cout<<"单位权中误差m0="<<m0<<endl; cout<<"测量精度:"<<endl; cout<<"Xs的精度为:"<<m0*sqrt(result1[0][0])<<endl; cout<<"Ys的精度为:"<<m0*sqrt(result1[1][1])<<endl; cout<<"Zs的精度为:"<<m0*sqrt(result1[2][2])<<endl; cout<<"q的精度为:"<<m0*sqrt(result1[3][3])<<endl; cout<<"w的精度为:"<<m0*sqrt(result1[4][4])<<endl; cout<<"k的精度为:"<<m0*sqrt(result1[5][5])<<endl;
九、程序运行结果 23像片运行结果
24号像片运行结果
10.实习心得
通过一周的实习，让我感觉得到了自己摄影测量知识的欠缺
和编程方面的无措。不过最后在同学们的帮助下中雨完成了这次实习 而在编程里也学习了好多以前没见过的函数，mult是矩阵相乘的函数、invers_矩阵求逆函数、transpose矩阵转换函数。这对以后的学习是大有裨益的。
通过此次实习，才知道书到用时方恨少，同时下定决心好好学习，非常感谢老师给我们提供此次实习。
-
数字摄影测量实习报告
数字摄影测量实习报告书学号20xx1000684班级序号11311205姓名舒超指导老师宋妍成绩中国地质大学武汉信息工程学院遥感科…
-
VirtuoZo NT全数字摄影测量系统实习报告
VirtuoZo全数字摄影测量系统实习报告一、实习目的通过本次实习，了解4D产品的生产过程，熟悉使用VirtuoZo全数字摄影测量…
-
数字摄影测量实习报告书
实习报告书专用纸1引言数字摄影测量是基于摄影测量的基本原理通过对所获取的数字数字化影像进行处理自动半自动提取被摄对象用数字方式表达…
-
全数字摄影测量实习报告
全数字摄影测量实习报告学号姓名班级专业课程名称数字摄影测量学指导老师20xx年1月一目的和要求了解全数字摄影测量的作业步骤二资料及…
-
数字摄影测量实习报告
目录1目的与要求12实习内容13操作步骤131新建测区设置测区参数文件132模型生成333模型定向与生成核线影像4331模型内定向…
-
摄影外业测量实习报告
组别：指导老师：姓名：班级：学号:黄河水利职业技术学院实习内容：航摄像片调绘、新增地物补测及选刺像控点三项实习时间：20xx年x月…
-
摄影测量实习报告
一、实习目的与要求本次实习是在摄影测量的教学基础上，理论实际相联系的动手操作实习，是我们在学习测量专业的一个重要的实习环节。一方面…
-
VirtuoZo NT全数字摄影测量系统实习报告
VirtuoZo全数字摄影测量系统实习报告一、实习目的通过本次实习，了解4D产品的生产过程，熟悉使用VirtuoZo全数字摄影测量…
-
摄影测量内业实习报告
摄影测量内业实习报告姓名xxx班级学号xxx指导老师XX20xx年1月1实习目的VirtuoZo从原始资料中间成果及最后产品等都是…
-
摄影测量实习报告
摄影测量学实习报告系别测绘与城市空间信息系专业姓名学号指导教师柏春岚杨明强河南城建学院20xx年11月26日目录前言1实习意义与目…
-
摄影测量-实习总结
黄河水利职业技术学院测绘工程系摄影测量实习总结班级：工程测量0905班姓名：任二朋学号：20xx020618时间：20xx年x月在…