OpenCV-11-背景提取

1 摘要

由于背景分割比较(Background Subtraction,也可以称为背景差分,Background Differencing)简单,并且在大多数场景中相机的位置也是相对固定的,因此在很多应用程序中它都是一个非常重要的图像处理操作,尤其是在视频安全程序中。在学习背景分割之前,首先需要了解背景模型。将学习到的背景模型(Background Model)和当前图像相比较,并去除匹配的背景部分,剩下的部分就被认为是新的前景对象。

当然背景是一个定义不明确的概念,这意味着其具体含义在不同程序中会不同。例如当你观察一个高速公路场景时,尽管车流是不停变换的,它也应当被认为是背景。通常情况下将场景中静止的或者周期移动的物体定义为背景(Background),剩余的部分定义为前景(Foreground)。另一种背景非静止的场景是以室外的树为背景,显然风会使得树来回摆动。静止和移动的背景通常分别在室内和室外的场景会遇到,我们需要找到一些工具能够处理这两种不同场景。

本章先会介绍传统的背景提取模型,并讨论其缺点,同时还会介绍一种快速背景提取方法,大多数情况下它都能够很好的处理室内光线改变不太强烈的静止背景。然后再介绍高级的场景模型,同时介绍“Code Book”方法,尽管他的效率相较于快速背景提取法更低,但是它能够很好的处理室内和室外场景,它还能够处理周期变化的背景(如风中摇摆的树)以及缓慢或者周期变化的光线环境。该方法还能在偶尔有前景物体移动的工程中学习背景模型。

接下来会讨论如何使用前面章节学到的连通分量在提取前景对象时执行一些清理任务。然后会比较快速和“Code Book”背景提取法的差异和适用场景。本章的结尾将会讨论OpenCV中关于两个现代背景提取算法的实现,这些算法适用了本章已经提及到的原理,但是它们进行了扩展和细节的优化,使得它们更适合于实际的应用程序。

2 背景提取的缺点

尽管这里提到的背景模型提取方法在处理简单场景的时候效果很不错,但是它们都基于一个假设的前提条件,即从统计学上看连续的视频帧中图像内不同位置像素的行为都是相互独立的,然而现实情况并非如此。本章即将介绍的背景模型算法都只考虑了单个位置像素的值在时间上的变化,并没有考虑它们和领域像素之间的关系。假设需要考虑相邻像素,我们可以设计一个更复杂的模型,比如考虑领域像素的亮度对某个像素的影响,即区分邻域像素明亮或者暗淡这两种情况。

此时对于每个位置的像素都训练两个模型,一个模型描述当前像素的邻域像素亮度高的情况,另一个模型描述当前像素的邻域像素亮度低的情况(考虑一个台灯背景,随着台灯亮度发生变化,则背景内单个像素强度会受到周围像素强度变化的影响)。但是这样对于每个像素都需要记录两个不同的值,会带来两倍的空间成本以及更多的计算成本。我们甚至可以将特定像素和邻域像素的关系由理想的高低亮度两种情况推广到更多的关系中,并生成更高维度直方图,通过几个步骤来实现更复杂的运算。显然空间和时间维度越复杂的模型需要更多的内存资源,更多的数据样本,以及更多的计算资源。

当计算机需要从数据中学习到一些东西时,通常需要克服的首要困难是获取到足够的数据。越复杂的模型越容易得到直观的结果,其表达能力更强,但是需要更多的训练数据,在本系列后续章节介绍机器学习时将会详细讨论这个议题。

由于复杂的模型会带来额外的成本,因此在实际应用中需要尽量避免使用这些复杂的模型。在遵循假设条件像素之间的行为是独立的前提下,通常会清理假正向(False Positive)像素从而提高算法效率。通常使用如腐蚀(cv::erode())、膨胀(cv::dilate())和浸水填充等图像处理操作来实现清理任务。在本系列的前面章节中介绍在包含噪声的图像中寻找大而紧凑的连通分量(Connected Components)时已经介绍过这些图像处理技术。在本章中我们还会用到连通分量的技巧,但是只考虑满足像素行为互相独立的假设前提条件。

3 场景建模

在正式学习场景建模之前,我们再回顾一下前景和背景的概念。如果观察一个停车场,此时一辆汽车驶入,则汽车就是一个新的前景物体,但是它一直都是前景吗?如果观察一片空地,一个垃圾桶从一个位置移动到了另外一个位置,则其移动前后的位置都属于新的前景对象,同样的它们会一直都是前景对象吗?如果对一个漆黑的房间建模,把等打开后房间会被照亮,是否整个房间都应该属于前景呢?为了回答这些问题,我们需要一个高层次的场景模型,从前景到背景分割为多层,并且需要一个基于时间变化的更新函数将静止的前景对象逐步降级为背景对象。同样检测全局场景是否发生变化,并在变化时训练新的模型。

一般来说场景模型分层为新的前景对象、老的前景对象和背景对象。同时可能需要一些运动检测方法,当物体移动时我们可以同时标记移动物体的正区域(Positive Aspect,即新的位置)和其负区域(Negative Aspect,即原来的位置)。

此时新的前景对象将会被放入前景对象层中,并标记为正或者负对象。如果前景对象在指定的时间内未移动,可以将其降级放入到老的前景对象层,此时算法统计其像素数据,当训练好的模型融入到背景模型中后,对象成为背景。

检测全局环境的改变如场景打开房间内的灯,可以使用全局帧差分技术。例如如果同时很多像素都发生了变化,我们可以认为是全局场景发生了改变,此时需要训练一个全新的背景模型。

3.1 像素切片

在讨论对像素变化建模之前,先了解如何观察图片中的像素在视频帧之间的变化。假设有一个视频拍摄的是随风摇摆的树,在某一时间伸手挡住了部分风景。则下图上半部分是视频文件60帧内抽取的三帧图形,而下半部分的折线图描述了60帧内上部分图片中黑色线条上所有像素的强度变化。我们将整棵树看做是背景,显然需要对这些波动进行建模。但是在这之前我们先额外花一点时间讨论如何对某条线上的像素进行采样,因为在创建特征和调试的时候这是一个非常有用的技巧。

图11-1 视频帧序列沿特定直线像素强度统计图

OpenCV提供了线迭代器(Line Iterator)用于对图像中的任意一条直线上的像素进行采样,其构造函数如下。

// image:需要采样的图片,可以是任意数据类型,可以包含任意数量通道
// pt1:直线的起点
// pt2:直线的终点
// connectivity:直线的连接方式,4邻域或者8邻域
// left_to_right:0表示从pt1采样至pt2,1表示根据x坐标值从左边的点采样至右边点点
cv::LineIterator::LineIterator(const cv::Mat& image,
                               cv::Point pt1, cv::Point pt2,
                               int connectivity = 8, int left_to_right = 0);

通过重载后的运算符++可以访问指定直线上的所有像素,通过重载后的运算符🌟可以获取某个像素的数据,但是其返回值是uchar🌟类型,需要我们自己做类型转换。

示例程序PixelOfARow读取了一个视频文件,并使用该函数对60个视频帧在指定的直线上的像素进行采样,并且得到了上面的统计折线图需要的数据,其核心代码如下。

int main(int argc, const char * argv[]) {
    // 读取视频文件
    cv::VideoCapture cap = cv::VideoCapture(argv[1]);
    // 创建用于播放视频的窗口
    cv::namedWindow(argv[0], cv::WINDOW_AUTOSIZE);
    
    // 准备写入视频帧某条直线所有像素BGR数据的文件
    std::ofstream bfile, gfile, rfile;
    bfile.open("blines.csv");
    gfile.open("glines.csv");
    rfile.open("rlines.csv");
    
    // 定义采样直线的两个端点
    cv::Point ptStart(10, 10), ptEnd(30, 30);
    // 定于接收视频数据的图片
    cv::Mat rawImage;
    cap >> rawImage;
    
    while (rawImage.data) {
        // 沿着指定直线采样
        cv::LineIterator it(rawImage, ptStart, ptEnd, 8);
        // 将数据写入文件
        for (int i = 0; i < it.count; ++i, ++it) {
            bfile << (int)(*it)[0] << ", ";
            gfile << (int)(*it)[1] << ", ";
            rfile << (int)(*it)[2] << ", ";
            // 将红色强度设为最大值
            (*it)[2] = 255;
        }
        
        // 显示图片
        cv::imshow(argv[0], rawImage);
        // 挂起程序等待用户事件,并以25FPS的速度播放视频
        cv::waitKey(1000/25);
        
        // 每帧的数据换行
        bfile << "\n"; gfile << "\n"; rfile << "\n";
        // 读取下一帧
        cap >> rawImage;
    }
    
    // 完成文件写入操作
    bfile << std::endl; gfile << std::endl; rfile << std::endl;
    bfile.close(); gfile.close(); rfile.close();
    return 0;
}

上面的代码使用了指针形式访问每个点的数据,另外一种实用的方式是创建一个合适数据类型的缓存,将所有数据放入缓存后再访问缓存内部的数据,代码如下。

cv::LineIterator it(rawImage, ptStart, ptEnd, 8);
vector<cv::Vec3b> buf(it.count);
for (int i = 0; i < it.count; i++, ++it) {
    buf[i] = &((const cv::Vec3b *)it);
};

接下来我们会对上图中的像素强度波动建模,我们会从简单模型介绍到复杂模型,这些模型的内存消耗和计算成本都会被限制在合理的范围内,从而使得它们可以处理实时视频数据。

3.2 帧差分

最简单的背景提取模型就是计算视频帧的差异,即用某一个视频帧减去时间轴上前面的某个视频帧(不强制要求计算连续视频帧的差异),得到的差值如果足够大,则认为这部分内容为前景对象。这个过程容易捕捉到移动物体的边缘。为了理解简单,考虑三个单通道视频帧frameTime1、frameTime2和frameForeground,可以通过如下代码计算前后两个视频帧之间的绝对差值得到前景图像。

cv::absdiff(frameTime1, frameTime2, frameForeground);

由于视频帧内总会存在噪声,帧间也总会存在波动,因此通常使用阈值函数处理上述代码得到的结果。如下面的代码将阈值设置为15,并且将高于该阈值的像素值设置为255。

cv::threshold(frameForeground, frameForeground, 15, 255, cv::THRESH_BINARY);

通过上述代码处理后得到的前景图像通常仍包含一些噪声,接下来可以使用前面章节学到的图像处理操作,如腐蚀操作cv::erode()或者是连通分量分析来清楚这些噪声。在处理彩色图像时,可以使用上述的操作分别处理三个颜色通道,最后再使用cv::max()合并三个通道结果得到最终的前景图像。

但是对于大多数应用程序而言,这种仅仅标记移动物体的模型过于简单。在更高效的背景提取模型中,(由于背景像素在不同帧之间存在波动,为了得到更准确对模型)我们需要分析场景中的像素值的平均值和平均偏差。帧差分的示例程序会在下文介绍完一些基本的背景知识后给出。

4 平均法

平均法学习所有像素的平均值和标准差(或者使用更快的平均偏差),并将其作为背景模型。在图11-1中我们统计了视频帧间每个像素的强度变化,而现在我们需要每个像素在视频帧间分布的平均值和平均偏差(这里统计的是两帧图像之间的差值,而不是和平均值的差值),统计结果如下图。在示例程序使用的视频文件中,某一时刻手从摄像头面前划过,其亮度比背景中的天空和树更低,在下图使用虚线表示。

图11-2 视频帧间像素强度平均值和平均偏差统计图

在上图所表示的背景提取模型中,对于需要比较的视频帧,位于阈值即粗黑线Mean±2✖️Avg_Diff内部的就是背景图像,外部的就是前景图像。尽管上图使用的阈值系数是2,但是实际应用时该系数还会调整。

4.1 平均偏差

示例程序AverageMethod使用平均偏差训练了背景提取模型,并标识出前景像素。程序首先定义了如下全局变量保存图像像素的统计数据。

// Global storage
// Float, 3-channel images needed by learning model.
cv::Mat previousFrame;
cv::Mat totalFrame, averageFrame, totalDifferenceFrame, averageDifferenceFrame;
cv::Mat upperLimitFrame, lowerLimitFrame;
cv::Mat forgroundFrame;
cv::Mat tmpFrame, tmpFrame2;

// Float, 1-channel images needed by learning model.
// Present each channel respectively.
std::vector<cv::Mat> channelFrames(3);
std::vector<cv::Mat> channelLowFrames(3);
std::vector<cv::Mat> channelHighFrames(3);

// Thresholds
// scaling the thesholds in backgroundDiff()
float highThreshFactor = 20.0;
float lowThreshFactor = 28.0;

// Counts number of images learned for averaging later
float learnedImageCount = 0;

接下来定义一个函数为这些需要使用到的矩阵对象分配内部数据的内存空间和默认值。需要注意这里背景模型使用到的累加器矩阵基本数据类型为cv::F32,当帧数不多时这没有什么问题,但是当帧数足够大时需要使用cv::F64避免数据溢出。

void AllocateImages(const cv::Size & size) {
    totalFrame = cv::Mat::zeros(size, CV_32FC3);
    averageFrame = cv::Mat::zeros(size, CV_32FC3);
    totalDifferenceFrame = cv::Mat::zeros(size, CV_32FC3);
    averageDifferenceFrame = cv::Mat::zeros(size, CV_32FC3);
    previousFrame = cv::Mat::zeros(size, CV_32FC3);
    upperLimitFrame = cv::Mat::zeros(size, CV_32FC3);
    lowerLimitFrame = cv::Mat::zeros(size, CV_32FC3);
    tmpFrame = cv::Mat::zeros(size, CV_32FC3);
    tmpFrame2 = cv::Mat::zeros(size, CV_32FC3);
    backgroundTempFrame = cv::Mat(size, CV_32FC1);
}

接下来定义一个函数来累加视频帧的像素强度,以及帧间偏差。注意这里使用的是计算速度稍快的帧间偏差,而不是标准差,在这个程序中它们的值相似,后面我们还会讨论到标准差。通常学习到帧数为30到1000,有时每秒视频数据中仅提取几帧,有时会提取全部可用的帧。

/// Learn the background statistics for one more frame
/// @param backgroundFrame backgroundFrame is a color sample of the background, 3-channel, 8u
void accumulateBackground(cv::Mat & backgroundFrame) {
    // nb. Not thread safe
    static int first = 1;
    // convert to float
    backgroundFrame.convertTo(tmpFrame, CV_32F);
    if (!first) {
        // 因为需要计算差分,因此首帧不处理
        totalFrame += tmpFrame;
        cv::absdiff(tmpFrame, previousFrame, tmpFrame2);
        totalDifferenceFrame += tmpFrame2;
        learnedImageCount += 1.0;
    }
    first = 0;
    previousFrame = tmpFrame;
}

当完成数据统计后,就可以使用均值的帧间绝对偏差均值(Frame-to-Frame Average Absolute Differences, FFAAD)计算上界和下界阈值,从而创建背景提取模型,其函数实现如下。

void updateHighThreshold(float scale) {
    upperLimitFrame = averageFrame + (averageDifferenceFrame * scale);
    cv::split(upperLimitFrame, channelHighFrames);
}

void updateLowThreshold(float scale) {
    lowerLimitFrame = averageFrame - (averageDifferenceFrame * scale);
    cv::split(lowerLimitFrame, channelLowFrames);
}

void createModelsfromStats() {
    averageFrame = totalFrame * (1.0 / learnedImageCount);
    averageDifferenceFrame = totalDifferenceFrame * (1.0 / learnedImageCount);
    
    // Make sure diff is always something
    averageDifferenceFrame += cv::Scalar(1.0, 1.0, 1.0);
    updateHighThreshold(highThreshFactor);
    updateLowThreshold(lowThreshFactor);
}

在上面的代码中,计算出背景阈值后,使用函数cv::split()将三个通道的值分别使用三个单通道矩阵存储,方便后续提取背景时逐通道比较。当背景模型训练完成后就可以用于分割位于两个阈值范围内对背景图像和剩下的前景图像,其函数实现如下。

// Byte, 1-channel image
cv::Mat backgroundTempFrame;
/// Create a binary: 0,255 mask where 255 (red) means foreground pixel
/// @param sourceFrame Input image, 3-channel, 8u
/// @param forgroundFrame Foreground image to be created, 1-channel 8u
void backgroundDiff(cv::Mat& sourceFrame, cv::Mat& forgroundFrame) {
    cv::Mat backgroundFrame = forgroundFrame;
    sourceFrame.convertTo(tmpFrame, CV_32F); // To float
    cv::split(tmpFrame, channelFrames);
    
    // Channel 1
    cv::inRange(channelFrames[0], channelLowFrames[0], channelHighFrames[0], backgroundFrame);
    // Channel 2
    cv::inRange(channelFrames[1], channelLowFrames[1], channelHighFrames[1], backgroundTempFrame);
    backgroundFrame = cv::min(backgroundFrame, backgroundTempFrame);
    // Channel 3
    cv::inRange(channelFrames[2], channelLowFrames[2], channelHighFrames[2], backgroundTempFrame);
    backgroundFrame = cv::min(backgroundFrame, backgroundTempFrame);

    // Finally, invert the results
    forgroundFrame = 255 - backgroundFrame;
}

上面的代码分别比较了视频帧的三个通道和背景提取模型的关系,并使用cv::inRange和一个临时矩阵变量backgroundTempFrame,这样位于阈值范围内的背景像素值就会被设置为255,前景像素值会被设置为0。最后合并三个通道的数据,假定任意通道被判定为前景,则该像素就是前景像素,当然你也可以使用逻辑异或运算符。最后反转背景矩阵就得到255表示前景的前景矩阵。

示例程序在主函数中组合了上述方法训练了一个背景模型,并对后面的视频帧提取前景像素将其标红。主函数的代码和一些其他细节这里就不再展示,可以前往链接的源码仓库查看。

这种简单的模型只适合不包含移动背景的场景,对包含如飘动窗帘、摇动树枝已经其他包含双或者多模态特征背景的场景效果并不好。使用该示例程序处理包含摇动树枝的背景训练得到的模型提取的前景效果如下图,其中左图中将原始图像中检测到是前景对象的像素红色通道强度设置为最大,右图是提取得到的背景像素二值图。

图11-3 使用平均法训练的背景模型效果图

4.2 累积均值

在使用平均法训练背景模型的过程中,我们计算了所有帧的像素强度从而求得平均值,实际上除了这种方式平均值还可以通过增量的方式计算,也就是说计算平均值可以不需要累积所有像素强度的和。但是在本章中我们并不严格区分增量计算的概念,这种通过累积所有像素强度和再处理样本数量的计算方式也称为增量计算,这种广义的增量计算方式效率会比狭义的增量计算方式(Nearly Incremental)稍快。

示例AverageMethod中计算平均值使用的方式是通过类cv::Mat重载的运算符+=累积所有帧的和,最后再除以帧的数量。而OpenCV还提供了如下累加函数。

// src:待累加的矩阵,1或者3通道,基础数据类型为U8或者F32
// dst:累加目标矩阵,基础数据类型为F32或者F64
// mask:蒙版矩阵,只有其中非零值对应的src像素值才会被累加到dst中
void accumulate(cv::InputArray src, cv::InputOutputArray dst,
                cv::InputArray mask = cv::noArray());

cv::Mat重载的运算符+=相比,该函数内部包含cv::Mat::convertTo()逻辑,并且它还支持使用蒙版矩阵对允许累加像素的位置精细控制。在训练背景模型的过程中,mask参数很有用,因为通常总有一些额外的信息会表明图像中的某一部分一定不是背景。例如在对高速公路或者其他颜色均一的场景训练背景模型时,通过颜色就可以判断出一些物体不属于背景。

另外一种平均值的计算方法是移动平均法(Running Average),其计算公式如下。其中acc(x, y)是像素P(x, y)的平均值,而a是一个常数,它可以被看作是控制先前帧的影响的衰减系数,其值越小,先前帧的影响衰减更快。

这种计算方式得到的值并不是严格的平均值,但是它允许每次计算都不需要依赖当前的累积量,从而实现增量计算。如使用一般法和运动平均法(这里令a = 0.5)计算样本[2, 3, 4]的平均值,则前者计算的结果为3,后者计算结果为3.25(0.5✖️2 + 0.5✖️3 = 2.5,0.5✖️2.5 + 0.5✖️4 = 3.25)。在这里运动平均法计算的结果更大是因为根据其定义,越新的点对计算结果的贡献越高。

OpenCV提供的运动平均值计算函数如下。

// src:当前样本矩阵,1或者3通道,基本数据类型为U8或者F32
// dst:运动平均值矩阵,基本数据类型为F32或者F64
// alpha:当前样本矩阵权重
// mask:蒙版矩阵,只有其中非零值对应的src像素值才会被计算
void accumulateWeighted(cv::InputArray src, cv::InputOutputArray dst,
                        double alpha, cv::InputArray mask = cv::noArray());

4.3 方差

下面介绍一个稍微复杂的模型,它通过平均值和方差建立高斯分布模型(Gaussian Model)来表示背景像素的变化行为。一维高斯模型使用一个期望(Mean)和方差(Variance)描述,其中期望表示背景像素强度的平均值,方差表示样本相对于期望的离散程度。D维高斯模型使用一个包含D个元素的向量表示期望,使用大小为D✖️D的协方差矩阵描述样本不同维度特征和期望的离散程度,以及不同维度之间的相关性。协方差到使用将在下一小节中计算。

在介绍使用方差法训练背景模型之前,先熟悉下方法的计算方式。在数学上,一个有限样本集合的方差计算公式如下。

其中N是样本个数,xi是第i个样本的值,减号后的表达式为集合所有样本的平均值。这个公式需要遍历集合两次,第一次计算平均值,第二次计算方差。因此通过简单的代数运算,上述公式可以转换为如下形式。

使用这种方法,通过一次便利计算所有样本的平方和以及和的平方就能计算出样本的方差。计算平方和最简单的方式是使用类似sqsum += I.mul(I)的方式,但是这种方式需要我们做额外的类型转换工作。此外,我们也可以使用如下OpenCV提供的平方和计算函数。

// src:当前样本矩阵,1或者3通道,基本数据类型为U8或者F32
// dst:平方和累加矩阵,基本数据类型为F32或者F64
// mask:蒙版矩阵,只有其中非零值对应的src像素值才会被计算
void accumulateSquare(cv::InputArray src, cv::InputOutputArray dst,
                      cv::InputArray mask = cv::noArray());

示例程序Variance使用方法法训练提取前景对象,其准备计算方差所需数据的代码如下。

// Store accumalative sum and square of correspoding pixels intensity in each video frame
cv::Mat sum, sqsum;
int trainedImageCount = 0;

/// Accumulate the information we need for our variance computation.
/// @param mat Image mat needed to accumulate into result.
void accumulateVariance(cv::Mat& mat) {
    if (sum.empty()) {
        sum = cv::Mat::zeros(mat.size(), CV_32FC(mat.channels()));
        sqsum = cv::Mat::zeros(mat.size(), CV_32FC(mat.channels()));
        trainedImageCount = 0;
    }
    cv::accumulate(mat, sum);
    cv::accumulateSquare(mat, sqsum);
    trainedImageCount++;
}

示例程序Variance中计算标准差的代码如下。

/// Compute standard deviation.
/// @param stDev Standard deviation that store the computed result.
void computeStdev(cv::Mat& stDev) {
    double scaleFactor = 1.0 / trainedImageCount;
    cv::sqrt(((scaleFactor * sqsum) - ((scaleFactor * scaleFactor) * sum.mul(sum))), stDev);
}

示例程序Variance的完整代码不在文中详细列出,可以参考链接。程序运行结果如下图,其中左图中将原始图像中检测到是前景对象的像素红色通道强度设置为最大,右图是提取得到的背景像素二值图。

图11-4 使用方差法训练的背景模型效果图

4.4 协方差

需要注意在示例程序Variance中并未使用到协方差数据,但是我们还是需要掌握该技巧。在训练背景模型时,当个通道的方差描述了背景像素在对应通道上强度的相似程度。而协方差描述的是背景像素各个通道强度之间的相互关系。例如当背景是大海时,海水的颜色红色通道的强度很低,它可以看作是由蓝色和绿色的混合色。并且这两个通道的强度会同时随着光照强度的改变而改变,即绿色通道的强度和蓝色通道的强度是成正比的。基于此可以判断对于新视频帧的某个像素,如果其蓝色通道增加,但是绿色通道并没有一起增加,则它就是前景像素。

下图是对不同光照条件下海水色的像素样本统计结果,这里只统计来蓝色通道(X)和绿色通道(Y)数据,其中图a使用方差来描述海水色背景区间,图b使用协方差来描述这个背景区间。另外,在图a中σxσy分布表示蓝色通道和绿色通道的方差,在图b中λ1λ2是协方差矩阵的两个特征值。可以明显看出使用协方差能够更好的描述海水的背景区域,其中囊括了更多的样本。

图11-5 使用方差和协方差法拟合样本集合

协方差矩阵的计算公式经过简单的代数运算后可以表示为如下形式。

在上述公式中,如果xy相等,即它们描述的是同一通道,则计算的结果将是该通道的方差。在d维空间上(如RGB颜色空间就是3维样本空间),通常使用d✖️d的协方差矩阵(Covariance Matrix)∑x,y来描述样本集的分布特征。该矩阵对对角线为每个维度的方差,而非对角线元素E(x, y)表示了xy维度的协方差。该矩阵是对称矩阵,即E(x, y) = E(y, x)

在本系列的前序文章中曾经介绍了函数cv::calcCovarMatrix(),它可以用于计算多个向量表示的单个集合数据,该集合由N个D维样本组成,计算结果为尺寸为D✖️D的协方差矩阵。但是在训练3通道彩色视频的背景模型时,我们需要计算M✖️N个集合数据,每个集合由FrameCount个3维样本组成,其中FrameCount是训练背景提取模型需要学习的视频帧数,3维表示的是RGB通道。或者对于每个协方差矩阵而言,至少需要计算其中的6个元素。

显然此时使用函数cv::calcCovarMatrix()并不是最好的选择,因为需要大量的内存空间来准备计算需要的数据,并且还需要重新组织数据结构。因此通常我们会结合上述公式和自定义的代码计算协方差。在此之前,先了解如下计算累积矩阵乘积函数。

// src1:当前样本矩阵1,1或者3通道,基本数据类型为U8或者F32、
// src2:当前样本矩阵2,1或者3通道,基本数据类型为U8或者F32
// dst:乘积累加矩阵,基本数据类型为F32或者F64
// mask:蒙版矩阵,只有其中非零值对应的src像素值才会被计算
void accumulateProduct(cv::InputArray src1, cv::InputArray src2,
                       cv::InputOutputArray dst, cv::InputArray mask = cv::noArray());

该函数不能自定义多通道输入矩阵在相乘时通道的对应关系,而我们需要计算的是不同通道之间的累积乘积,因此需要先使用函数cv::Split()将输入图像分割为单通道矩阵。准备计算协方差所需数据的相关代码如下。

// Store accumulative value of pixels intensity in each channel during model learning
std::vector<cv::Mat> channelSums(3);
// Store accumulative dot value of pixels intensity in each channel pair during model learning
// to be used calculate coavariance
std::vector<cv::Mat> xySums(6);
// Temporary mat use as data conationer of each channel pixels.
std::vector<cv::Mat> acCovarianceChannels(3);
void accumulateCovariance(cv::Mat& mat) {
    // Allocation operation should actually implemented at the begining of application,
    // and this form just show as example presentation.
    if (sum.empty()) {
        trainedImageCount = 0;
        for (int i = 0; i < 3; i++) {
            // the r, g, and b sums
            channelSums[i] = cv::Mat::zeros(mat.size(), CV_32FC1);
        }
        for (int i = 0; i < 6; i++) {
            // the rr, rg, rb, gg, gb, and bb elements
            xySums[i] = cv::Mat::zeros(mat.size(), CV_32FC1);
        }
    }
    
    cv::split(mat, acCovarianceChannels);
    for (int i = 0; i < 3; i++) {
        cv::accumulate(acCovarianceChannels[i], channelSums[I]);
    }
    
    int n = 0;
    // Coompute the accumulation of rr, rg, rb, gg, gb, and bb channel dot value.
    for (int i = 0; i < 3; i++) {
        for (int j = i; j < 3; j++ ) {
            n++;
            cv::accumulateProduct(acCovarianceChannels[i], acCovarianceChannels[j], xySums[n]);
        }
    }
    
    trainedImageCount++;
}

计算协方差的示例代码如下。

/// Compute covariance.
/// @param covariance Covariance mat  that store the computed result.
/// A six-channel array, channels store the covariance of rr, rg, rb, gg, gb, and bb channels.
void computeCoariance(cv::Mat& covariance) {
    double scaleFactor = 1.0 / trainedImageCount;
    
    // Compute the covariance of rr, rg, rb, gg, gb, and bb channels.
    int n = 0;
    for (int i = 0; i < 3; i++) {
        for (int j = i; j < 3; j++) {
            n++;
            xySums[n] = (scaleFactor * xySums[n]) - 
                        ((scaleFactor * scaleFactor) * channelSums[i].mul(channelSums[j]));
        }
    }
    
    // Reassemble the six individual elements into a six-channel array
    cv::merge(xySums, covariance);
}
4.4.1 背景像素评判方法

尽管在示例程序Variance中我们使用方差法提取了背景模型,但是判断背景像素使用的方法还是相对简单。在使用方差法提取背景模型时,对每个通道都计算互相独立的高斯分布模型,由于各个通道的方差通常不同,因此通常计算z值(z-score)描述特定像素属于特定分布的概率,在这里表示特定像素属于背景的概率。单个通道的z值计算公式为z = (x - x’) / σ,其中x是特定像素在某一通道的强度,x’σ分别是背景模型中该通道的均值和标准差。多通道的z值计算公式如下。

在使用协方差法提取背景模型时,通常计算马氏距离(Mahalanobis)描述特定像素和期望值的距离,它本质上是以恒定概率等高线测量点到均值的距离。在这里对于仅使用方差提取背景模型的场景,由于各个维度数据不相关,因此协方差矩阵为对角矩阵,因此马氏距离公式退化就是z值计算公式。参考图11-5,图a和图b中均值左上角同一个样本分别具有更小和更大的马氏距离。OpenCV中计算两个点的马氏距离的函数原型如下。

// 返回值:两点之间的马氏距离,数据类型为F64
// vec1:待计算的点A,N✖️1或1✖️N矩阵
// vec1:待计算的点B,N✖️1或1✖️N矩阵
// icovar:样本分布的逆协方差矩阵,N✖️N矩阵
double cv::Mahalanobis(cv::InputArray vec1, cv::InputArray vec2, cv::InputArray icovar);

这里使用逆协方差矩阵的原因是马氏距离内部需要使用到逆协方差矩阵,而通常我们会使用同一个协方差矩阵计算多组向量的马氏距离,而逆协方差矩阵计算成本较高。因此我们应当只计算一次逆协方差矩阵,并使用该矩阵多次调用上述函数。

在背景提取时使用该函数计算马氏距离并不便利,因为我们需要为每个像素调用一次该函数。遗憾的是在OpenCV中并没有该函数的矩阵版本。因此在实际工作中,训练模型是对于每一个像素都需要计算其协方差矩阵,再计算并存储其逆矩阵。在背景分割时,则需要遍历图像中的每个像素,使用对应的逆协方差矩阵和函数Mahalanobis计算它们和均值的马氏距离。

5 码书法

很多时候背景会包含复杂移动物体,如风中挥舞的树木,旋转的扇叶以及飘动的窗帘。另外这些复杂模型也可能包含光线的变,如飘动的云以及窗口或者门的开关都会导致光线发生变化。

一个处理这种问题的好的方式是基于时间序列为每个或者每组像素建立模型,即每帧图像独立建立模型,这将会消耗大量的内存。即如果为帧率为30FPS的视频文件截取2秒时长的资源训练背景模型,每个像素都需要60个样本。最后对于每个像素,通过60个不同的权重值对这些样本计算得到最终的模型。通常模型训练的时间会超过2秒,这种巨大的内存消耗在当前的硬件环境下是不切实际的。

为了接近自适应滤镜的性能,我们从视频压缩技术中获得灵感,尝试使用YUV码书(CodeBook)表示背景中的显著状态。此外背景分割的后期工作还允许任意的相机运动,以及使用均值位移算法进行动态建模,这里不详细介绍。最简单的方式是比较观察到的像素值和该像素在前一刻的值,如果它们很接近,则认为这是颜色扰动,如果相隔较远,则认为它是该像素作为背景的一组新的颜色。训练得到的结果可以理解为在RGB颜色空间上的一组浮动的颜色气泡,每个气泡都表示一组作为背景的颜色值。

在实际应用中选用RGB颜色空间并不是最优的,使用包含亮度轴的颜色空间总能得到更好的效果,例如YUV颜色空间。尽管YUV是最佳选择,但是如HSV等颜色空间也能获得类似的效果。这样做的原因是背景颜色通常只在亮度上发生变化,而不是“颜色”。

接下来应该考虑如何对每个颜色气泡建模,可以使用前文介绍的方式使用高斯簇,即使用期望和协方差概括每个气泡。事实证明计算使用最简单的方式,直接在每个轴上使用阈值来概括每个气泡也能够取得很好的效果。这种方式使用的内存成本更低,并且判断一个新观察到的像素是否位于学习到的颜色气泡,即颜色盒子内部的计算成本也会更低。

码书法的示意图如下,码书由一系列盒子组成,在学习的过程中这些盒子会不断增长以覆盖常见的颜色。图中上半部分描述了某个像素点亮度值随时间的变化曲线,下半部分演示了盒子的生长过程,首先只覆盖一小部分亮度区间,随着时间的变化,接近的亮度值会使盒子长大,而距离相差较大的亮度值会长出新的盒子。

图11-6 码书法示意图

在训练背景模型时,使用到的码书是由包含3个通道的盒子组成,这些通道对应了像素的不同颜色通道。下图是取六个不同的像素点学习到的亮度码书,当然这里只是示意,实际工作中会对每个像素学习到一个码书。这种方法能够很好的处理颜色剧烈变化的背景,如风中摇晃的树这一场景,某个背景像素可能是不同树叶的颜色,也可能是树后蓝天白云的颜色。通过这种更精确的方法,我们能够检测出颜色在这些背景颜色之间的前景物体。和图11-2平均法不能很好的区分手的颜色相比较,码书法具有更好的性能。

图11-7 使用码书法训练背景模型示意图

在使用该方法实际训练背景模型时,会在三个颜色轴上各计算出两个背景阈值(maxmin),如果新观察到的像素值位于背景阈值之外,但是位于训练阈值(learnHighlearnLow)之内,则码书向该值扩张,并更新背景阈值。如果新观察到的像素值位于背景阈值之外,也位于训练阈值之外,则新建一个盒子。挡模型训练完成后,在提取背景时还允许设置一个误差阈值maxModminMod,即如果待测试像素不在背景阈值内,但是在误差阈值内(max+maxModmin-minMod)内也会被判定为背景像素,通常误差阈值直接设置为0

这里不讨论使用云台摄像机(pan-tilt camera)拍摄大场景的情况,此时我们应当考虑学习到的模型和云台摄像机的角度。

5.1 结构

示例程序CodeBook使用码书法训练背景模型,并提取了前景对象。其中码书的结构定义如下。

// The variable t counts the number of points we’ve accumulated since the start or 
// the last clear operation.
// You need one of these for each pixel in the video image (rowXcol)
// 该类继承于标准模版库的向量,这样做效率比较低,实际工作中可以直接使用C语言风格的数组
class CodeBook : public std::vector<CodeElement> {
    public:
    // Count of every image learned on
    // 训练码书使用的总图片数
    int t;
    
    // count every access
    CodeBook() {
        t = 0;
    }
    
    // Default is an empty book
    // Construct book of size n
    CodeBook(int n) : std::vector<CodeElement>(n) {
        t = 0;
    }
};

码书是由很多个盒子构成的数据结构,其中每个盒子可以理解是一个背景像素的气泡,这些盒子的数据结构定义如下。CHANNELS通常为1或者3,分别用于处理只包含亮度通道的图像和YUV等亮度颜色通道。对于一个新的像素P,如果其任意一个通道的值p[i]不在区间[learnLow[i], learnHigh[i]]内,都会创建一个新的码书元素实例。stale记录了码书元素最常的不活跃帧数,可以在模型训练中定期清理掉常时间未更新的码书元素,因为它们可能是噪声,即偶尔飘过的前景物体。

// Here’s how the actual codebook elements are described:
class CodeElement {
    public:
        // High side threshold for learning
        unsigned char learnHigh[CHANNELS];
        // Low side threshold for learning
        unsigned char learnLow[CHANNELS];
        // High side of box boundary
        unsigned char max[CHANNELS];
        // Low side of box boundary
        unsigned char min[CHANNELS];
        // Allow us to kill stale entries
        // 最近更新对应的模型训练索引,即视频帧索引
        int t_last_update;
        // max negative run (longest period of inactivity)
        // 最大空闲时间,即在模型训练时最大经历多少帧仍未被更新
        int stale;

        CodeElement() {
            for (int i = 0; i < CHANNELS; i++) {
                learnHigh[i] = learnLow[i] = max[i] = min[i] = 0;
            }
            t_last_update = stale = 0;
        }
        
        CodeElement& operator = (const CodeElement& ce) {
            for (int i = 0; i < CHANNELS; i++) {
                learnHigh[i] = ce.learnHigh[i];
                learnLow[i] = ce.learnLow[i];
                min[i] = ce.min[i];
                max[i] = ce.max[i];
            }
            t_last_update = ce.t_last_update;
            stale = ce.stale;
            return *this;
        }
        
        CodeElement(const CodeElement& ce) {
            *this = ce;
        }
};

5.2 训练背景模型

对于待训练的图像而言,每个像素都需要学习到一个码书模型。即对于训练图片集合中的每一幅图片,其中的每一个像素都需要调用一次如下函数训练码书模型。

/// Updates the codebook entry with a new data point
/// Note: cbBounds must be of length equal to numChannels
/// Return CodeBook index
/// @param p incoming YUV pixel
/// @param c CodeBook for the pixel
/// @param cbBounds Bounds for codebook (usually: {10,10,10})
/// @param numChannels Number of color channels we're learning
int updateCodebook(const cv::Vec3b& p, CodeBook& c, int* cbBounds, int numChannels) {
    if (c.size() == 0) {
        c.t = 0;
    }
    // Record learning event
    c.t += 1;
    // SET HIGH AND LOW BOUNDS
    // 计算以当前像素为中心的盒子上下边界值
    unsigned int high[3], low[3];
    for (int i = 0; i < numChannels; i++) {
        high[i] = p[i] + *(cbBounds + i);
        if (high[i] > 255) {
            high[i] = 255;
        }
        low[i] = p[i] - *(cbBounds + i);
        if (low[i] < 0) {
            low[i] = 0;
        }
    }

    // SEE IF THIS FITS AN EXISTING CODEWORD
    int matchChannel;
    bool elementIsExist = false;
    for (int i = 0; i < c.size(); i++) {
        matchChannel = 0;
        for (int j = 0; j < numChannels; j++) {
            // Found an entry for this channel
            if ((c[i].learnLow[j] <= p[j]) && (p[j] <= c[i].learnHigh[j])) {
                matchChannel++;
            }
        }

        // If an entry was found
        if (matchChannel == numChannels) {
            elementIsExist = true;
            c[i].t_last_update = c.t;

            // adjust this codeword for the first channel
            for (int j = 0; j < numChannels; j++) {
                if (c[i].max[j] < p[j]) {
                    c[i].max[j] = p[j];
                } else if (c[i].min[j] > p[j]) {
                    c[i].min[j] = p[j];
                }

                // SLOWLY ADJUST LEARNING BOUNDS
                // 如果匹配的码书元素学习阈值位于当前像素为中心的盒子边界值外,则逐渐扩张匹配到的码书元素学习阈值
                if (c[i].learnHigh[j] < high[j]) {
                    c[i].learnHigh[j] += 1;
                }
                if (c[i].learnLow[j] > low[j]) {
                    c[i].learnLow[j] -= 1;
                }
            }
            break;
        }
    }

    // OVERHEAD TO TRACK POTENTIAL STALE ENTRIES
    for (int i = 0; i < c.size(); i++) {
        // Track which codebook entries are going stale:
        int negRun = c.t - c[i].t_last_update;
        if (c[i].stale < negRun) {
            c[i].stale = negRun;
        }
    }

    // ENTER A NEW CODEWORD IF NEEDED
    if (!elementIsExist) {
        // if no existing codeword found, make one
        CodeElement ce;
        for (int i = 0; i < numChannels; i++) {
            ce.learnHigh[i] = high[I];
            ce.learnLow[i] = low[I];
            ce.max[i] = p[i];
            ce.min[i] = p[i];
        }

        ce.t_last_update = c.t;
        ce.stale = 0;
        c.push_back(ce);
    }

    return int(c.size());
}

5.3 清除移动前景

在背景学习的过程中可能图像中偶尔会出现移动的前景物体,在学习到的码书中,这种情况的具体现象就是会得到一些很长时间都不会更新的码书元素,可以通过周期性调用如下函数移除这些噪声。在我们的示例程序中是在50帧图片训练完成后调用概函数。

// During learning, after you've learned for some period of time,
// periodically call this to clear out stale codebook entries
// 长时间未被更新的码书元素可以认为是背景中偶尔出现的前景物体,这些元素需要删除
/// return number of entries cleared
/// @param c Codebook to clean up
int clearStaleEntries(CodeBook &c) {
    // 设置一个不活跃的时间阈值
    int staleThresh = c.t >> 1;
    // 标记数组
    int *keep = new int[c.size()];
    int keepCnt = 0;
    
    // SEE WHICH CODEBOOK ENTRIES ARE TOO STALE
    for (int i = 0; i < c.size(); i++) {
        if (c[i].stale > staleThresh) {
            // Mark for destruction
            keep[i] = 0;
        } else {
            // Mark to keep
            keep[i] = 1;
            keepCnt += 1;
        }
    }
    
    // move the entries we want to keep to the front of the vector and then
    // truncate to the correct length once all of the good stuff is saved.
    int k = 0;
    int numCleared = 0;
    for (int i = 0; i < c.size(); i++) {
        if (keep[i]) {
            c[k] = c[I];
            // We have to refresh these entries for next clearStale
            c[k].t_last_update = 0;
            k++;
        } else {
            numCleared++;
        }
    }
    c.resize(keepCnt);
    delete[] keep;
    
    return numCleared;
}

5.4 提取前景对象

当通过上述函数训练好背景模型后,就可以通过如下函数提取图片中的前景对象。

/// Given a pixel and a codebook, determine whether the pixel is covered by the codebook
/// NOTES:
/// minMod and maxMod must have length numChannels,
/// e.g. 3 channels => minMod[3], maxMod[3]. There is one min and one max threshold per channel.
///
/// return 0 => background, 255 => foreground
/// @param p Pixel (YUV)
/// @param c Codebook
/// @param numChannels Number of channels we are testing
/// @param minMod_ Add this (possibly negative) number onto max level when determining 
/// whether new pixel is foreground
/// @param maxMod_ Subtract this (possibly negative) number from min level when
/// determining whether new pixel is foreground
uchar backgroundDiff(const cv::Vec3b& p, CodeBook& c, int numChannels, int* minMod_, int* maxMod_) {
    int matchChannel;
    
    // SEE IF THIS FITS AN EXISTING CODEWORD
    int I;
    for (i = 0; i < c.size(); i++) {
        matchChannel = 0;
        for (int n = 0; n < numChannels; n++) {
            if((c[i].min[n] - minMod_[n] <= p[n]) && (p[n] <= c[i].max[n] + maxMod_[n])) {
                matchChannel++; // Found an entry for this channel
            } else {
                break;
            }
        }
        // Found an entry that matched all channels
        if (matchChannel == numChannels) {
            break;
        }
    }
    
    // No match with codebook => foreground
    if (i >= c.size()) {
        return 255;
    }
    // Else background
    return 0;
}

示例程序CodeBook的运行效果如下图。

图11-8 使用码书法提取前景对象效果图

总的来说使用码书法分离前景对象需要如下几步。

  • 调用函数updateCodebook()从几秒或者几分的视频片段中训练得到一个基本的背景模型。
  • 调用函数clearStaleEntries()清除长时间未更新的码书元素,消除噪声。
  • 调整允许误差minMod和maxMod优化前景分割的效果。
  • 得到前景提取的最终模型(维护一个前文提到的高抽象层的场景模型)。
  • 使用函数backgroundDiff()和学习到的场景模型分割前景对象。
  • 周期性更新学习到的背景像素(更新场景模型)。
  • 调用函数clearStaleEntries()以更低的频率清楚长时间未更新的码书元素。

5.5 码书法的一些想法

总的来说,码书法在处理很多场景时都有相当不错的效果,它在训练和运行模型时都有较高的效率。但是它不能处理光照模式变化的场景,如清晨、中午和晚上不同的光照条件,或者是有人在室内打开或关上灯。我们可以使用几个不同的码书模型来处理这种场景,每种不同的光照模式对应一个码书模型,然后根据光照模式决定启用其中某个模型。

6 连通分量分析

6.1 联通分量分析与码书法

通过联通分量可以进一步处理提取到的前景图像,消除其内部噪声。这里首先使用本系列文章中滤波和卷积章节介绍到的图像形态学操作对图像预处理,具体方法是使用开操作移除高亮噪声,再使用闭操作填充小块阴影区域。接下来找到处理后图像中长度位于某个阈值之上的轮廓,或者是最大的轮廓,并计算一些诸如联通分量中心点以及围绕矩形等统计信息。

示例程序ConnectedComponents在示例程序CodeBook的基础上使用了联通分量分析的方式消除了前景对象中的噪声,其中联通分量分析的函数实现如下。

/// 寻找联通分量
/// @param mask Is a grayscale (8-bit depth) "raw" mask image that will be cleaned up
/// @param poly1_hull0 If set, approximate connected component by (DEFAULT: 1) polygon,
/// or else convex hull (0)
/// @param perimScale Len = (width+height)/perimScale. If contour len < this, delete that
/// contour (DEFAULT: 4)
/// @param bbs Ref to bounding box rectangle return vector
/// @param centers Ref to contour centers return vector
void findConnectedComponents(cv::Mat& mask, int poly1_hull0, float perimScale, 
                             std::vector<cv::Rect>& bbs, std::vector<cv::Point>& centers) {
    bbs.clear();
    centers.clear();
    
    // 1 使用基础形态学操作去除噪点
    // CLEAN UP RAW MASK
    cv::morphologyEx(mask, mask, cv::MORPH_OPEN, cv::Mat(), cv::Point(-1,-1), CVCLOSE_ITR);
    cv::morphologyEx(mask, mask, cv::MORPH_CLOSE, cv::Mat(), cv::Point(-1,-1), CVCLOSE_ITR);
    
    // 2 提取高于阈值点有效轮廓,并根据参数设置使用DP法或者轮廓突包法计算近似多边形,减少算法复杂度
    // FIND CONTOURS AROUND ONLY BIGGER REGIONS
    // all contours found
    std::vector<std::vector<cv::Point>> contours_all;
    std::vector<std::vector<cv::Point>> contours;

    // just the ones we want to keep
    cv::findContours(mask, contours_all, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_SIMPLE);

    std::vector<std::vector<cv::Point>>::iterator call_i = contours_all.begin();
    for (; call_i != contours_all.end(); ++call_i) {
        // length of this contour
        int len = cv::arcLength(*call_i, true);
        
        // length threshold a fraction of image perimeter
        double q = (mask.rows + mask.cols) / DP_EPSILON_DENOMINATOR;
        // If the contour is long enough to keep...
        if (len >= q) {
            std::vector<cv::Point> c_new;
            if (poly1_hull0) {
                // If the caller wants results as reduced polygons...
                cv::approxPolyDP(*call_i, c_new, len/200.0, true);
            } else {
                // Convex Hull of the segmentation
                cv::convexHull(*call_i, c_new);
            }
            contours.push_back(c_new);
        }
    }

    // Just some convenience variables
    const cv::Scalar CVX_WHITE(0xff, 0xff, 0xff);
    const cv::Scalar CVX_BLACK(0x00, 0x00, 0x00);
    
    // 3 计算有效轮廓的几何中心以及包裹矩形
    // CALC CENTER OF MASS AND/OR BOUNDING RECTANGLES
    int idx = 0;
    cv::Moments moments;
    cv::Mat scratch = mask.clone();
    std::vector<std::vector<cv::Point>>::iterator c_i = contours.begin();
    for (; c_i != contours.end(); c_i++, idx++) {
        cv::drawContours(scratch, contours, idx, CVX_WHITE, cv::FILLED);
        
        // Find the center of each contour
        moments = cv::moments(scratch, true);
        cv::Point p;
        p.x = (int)(moments.m10 / moments.m00);
        p.y = (int)(moments.m01 / moments.m00);
        centers.push_back(p);
        bbs.push_back(cv::boundingRect(*c_i));
        scratch.setTo(0);
    }

    // 4 绘制所有有效轮廓
    // PAINT THE FOUND REGIONS BACK INTO THE IMAGE
    mask.setTo(0);
    cv::drawContours(mask, contours, -1, CVX_WHITE, cv::FILLED);
}

该函数依赖两个宏定义的变量DP_EPSILON_DENOMINATORCVCLOSE_ITR,其中DP_EPSILON_DENOMINATOR控制使用DP算法计算近似轮廓时,需要得到的轮廓简单程度,该值越大,得到的轮廓越简单,算法开销越低。CVCLOSE_ITR用于控制形态学操作的迭代次数,如在开操作中先执行CVCLOSE_ITR次腐蚀操作,再执行CVCLOSE_ITR次膨胀操作。该值的设置取决于图像的分辨率,如果图像的分辨率很高,该值设置为1时可能并不能得到很好的效果。当然该值越大,算法的开销也越高,但是能够清除的噪声块区域也越大。

在上述代码中使用到了OpenCV提供的函数cv::findContours()寻找较大面积的联通分量,在OpenCV3版本以后还可以在该函数之前调用函数cv::connectedComponentsWithStats()标记联通分量,并删除其中面积较少者。示例程序ConnectedComponents运行结果如下图。

图11-9 使用码书法提取前景对象效果图

6.2 联通分量分析与帧间差分法

除了使用码书法外,还可以使用前文提到的帧间差分法,它是相当简单的一种背景提取方法。它通过计算时间轴上不同的两帧画面,使用较晚的图片减去较较早的图片计算其差值,再使用阈值处理得到的结果,从而提取前景对象。一般情况下连续的视频帧很相似,除非场景中包含快速移动的前景对象。但是相邻的视频帧相减的结果也会得到很多噪声,因此实际使用帧间差分法需要处理的主要问题是如何清除这些噪声。

为了更好的理解这些噪声,首先考虑如下图的两个不包含前景对象的视频帧差分结果。其中左上的图片是某个时刻的视频帧,而右上角的图片是上一时刻的视频帧。它们直接相减的结果使用15为阈值处理后得到左下的图片,可以明显看到移动的树叶产生了大量的噪声。连通分量分析法能够很好的处理这些噪声,处理结果如右下图。噪声信号是空间无关的,它们通常无规律的分布在图像中,因此通过设置阈值删除小的联通分量能够达到清理噪声的目的。

图11-10 联通分量分析不包含前景对象的帧间差分结果

联通分量分析的阈值通常需要调整为对于无前景对象的帧间差分结果其给出的处理图是不包含任何信号的,即全是黑色部分。通常这样处理后能够得到较好的前景对象分割效果。

接下来考虑包含前景对象的场景,当手从左向右晃动时,获取到当前视频帧如下左上图,以及上一时刻的视频帧如下右图。和前面的例子类似,计算它们的差值并通过阈值处理得到如下左下图,最后通过联通分量分析得到右下图。

图11-11 联通分量分析帧间差分结果

可以明显看到帧间差分法的一个缺点,就是不能够明确的区分出前景物体移动前的位置和物体当前的位置,并且在移动前后位置重叠的地方通常会形成了一个空洞,因为相同的肤色相减结果为0,或者说位于阈值之下。

7 方差法及码书法的比较

本章主要介绍了两种背景模型训练方法,分别是平均法以及其包含方差法等变体,和码书法。决定使用何种方法等最好方式就是暴力的将所有可用的方法都尝试一次。这里继续使用窗外飘动树叶的视频为例,在该视频中除了晃动的树外,还有来自于右侧建筑物和左侧内墙壁折射出的光,这使得场景建模更具挑战性。

分别使用平均法和码书法训练背景模型,提取到的前景对象如下图。上面两幅图分别是使用平均法提取到的前景图像,和使用联通分量分析后得到的结果。下面两幅则是使用码书法提取得到的前景图像以及使用联通分量分析后得到的结果。可以明显看到平均法得到了一个更稀疏的结果,并且将识别到的前景对象切割为两部分。因为平均差分建立的背景模型会将部分手的颜色值包含在模型内部。而码书法能够更准确的对树叶和树枝晃动带来的背景像素波动建模。

图11-12 平均法和码书法提取前景对象的比较

8 OpenCV中背景提取的封装

前文介绍了背景提取的一些基本方法,OpenCV在这些方法上封装了一系列工具类,它们包含更复杂的逻辑,也能更好的完成背景提取任务。OpenCV首先定义了一个背景提取的接口类,目前为止它有两个实现,在以后可能会有更多。这两个实现都基于高斯混合模型(Mixture of Gaussians, MOG),它本质上使用的是前文介绍的简单背景建模方法的统计学原理,如累积平均值、方差和协方差,再辅以码书法的思想。两个算法都是21世纪提出,能够应对大多数日常实际任务。

8.1 接口基类

接口基类cv::BackgroundSubtractor的定义如下。该基类只定义了两个方法,第一个用于分析一张新的图像,并计算出其对应的前景图像。第二个函数用于获取背景图像。看上去该类并未定义任何累加图像用于训练背景模型的方法,而只定义了提取前景对象的方法,这似乎是设计上的一个漏洞。实际上方法apply既用于学习背景模型,又用于提取前景图像。学术界有一个共识,那就是任意背景提取算法都应当具备连续学习的能力。这样要求的原因有很多,如停车场内车子停好一段时间后就应当被认为是背景。类似的例子还很多,因此该基类在设计时就不再区分训练模式和工作模式了。

class cv::BackgroundSubtractor {
public:
    virtual ~BackgroundSubtractor();
  
    virtual void apply() (cv::InputArray image, cv::OutputArray fgmask, double learningRate = -1);

    virtual void getBackgroundImage(cv::OutputArray backgroundImage) const;
};

8.2 KaewTraKuPong和Bowden方法

KaewTraKuPong和Bowden(KB)算法在背景分割领域带来了一些处理真实挑战的新能力。它们分别是多重模型、连续在线训练、两个提高初始化性能的分离式自动训练模式、显示阴影检测和阴影去除能力。这些能力对于用户几乎都是不可见的,但是不出意料的,OpenCV提供了一些参数用于控制该算法的实现。它们是历史参数(History)、高斯模型的数量(Number of Gaussian Mixtures)、背景比率(Background Ratio)和噪声强度(Noise Strength)。

第一个参数history用于指定算法初始化阶段的帧数,默认值为200。尽管现代背景提取算法不区分训练模式和工作模式,但是它们仍然需要一个初始化过程。

第二个参数高斯模型数量nmixtures控制了背景建模时每个像素需要的高斯模型数,默认值为5。

对于模型中的每个高斯分量都会包含一个权重值,表示该高斯分量的贡献像素观察值数量在该像素的所有观察值数量里面的比例。由于有时场景内会包含移动的前景物体,因此它们会生成权重较低的高斯分量。根据权重对这些高斯分量进行排序,其权重累加值称为背景比例backgroundRatio,通过特定参数控制,期默认值为0.7。即对于某个像素的背景模型包含5个高斯分量,其权重分别为0.4、0.25、0.2、0.1、0.05,则当背景比率为默认值0.7时,只有前3个模型才会被认为是背景。

最后一个参数噪声强度noiseSigma控制了新高斯分量创建时其不确定度。当某个像素观察到的新值无法被已有的模型解释时,会生成新的模型。如果此时已有的高斯分量未达到最大值就会创建新的高斯分量,如果已经达到最大值,则会循环使用最没有价值的高斯分量。当增加噪声强度时,会提高高斯分量能够表示的信息,即能够包含更多的值。当然这其中的权衡是高噪声强度会使高斯模型可能能够表示未被观察到的值。噪声强度的默认值为15,单位为像数强度,像素强度取值区间为[0,255]。

KB算法的实现是类cv::BackgroundSubtractorMOG,除了使用继承自父类对构造函数外,在opencv_contrib库的bgsegm模块中定义了如下构造函数来创建该类的实例。

cv::Ptr<cv::bgsegm::BackgroundSubtractorMOG> cv::bgsegm::createBackgroundSubtractorMOG(
    int history = 200, int nmixtures = 5, double backgroundRatio = 0.7, double noiseSigma = 0);

8.3 Zivkovic 方法

和KB方法类似,Zivkovic方法也适用高斯混合模型(Gaussian Mixture Model)对像素的颜色分布进行建模。它和KB法最大的区域是它会基于具体的统计分布动态调节每个像素建模适用的高斯分量个数,而KB法是固定的。这样也会带来一个缺点,即使用的高斯分量越多,在更新和比较时消耗的计算资源也越多,当然它也会使得训练得到的模型更准确。

该算法引入了一些新的参数,但是总的来看它所有参数中只有两个非常重要,其他的都可以直接使用默认值。这两个重要的参数分别是历史参数(history),也叫衰减参数(Decay Papameter),和方差阈值(Variance Threshold)。

其中参数history控制了对于特定像素新的观察值的过期时间,在这个时间内它都会对背景模型有所贡献,当超过这个时间则不再考虑,其默认值是500帧。实际上算法内部是以指数形式计算每个观察值的影响系数的,即参数history为500时,像素的衰减系数a - 1 / 500 = 0.002,则新观察值的影响系数计算公式为(1 - a)^t。

参数方差阈值varThreshold决定了像素的一个新观察值被判断为一个已经存在的高斯混合模型的置信水平,其单位是平方马氏距离。即当你需要设置距离分布中心3个单位马氏距离的观察值都应该认为是该分布的一部分,大于此值的不属于分布时,改参数需要设置为9,其默认值为4✖️4 = 16。

Zivkovic算法的实现类为cv::BackgroundSubtractorMOG2,其构造函数如下,在构造函数中可以设置3个重要的参数,其他参数可以在构建实例后设置。

Ptr<BackgroundSubtractorMOG2> createBackgroundSubtractorMOG2(
    int history = 500, double varThreshold  = 16, bool detectShadows = true);

参数detectShadows用于开启阴影检测和剔除,当其设置为true时,算法在检测到某个像素观察值不属于明显的背景时,还会考虑它是否是颜色变暗的背景。如果检测到属于这种情况,则会使用特殊的值标记这些像素,从而和前景和背景像素区分。

另外你还可以通过以下参数更个性化的定制算法的实现逻辑。

class cv::BackgroundSubtractorMOG2 {
...
public:
...
    // 每个像素使用的高斯混合模型最小数量,默认值为5
    // 该值越高,训练的模型越准确,但是计算成本越高
    int getNMixtures() const;
    void setNMixtures(int nmixtures);

    // 背景比率,其含义和KB算法内相同,默认值为0.9
    double getBackgroundRatio() const;
    void setBackgroundRatio(double backgroundRatio);

    // 高斯混合模型分量的初始方差,默认值为15
    // 该值含义和KB算法内的参数noiseSigma相似
    double getVarInit() const;
    void setVarInit(double varInit) const;

    // 高斯混合模型分量的最小允许方差,默认值为4
    double getVarMin() const;
    void setVarMin(double varMin);

    // 高斯混合模型分量的最大允许方差,默认值为75
    double getVarMax() const;
    void setVarMax(double varMax);

    // 复杂度缩减阈值,即生成新的高斯分量需要的最低像本数最,默认值为0.05
    // 当该值设置为0时将会显著简化算法,但是得到的结果质量也相应降低
    double getComplexityReductionThreshold() const;
    void setComplexityReductionThreshold(double CT);

    // 是否检测阴影
    bool getDetectShadows() const;
    void setDetectShadows(bool detectShadows);

    // 在输出的蒙版中表示阴影像素的值,默认值为127
    int getShadowValue() const;
    void setShadowValue(int value);

    // 阴影阈值,可以认为是某个颜色被认为是已有模型表示颜色的阴影时,其亮度百分比值下限,其默认值为0.5
    // 即当该值设置为0.6时,表示如果待比较多颜色和已有的模型具有相同的颜色,
    // 但是其亮度值介于0.6到1.0倍时,会被标记为已有模型的阴影
    double getShadowThreshold() const;
    void setShadowThreshold(double shadowThreshold);
...
};

示例程序Subtractor使用Zivkovic方法提取了如下场景中的前景对象,这里只是简单的演示,为了获得更好的效果,你需要自行优化算法的参数。

图11-13 Zivkovic法提取前景对象效果图

9 小结

本章主要介绍了背景分割技术,该技术在工业自动化、安全及机器人领域都扮演着重要的角色。在介绍完背景分割概念后,使用两个模型论述了背景分割技术是如何基于简单的统计学方法完成的。接下来演示了如何使用联通分量分析优化背景分割的结果,并比较了两个基本的背景分割方法效果。

随后介绍了更高级的背景分割方法,它们通过OpenCV库内的工具类提供。它们基于我们在前半部分文章中介绍的那些背景分割基本方法,但是包含更复杂的逻辑,具有更好的效果,更适合处理充满挑战的真实背景分割案例。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,271评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,275评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,151评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,550评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,553评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,559评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,924评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,580评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,826评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,363评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,940评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2 342