OpenCV图像处理——积分图像计算（C++/Python）

概述

积分图像是一种高效的图像处理技术，最初由Crow在1984年提出，目的是为了提高多尺度透视投影的渲染速度。它通过构建一个积分图，使得图像中任意矩形区域的像素和能够在常数时间内快速计算出来，极大地减少了在图像模糊、边缘提取、对象检测等操作中的计算量，提高了计算速度。

积分图像的构建基于一个简单但有效的概念：图像中的每个点都存储了从图像左上角到该点的区域的像素值之和。这意味着，一旦积分图像构建完成，就可以通过简单的查表和有限次运算来快速获取任何区域的像素和，从而加速了多种图像处理算法的执行。

在实际应用中，积分图像被广泛用于多种图像处理算法中。例如，在Viola-Jones的对象检测框架中，积分图像被用来加速Haar特征的计算。此外，SURF特征提取算法也利用积分图像来快速计算图像特征，从而保持了尺度不变性和旋转不变性，同时对光照变化和放射变化具有很强的鲁棒性。

除了在特征提取中的应用，积分图像也用于二值图像分析，如图像的腐蚀和膨胀操作，以及图像相似相关性NCC（归一化互相关）的计算中。这些应用展示了积分图像在图像处理领域的多样性和实用性。

积分图（Integral Image）的定义：取图像左上侧的全部像素计算累加和，并用这个累加和替换图像中的每一个像素，使用这种方式得到的图像称为积分图像。

公式：

积分图又称总和面积表（summed area table，简称SAT），是一个快速且有效的对一个网格的矩形子区域中计算和的数据结构和算法。

积分图可以只遍历一次图像即可有效的计算出来，其通常被用来加速计算过程。一旦积分图计算完毕，对任意矩形区域的和的计算就可以在常数时间内（一次加法，两次减法）完成。如下图中，阴影矩形区域的值：

其中公式如下：
图像积分图建立与查找，在积分图像(Integral Image - ii)上任意位置(x, y)处的ii(x, y)表示该点左上角所有像素之和， 其中(x,y)是图像像素点坐标。

API

integral(
InputArray src, // 输入图像
OutputArray sum, // 和表
OutputArray sqsum, // 平方和表
OutputArray tilted, // 瓦块和表
int sdepth = -1, // 和表数据深度常见CV_32S
int sqdepth = -1 // 平方和表数据深度 常见 CV_32F
)

C++ 实现

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

using namespace std;
using namespace cv;

void blur_demo(Mat &image, Mat &sum);
void edge_demo(Mat &image, Mat &sum);
int getblockSum(Mat &sum, int x1, int y1, int x2, int y2, int i);

/*
 * 图像积分图算法
 */
int main() {
    Mat src = imread("../images/test.png");
    if (src.empty()) {
        cout << "could not load image.." << endl;
    }
    imshow("input", src);

    // 计算积分图
    Mat sum, sqrsum;
    integral(src, sum, sqrsum);

    /*
     * 积分图应用
    */
    int type = 0;
    // 模糊应用
    blur_demo(src, sum);
    // 边缘检测
    edge_demo(src, sum);

    waitKey(0);
    return 0;
}

void blur_demo(Mat &image, Mat &sum) {
    int w = image.cols;
    int h = image.rows;
    Mat result = Mat::zeros(image.size(), image.type());
    int x2 = 0, y2 = 0;
    int x1 = 0, y1 = 0;
    int ksize = 5;
    int radius = ksize / 2;
    int ch = image.channels();
    int cx = 0, cy = 0;
    for (int row = 0; row < h + radius; row++) {
        y2 = (row + 1)>h ? h : (row + 1);
        y1 = (row - ksize) < 0 ? 0 : (row - ksize);
        for (int col = 0; col < w + radius; col++) {
            x2 = (col + 1)>w ? w : (col + 1);
            x1 = (col - ksize) < 0 ? 0 : (col - ksize);
            cx = (col - radius) < 0 ? 0 : col - radius;
            cy = (row - radius) < 0 ? 0 : row - radius;
            int num = (x2 - x1)*(y2 - y1);
            for (int i = 0; i < ch; i++) {
                // 积分图查找和表，计算卷积
                int s = getblockSum(sum, x1, y1, x2, y2, i);
                result.at<Vec3b>(cy, cx)[i] = saturate_cast<uchar>(s / num);
            }
        }
    }
    imshow("blur_demo", result);
}
/**
* 3x3 sobel 垂直边缘检测演示
*/
void edge_demo(Mat &image, Mat &sum) {
    int w = image.cols;
    int h = image.rows;
    Mat result = Mat::zeros(image.size(), CV_32SC3);
    int x2 = 0, y2 = 0;
    int x1 = 0, y1 = 0;
    int ksize = 3; // 算子大小，可以修改，越大边缘效应越明显
    int radius = ksize / 2;
    int ch = image.channels();
    int cx = 0, cy = 0;
    for (int row = 0; row < h + radius; row++) {
        y2 = (row + 1)>h ? h : (row + 1);
        y1 = (row - ksize) < 0 ? 0 : (row - ksize);
        for (int col = 0; col < w + radius; col++) {
            x2 = (col + 1)>w ? w : (col + 1);
            x1 = (col - ksize) < 0 ? 0 : (col - ksize);
            cx = (col - radius) < 0 ? 0 : col - radius;
            cy = (row - radius) < 0 ? 0 : row - radius;
            int num = (x2 - x1)*(y2 - y1);
            for (int i = 0; i < ch; i++) {
                // 积分图查找和表，计算卷积
                int s1 = getblockSum(sum, x1, y1, cx, y2, i);
                int s2 = getblockSum(sum, cx, y1, x2, y2, i);
                result.at<Vec3i>(cy, cx)[i] = saturate_cast<int>(s2 - s1);
            }
        }
    }
    Mat dst, gray;
    convertScaleAbs(result, dst);
    normalize(dst, dst, 0, 255, NORM_MINMAX);
    cvtColor(dst, gray, COLOR_BGR2GRAY);
    imshow("edge_demo", gray);
}
int getblockSum(Mat &sum, int x1, int y1, int x2, int y2, int i) {
    int tl = sum.at<Vec3i>(y1, x1)[i];
    int tr = sum.at<Vec3i>(y2, x1)[i];
    int bl = sum.at<Vec3i>(y1, x2)[i];
    int br = sum.at<Vec3i>(y2, x2)[i];
    int s = (br - bl - tr + tl);
    return s;
}

Python代码实现

import cv2 as cv
import numpy as np


def get_block_sum(ii, x1, y1, x2, y2, index):
    tl = ii[y1, x1][index]
    tr = ii[y2, x1][index]
    bl = ii[y1, x2][index]
    br = ii[y2, x2][index]
    s = (br - bl - tr + tl)
    return s


def blur_demo(image, ii):
    h, w, dims = image.shape
    result = np.zeros(image.shape, image.dtype)
    ksize = 15
    radius = ksize // 2
    for row in range(0, h + radius, 1):
        y2 = h if (row + 1)> h else (row + 1)
        y1 = 0 if (row - ksize) < 0 else (row - ksize)
        for col in range(0, w + radius, 1):
            x2 = w if (col + 1)>w else (col + 1)
            x1 = 0 if (col - ksize) < 0 else (col - ksize)
            cx = 0 if (col - radius) < 0 else (col - radius)
            cy = 0 if (row - radius) < 0 else (row - radius)
            num = (x2 - x1)*(y2 - y1)
            for i in range(0, 3, 1):
                s = get_block_sum(ii, x1, y1, x2, y2, i)
                result[cy, cx][i] = s // num

    cv.imshow("integral fast blur", result)
    cv.imwrite("D:/result.png", result)


src = cv.imread("D:/images/test1.png")
cv.namedWindow("input", cv.WINDOW_AUTOSIZE)
cv.imshow("input", src)
sum_table = cv.integral(src, sdepth=cv.CV_32S)
blur_demo(src, sum_table)

cv.waitKey(0)
cv.destroyAllWindows()