基于灰度投影法的图片偏移测量

目录

前言

一、灰度投影法简介 

二、代码实现

总结

---

前言

前几天导师给了个需求，要求给出经过运动放大后的视频和原视频对应帧图片的整体偏移量，上网搜索，似乎光流法非常流行，但是还有一种灰度投影法看大家评论也很好，所以就在这里写下自己对该算法的理解和感悟，并且附上自己的代码，有劳各位指正，也希望能对大家有所帮助，共同进步！

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、灰度投影法简介 

如何定义图片整体的位移？最理想的情况是，图片上的每个像素点都统一的变化到另一个位置（类似于理论力学中刚体的平动），不妨假设图片整体向左移动dw，向上移动dh，那么对任意（h,w）∈图片像素范围，新图片的（h,w）像素点和原图片的（h+dh,w+dw）应该完全一致。也即新图片第h行和原图片第h+dh行应该完全一样，新图片第w列和原图片第w+dw列应该完全一样。

让我们把问题变得稍微复杂一些，图片上可能各行各列的像素点位移不一致，尽管他们在大方向上具有同样的位移趋势，如何利用第一段里的观察到的现象，去度量这种位移呢？我们给出行灰度投影与列灰度投影的公式如下：

 

上式中，Gk(i,j)代表第k帧图像在(i,j)处的像素灰度值，Gk(i)代表了图像第i行的像素灰度值，Gk(j)代表图像第j列的像素灰度值。

---

下面所有公式，我们都仅给出行的例子，列的公式可简单类比不再赘述。

下面计算行灰度投影值，假设图像大小为M*N：

Mean代表第k帧图像的行灰度均值，rProj(i)为第i行的行灰度投影。个人猜测，这是利用了线性代数对于投影空间的定义，Gk(i)为向量，而Mean为该向量在投影空间的补空间分量，两者相减得出投影量。当然上述论述我暂时没有想到什么数学依据，权当一个记忆方式吧。

现在我们有了图像行与列的灰度投影，他们能够反映当前行与列的特征，然而我们现在需要找到一种度量方式，他至少要满足以下的需求：一、他能够反映图像整体对于在h和w方向上的变化情况；二、他的极值点能被拿来当作，或用来求得图像整体的偏移量，且偏移量应具有正负。

基于此，我们给出如下的计算公式：

式中，rProj和rProj0分别代表你想判断的两个图片的行灰度投影。可以看出，这个式子对j进行了求和，并且把图片上h-m个（m的含义之后解释，稍安勿躁）间距为i-m/2的行的灰度投影值的差做了平方，再求和。求和这个操作可以在一定程度上满足我们提出的需求一，那么，他的极小值点又有什么意义呢？

假定i0是上述公式的极小值点，定性地讲，整体上看，图像上间隔为i0-m/2的行之间灰度投影值差距很小，我们有理由认为，i0-m/2就是一个对于图像h方向上整体偏移量的很好的估计。所以，我们就认为图像的h方向位移是i0-m/2。

更严谨地数学证明可以去参考一下互相关系数，限于本人知识水平暂无法给出证明，留待日后有时间再做吧！

我们把0≤i≤m带入i-m/2，可以看出偏移量的取值范围是[-m/2,m/2]，这满足了我们对于偏移量应具有正负的要求。而m则代表了本算法对于偏移量的搜索范围。

至此，本算法原理简介完毕

二、代码实现

读取两个一一对应的图片序列，利用灰度投影算法计算出对应帧的位移后存入txt文件中：

import numpy as np
import cv2
import matplotlib.pyplot as plt
import os

#返回灰度投影值
def calculate(img):
    h,w=img.shape

    #先做行向量的投影
    rsum=np.sum(img,axis=1)
    rsum_mean=np.sum(rsum)/h
    rProj=rsum-rsum_mean
    
    #再做列向量的投影
    csum=np.sum(img,axis=0)
    csum_mean=np.sum(csum)/w
    cProj=csum-csum_mean
    
    return rProj,cProj

#利用灰度投影计算出位移
def compare(rProj,rProj0,cProj,cProj0,img,idx):
    h,w=img.shape
    
    R=np.zeros(idx)
    C=np.zeros(idx)
    
    for i in range(idx):
        for j in range(h-idx):
            R[i]=R[i]+(rProj[i+j]-rProj0[int(idx/2)+j])*(rProj[i+j]-rProj0[int(idx/2)+j])
    dh=returnindex(R)-idx/2
    # print(R)
    
    for i in range(idx):
        for j in range(w-idx):
            C[i]=C[i]+(cProj[i+j]-cProj0[int(idx/2)+j])*(cProj[i+j]-cProj0[int(idx/2)+j])
    dw=returnindex(C)-idx/2
    # print(C)
    
    return dh,dw

#这里我重新定义了一个去最值函数，返回数组中最大值对应的最大下标
#不太清楚np.argmax()是否能有同样的效果
def returnindex(A):
    temp=A[0]
    global index
    for i in range(len(A)):
        if A[i]<=temp:
            temp=A[i]
            index=i
    return index

img0=cv2.imread("···")#这里改成你要读取的第一个图片序列
img0=cv2.cvtColor(img0,cv2.COLOR_BGR2GRAY)

img=cv2.imread("···")#这里改成你要读取的第二个图片序列
img=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

rProj,cProj=calculate(img)
rProj0,cProj0=calculate(img0)
dh,dw=compare(rProj,rProj0,cProj,cProj0,img,500)#这里的500就是算法介绍中的m值，是搜索范围

length=(dh*dh+dw*dw)**0.5
print(dh)
print(dw)
print((dh*dh+dw*dw)**0.5)

#这之后的代码是实际操作过程，因为我在实验室服务器上进行操作
#所以我将水平位移，竖直位移和总位移分别按行记在了三个txt文件中
#要用的时候我就读文件，然后用python画图

f1=open('···','w')
f2=open('···','w')
f3=open('···','w')
root_img='···'
root_changeimg='···'
filelist=os.listdir(root_img)
filelist.sort()
for file in filelist:
    
    img0=cv2.imread(root_img+file)
    img0=cv2.cvtColor(img0,cv2.COLOR_BGR2GRAY)
    img=cv2.imread(root_changeimg+file)
    img=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    rProj,cProj,rsum,csum=calculate(img)
    rProj0,cProj0,rsum0,csum0=calculate(img0)
    dh,dw=compare(rProj,rProj0,cProj,cProj0,img,500)
    height_list.append(dh)
    weight_list.append(dw)
    length=(dh*dh+dw*dw)**0.5
    total_list.append(length)
    f1.write(str(length)+'\n')
    f2.write(str(dh)+'\n')
    f3.write(str(dw)+'\n')
    print((dh*dh+dw*dw)**0.5)
    

f1.close()
f2.close()
f3.close()

---

总结

灰度投影法是一个相对来说比较容易理解，而且较易实现的算法。本文中是它最粗浅的实现，各位也可以根据自己的需求做出自己的改进，相信能取得更好的效果。有时间的话我也把光流法看一下，看能不能再水一篇文章