master

分支 (1)

管理

管理

master

opencv
/
_8depth.py

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 19 19:31:12 2018

@author: lenovo
"""
'''
深度估计、双目匹配
'''

import numpy as np
import cv2

'''
1、使用普通摄像机进行深度估计
真实场景的双目立体匹配（Stereo Matching）获取深度图详解：https://www.cnblogs.com/riddick/p/8486223.html
双目立体匹配一直是双目视觉的研究热点，双目相机拍摄同一场景的左、右两幅视点图像，运用立体匹配匹配算法获取视
差图，进而获取深度图。而深度图的应用范围非常广泛，由于其能够记录场景中物体距离摄像机的距离，可以用以测量、
三维重建、以及虚拟视点的合成等。

https://blog.csdn.net/qq_37059483/article/details/79481014
https://blog.csdn.net/wangxiaokun671903/article/details/37935113
https://blog.csdn.net/ccblogger/article/details/72900316
https://blog.csdn.net/qq_37059483/article/details/79481014

双目立体视觉匹配算法-----SAD匹配算法、BM算法、SGBM算法、GC算法：https://blog.csdn.net/liulina603/article/details/53302168
'''

'''
什么是视差，假设两个相机平行，相距d，且两个相机成像大小一样
假设A相机拍出图片为imgL
假设B相机拍出图片为imgR

由于两相机平行，imgR可以看做imgL的平移，则同一物体在imgL的成像位置和imgR的成像的位置差就是视差
'''


'''
1、双目立体匹配SAD算法 (Sum of absolute differences)
https://www.cnblogs.com/Crazy-Dog123/articles/5041950.html
http://www.cnblogs.com/Crazy-Dog123/articles/5043864.html
'''

def sad(imgL,imgR):
    '''
    第一步：设定一个窗口D对图像进行邻域求和，然后遍历整个图像，这样，我们就得到了两个计算了邻域值的两个
    矩阵，里面保存了矩阵的灰度值求和
    输入：两个m*n的灰度图像
    输出：两个m*n的灰度矩阵，这两个矩阵分别保存了两个图像在D邻域下的灰度和
    算法： 用D窗口实现对每一个像素的滚动操作，最终生成矩阵L1和R1
    '''
    #高、宽
    m = imgL.shape[0]
    n =  imgL.shape[1]

    #滑动窗口的大小
    dsize = 8
    #设置最大的视差，30个像素
    S = 30

    #窗口大小的一半
    N = int(dsize/2)

    sumL = np.zeros((m,n))
    sumR = np.zeros((m,n))

    for i in range(N,m-N):       #[N,m-N)
        for j in range(N,n-N):
            #对D窗口内的像素值进行求和
            sumL[i,j] = np.sum(imgL[i-N:i+N,j-N:j+N])
            sumR[i,j] = np.sum(imgR[i-N:i+N,j-N:j+N])

    '''
    第二步：搜索和匹配，在用imgL的矩阵值剪掉imgR的矩阵值，并求取绝对值，在视差偏移量最大的范围内通过比较
    矩阵imgL某个值与范围内imgR的值，并记录下最小绝对值与偏移量
    这里最小的便宜量对应了偏移量就是对应了视差图的灰度值，这样生成视差图矩阵S1
    输入：imgL矩阵，imgR矩阵
    输出：disparity矩阵
    算法：在视差偏移量允许范围内，找到最佳匹配点，可以等效于求取最小代价函数scale，找到最小的数值点便
    确定了偏移量
    '''
    #视差图
    disparity = np.zeros((m,n))

    for i in range(N,m-N):
        for j in range(N,n-N-S):
            sad_sum = 1e6
            scale = 0.0
            dis = 0
            #遍历视差
            for k in range(0,S):
                #做差  水平坐差，主要是因为相机经过标定校正后，匹配点都在同意行上
                scale = np.abs(sumL[i,j] - sumR[i,j+k])
                if scale < sad_sum:
                    #寻找最小的匹配代价函数
                    sad_sum = scale
                    #记录最小的能量位置
                    dis = k
            #保存
            disparity[i,j] = dis
    '''
    第三步：后期处理和显示
    '''
    disparity = disparity*255/S
    dispmap = np.uint8(disparity)
    cv2.imshow('sad',dispmap)


from operator import itemgetter, attrgetter

def sad_improve(imgL,imgR):
    '''
    sad函数中两个图如果直接求和会造成两个图的数值过分大，数值太大一方面会影响运算的速度，
    另外一方面影响求和的精度，所以建议先实现图像的相减操作，然后再求和，同时保存记录好偏移量，
    这样数值就不会很大了
    '''
    '''


    输入：两幅图像，一幅Left-Image，一幅Right-Image

    对左图，依次扫描，选定一个锚点：

    （1）构造一个小窗口,类似于卷积核；
    （2）用窗口覆盖左边的图像,选择出窗口覆盖区域内的所有像素点；
    （3）同样用窗口覆盖右边的图像并选择出覆盖区域的像素点；
    （4）左边覆盖区域减去右边覆盖区域，并求出所有像素点灰度差的绝对值之和；
    （5）移动右边图像的窗口，重复（3）-（4）的处理（这里有个搜索范围,超过这个范围跳出）；
    （6）找到这个范围内SAD值最小的窗口,即找到了左图锚点的最佳匹配的像素块。
    '''
    #高和宽
    H = imgL.shape[0]
    W =  imgL.shape[1]

    #转换为灰度图
    if imgL.shape[2] == 3:
        imgL = cv2.cvtColor(imgL,cv2.COLOR_BGR2GRAY)

    if imgR.shape[2] == 3:
        imgR = cv2.cvtColor(imgR,cv2.COLOR_BGR2GRAY)

    imgL = np.double(imgL)
    imgR = np.double(imgR)

    #滑动窗口的大小的一半
    N = 9
    #最大视差
    D = 20

    #计算右图减去左图，相减产生D个矩形放到imgDiff中
    imgDiff = np.zeros((H,W,D))
    e = np.zeros((H,W))


    for i in range(D):
        print(i)
        #右图 - 左图
        e[:,0:W-i] = np.abs(imgR[:,0:W-i] - imgL[:,i:W])

        #计算窗口内的和
        e2 = np.zeros((H,W))
        for y in range(N,H-N):
            for x in range(N,W-N):
                #对D窗口内的像素值进行求和
                e2[y,x] = np.sum(e[y-N:y+N,x-N:x+N])
        imgDiff[:,:,i] = e2

    #找到最小的视差，得到dispMap
    dispMap = np.zeros((H,W))
    for y in range(H):
        for x in range(W):
            sort_data = sorted(enumerate(list(imgDiff[y,x,:])),key=itemgetter(1))
            #if np.abs(sort_data[0][1]-sort_data[1][1]) > 10:
            dispMap[y,x] = sort_data[0][0]

    dispMap = dispMap*255/D
    dispmap = np.uint8(dispMap)
    cv2.imshow('sad_improve',dispmap)


'''
2、双目匹配值SGBM算法 (Semi-global matching) 半全局匹配算法
http://www.cnblogs.com/Crazy-Dog123/articles/5057179.html
'''

min_disp = 16
num_disp = 192 - min_disp
def update(stereo,val = 0):
    stereo.setBlockSize(cv2.getTrackbarPos('window_size','disparity'))
    stereo.setUniquenessRatio(cv2.getTrackbarPos('uniquenessRatio','disparity'))
    stereo.setSpeckleWindowSize(cv2.getTrackbarPos('speckleWindowSize','disparity'))
    stereo.setSpeckleRange(cv2.getTrackbarPos('speckleRange','disparity'))
    stereo.setDisp12MaxDiff(cv2.getTrackbarPos('disp12MaxDiff','disparity'))
    print('computing disparity...')

    #得到一个视差图
    disp = stereo.compute(imgL,imgR).astype(np.float32)/16.0

    #cv2.imshow('left',imgL)  归一化到0-1之间
    cv2.imshow('disparity',(disp - min_disp)/num_disp)

def sgbm(imgL,imgR):
    window_size = 2
    #blockSize = window_size
    uniquenessRatio = 1
    speckleRange = 3
    speckleWindowSize =3
    disp12MaxDiff = 200
    P1 = 600
    P2 = 2400

    #创建几个跟踪条来调整算法参数，然后调用update参数
    cv2.namedWindow('disparity')
    #cv2.createTrackbar('speckleRange','disparity',speckleRange,50,update)
    #cv2.createTrackbar('window_size','disparity',window_size,21,update)
    #cv2.createTrackbar('speckleWindowSize','disparity',speckleWindowSize,200,update)
    #cv2.createTrackbar('uniquenessRatio','disparity',uniquenessRatio,50,update)
    #cv2.createTrackbar('disp12MaxDiff','disparity',disp12MaxDiff,250,update)

    #创建一个StereoSGBG实例，这是一个计算视差图的算法
    stereo = cv2.StereoSGBM_create(
            minDisparity = min_disp,        #最小的视差，一般设置为0
            numDisparities = num_disp,      #最大的视差，这个需要自己定义，这个数值比0大，而且要被16整除
            blockSize = window_size,        #滑动窗口的大小，必须要大于1，一般取3-11之间
            uniquenessRatio = uniquenessRatio,  #最好的代价方程值“赢了”第二好的代价方程值的概率，通常设置为5-15之间效果达到最佳
            speckleRange = speckleRange,        #相邻像素点的视差值浮动范围，通常设置为1-2就好了，这个系数会被乘以16输入到程序中
            speckleWindowSize = speckleWindowSize,  #针对散斑滤波的窗口大小，如果设置为0则不允许散斑滤波，否则设置为50-200之间
            disp12MaxDiff = disp12MaxDiff,      #视差图的像素点检查，设置为1
            P1 = P1,   #惩罚系数，用于控制图像的平滑，具体需要看这个函数的算法，很多博客有介绍P1和P2的惩罚系数
            P2 = P2)

    update(stereo)


if __name__ == '__main__':
    #加载两张图
    imgL = cv2.imread('image/11.bmp')
    #imgL = cv2.resize(imgL,None,fx=0.2,fy=0.2,interpolation=cv2.INTER_CUBIC)
    imgR = cv2.imread('image/22.bmp')
    #imgR = cv2.resize(imgR,None,fx=0.2,fy=0.2,interpolation=cv2.INTER_CUBIC)
    #sad(imgL,imgR)
    #sad_improve(imgL,imgR)
    sgbm(imgL,imgR)

    #3X3
    R = [[9.9998525621483181e-01, 3.5314113921735158e-03,
       4.1251044273392700e-03], [-3.5020112874907588e-03,
       9.9996857200336808e-01, -7.1127296087366331e-03],
       [-4.1500927579409897e-03, 7.0981785779126857e-03,
       9.9996619572412371e-01 ]]

    #3X3
    R1 = [[9.9995544861535535e-01, -7.5529244715175194e-03,
       5.6616354872873727e-03],[ 7.5729192221467706e-03,
       9.9996513477129034e-01, -3.5185417846053421e-03],
       [-5.6348628127222394e-03, 3.5612601369073723e-03,
       9.9997778262685366e-01]]

    #3X3
    R2 = [[9.9993738784113484e-01, -1.1094812048627922e-02,
       1.4579242275623966e-03], [1.1089581485058921e-02,
       9.9993221385813869e-01, 3.5480799382088197e-03],
       [-1.4971906805518226e-03, -3.5316900157434700e-03,
       9.9999264276578492e-01]]

    #3X4
    P1 = [[4.4828811368916126e+02, 0., 3.3378979110717773e+02, 0.],[ 0.,
       4.4828811368916126e+02, 2.4358526802062988e+02, 0.],[ 0., 0., 1.,
       0.]]

    P2 = [ [4.4828811368916126e+02, 0., 3.3378979110717773e+02,
       -2.9831681553297447e+04], [0., 4.4828811368916126e+02,
       2.4358526802062988e+02, 0.], [0., 0., 1., 0. ]]

    T = [[-6.6541612004454194e+01],[ 7.3831290596785981e-01],
       [-9.7018702833071249e-02]]


    #进行验证
    P_2 = np.dot(R,P1) + T
    print(P_2)

    cv2.waitKey()
    cv2.destroyAllWindows()