摄影测量手动提取立体像对同名点像素，相对定向解算，特征点提取，影像匹配同名点的python实现
-1 2023-6-2更
很高兴有那么多人看，可以说非常意外了。
这篇博客是本人一年前还在学校的时候写的，当时是结课作业，没想太多就照着书上的算法实现了一遍，一个礼拜的时间，算是体验了一下造轮子的痛苦……
写这些代码的时候，我还是非常稚嫩的，所以有很多地方可能写的不对，也来不及去验证结果的准确性了，后面和我的研究方向也并不契合，所以就一直搁置在这了。
很久没看，忘记源码要积分……对不起，我已经修改成0积分了下载了，有需要的可以下载下来看看，如果不嫌弃如此粗糙的代码的话……
已经工作快一年了，很早就想着做一个自己的网站，后面要写博客文章什么的，也会放到自己网站上去，准备用gitpage做个简单的静态网站，而且现在工作也没做前端了……真是造化弄人。
咳咳，附上地址，龟速开发中……
gitpage
0 写在前面
- 摄影测量课程布置的课程大作业，主要要求是：
独立完成一种相对定向的算法程序设计与开发。功能包括：
（1）打开立体像对，手工屏幕量测像点坐标，转换成像平面坐标（o-xy），取得5个以上同名像点坐标并保存；
（2）进行相对定向计算；
（3）统计观测值的单位权中误差，并实施一种剔除粗差的处理。
（4）采用任意一种匹配算法，通过影像匹配取得同名像点，保存6-10个质量可靠、位置合理的同名像点坐标。重复上面的（2）（3）步骤；
（5）比较手工量测和计算机自动匹配同名像点，两者对相对定向结果有无差别。 - 使用的python库
numpy 操作数组和矩阵
cv2 操作图像
pandas 读写文件 - 文中具体公式原理参照：摄影测量学第二版武汉大学出版社、误差理论与测量平差基础第三版武汉大学出版社。同时，代码也写了注释
- 所有的文件代码都会上传，下载地址在文章最后
1 进行立体相对的手工屏幕量测
- 使用OpenCV，定义一个供setMouseCallback使用的回调函数mouse，这个回调函数在捕获到鼠标左键点击事件(cv2.EVENT_LBUTTONDOWN)时，获取点击的像素点坐标，并绘制一个实心的圆且显示出坐标。
具体实现参照：
python+OpenCV交互图片获取鼠标单击坐标点（解决大图片显示不全问题） - 获取得到5对以上的同名点坐标，我选取了9对：
uvs1 = [[275, 1566, 1], [3205, 1979, 2], [5729, 1748, 4], [6075, 3166, 5], [3933, 3189, 6], [566, 3121, 7], [326, 4469, 8], [3499, 4366, 9], [6110, 4535, 10]] #左片采样点像素坐标,其中3号点位为错误点
uvs2 = [[676, 315, 1], [3513, 521, 2], [6061, 56, 3], [6502, 1429, 4], [4325, 1659, 5], [1093, 1840, 6], [975, 3127, 7], [3979, 2832, 8], [6633, 2836, 9]] #右片采样点像素坐标[y,x,点号]
2 像素坐标与像平面、像空间、像空间辅助坐标的转换
- 像素点坐标转像平面坐标
坐标转换参照：
计算机视觉：相机成像原理：世界坐标系、相机坐标系、图像坐标系、像素坐标系之间的转换
xArrs1 = [] #存放得到的左片像点平面坐标[-x,-y,1] 因为左片在下，右片在上
xArrs2 = [] #存放得到的右片像点平面坐标
# 像素点坐标转像平面坐标
def pixelToImgPoint():
senserSize = [35.9, 24] #传感器尺寸 宽*高（mm）
pixelSize = [7360,4912] #像素尺寸 宽*高 （像素）
u0 = pixelSize[0] / 2 #中心点像素横坐标
v0 = pixelSize[1] / 2 #中心点像素横坐标
dx = senserSize[0] / pixelSize[0] #x轴1像素实际大小(mm)
dy= senserSize[1] / pixelSize[1] #y轴1像素实际大小
f = 35 #焦距（mm）
uvs1 = [[275, 1566, 1], [3205, 1979, 2], [5729, 1748, 4], [6075, 3166, 5], [3933, 3189, 6], [566, 3121, 7], [326, 4469, 8], [3499, 4366, 9], [6110, 4535, 10]] #左片采样点像素坐标,其中3号点位为错误点
uvs2 = [[676, 315, 1], [3513, 521, 2], [6061, 56, 3], [6502, 1429, 4], [4325, 1659, 5], [1093, 1840, 6], [975, 3127, 7], [3979, 2832, 8], [6633, 2836, 9]] #右片采样点像素坐标[y,x,点号]
A = np.array([[1/dx, 0, u0], [0, 1/dy, v0], [0, 0, 1]])
for arr1 in uvs1:
b = np.transpose(np.array([[arr1[0], arr1[1], 1]]))
x = np.linalg.solve(A, b)
xArrs1.append(x)
for arr2 in uvs2:
b = np.transpose(np.array([[arr2[0], arr2[1], 1]]))
x = np.linalg.solve(A, b)
xArrs2.append(x)
- 像平面转像空间辅助坐标
spaceArr1 = [] # 左像空间辅助坐标集合
spaceArr2 = [] # 右像空间辅助坐标集合
# 计算像空间辅助坐标
def computerPixelAssist(r, w, k):
f = 35.0 # 焦距（mm）
global spaceArr1,spaceArr2
spaceArr1, spaceArr2 = [], [] # 清空数组
if len(xArrs1) ==0 | len(xArrs2) == 0:
return
# 左相片做换为像空坐标
for arr1 in xArrs1:
# 交换x y轴 再取负
t = np.copy(arr1[0])
arr1[0] = -arr1[1]
arr1[1] = -t
arr1[2] = -f
spaceArr1.append(arr1)
#计算右片旋转矩阵
# r, w, k = 0, 0, 0
a1 = np.cos(r)* np.cos(k) - np.sin(r) * np.sin(w) * np.sin(k)
a2 = -np.cos(r)* np.sin(k) - np.sin(r) * np.sin(w) * np.cos(k)
a3 = -np.sin(r)* np.cos(w)
b1 = np.cos(w)* np.sin(k)
b2 = np.cos(w)* np.cos(k)
b3 = -np.sin(w)
c1 = np.sin(r) * np.cos(k) + np.cos(r) * np.sin(w) * np.sin(k)
c2 = -np.sin(r) * np.sin(k) + np.cos(r) * np.sin(w) * np.cos(k)
c3 = np.cos(r) * np.cos(w)
R = np.array([[a1,b1,c1], [a2,b2,c2], [a3,b3,c3]]) #旋转矩阵
# 计算右片像空间辅助坐标
for arr2 in xArrs2:
# 交换x y轴 再取负
t = np.copy(arr2[0])
arr2[0] = -arr2[1]
arr2[1] = -t
arr2[2] = -f
b = arr2
x = np.linalg.solve(R, b)
spaceArr2.append(x)
3 粗差检测、计算各个相对定向参数、相对定向解算、迭代计算、精度评定
采用连续法相对定向
相对定向、共面方程原理可参照：
摄测与CV 3: 共面方程，相对定向，绝对定向
- 粗差处理 数据探测法假设检验粗差，进行了粗差检验，并没有处理只是输出第几对像点存在粗差，粗差改正可后续添加
注：数据探测法假设检验参照误差理论与测量平差基础
def findCrossError():
global X
if (len(A) == 0) | (len(V) == 0):
return
q = 1 # 单位权中误差
Q = np.eye(9,9) # 9*9的协因数阵
Qvv = Q - (A * (A.T * A).I * A.T)
t = 0.05 # 正态分布的阿尔法值
ut = 1.96 # 正态分布表里的u0.025的限值
i = 0
for arr in V:
u = arr / (q * np.sqrt(Qvv[i,i]))
if np.abs(u) > ut:
print("第{}个点存在粗差,粗差为：{}".format(i,arr))
# 粗差处理
# X = np.zeros(shape=(5,1))# 舍去改正值
# print("已改正...")
else:
if i == len(V)-1: # 循环结束
print("未发现粗差...")
i += 1
- 各相对定向参数结算
解算参数和过程参照摄影测量学
X = 0 # A矩阵 L由观测值计算得到的q值 V五个参数的改正值 X真值
V = 0 # 存放改正值
A = 0 # 系数矩阵
def computerNewValue():
global X, V, A
L = 0
i = 0
while i < len(QArr):
x2, y2, z2, n1 = spaceArr2[i][0], spaceArr2[i][1], spaceArr2[i][2], N1Arr[i]
a1 = - (x2 * y2 / z2) * n1
a2 = - (z2 + y2 * y2 / z2) * n1
a3 = x2 * n1
a4 = Bx
a5 = - (y2 / z2) * Bx
a = np.array([[a1[0], a2[0], a3[0], a4[0], a5[0]]])
if i == 0:
A = a
else:
A = np.row_stack((A, a))
i += 1
# print(len(A))
A = np.mat(A) # 转换成矩阵
L = np.transpose(np.array([QArr]))
L = np.mat(L)
X = (A.T.dot(A)).I.dot(A.T).dot(L) # A.T求转置 A.I求逆
V = (A.dot(X) - L)
- 入口函数、迭代计算和精度评定
由于得到的未知数改正值有时无法收敛到限值，未找到原因，故没有以改正值为迭代的收敛条件，而是以迭代次数为收敛条件
if __name__ == '__main__':
r, w, k = 0, 0, 0 # 起始值
u, v = 0, 0
# V = 0 #存放差值
rDiffi, wDiffi, kDiffi, uDiffi, vDiffi = 0,0,0,0,0 #各差值
pixelToImgPoint() # 像素转像平面坐标
def computer():
global r,w,k,u,v,V
global rDiffi,wDiffi,kDiffi,uDiffi,vDiffi
computerPixelAssist(r, w, k) # 转像空间辅助坐标
computerParameter(u, v) # 计算各参数
computerNewValue() #计算新值
# 计算差值
# V = X - np.mat([[r,w,k,u,v]]).T
# print(X)
rDiffi = X[0, 0] # 各个差值
wDiffi = X[1, 0]
kDiffi = X[2, 0]
uDiffi = X[3, 0]
vDiffi = X[4, 0]
computer()
findCrossError()
# 迭代计算，判断是否满足要求
limit = 0.08 # 限值
flag = True # 是否继续迭代的指示
num = 1 # 迭代次数
# r, w, k, u, v = V[0, 0], V[1, 0], V[2, 0], V[3, 0], V[4, 0]
# computer()
while num < 100: #迭代多少次停止
# if (np.abs(rDiffi) > limit) | (np.abs(wDiffi) > limit) | (np.abs(kDiffi) > limit) | (np.abs(uDiffi) > limit)| (np.abs(vDiffi) > limit): # 超出限值
print("已迭代：%d次..."%num)
r,w,k,u,v = rDiffi+r,wDiffi+w,kDiffi+k,uDiffi+u,vDiffi+ v # 参数改正
precision = np.sqrt((V.T.dot(V))/4)[0,0]
print("精度：{}毫米".format(precision))
computer()
num += 1
print("已迭代：%d次..." % num)
# else:
# flag = False
# break
print("================================================")
print("改正值：r:{}, w:{}, k:{}, u:{}, v:{}".format(rDiffi, wDiffi, kDiffi, uDiffi, vDiffi))
print("最终精度：{}".format(precision))
print("改正后的值：r:{}, w:{}, k:{}, u:{}, v:{}".format(r,w,k,u,v))
运行输出结果如下：
迭代次数、每次迭代后的精度、结果的改正值、结果的精度、相对定向的5参数
到此相对定向解算完成
4 特征点提取
采用Moravex算子进行特征点提取，提取过程如下：
提取左影像的特征点、保存左影像特征点像素坐标、提取部分特征点进项右影像匹配、保存左右影像同名特征点
- 兴趣值计算
兴趣值计算比较耗时，需要计算接近7360*4912个的像素的兴趣值
# 计算感兴趣值 mat:灰度矩阵 fullSize窗口大小 x中心像元行号 y列号
def getInterestValue(mat, fullSize, x, y):
halfSize = fullSize // 2
i = 0
v1, v2,v3,v4 = 0,0,0,0
while i < halfSize * 2 :
v1 += np.power((int(mat[x-halfSize+i,y]) - int(mat[x-halfSize+i+1,y])),2)
v2 += np.power((int(mat[x-halfSize+i,y-halfSize+i]) - int(mat[x-halfSize+i+1,y-halfSize+i+1])),2)
v3 += np.power((int(mat[x,y-halfSize+i]) - int(mat[x,y-halfSize+i+1])),2)
v4 += np.power((int(mat[x-halfSize+i,y-halfSize-i]) - int(mat[x-halfSize+i+1,y-halfSize-i-1])),2)
i += 1
v = [v1, v2, v3, v4]
interest = np.min(v)
return interest
inputimagepath = "./image/left.JPG"
def gray_cvt(inputimagepath): # 入口函数
img = cv2.imread(inputimagepath)
gray = cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)# 灰度化 gray_cvt_image为灰度矩阵宽*高
# cv2.namedWindow('gray') # 这行没啥用 控制显示图片窗口的名字
h, w = gray.shape
# Moravex算子
grayMat = np.mat(gray)
interests = np.zeros((h,w)) #兴趣值数组
fullSize = 9 #窗口大小
inhibitionsize = 17 #抑制局部最大的窗口大小
halfSize = fullSize // 2
halfInhiSize = inhibitionsize // 2
i, j = halfSize, halfSize
sum = 0
while i < h-halfSize: # 行列数搞错了，计算结果出错 高度为行数
print("兴趣值计算到{}行".format(i+1))
j = halfSize # 找半天这个bug 还是太菜了
while j < w-halfSize:
interest = getInterestValue(grayMat,fullSize,i,j)
interests[i][j] = interest #bug 若直接拷贝原数组会部分赋值失败，可能是像素超过255便不赋值
sum += interest
j +=1
i+=1
print("所有兴趣值计算完毕...")
- 确定阈值、阈值过滤
滤值是由所有兴趣值的均值确定的
注：以下代码在在主函数内
i, j = halfSize, halfSize
mean = sum / ((w-fullSize )* ( h-fullSize)) #阈值
print("阈值{}",format(mean))
while i < h-halfSize:
print("阈值过滤到{}行".format(i+1))
j = halfSize
while j < w-halfSize:
value = interests[i][j]
if value < mean:
interests[i][j] = 0
j +=1
i+=1
print("阈值过滤完毕...")
- 抑制局部最大
# 抑制局部最大
i, j = halfSize, halfSize
while i < h-halfInhiSize:
print("抑制局部最大到{}行".format(i+1))
j = halfInhiSize
while j < w-halfInhiSize:
max = interests[i][j]
m, n = 0, 0
while m < inhibitionsize:
n = 0
while n < inhibitionsize:
value = interests[i - halfInhiSize + m][j - halfInhiSize + n]
if value < max:
interests[i - halfInhiSize + m][j - halfInhiSize + n] = 0
n += 1
m += 1
j +=1
i+=1
print("抑制局部非最大完毕...")
- 最后保存结果，结果为特征点的像素坐标和兴趣值
# 保存结果
i, j = halfSize, halfSize
xs, ys, grayValue = [],[],[]
while i < h - halfSize:
j = halfSize
while j < w - halfSize:
value = interests[i][j]
if value > 0:
xs.append(j)
ys.append(x)
grayValue.append(value)
j += 1
i += 1
# 字典中的key值即为csv中列名
dataframe = pd.DataFrame({'X坐标': xs, 'Y坐标': ys ,'兴趣值': grayValue})
# 将DataFrame存储为csv,index表示是否显示行名，default=True
dataframe.to_csv("featurePoints.csv", index=False, sep=',')
print("所有特征点写入完毕...")
cv2.waitKey()#等待操作
cv2.destroyAllWindows()#关闭显示图像的窗口
- 运行结果，featurePoints.csv文件：
注：运行时没有写入兴趣值
5 影像匹配
相关系数匹配特征点
思路是：读取左影像部分特征点，获得特征点特定大小（这里选取5*5像素）的窗口（目标窗口），获得右影像灰度矩阵（搜素窗口），右影像矩阵划分成目标窗口大小的矩阵，选定目标窗口逐个对搜素窗口划分后的小窗口进行系数计算，求出最大的系数，该系数对应的小窗口的中心像素即为匹配得到的右影像特征点坐标，最后写入文件
- 读取影像函数
def readImage(imageName):
img = cv2.imread(imageName)
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) # 灰度化 gray_cvt_image为灰度矩阵宽*高
# h, w = gray.shape
cv2.waitKey() # 等待操作
cv2.destroyAllWindows() # 关闭显示图像的窗口
return gray
- 读写文件函数
def readFile():
xys = []
df = pd.read_csv(CSV_FILE_PATH)
for i in range(10000,100000,5000): #间隔取点
x = df["X坐标"][i]
y = df["Y坐标"][i]
xys.append([x,y])
return xys
def writeFile(xs1,ys1,xs2,ys2):
dataframe = pd.DataFrame({'左影像X坐标': xs1, '左影像Y坐标': ys1,'匹配的右影像X坐标': xs2, '匹配的右影像Y坐标': ys2})
dataframe.to_csv("featurePoints_matching.csv", index=False, sep=',')
print("所有匹配特征点写入完毕...")
- 划分数组
# 划分数组
def divideArr(arr, rows, cloumns):
h, w = arr.shape
rowsNum, colNum = [],[]
i, j = 0,0
numRow = h // rows #分成几等行
numClo = w // cloumns #分成几等列
sigleXYS = [] # 单个搜索区窗口列表
# 拆分行
while i < numRow * rows:
i += rows
rowsNum.append(i)
sigleX = np.split(arr, rowsNum, axis=0)
finallymat = sigleX[len(sigleX)- 1]
gao1 ,kuan1 = finallymat.shape
if gao1 < rows:
del sigleX[len(sigleX)- 1]
t = 0
while t < len(sigleX):
# 拆分列
while j < numClo * cloumns:
j += cloumns
colNum.append(j)
sigleY = np.split(sigleX[t], colNum, axis=1)
finallymat2 = sigleY[len(sigleY) - 1]
gao1, kuan1 = finallymat2.shape
if kuan1 < cloumns:
del sigleY[len(sigleY) - 1]
sigleXYS.append(sigleY)
t += 1
return sigleXYS
4.求相关系数
# 求相关系数
def matching(mat1, mat2):
m,n = mat1.shape
i,j = 0,0
fz1,fz2,fz3,fm1,fm2 = 0,0,0,0,0
while i <len(mat1):
j = 0
while j < len(mat2):
a,b = int(mat1[i][j]), int(mat2[i][j])
fz1 += a * b
fz2 += a
fz3 += b
fm1 += a *a
fm2 += b *b
j += 1
i += 1
fzp = fz1 - (1/(m* n)) *fz2 *fz3
fmp1 = fm1 - (1/(m*n))* fz2 *fz2
fmp2 = fm2 - (1/(m*n)* fz3 *fz3)
if fmp1* fmp2 == 0: #除数太小 python会看成0 报错
return 0
p = fzp / (np.sqrt(fmp1* fmp2)) #相关系数
return p
- 入口函数
def computerValue(g,g1):
h, w = g1.shape
m, n,k,l =5, 5,h, w # m n为目标窗口大小 k l 为搜索影像区域
r,c = m // 2, n//2
i,j = 0,0
sigleXYS = divideArr(g1,m, n) #982*1472
finalCol = len(sigleXYS)
finalRow = len(sigleXYS[0])
xys = readFile()
if len(xys) == 0:
print("读取文件失败...")
pipeiNum = 0
xbefore, ybefore, xmatch, ymatch = [], [],[],[] # 匹配的坐标
for xy in xys:
x,y = xy[0],xy[1]
xbefore.append(x)
ybefore.append(y)
E = range(x-r, x+r+1) # 定义列数
F = range(y-c, y+c+1) # 定义行数
goalMat = g[F] # 取出需要的行
goalMat = goalMat[:,E] # 取出需要的列 X_3为目标矩阵
ps = [] # 相关系数集合
i = 0
while i < finalCol:
j = 0
while j < finalRow:
sigleMat = sigleXYS[i][j]
p = matching(goalMat,sigleMat)
ps.append(p)
j +=1
i += 1
pmax = np.max(ps)
maxIndex = ps.index(pmax)
goalRow = (maxIndex) // finalCol
goalCol = (maxIndex) % finalCol
goalRow = goalRow * m - c
goalCol = goalCol * n - r
xmatch.append(goalRow)
ymatch.append(goalCol)
pipeiNum += 1
print("已匹配:{}个像点".format(pipeiNum))
writeFile(xbefore,ybefore,xmatch,ymatch) #写入文件
- featurePoints_matching.csv内容
6 总结
- 相对定向的收敛条件无法满足，可能的原因是：1.相对定向计算错误，本人检查了几遍，没有发现错误之处。2. 手动选取的同名像点坐标不准确。3.应该不是粗差的问题，假设检验后提示未发现粗差。
- 特征点的计算出的阈值太小，导致提取的特征点非常多，可按照经验值重新设置。
- 影像匹配的目标窗口和搜索窗口有时候并没有重叠，导致匹配出的特征点不准确，同时，影像匹配后的同名点存在不准确的情况，有待进一步改正。
- 计算效率和准确度都没有达到预期效果，不能算是成功完成，只能算是学习的过程。
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
