接OpenCV计算机视觉整理

图像的分割与修复

图像分割的基本概念

• 什么是图像分割

将前景物体从背景中分离出来。

• 图像分割的方法

之前的很多方法都是图像分割的前置步骤，比如腐蚀、膨胀、二值化等等。图像分割方法又分为传统的图像分割方法和基于深度学习的图像分割方法。

• 传统的图像分割方法

1. 分水岭法
2. GrabCut方法
3. MeanShift法
4. 背景扣除

• 分水岭法原理

上图表示图像有一定的梯度，0代表黑色，代表比较低洼的地方，白色是255，代表一个峰点。当我们使用不同的颜色向低洼处灌水。

当两个低洼处的水要打通的时候，此时这个颜色就会产生一个冲突。在冲突点的地方设置一个边界点，这样就将不同的区域给分割开来。但它也是有问题的。

当图像中存在过多的极小区域而产生许多小的集水盆，但实际上我们真正想要的是一个大块，视为一个整体。但通过分水岭法就可能分成很多的小块。对于传统的分水岭法可能会把一张图切割的很碎，不利于我们后面的处理。

但幸运的是，OpenCV的分水岭法可以将一大块分成一整块。

• 分水岭法的处理步骤

1. 标记背景
2. 标记前景
3. 标记未知域
4. 进行分割

现在我们要对这张图片进行分割

import cv2
import numpy as np

if __name__ == "__main__":

    cv2.namedWindow('img', cv2.WINDOW_NORMAL)
    img = cv2.imread("/Users/admin/Documents/9527.png")
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    # 二值化，cv2.THRESH_OTSU表示自适应阈值
    ret, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))
    # 进行两次开运算
    open = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)
    # 进行一次膨胀
    bg = cv2.dilate(open, kernel, iterations=1)
    # 获取前景物体，计算非0值到离它距离最近的0值之间的距离
    # cv2.DIST_L2为距离的计算方式，这里为欧式距离
    # maskSize:扫描时的卷积核大小，L1时用3*3，L2时用5*5
    dist = cv2.distanceTransform(open, cv2.DIST_L2, 5)
    ret, fg = cv2.threshold(dist, 0.7 * dist.max(), 255, cv2.THRESH_BINARY)
    fg = np.uint8(fg)
    # 获取未知区域
    unknow = cv2.subtract(bg, fg)
    # 创建连通域，求所有非0相同像素的连通域，它还有一个参数为连通周围的4个还是
    # 8个，默认为8个
    ret, marker = cv2.connectedComponents(fg)
    marker += 1
    marker[unknow == 255] = 0
    # 进行图像分割
    result = cv2.watershed(img, marker)
    img[result == -1] = [0, 0, 255]
    while True:
        cv2.imshow('img', img)
        key = cv2.waitKey()
        if key & 0xFF == ord('q'):
            break
    cv2.destroyAllWindows()

运行结果

GrabCut方法

通过交互的方式获得前景物体。当我们有一张图片的时候可以通过键盘或者鼠标对其进行指定，通过框住我们需要分离的物品，就可以将物品进行分离。

• GrabCut原理

1. 用户指定前景的大体区域，剩下的为背景区域。
2. 用户还可以明确指定某些地方为前景或背景。
3. GrabCut采用分段迭代的方法分析前景物体形成模型树。
4. 最后根据权重决定某个像素是前景还是背景。

现在我们使用

来进行图像的分割

import cv2
import numpy as np

class App:

    def __init__(self, img):
        self.img = img
        self.img2 = img.copy()
        self.rect = (0, 0, 0, 0)
        self.flag_rectangle = False
        self.startx = 0
        self.starty = 0
        self.mask = np.zeros(self.img.shape[:2], np.uint8)
        self.output = np.zeros(self.img.shape, np.uint8)

    def onmouse(self, event, x, y, flags, param):
        if event == cv2.EVENT_LBUTTONDOWN:
            self.flag_rectangle = True
            self.startx = x
            self.starty = y
        elif event == cv2.EVENT_LBUTTONUP:
            self.flag_rectangle = False
            cv2.rectangle(self.img, (self.startx, self.starty),
                          (x, y), (0, 0, 255), 5)
            self.rect = (min(self.startx, x), min(self.starty, y),
                         abs(self.startx - x), abs(self.starty - y))
        elif event == cv2.EVENT_MOUSEMOVE:
            if self.flag_rectangle:
                # 每次绘制都在新的图像中绘制
                self.img = self.img2.copy()
                cv2.rectangle(self.img, (self.startx, self.starty),
                                (x, y), (255, 0, 0), 5)

    def run(self):
        cv2.setMouseCallback('img', self.onmouse)
        while True:
            cv2.imshow('img', self.img)
            cv2.imshow('output', self.output)
            key = cv2.waitKey(100)
            if key & 0xFF == ord('q'):
                break
            if key & 0xFF == ord('g'):
                bgdmodel = np.zeros((1, 65), np.float64)
                fgdmodel = np.zeros((1, 65), np.float64)
                # 进行grabCut分割
                # self.mask是分割之后产生的掩码，BGD背景0，FGD前景1,
                # PR_BGD可能是背景2，PR_FGD可能是前景3
                # self.rect为我们选取的区域，
                # bgdmodel, fgdmodel是两个固定值
                # 1是迭代次数，
                # cv2.GC_INIT_WITH_RECT初始化时分离前后景使用的mode
                # 还有一个cv2.GC_INIT_WITH_MASK在后续迭代中使用的mode
                cv2.grabCut(self.img2, self.mask, self.rect,
                            bgdmodel, fgdmodel, 1,
                            cv2.GC_INIT_WITH_RECT)
            mask2 = np.where((self.mask == 1) | (self.mask == 3), 255, 0).astype('uint8')
            self.output = cv2.bitwise_and(self.img2, self.img2, mask=mask2)
        cv2.destroyAllWindows()

if __name__ == "__main__":

    cv2.namedWindow('img', cv2.WINDOW_NORMAL)
    cv2.namedWindow('output', cv2.WINDOW_NORMAL)
    img = cv2.imread("/Users/admin/Documents/111.jpeg")
    app = App(img)
    app.run()

运行结果，当我们绘制区域，并按下g键时

meanshift图像分割

• MeanShift原理

1. 严格来说该方法并不是用来对图像分割的，而是在色彩层面的平滑滤波。
2. 它会中和色彩分布相近的颜色，平滑色彩细节，侵蚀掉面积较小的颜色区域。
3. 它以图像上任一点P为圆心，半径为sp，色彩幅值(色彩相近的范围)为sr进行不断迭代。

我们以这张图为例来看看效果

import cv2

if __name__ == "__main__":

    cv2.namedWindow('img', cv2.WINDOW_NORMAL)
    img = cv2.imread("/Users/admin/Documents/10086.png")
    # 对图像进行meanshift处理，20为半径，30为幅值
    mean_img = cv2.pyrMeanShiftFiltering(img, 20, 30)
    while True:
        cv2.imshow('img', mean_img)
        key = cv2.waitKey()
        if key & 0xFF == ord('q'):
            break
    cv2.destroyAllWindows()

运行结果

这样我们就可以看到树叶和部分红花就给同化了，变成了一个抽象画。不过这样并不能完成图像分割，我们还要做进一步处理。

import cv2

if __name__ == "__main__":

    cv2.namedWindow('img', cv2.WINDOW_NORMAL)
    img = cv2.imread("/Users/admin/Documents/10086.png")
    # 对图像进行meanshift处理，20为半径，30为幅值
    mean_img = cv2.pyrMeanShiftFiltering(img, 20, 30)
    # 对meanshift处理后的图像边缘
    imgcanny = cv2.Canny(mean_img, 150, 300)
    # 获取边缘图像的最外层轮廓
    contours, _ = cv2.findContours(imgcanny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    # 在愿图像中绘制轮廓
    cv2.drawContours(img, contours, -1, (0, 0, 255), 3)
    while True:
        cv2.imshow('img', img)
        key = cv2.waitKey()
        if key & 0xFF == ord('q'):
            break
    cv2.destroyAllWindows()

运行结果

图像修复

这里我们需要准备两张图，一张是带有失真的图片，如

另外一张为与上图中失真部分一模一样的灰度图

现在就可以对失真的图进行修复

import cv2

if __name__ == "__main__":

    cv2.namedWindow('img', cv2.WINDOW_NORMAL)
    img = cv2.imread("/Users/admin/Documents/2122.png")
    mask = cv2.imread("/Users/admin/Documents/2121.png", 0)
    # 图像修复，5为每个点的圆型邻域半径
    # cv2.INPAINT_TELEA计算破损周围像素点的加权平均值
    # 当作它恢复的像素值
    # 还有一个cv2.INPAINT_NS算法较为复杂
    dst = cv2.inpaint(img, mask, 5, cv2.INPAINT_TELEA)
    while True:
        cv2.imshow('img', dst)
        key = cv2.waitKey()
        if key & 0xFF == ord('q'):
            break
    cv2.destroyAllWindows()

运行结果

这是一种比较传统的图像修复技术，而且对抠图的技术要求比较高。只能做一些简单的修复，要做比较大块的修复是不可能的。

标志码

在opencv中有一种标志码如下

它可以用来定位，相机标定以及可以得到该码在相机空间坐标系中的位置。也可以当作棋盘格标定板使用

要使用最新版本的标志码，我们需要做如下操作

pip uninstall opencv-contrib-python
pip uninstall opencv-python
pip install opencv-contrib-python

• 我们来看看如何生成一张简单的标志码

import cv2
import cv2.aruco as aruco

if __name__ == '__main__':

    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)

    # 从标记创建图像
    # 第二个参数是ID号
    # 最后一个参数是总图像大小200*200
    img = aruco.generateImageMarker(aruco_dict, 1, 200)

    cv2.imshow('frame', img)

    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果

这个样子的标志码我们是无法使用扫码工具去扫描的，它得需要一个白边，所以我们需要给该标志码添加白边。

import numpy as np
import cv2
import cv2.aruco as aruco

if __name__ == '__main__':

    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)

    # 添加白色背景
    canvas = (np.ones((300, 300)) * 255).astype(np.uint8)
    # 从标记创建图像
    # 第二个参数是ID号
    # 最后一个参数是总图像大小
    img = aruco.generateImageMarker(aruco_dict, 1, 200)
    canvas[50:250, 50:250] = img

    cv2.imshow('frame', canvas)

    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果

• 生成标志码阵列

import cv2
import cv2.aruco as aruco
import numpy as np

if __name__ == '__main__':

    # 创建网格板，这是一组Aruco标记
    # 以下调用得到一个3宽X4高的标记板
    markersX = 3
    markersY = 4
    gridboard = aruco.GridBoard(
            (markersX, markersY),
            markerLength=0.1,
            markerSeparation=0.01,
            dictionary=aruco.getPredefinedDictionary(aruco.DICT_5X5_50))

    # 从网格板创建图像
    img = gridboard.generateImage(outSize=(988, 1400))
    canvas = (np.ones((1500, 1088)) * 255).astype(np.uint8)
    canvas[50:1450, 50:1038] = img

    cv2.imshow('Gridboard', canvas)

    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果

• 生成棋盘格标定板

import cv2
import cv2.aruco as aruco

if __name__ == '__main__':

    # 创建ChArUco棋盘，这是棋盘设置中的一组Aruco标记
    # 用于校准
    # 以下调用将得到一个由5宽X7高的瓷砖组成的ChArUco板
    squaresX = 6
    squaresY = 8
    gridboard = aruco.CharucoBoard(
            (squaresX, squaresY),
            squareLength=0.04,
            markerLength=0.02,
            dictionary=aruco.getPredefinedDictionary(aruco.DICT_5X5_50))

    # 从网格板创建图像
    img = gridboard.generateImage(outSize=(988, 1400))

    cv2.imshow('Gridboard', img)

    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果

由于棋盘格标定板检测的是棋盘格内角点，所以这里无需给其添加白边。具体标定方式可以参考双目视觉理论篇 中的基于 OpenCV 的相机标定。

• 位姿估计

位姿估计是画出图像位于真实3D空间的x,y,z坐标轴，我们的原始图像如下

首先我们需要先检测出该图像中的标志，由于上图中标志图中上下两行没有白边，故只能检测到中间一行的标志

import cv2
import cv2.aruco as aruco

if __name__ == '__main__':

    img = cv2.imread('IMG_20230222_212055.jpg')
    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)
    # 标记检测器参数
    parameters = aruco.DetectorParameters()
    # 检测标记码的标记角坐标和标记id
    corners, ids, _ = aruco.detectMarkers(img, aruco_dict, parameters=parameters)
    # 画出标记角
    aruco.drawDetectedMarkers(img, corners, ids)
    cv2.imshow('img', img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果，这里我们只画出识别到的部分标记

识别到的标志的中心会有一个id号，如id=1，id=4；并且每一个被识别到的标志会有一个绿色的边框。

import cv2
import cv2.aruco as aruco
import numpy as np
import imutils

if __name__ == '__main__':

    img = cv2.imread('IMG_20230222_212055.jpg')
    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)
    # 标记检测器参数
    parameters = aruco.DetectorParameters()
    # 相机内参矩阵
    camera_matrix = np.array([ [3.3811614662976608e+03, 0., 2.0318859679474272e+03], [0.,
           3.3833367964215813e+03, 1.5347124600270040e+03], [0., 0., 1.] ])
    # 相机畸变系数
    dist_coeff = np.array([ 1.4295776727464021e-01, -5.4377825900331767e-01,
           8.0729484571478852e-04, -2.3554260458140577e-03,
           6.3943687964321116e-01 ])
    # 检测标记码的标记角坐标(上下左右四个点的坐标)和标记id
    corners, ids, _ = aruco.detectMarkers(img, aruco_dict, parameters=parameters)
    # 获取相机坐标系下的旋转向量和平移向量
    rvecs, tvecs, _ = aruco.estimatePoseSingleMarkers(corners, 0.05, camera_matrix, dist_coeff)
    # 画出标记角
    aruco.drawDetectedMarkers(img, corners, ids)

    for i in range(rvecs.shape[0]):
        # 绘制3D坐标轴
        cv2.drawFrameAxes(img, camera_matrix, dist_coeff, rvecs[i], tvecs[i], 0.1)
        # 罗德里格斯变换，由旋转向量转换为旋转矩阵
        R, _ = cv2.Rodrigues(rvecs[i])
        # 获取平移向量
        T = tvecs[i]
        # 将旋转矩阵与平移向量组合成位姿矩阵
        pose_mat = np.concatenate((R, T.T), axis=1)
        # 将位姿矩阵扩大一个维度，增加一行
        pose_mat = np.concatenate((pose_mat, np.array([[0, 0, 0, 1]])), axis=0)
        # 获取该矩阵的逆
        pose_mat_inv = np.linalg.inv(pose_mat)
        # 获取标记的位姿
        marker_pos = np.dot(pose_mat, np.array([0, 0, 0, 1]))
        t = tvecs[i]
        r = rvecs[i]
        # 将3D标记中心点坐标转换成2D图像坐标
        res = cv2.projectPoints(np.array([[0.0, 0.0, 0.0]]), r, t, camera_matrix, dist_coeff)
        x = int(res[0][0][0][0])
        y = int(res[0][0][0][1])
        # 在标记中心上画一个红色的圆
        cv2.circle(img, (x, y), 50, (0, 0, 255), -1)
        print(marker_pos[0:3])
    img_resize = imutils.resize(img, width=1024)
    cv2.imshow('', img_resize)
    cv2.waitKey()

运行结果

• 计算两点的距离

这里依然以这张图像作为原始图像，以第二排第一个标志作为我们的基准，所以我们需要知道该基准真实的边长，经过测量，我们知道该基准边长为50mm，先进行透视变换后，再对其他物体的像素长度来计算其真实距离。

import cv2
import cv2.aruco as aruco
import numpy as np

if __name__ == '__main__':

    # 获取需要校正的图像
    img = cv2.imread('IMG_20230222_212055.jpg')
    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)
    # 标记检测器参数
    parameters = aruco.DetectorParameters()
    # 相机内参矩阵
    camera_matrix = np.array([[3.3811614662976608e+03, 0., 2.0318859679474272e+03],
                              [0., 3.3833367964215813e+03, 1.5347124600270040e+03], [0., 0., 1.]])
    # 相机畸变系数
    dist_coeff = np.array([1.4295776727464021e-01, -5.4377825900331767e-01,
                           8.0729484571478852e-04, -2.3554260458140577e-03,
                           6.3943687964321116e-01])

    h, w = img.shape[:2]
    # 畸变去除，获取新的相机内参矩阵
    newcameramtx, roi = cv2.getOptimalNewCameraMatrix(camera_matrix, dist_coeff, (w, h), 1, (w, h))
    # 对图像进行畸变校正，得到校正后的新图像，该图像完全符合小孔成像的比例关系
    img_undist = cv2.undistort(img, camera_matrix, dist_coeff, None, newcameramtx)
    # 检测标记码的标记角坐标(上下左右四个点的坐标)和标记id
    corners, ids, _ = aruco.detectMarkers(img_undist, aruco_dict, parameters=parameters)
    # 获取相机坐标系下的旋转向量和平移向量
    rvecs, tvecs, _ = aruco.estimatePoseSingleMarkers(corners, 0.05, camera_matrix, dist_coeff)
    # 画出标记角
    aruco.drawDetectedMarkers(img_undist, corners, ids)
    # 获取第一个标记的四个点的坐标作为透视变换的源区域
    corner0 = corners[0][0]
    # 透视变换的目标区域
    points2 = np.float32([[500, 1000], [500, 500], [1000, 500], [1000, 1000]])
    # 无畸变的畸变系数
    dist_coeff_undist = np.float32([0, 0, 0, 0, 0])
    # 将3D标记中心点坐标转换成2D图像坐标
    # 对于res来说，它是以屏幕中的图像为相机坐标中的点，故它没有厚度
    # 对于res2来说，可能是实物板，它是有厚度的，我们需要减去这个厚度，故为-0.003，单位米
    # res2 = cv2.projectPoints(np.array([[0.0, 0.0, -0.003]]), rvecs[0], tvecs[0], newcameramtx, dist_coeff_undist)
    res = cv2.projectPoints(np.array([[0.0, 0.0, 0.0]]), rvecs[0], tvecs[0], newcameramtx, dist_coeff_undist)

    # x = int(res2[0][0][0][0])
    # y = int(res2[0][0][0][1])
    # cv2.circle(img_undist, (x, y), 5, (0, 0, 255), -1)
    x1 = int(res[0][0][0][0])
    y1 = int(res[0][0][0][1])
    cv2.circle(img_undist, (x1, y1), 5, (0, 255, 0), -1)
    # 获取透视变换转换矩阵
    M = cv2.getPerspectiveTransform(corner0, points2)
    # 进行透视变换，只有进行了透视变换，我们才可以去真正计算两点的距离和区域面积
    Affine_img = cv2.warpPerspective(img_undist, M, (img_undist.shape[1], img_undist.shape[0]))
    cv2.imshow('img', Affine_img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

运行结果

现在我们来计算第二排第二个标志的边长。

import cv2
import cv2.aruco as aruco
import numpy as np

if __name__ == '__main__':

    # 获取需要校正的图像
    img = cv2.imread('IMG_20230222_212055.jpg')
    # 选择aruco标记的类型（大小）
    aruco_dict = aruco.getPredefinedDictionary(aruco.DICT_5X5_50)
    # 标记检测器参数
    parameters = aruco.DetectorParameters()
    # 相机内参矩阵
    camera_matrix = np.array([[3.3811614662976608e+03, 0., 2.0318859679474272e+03],
                              [0., 3.3833367964215813e+03, 1.5347124600270040e+03], [0., 0., 1.]])
    # 相机畸变系数
    dist_coeff = np.array([1.4295776727464021e-01, -5.4377825900331767e-01,
                           8.0729484571478852e-04, -2.3554260458140577e-03,
                           6.3943687964321116e-01])

    h, w = img.shape[:2]
    # 畸变去除，获取新的相机内参矩阵
    newcameramtx, roi = cv2.getOptimalNewCameraMatrix(camera_matrix, dist_coeff, (w, h), 1, (w, h))
    # 对图像进行畸变校正，得到校正后的新图像，该图像完全符合小孔成像的比例关系
    img_undist = cv2.undistort(img, camera_matrix, dist_coeff, None, newcameramtx)
    # 检测标记码的标记角坐标(上下左右四个点的坐标)和标记id
    corners, ids, _ = aruco.detectMarkers(img_undist, aruco_dict, parameters=parameters)
    # 获取相机坐标系下的旋转向量和平移向量
    rvecs, tvecs, _ = aruco.estimatePoseSingleMarkers(corners, 0.05, camera_matrix, dist_coeff)
    # 画出标记角
    aruco.drawDetectedMarkers(img_undist, corners, ids)
    # 获取第一个标记的四个点的坐标作为透视变换的源区域
    corner0 = corners[0][0]
    # 透视变换的目标区域
    points2 = np.float32([[500, 1000], [500, 500], [1000, 500], [1000, 1000]])
    # 无畸变的畸变系数
    dist_coeff_undist = np.float32([0, 0, 0, 0, 0])
    # 将3D标记中心点坐标转换成2D图像坐标
    # 对于res来说，它是以屏幕中的图像为相机坐标中的点，故它没有厚度
    # 对于res2来说，可能是实物板，它是有厚度的，我们需要减去这个厚度，故为-0.003，单位米
    # res2 = cv2.projectPoints(np.array([[0.0, 0.0, -0.003]]), rvecs[0], tvecs[0], newcameramtx, dist_coeff_undist)
    res = cv2.projectPoints(np.array([[0.0, 0.0, 0.0]]), rvecs[0], tvecs[0], newcameramtx, dist_coeff_undist)

    # x = int(res2[0][0][0][0])
    # y = int(res2[0][0][0][1])
    # cv2.circle(img_undist, (x, y), 5, (0, 0, 255), -1)
    x1 = int(res[0][0][0][0])
    y1 = int(res[0][0][0][1])
    cv2.circle(img_undist, (x1, y1), 5, (0, 255, 0), -1)
    # 获取透视变换转换矩阵
    M = cv2.getPerspectiveTransform(corner0, points2)
    # 进行透视变换，只有进行了透视变换，我们才可以去真正计算两点的距离和区域面积
    Affine_img = cv2.warpPerspective(img_undist, M, (img_undist.shape[1], img_undist.shape[0]))

    corners, ids, _ = aruco.detectMarkers(Affine_img, aruco_dict, parameters=parameters)
    # 获取第二个标志的一条边长，单位像素
    print(abs(corners[1][0][0][0] - corners[1][0][2][0]))

运行结果

504.0

C++版本的OpenCV

下载完OpenCV并编译完成之后(安装方式可以参考模型部署篇 之编译 tensorrtx/yolov5)，在CMakeLists.txt中添加

add_executable(untitled13 opencv.cpp)
find_package(OpenCV REQUIRED)
include_directories(${OpenCV_INCLUDE_DIRS})target_link_libraries(untitled13 ${OpenCV_LIBS})

Makefile方式编译，内容如下

EXE=op

#SUBDIR=src object
PKGS=opencv4
CFLAGS= `pkg-config --cflags $(PKGS)`
LIBS= `pkg-config --libs $(PKGS)`

CXX_OBJECTS := $(patsubst %.cpp,%.o,$(shell find . -name "*.cpp"))
DEP_FILES  =$(patsubst  %.o,  %.d, $(CXX_OBJECTS))

$(EXE): $(CXX_OBJECTS)
	$(CXX)  $(CXX_OBJECTS) -o $(EXE) $(LIBS)
	
%.o: %.cpp
	$(CXX) -c -o $@ $(CFLAGS) $<

clean: 
	rm  -rf  $(CXX_OBJECTS)  $(DEP_FILES)  $(EXE)

test:
	echo $(CXX_OBJECTS)

读取图片

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    Mat img = imread("/Users/admin/Documents/1234.jpeg");
    imshow("img",img);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

打开一个视频文件

#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main(int argc, char* argv[]) {
    namedWindow("video", WINDOW_NORMAL);
    VideoCapture capture("/Users/admin/Documents/Video_20230318115712706_20230318_125022_0.mp4");
    while (capture.isOpened()) {
        Mat frame;
        capture >> frame;
        imshow("video", frame);
        int key = waitKey(10);
        if (key == 113) {
            break;
        }
    }
    capture.release();
    destroyAllWindows();
    return 0;
}

摄像头视频采集

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("video",WINDOW_NORMAL);
    resizeWindow("video",640,640);
    VideoCapture capture(0);
    while (capture.isOpened()) {
        Mat frame;
        capture >> frame;
        imshow("video",frame);
        int key = waitKey(30);
//        按下q键关闭窗口,q的acsii值为113
        if (key == 113) {
            break;
        }
    }
    capture.release();
    destroyAllWindows();
    return 0;
}

运行结果

获取图像像素指针

矩阵的掩膜操作来重新计算每个像素的像素值，掩膜(mask)也被称为kernel。

通过掩膜操作实现图像对比度提高

这里I(i,j)代表图像中的某一个像素。

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/1234.jpeg");
    int cols = (img.cols - 1) * img.channels();
    int offset = img.channels();
    int rows = img.rows;
    Mat dst = Mat::zeros(img.size(),img.type());
    for (int row = 1; row < (rows - 1); ++row) {
//        获取当前行，上一行，下一行
        const uchar* current = img.ptr<uchar>(row);
        const uchar* pre = img.ptr<uchar>(row - 1);
        const uchar* next = img.ptr<uchar>(row + 1);
        uchar* output = dst.ptr<uchar>(row);
        for (int col = offset; col < cols; ++col) {
            output[col] = 5 * current[col] - (current[col - offset] + current[col + offset]
                    + pre[col] + next[col]);
        }
    }
    imshow("img",dst);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

其实这些功能都是我们自己去实现的，但是opencv提供了一个算子，可以直接完成这一步。

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/1234.jpeg");
    Mat dst = Mat::zeros(img.size(),img.type());
    Mat kernel = (Mat_<char>(3,3) << 0,-1,0,-1,5,-1,0,-1,0);
    filter2D(img,dst,img.depth(),kernel);
    imshow("img",dst);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

这里我们将Python版本的filter2D功能做相同的实现

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/111.jpeg");
    Mat dst = Mat::zeros(img.size(),img.type());
    Mat kernel = Mat::ones(5,5,CV_32F) / 25;
    filter2D(img,dst,-1,kernel);
    imshow("img",dst);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

图像色彩空间转换

转换成灰度图像

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/111.jpeg");
    Mat gray;
    cvtColor(img,gray,COLOR_BGR2GRAY);
    imshow("img",gray);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

转化为hsv图像

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/111.jpeg");
    Mat hsv;
    cvtColor(img,hsv,COLOR_BGR2HSV);
    imshow("img",hsv);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

创建一张纯蓝色图像

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat dst = Mat(1080,1920,CV_8UC3);
    dst = Scalar(255,0,0);
    imshow("img",dst);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

深拷贝、修改图像区域颜色

#include <opencv2/opencv.hpp>
//std指标准库
using namespace std;
using namespace cv;

int main(int argc,char** argv) {
    namedWindow("img",WINDOW_NORMAL);
    Mat img = imread("/Users/admin/Documents/111.jpeg");
    Mat dst = img.clone();
    for (int i = 10; i <100; ++i) {
        for (int j = 10; j < 100; ++j) {
            dst.at<Vec3b>(i,j)[0] = 0;
            dst.at<Vec3b>(i,j)[1] = 0;
            dst.at<Vec3b>(i,j)[2] = 255;
        }
    }
    imshow("img",img);
    imshow("dst",dst);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

运行结果

判断一个点是否在一个多边形内

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>

using namespace cv;
using namespace std;

/** @function main */
int main( int argc, char** argv )
{
  /// 创建一个全黑的图像
  const int r = 100;
  Mat src = Mat::zeros( Size( 4*r, 4*r ), CV_8UC1 );

  /// 创建点序列以生成轮廓
  vector<Point2f> vert(6);

  vert[0] = Point( 1.5*r, 1.34*r );
  vert[1] = Point( 1*r, 2*r );
  vert[2] = Point( 1.5*r, 2.866*r );
  vert[3] = Point( 2.5*r, 2.866*r );
  vert[4] = Point( 3*r, 2*r );
  vert[5] = Point( 2.5*r, 1.34*r );

  /// 连接所有的点序列，生成一个六边形
  for( int j = 0; j < 6; j++ )
     { line( src, vert[j],  vert[(j+1)%6], Scalar( 255 ), 3, 8 ); }

  vector<vector<Point> > contours; vector<Vec4i> hierarchy;
  /// 拷贝原图
  Mat src_copy = src.clone();
  /// 获取轮廓contours
  findContours( src_copy, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE);

  /// Calculate the distances to the contour
  Mat raw_dist( src.size(), CV_32FC1 );
  ///计算原图中每一个像素点到轮廓(多边形)的距离，并记录在raw_dist中
  ///如果pointPolygonTest的第三个参数为true，返回距离值，如果为正表示在多边形内部，如果为负表示在多边形外部，0表示在多边形上
  ///如果为false，返回-1，0，1。-1表示在多边形外，0表示在多边形轮廓上，1表示在多边形内
  for( int j = 0; j < src.rows; j++ )
     { for( int i = 0; i < src.cols; i++ )
          { raw_dist.at<float>(j,i) = pointPolygonTest( contours[0], Point2f(i,j), true ); }
     }

  double minVal; double maxVal;
  minMaxLoc( raw_dist, &minVal, &maxVal, 0, 0, Mat() );
  minVal = abs(minVal); maxVal = abs(maxVal);

  /// 以图形方式描述距离
  Mat drawing = Mat::zeros( src.size(), CV_8UC3 );
  /// 以距离为依据，采用不同的像素值来填充描述图像,分为轮廓内的，轮廓外的，轮廓上的
  for( int j = 0; j < src.rows; j++ )
     { for( int i = 0; i < src.cols; i++ )
          {
            if( raw_dist.at<float>(j,i) < 0 )
              { drawing.at<Vec3b>(j,i)[0] = 255 - (int) abs(raw_dist.at<float>(j,i))*255/minVal; }
            else if( raw_dist.at<float>(j,i) > 0 )
              { drawing.at<Vec3b>(j,i)[2] = 255 - (int) raw_dist.at<float>(j,i)*255/maxVal; }
            else
              { drawing.at<Vec3b>(j,i)[0] = 255; drawing.at<Vec3b>(j,i)[1] = 255; drawing.at<Vec3b>(j,i)[2] = 255; }
          }
     }

  /// Create Window and show your results
  char* source_window = "Source";
  namedWindow( source_window, CV_MINOR_VERSION);
  imshow( source_window, src );
  namedWindow( "Distance", CV_MINOR_VERSION);
  imshow( "Distance", drawing );

  waitKey(0);
  return(0);
}

运行结果

ROI区域获取、二值化、腐蚀、膨胀、轮廓查找、轮廓外接矩形、周长、面积、屏幕写字

#include <opencv2/opencv.hpp>
#include <iostream>
#include <sstream>

using namespace std;
using namespace cv;

int main(int argc, char* argv[]) {
    namedWindow("video", WINDOW_NORMAL);
    VideoCapture capture("/Users/admin/Downloads/1234.mpeg");
    RotatedRect rect;
    Point2f P[4];
    int i = 0;
    double max_area;
    int show_area;
    while (capture.isOpened()) {
        Mat frame;
        capture >> frame;
        //获取激光区域
        Mat roi = frame(Rect(76, 320, 1108, 250));
        Mat gray, thresh, er1, er2, dil, element1, element2, element3;
        element1 = getStructuringElement(MORPH_RECT, Size(9, 9));
        element2 = getStructuringElement(MORPH_RECT, Size(5, 5));
        element3 = getStructuringElement(MORPH_RECT, Size(3, 3));
        //将该区域转化为灰度模式
        cvtColor(roi, gray, COLOR_BGR2GRAY);
        //二值化
        threshold(gray, thresh, 20, 255, THRESH_BINARY);
        //第一次腐蚀
        erode(thresh, er1, element1);
        //对部分区域遮盖
        Mat mask = Mat::zeros(Size(30, 230), CV_8UC1);
        Mat roi_mask = er1(Rect(0, 20, 30, 230));
        mask.copyTo(roi_mask);
        //第二次腐蚀
        erode(er1, er2, element2);
        //膨胀
        dilate(er2, dil, element3);
        //第一次轮廓查找
        vector<vector<Point>> contours;
        vector<Vec4i> hierarchy;
        findContours(dil, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE);
        Scalar color(255, 255, 255);
        //遍历所有的轮廓
        for (int i = 0; i < contours.size(); ++i) {
            double len = arcLength(contours[i], true);
            //如果该轮廓周长足够大
            if (len > 2500) {
                //获取该轮廓的最小外接矩形
                rect = minAreaRect(contours[i]);
                rect.points(P);
                //画出外接矩形
//                for (int j = 0; j <= 3; j++) {
//                    line(dil, P[j], P[(j + 1) % 4], color, 2);
//                }
                line(dil, P[1], P[2], color, 4);
                circle(dil, P[2], 25, color, -1);
            } else {
//                for (int j = 0; j <= 3; j++) {
//                    line(dil, P[j], P[(j + 1) % 4], color, 2);
//                }
                line(dil, P[1], P[2], color, 4);
                circle(dil, P[2], 25, color, -1);
            }
        }
        //第二次轮廓查找
        findContours(dil, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE);
        //画出所有轮廓
        drawContours(roi, contours, -1, Scalar(0,0,255), 2, 8, hierarchy);
        //获取所有轮廓的最大面积
        double area = 0;
        for (int i = 0; i < contours.size(); ++i) {
            double temp_area = contourArea(contours[i]);
            if (temp_area > area) {
                area = temp_area;
            }
        }
        //设置基准最大面积
        if (i == 0) {
            max_area = area;
            i++;
        }
        stringstream ss;
        if (area * 100 / max_area > 100) {
            putText(roi, "0%", Point(10, 230), FONT_HERSHEY_SIMPLEX, 1, Scalar(0,255,255), 2);
            cout << "0%" << endl;
        } else if (area * 100 / max_area > 40) {
            show_area = int((1.0 - area / max_area) * 100);
            ss << show_area << "%";
            putText(roi, ss.str(), Point(10, 230), FONT_HERSHEY_SIMPLEX, 1, Scalar(0,255,255), 2);
            cout << show_area << "%" << endl;
        } else {
            ss << show_area << "%";
            putText(roi, ss.str(), Point(10, 230), FONT_HERSHEY_SIMPLEX, 1, Scalar(0,255,255), 2);
            cout << show_area << "%" << endl;
        }
        imshow("video", roi);
//        cout << frame.size() << endl;
//        if (i == 0) {
//            imwrite("/Users/admin/Downloads/1234.jpg", er);
//        }
//        i++;
        int key = waitKey(1);
        if (key == 113) {
            break;
        }
    }
    capture.release();
    destroyAllWindows();
    return 0;
}

区域消除与区域连接(部分代码)

由于有数学计算，需要增加头文件

#include <cmath>
#include <opencv2/core.hpp>
#include <opencv2/videoio.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/imgproc/types_c.h>

主要部分代码如下

    Mat thresh, mor, dil, element, element1;
    element = getStructuringElement(MORPH_RECT, Size(5, 5));
    element1 = getStructuringElement(MORPH_RECT, Size(5, 5));
	//二值化
    threshold(srcImage, thresh, 5, 255, THRESH_BINARY);
	//开运算消除杂点
    morphologyEx(thresh, mor, MORPH_OPEN, element);
	//膨胀
    dilate(mor, dil, element1);
    vector<vector<Point>> contours;
	//需要消除掉的区域
    vector<vector<Point>> blacks;
	//保留下来的区域
    vector<vector<Point>> boundingPoints;
    vector<Vec4i> hierarchy;
	//第一次查找区域
    findContours(dil, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE);
    Scalar color(255, 255, 255);
    for (int i = 0; i < contours.size(); ++i) {
         //如果该区域的长度小于200则抹除
         double len = arcLength(contours[i], true);
         if (len < 200) {
             blacks.push_back(contours[i]);
         }
    }
	//在第一次查找中需要的消除区域进行黑色填充
    fillPoly(dil, blacks, Scalar(0, 0 ,0));
    //两点连接的最小距离
    float minD;
    //两个区域点的索引
    int mink, minl；
    //距离阈值
    float T = 600;
	//第二次区域查找
    findContours(dil, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE);
	//遍历所有的区域
    for (int i = 0; i < contours.size(); ++i) {
		//获取该区域所有的点
		vector<Point> p1 = contours[i];
		//遍历除了本区域的其他区域
		for (int j = i + 1; j < contours.size(); ++j) {
			//获取其他区域的所有的点
			vector<Point> p2 = contours[j];
			for (int k = 0; k < p1.size(); k++) {
                for (int l = 0; l < p2.size(); l++) {
                    //计算两个区域所有点的距离，找出两个区域中最短的两点的距离
                    float d = sqrt(pow(float(p1[k].x - p2[l].x), 2) + pow(float(p1[k].y - p2[l].y), 2));
                    if (k == 0 && l == 0) {
                        minD = d;
                        mink = k;
                        minl = l;
                    }
                    if (d < minD) {
                        minD = d;
                        mink = k;
                        minl = l;
                    }
                }
            }
            //如果该最短距离小于阈值则进行连线
            if (minD < T) {
                line(dil, p1[mink], p2[minl], color, 8, LINE_AA, 0);
            }
		}	
    }

判断视频运动还是静止

#include <opencv2/core.hpp>
#include <opencv2/videoio.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/imgproc/types_c.h>
#include <iostream>
#include <vector>

using namespace cv;
using namespace std;

/*
 * 获取两幅图像差值的均值
 */
int bMean(Mat src)
{
    int counter = 0;
    int mean = 0;
    Mat_<uchar>::iterator it = src.begin<uchar>();
    Mat_<uchar>::iterator itend = src.end<uchar>();
    for (; it != itend; ++it) {
        if ((*it) > 0) { 
			mean += *it;
            counter++;
		}
    }
    return mean / counter;
}

int main() {
    VideoCapture capture("rtsp://address");
    Mat gry, prev_img;
    bool first_frame = true;
    time_t start = time(NULL);
    while(capture.isOpened()) {
		Mat img;
		time_t end = time(NULL);
		if (end - start > 2) {
			start = time(NULL);
			capture >> img;
			if (!img.empty()) {
				cvtColor(img, gry, CV_BGR2GRAY);
			} else {
				capture.open("rtsp://address");
				continue;
			}
			if (first_frame) {
				prev_img = gry.clone();
				first_frame = false;
				continue;
			}
			int col = gry.cols;
			int row = gry.rows;
			Mat diff = Mat(row, col, CV_8UC1);
			//比较两幅图像，并将差值放入diff
			absdiff(prev_img, gry, diff);
			int change = bMean(diff);
			cout << change << endl;
			if (change > 1) {
			cout << "视频运动中" << endl;
			} else {
			cout << "视频静止或者无煤运行中" << endl;
			}
			prev_img = gry.clone();
		}
    }
    return 0;
}

Mat矩阵

• Mat对象中的属性

1. data：uchar类型的指针，指向Mat数据矩阵的首地址
2. depth：度量像素每一个通道的精度，它本身与图像的通道数无关！depth数值越大，精度越高。eum{CV_8U=0，CV_8S=1，CV_16S=3，CV_32S=4，CV_32F=5，CV_64F=6}
3. elemSize：表示矩阵中每一个元素的数据大小（结合通道数，单位：字节）。如Mat的数据类型是CV_8UC1，那么elemSize==1；CV_8UC3或CV_8SC3，那么elemSize==3
4. elemSize1：表示矩阵中每一个元素单通道下数据大小（单位：字节），所以有：elemSize1==elemSize/channels。
5. step：矩阵的有效宽度。是累计了一行中所有元素、所有通道、所有通道的elemSize1之后的值。step==cols*elemSize1*channels。

消除反光

上图是整体的流程，对于输入的图像，我们会先通过一个卷积核来分别获取明亮通道的图像。

• 获取明暗通道先验

先验估计暗通道

\(I^{dark}(x)=\min_{c∈[r,g,b]}(\min_{y∈Ω(x)}(I^c(y)))\)

这里\(I^{dark}(x)\)表示低光图像，\(I^c(y)\)表示原始图像的单通道，Ω(x)表示局部区域。

先验估计明亮通道

\(I^{bright}(x)=\max_{c∈[r,g,b]}(\max_{y∈Ω(x)}(I^c(y)))\)

Python代码

import cv2
import numpy as np

def get_illumination_channel(I, w):
    M, N, _ = I.shape
    # 通道填充
    padded = np.pad(I, ((int(w/2), int(w/2)), (int(w/2), int(w/2)), (0, 0)), 'edge')
    darkch = np.zeros((M, N))
    brightch = np.zeros((M, N))
    for i, j in np.ndindex(darkch.shape):
        darkch[i, j] = np.min(padded[i:i + w, j:j + w, :]) # dark channel
        brightch[i, j] = np.max(padded[i:i + w, j:j + w, :]) # bright channel

    return darkch, brightch

# 读取图像
image = cv2.imread('/Users/admin/Documents/WechatIMG1287.jpg')
d, b = get_illumination_channel(image, 5)
b = b.astype(np.uint8)
d = d.astype(np.uint8)

# 显示结果
cv2.imshow('Result', b)
cv2.imshow('Result', d)
cv2.waitKey(0)
cv2.destroyAllWindows()

输入图像为

暗通道图像为

亮通道图像为

C++代码

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

pair<Mat, Mat> get_illumination_channel(Mat I, float w) {
    int N = I.size[0];
    int M = I.size[1];
    Mat darkch = Mat::zeros(Size(M, N), CV_32FC1);
    Mat brightch = Mat::zeros(Size(M, N), CV_32FC1);

    int padding = int(w/2);
    // 通道填充
    Mat padded = Mat::zeros(Size(M + 2*padding, N + 2*padding), CV_32FC3);

    for (int i = padding; i < padding + M; i++) {
        for (int j = padding; j < padding + N; j++) {
            padded.at<Vec3f>(j, i).val[0] = (float)I.at<Vec3b>(j-padding, i-padding).val[0] / 255;
            padded.at<Vec3f>(j, i).val[1] = (float)I.at<Vec3b>(j-padding, i-padding).val[1] / 255;
            padded.at<Vec3f>(j, i).val[2] = (float)I.at<Vec3b>(j-padding, i-padding).val[2] / 255;
        }
    }

    for (int i = 0; i < darkch.size[1]; i++) {
        int col_up, row_up;
        col_up = int(i+w);
        for (int j = 0; j < darkch.size[0]; j++) {
            double minVal, maxVal;
            row_up = int(j + w);
            minMaxLoc(padded.colRange(i, col_up).rowRange(j, row_up), &minVal, &maxVal);
            darkch.at<float>(j,i) = minVal; //dark channel
            brightch.at<float>(j,i) = maxVal; //bright channel
        }
    }

    return make_pair(darkch, brightch);
}

int main(int argc, char* argv[]) {
    Mat src = imread("/Users/admin/Documents/WechatIMG1287.jpg");
    pair<Mat, Mat> result = get_illumination_channel(src, 5.0);
    imshow("result", result.second);
    imshow("result", result.first);
    waitKey(0);
    destroyAllWindows();
    return 0;
}

• 计算全局大气光照

计算全局大气光照，是通过上面获得的明亮通道取前10%强度的平均值来计算的。取10%的值是为了确保一个小的异常不会对其产生很大的影响。

Python代码

def get_atmosphere(I, brightch, p=0.1):
    M, N = brightch.shape
    flatI = I.reshape(M * N, 3)  # 扁平化原始图像
    flatbright = brightch.ravel()  # 扁平化明亮通道图像
    searchidx = (-flatbright).argsort()[:int(M * N * p)]  # 获取前10%的亮度值(从大到小)
    A = np.mean(flatI.take(searchidx, axis=0), dtype=np.float64, axis=0)  # 对原始图像的三通道取平均值
    return A

A = get_atmosphere(image, b)
print(A)

运行结果

[240.43043996 231.86428468 216.26792137]

C++代码

Mat get_atmosphere(Mat I, Mat brightch, float p=0.1) {
    int N = brightch.size[0];
    int M = brightch.size[1];

    // 扁平化原始图像
    Mat flatI(Size(1, N * M), CV_8UC3);
    // 扁平化明亮通道图像(包含数值和索引)
    vector<pair<float, int>> flatBright;

    for (int i=0; i < M; i++) {
        for (int j=0; j < N; j++) {
            int index = i*N + j;
            flatI.at<Vec3b>(index, 0).val[0] = I.at<Vec3b>(j, i).val[0];
            flatI.at<Vec3b>(index, 0).val[1] = I.at<Vec3b>(j, i).val[1];
            flatI.at<Vec3b>(index, 0).val[2] = I.at<Vec3b>(j, i).val[2];

            flatBright.push_back(make_pair(-brightch.at<float>(j, i), index));
        }
    }


    // 对扁平化后的明亮通道图像进行排序
    sort(flatBright.begin(), flatBright.end());

    Mat A = Mat::zeros(Size(1, 3), CV_32FC1);

    for (int k=0; k < int(M*N*p); k++) {
//        获取索引值
        int sindex = flatBright[k].second;
//        获取原始图像相应索引上的图像值
        A.at<float>(0, 0) = A.at<float>(0, 0) + (float)flatI.at<cv::Vec3b>(sindex, 0).val[0];
        A.at<float>(1, 0) = A.at<float>(1, 0) + (float)flatI.at<cv::Vec3b>(sindex, 0).val[1];
        A.at<float>(2, 0) = A.at<float>(2, 0) + (float)flatI.at<cv::Vec3b>(sindex, 0).val[2];
    }
//    取10%
    A = A / int(M * N * p);

    return A;
}

Mat A = get_atmosphere(src, result.second);
cout << A << endl;

• 查找初始透射图

投射图描述了未被散射并到达相机的光部分，从以下等式从明亮通道先验估计得到

\(t^{bright}={I^{bright}(x)-A^c\over 255-A^c}\)

\(A^c\)是大气光局部区域的最大值

Python代码

def get_initial_transmission(A, brightch):
    Ac = np.max(A)
    init_t = (brightch - Ac) / (1. - Ac)  # 初始透射图
    return (init_t - np.min(init_t)) / (np.max(init_t) - np.min(init_t))  # 归一化初始透射图

t = get_initial_transmission(A, b)
# 显示结果
cv2.imshow('Result', t)
cv2.waitKey(0)
cv2.destroyAllWindows()

运行结果

C++代码

Mat get_initial_transmission(Mat A, Mat brightch) {
    double A_n, A_x, minVal, maxVal;
//    查找A中的最大值A_x，最小值无视
    minMaxLoc(A, &A_n, &A_x);
    Mat init_t(brightch.size(), CV_32FC1);
    init_t = brightch.clone();
    // 初始透射图
    init_t = (init_t - A_x) / (1.0 - A_x);
    minMaxLoc(init_t, &minVal, &maxVal);
    // 归一化初始透射图
    init_t = (init_t - minVal) / (maxVal - minVal);

    return init_t;
}

Mat t = get_initial_transmission(A, result.second);
imshow("result", t);
waitKey(0);
destroyAllWindows();

• 使用暗通道校正透射图

通过暗通道先验计算投射图，并计算先验之间的差异。此计算是为了纠正从明亮通道先验获得的潜在错误透射估计。

Python代码

def get_corrected_transmission(I, A, darkch, brightch, init_t, alpha, omega, w):
    im = np.empty(I.shape, I.dtype)
    for ind in range(0, 3):
        im[:, :, ind] = I[:, :, ind] / A[ind]  # 将像素值除以大气光
    dark_c, _ = get_illumination_channel(im, w)  # 暗通道透射图
    dark_t = 1 - omega * dark_c  # 修正暗通道透射图
    corrected_t = init_t  # 用初始透射图初始化校正透射图
    diffch = brightch - darkch  # 投射图之间的差

    for i in range(diffch.shape[0]):
        for j in range(diffch.shape[1]):
            if(diffch[i, j] < alpha):
                corrected_t[i, j] = dark_t[i, j] * init_t[i, j]

    return np.abs(corrected_t)

ct = get_corrected_transmission(image, A, d, b, t, 0.4, 0.75, 5)
# 显示结果
cv2.imshow('Result', ct)
cv2.waitKey(0)
cv2.destroyAllWindows()

运行结果

C++代码

Mat get_corrected_transmission(Mat I, Mat A, Mat darkch, Mat brightch, Mat init_t, float alpha, float omega, int w) {
    Mat im3(I.size(), CV_32FC3);
    //将像素值除以大气光
    for (int i=0; i < I.size[1]; i++) {
        for (int j=0; j < I.size[0]; j++) {
            im3.at<Vec3f>(j, i).val[0] = (float)I.at<Vec3b>(j, i).val[0] / A.at<float>(0, 0);
            im3.at<Vec3f>(j, i).val[1] = (float)I.at<Vec3b>(j, i).val[1] / A.at<float>(1, 0);
            im3.at<Vec3f>(j, i).val[2] = (float)I.at<Vec3b>(j, i).val[2] / A.at<float>(2, 0);
        }
    }

    Mat dark_c, dark_t, diffch;

    pair<Mat, Mat> illuminate_channels = get_illumination_channel(im3, w);
    // 暗通道投射图
    dark_c = illuminate_channels.first;
    // 修正暗通道透射图
    dark_t = 1 - omega * dark_c;
    Mat corrected_t = init_t;
    diffch = brightch - darkch; //投射图之间的差

    for (int i=0; i < diffch.size[1]; i++) {
        for (int j=0; j < diffch.size[0]; j++) {
            if (diffch.at<float>(j, i) < alpha) {
                // 用初始透射图初始化校正透射图
                corrected_t.at<float>(j, i) = abs(dark_t.at<float>(j, i)*init_t.at<float>(j, i));
            }
        }
    }

    return corrected_t;
}

Mat ct = get_corrected_transmission(src, A, result.first, result.second, t, 0.4, 0.75, 5);
imshow("result", ct);
waitKey(0);
destroyAllWindows();

• 计算结果图像

Python代码

def get_final_image(I, A, refined_t, tmin):
    refined_t_broadcasted = np.broadcast_to(refined_t[:, :, None], (refined_t.shape[0], refined_t.shape[1], 3))  # 将2D精制图的通道复制到3个通道
    J = (I - A) / (np.where(refined_t_broadcasted < tmin, tmin, refined_t_broadcasted)) + A  # 得到最终结果

    return (J - np.min(J)) / (np.max(J) - np.min(J))  # 归一化图像

f = get_final_image(image, A, ct, 0.6)
# 显示结果
cv2.imshow('Result', f)
cv2.waitKey(0)
cv2.destroyAllWindows()

运行结果

C++代码

Mat get_final_image(Mat I, Mat A, Mat refined_t, float tmin) {
    Mat J(I.size(), CV_32FC3);
    for (int i = 0; i < refined_t.size[1]; i++) {
        for (int j = 0; j < refined_t.size[0]; j++) {
            float temp = refined_t.at<float>(j, i);
            if (temp < tmin) {
                temp = tmin;
            }
            // finding result
            J.at<Vec3f>(j, i).val[0] = (I.at<Vec3b>(j, i).val[0] - A.at<float>(0,0)) / temp + A.at<float>(0,0);
            J.at<Vec3f>(j, i).val[1] = (I.at<Vec3b>(j, i).val[1] - A.at<float>(1,0)) / temp + A.at<float>(1,0);
            J.at<Vec3f>(j, i).val[2] = (I.at<Vec3b>(j, i).val[2] - A.at<float>(2,0)) / temp + A.at<float>(2,0);
        }
    }
    double minVal, maxVal;
    minMaxLoc(J, &minVal, &maxVal);

    // 归一化
    for (int i = 0; i < J.size[1]; i++) {
        for (int j = 0; j < J.size[0]; j++) {
            J.at<Vec3f>(j, i).val[0] = (J.at<Vec3f>(j, i).val[0] - minVal) / (maxVal - minVal);
            J.at<Vec3f>(j, i).val[1] = (J.at<Vec3f>(j, i).val[1] - minVal) / (maxVal - minVal);
            J.at<Vec3f>(j, i).val[2] = (J.at<Vec3f>(j, i).val[2] - minVal) / (maxVal - minVal);
        }
    }

    return J;
}

Mat f = get_final_image(src, A, ct, 0.6);
imshow("result", f);
waitKey(0);
destroyAllWindows();