一、原理

用GMM对背景建模的基本思想是把每一个像素点所呈现的颜色用K个高斯分布的叠加来表示，通常K取3-5之间。将像素点所呈现的颜色X认为是随机变量，则在每时刻t=1,…,T所得到视频帧图像的像素值只是随机变量X的采样值（观值）。

在进行前景检测前，先对背景进行训练，对每一帧图像中每个背景采用一个混合高斯模型进行模拟，背景一旦提取出来，前景的检测就简单了，检查像素是否与背景的高斯模型匹配，匹配是背景，不匹配就是前景。

所以关键就是混合高斯背景模型的建立。

GMM之所以能够将前景和背景分开是基于如下两点事实的：

(1)在长期观测的场景中，背景占大多数时间，更多的数据是支持背景分布的

(2)即使是相对颜色一致的运动物体也会比背景产生更多变化，况且一般情况下物体都是带有不同颜色的。

二、算法流程

1. 参数初始化：定义对每个像素建立几个单高斯组成混合高斯模型，一般去3到5，我取4，变量定义为C=4。另外每个高斯分布需要用三个变量表示，分别是权重w，均值mean和方差pixel_sd。对图像的每个像素都需要建立C个高斯分布，所以这三个高斯分布参数变量是三个height*width*C大小的三维矩阵。初始化时，令mean在像素取值范围内去随机数，w取1/C，pixel_sd人工进行初始化，我取值为9。另外几个参数后面提到再详细说明。
2. 一次读取一幅图像（假设像素是单通道，多通道的图像分别计算），对每个像素计算与本像素各个高斯的均值的差的绝对值，然后计算这个差值是否小于对应高斯的方差的D倍，这个是为了容忍噪音的，否则当北京训练趋于稳定的时候，方差会很小，很容易就与当前高斯模型不匹配。D通常定义为2.5。对于匹配的高斯模型参数进行如下调整（其中p是alpha/w，alpha是人为定义的学习速率，用于更新权重，p用于更新高斯模型，除以w的作用是使得权重越大的模型越稳定，而权重过小的新定义的模型可以更快的收敛）：
   

对于所有的高斯模型的权重进行如下更新：

另一种参考文章里的更好的更新方式：

M是match的意思，也就是说如果当前像素与对应高斯匹配，M=1，否则M=0。

3. 我们发现进行上述更新后，匹配到的高斯权重会变大，没有被匹配的高斯权重会变变小，所以更新完成后，权重值和可能不唯一，所以要进行权重的归一化。
4. 然后对每个像素的所有高斯进行rank的计算，rank=w/pixel_sd，也就是说rank分高意味着权重大、方差小。如果某个像素点的所有高斯都没有被匹配，就说明那个像素的模型建立得不好，就把rank最小的高斯模型换成一个匹配当前像素值的高斯，即mean设为当前高斯像素值，方差用人工定义的值初始化，w不变。
5. 然后进行北京的选择，这里需要定义一个阈值thresh，表示所有高斯模型中有多少比例构成北京，按照rank顺序将w进行累加，知道前B个高斯的权重满足thresh，就完成北京的选择。

三、实验

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

clear all
close all
clc
% ----------------------- 使用第一帧获取图片大小 -----------------------

addpath('./data');
img_dir = dir('./data/*.jpg');
fr = imread(img_dir(1).name);           % 读第一幅图片
[height,width] = size(fr(:,:,3)); 

% ----------------------- 变量初始化 -----------------------------------

tolerance = 1;
C_max = 4;                       % 组成混合高斯模型的最多单高斯数目
D = 2.5;                     % 阀值（一般为2.5个标准差）
alpha = 0.01;               % 学习速率（一般为0到1）
thresh = 0.6;               % 前景阀值（0.25到0.47）
sd_init = 9;                 % 初始化标准差

% -------------------------初始化均值,权值和方差----------------------------
 
pixel_depth = 8;                        %  像素深度为8位
pixel_range = 2^pixel_depth - 1;         %  灰度范围

w=ones(height, width, C_max)/C_max;
pixel_sd=ones(height, width, C_max)*sd_init;
mean=rand(height, width, C_max)*pixel_range;
bg=zeros(height, width,3);
% -------------------------用于调试-------------------------
mean_history=zeros(C_max, 200);
pixel_sd_history=zeros(C_max, 200);
w_history=zeros(C_max, 200);
sort_history = zeros(C_max, 200);
match_history = zeros(C_max, 200);
%--------------------- 处理每一帧 -----------------------------------
for c =1:3
    for n = 1:150,


        % ----------------- diff ------------------
        fr = imread(img_dir(n).name);       

        % 把图像赋值C_max层，便于后面的矢量化操作
        pixel=zeros(height,width,C_max);
        for i = 1:C_max
            pixel(:,:,i) = fr(:,:,c);
        end

        % 计算像素差的绝对值
        u_diff = abs(double(pixel) - double(mean));



        % 获得针对每个像素各个高斯分布的c_max层的match矩阵，为1的元素是匹配到的高斯成分
        match = (u_diff <= (D*pixel_sd+tolerance)) ;
        for i = 1:C_max
            match_history(i,n) = sum(sum(match(:,:,i)))/(576*720);
        end
        % 根据match更新各个高斯分布参数
        p = ones(height,width,C_max)*alpha./w;
        w = w*(1-alpha)+match*alpha;
        delta = pixel.*match + (1-match).*mean;
        mean = (1-p).*mean + p.*delta;
        delta = sqrt((1-p).*pixel_sd.^2 + p.*(double(pixel) - mean).^2);
        pixel_sd = match.*delta +(1-match).*pixel_sd;

        % 获得针对每个像素各个高斯分布的c_max层的match矩阵
        % 为1的元素是说明对应位置的像素没有任何一个高斯成分被匹配到，需要把评分最小的高斯成分进行替换
        unmatch=sum(match,3)==0;


        % 计算排序的rank矩阵
        rank=w./pixel_sd;
        [min_rank,min_rank_index] = min(rank,[],3);

        % 根据min_rank_index生成与w和pixel_sd大小一致的索引old_index，索引矩阵素为1的位置是需要被替换掉的高斯分量
        old_index=ones(height,width,C_max);
        for i =1:C_max
            old_index(:,:,i) = ((old_index(:,:,i)*i).*unmatch) == min_rank_index;   
        end

        % 对需要替换的高斯分量进行替换
        mean = (1-old_index).*mean + old_index.*double(pixel);
        pixel_sd=(1-old_index).*pixel_sd + old_index.*(ones(height,width,C_max)*sd_init);

        % 权值归一化   
        sum_w=sum(w,3);
        for i=1:C_max
            w(:,:,i) = w(:,:,i)./sum_w;
        end


        rank=real(w./pixel_sd);    % 不知道为什么算出来时虚数，但是虚部为零，用real处理一下
        [sort_rank, sort_rank_index] = sort(rank, 3, 'descend');
        background = zeros(height,width);
        now_thresh = zeros(height,width);
        bg_match = zeros(height,width,C_max);
        for i = 1:C_max
            [x,y]=meshgrid(1:width,1:height);
            now_index = (sort_rank_index(:,:,i)-1)*height*width+(x-1)*height+y;
            bg_match(:,:,i) = now_thresh<=thresh;
            now_thresh = now_thresh + w(now_index);
            background = background + mean(now_index).*w(now_index).*bg_match(:,:,i);
            sort_history(i,n) = sum(sum(sort_rank_index(:,:,i)));
        end

        bg_match = (1-match).*(1-bg_match);
        bg_match = sum(bg_match,3)>0;



        % -------------统计---------------
        for i = 1:C_max
            w_history(i,n) = sum(sum(w(:,:,i)))/(576*720);
            mean_history(i,n) = sum(sum(mean(:,:,i)))/(576*720);
            pixel_sd_history(i,n) = sum(sum(pixel_sd(:,:,i)))/(576*720);

        end

        
        imshow(uint8(background))
        drawnow
        n

    end
    bg(:,:,c)=background;
end
imwrite(uint8(bg),'bg2.jpg')

视频截图：

选择150张图片进行背景建模，分通道提取单色背景，然后将三个颜色进行叠加，背景如下：

四、调参

1. 矢量化编程

开始的时候使用for循环对每张图片每个像素进行上述计算，在matlab上一次处理200张图片需要2个小时左右，这使得尝试新参数变得很困难，因为一个下午满打满算也只能尝试两三组参数。

后来在深度学习相关资料上看到了矢量化编程的概念，就是在python或者matlab等类似较高层语言上进行数学计算的时候，尽量不要使用for循环，要利用内建的适量运算来达到并行，然后按照矢量化的方法进行重写之后，一次运行需要5到10分钟左右。

具体方法是一幅图像看成一个height*width*C_max的图像，然后直接对这个三维数组进行运算，如果涉及到匹配与不匹配的运算不同的情况，添加一个同样大小的match数组来表明当前位时候需要计算。

2. 关于调参
   

如上图示第一个高斯配匹配的命中率，其他几个高斯分布大致如此不过第一个的命中率最高。关于阈值，基本上最大的高斯分布的阈值在0.3-0.5左右，如果阈值过小，会只取到一个高斯，阈值过大又会取到不属于背景的部分。最后经过尝试thresh=0.6

参考：前景检测算法GMM

一个简单的问题

问从一个横边长为m，竖边长为n方格棋盘的一个角走到对角点的最短路径共有多少种方法？这个问题应该会有很多人见到过，想明白的话这也是一个很简单的问题，但是如果在这个问题上进行一点小小的改变，这个问题也是可以很有意思的。

level 0

方便起见就不失一般性的假设棋盘大小是4*3：

那么问题就是从左下角到右上角的路经数目有多少？
一种思路就是用动态规划的思想，到格子（i，j）处的只有两种方法：从（i-1，j）向右走一步或者从（i，j-1）处向上走一步。用N（i，j）表示从起点到（i，j）的路径数，则：

N（i，j）=N（i-1，j）+N（i，j-1）

所以这个3*4的棋盘问题可解：

但是说的高大上是动态规划算法，说的不好听了不就是傻算嘛……换一个角度考虑就能转化成一个高中就能瞬间写出答案的问题，其实就是一种选择的问题。因为在每个路口都要做出选择，是向上还是向右走。就好像我高中骑车上学一样。假设我家在左下角（0，0），学校在右上角（4，3）。那么我可以先向右骑一段路，到路口（1，0），我可以选择向东直行到（2，0）或者左拐到（1， 1）处，每一个路口都要做出选择，但是全程我必须也只能向北走三段路，同理向东走四段路，即在全程长为7的路程中选择三段路向北走，四段路向东走。
例如（U U R U R R R ）或者（R U R U R U R）或者（R R R R U U U）都是可能的路径。（U表示向上走，R表示向右走）。
所以总路程数就是：

C(7,3) = C(7,4) = 35

同理从（0，0）点走到（m，n）处的路线数就是C（m+n，n）或者C（m+n，m）。

level 1

好！小棋盘升级了！现在我们的问题改变了，在一个mm正方形棋盘上从角点（0，0）走到对角线上的最短路线有多少种？

升了一级也没感觉怎么样嘛！每一个路口有2种向上向右两种走法，走m步，共$2^m$种路线。就拿这个99的图来说，答案就是$2^9$呗！so easy！我想说的是，除了这一个思路，用level 0的想法解这道题。
从上往下，到达第一个点的路径数目是C（9，0），到达第二个点的路径数目是C（9，1），依次类推，到达第i个点的路径数目是C（9，i-1），这里 i 的取值区间是大于等于0，小于等于9。也就是说，总共的路径数是：

$\sum_{i=0}^9C(9,i)$

这十个值分别是1,9,36,84,126,126,84,36,9,1，而如果把这个计算过程写下来的话，我们会发现这不就是杨辉三角嘛！

—（图片来自百度，侵删）—

Read More

参数问题

（自己挖的坑，跪着也要写完！(╯‵□′)╯︵┻━┻）

理论知识

前面说了HMM的解码和评估，但是这都是建立在我们有了HMM的正确参数的前提条件下。然而如果我们真的想解决实际问题，比如动作识别，语音识别等问题时，很难有最优的参数用来预设值。如果我们有先验知识，可以在初始化的时候使参数向最优值靠近一些，如果没有的时候，甚至我们只能用将概率1均分来赋初值。那么这种时候怎么办，连模型参数都没有，怎么评估，怎么解码？

在第一篇文章中，讲前向算法的地方，我们假设，我们统计了一学期的食堂后厨的轮班记录，也就是说我们有一个小间谍，使得我可以掌握隐藏状态的序列，那么整个模型的状态图就是完全可见的了，只要数据够多，我只需要用最大似然估计去统计频率就可以去无限的逼近真实参数了。

现在问题就是我们无法拿到内部状态序列的数据，无法求HMM的参数Pi，A，B，也就无法进行前向变量、后向变量和维特比算法的计算。这个问题看似无法解决，但是Baum大神还是找到了一个解决这个问题的方法，这类算法都属于一类叫做EM算法的范畴。
我对这类算法的浅显的理解是，在数据量足够的情况下通过一次次的递归来笔记呢模型中的特定隐含参数的真实值。首先随意对内部参数赋值，然后得到一种解释输出的分布，然后对新得到的分布求内部参数的期望，得到新的内部参数值，在重新对可能的分布进行估值，以此类推。我的解释可能不到位，引用HMM学习最佳范例中七中的解释：

直观地理解EM算法，它也可被看作为一个逐次逼近算法：事先并不知道模型的参数，可以随机的选择一套参数或者事先粗略地给定某个初始参数λ0 ，确定出对应于这组参数的最可能的状态，计算每个训练样本的可能结果的概率，在当前的状态下再由样本对参数修正，重新估计参数λ ，并在新的参数下重新确定模型的状态，这样，通过多次的迭代，循环直至某个收敛条件满足为止，就可以使得模型的参数逐渐逼近真实参数。

对EM算法的一个简单实现是聚类算法中的K-means算法，实现起来比较简单，可以直观的理解算法背后的EM思想。

那么在Baum-Walch算法中怎么实现了EM思想求内部参数呢，这个算法中定义了两个关键变量，是这个算法的重中之重，符号是\(\xi\)和\(\gamma\)，我在后面就称之为gamma变量和xi变量。

首先解释gamma变量，它的计算方法是：

用图来表示就是这个意思：

在前面说过，前向变量表示的是在t时刻以状态Si结束时总的概率（给定观测序列\(O_{1…t}\)），也就是计算的是图中的红色框部分。后向变量的意思就是假定时间t时的内部状态是Si，模型输出\(O_{t+1…T}\)的概率，也就是计算上图中的蓝色圆圈中的部分。那么在此基础上，gamma变量的分子的意思就好理解了，它表示的是，在输出完整的观测序列\(O_{1…T}\)的条件下，t时刻的内部状态是Si的概率。而分母就是在输出完整的观测序列的情况下，所有可能的内部状态的和，也就是我们在第一篇文章中说的，整个序列的概率。所以gamma变量的意义就是t时刻出现内部状态Si的概率，这个概率在t时刻是归一的，就是说t时刻的所有gamma变量相加为一。

然后xi变量比gamma变量稍微复杂一点，但是基本思想是一样的，表示的是：由t时刻到t+1时刻出现内部状态Si到Sj的这种转移的概率，公式如下：

Read More

第二个问题：Viterbi算法

土方法

现在来到了第二个问题，我吃了三天食堂都是红烧肉，我想知道这三天的大厨最有可能是谁，怎么办？直观的思路还是使用枚举法（没办法啊。。。我总是想出这么简单naive的小白算法）。但是前面已经讨论，这种算法无论是从效率上还是逼格上都太低，因此我们考虑其他的方法。

还算好但是不常用的方法

我们上一篇文章提到了我们可以用前向变量和后向变量搭配使用求解t时刻内部变量是\(S_i\)且输出观测序列O的概率，那么算出t时刻全部的内部状态的各自的概率，那么概率最大的那个就是t时刻最有可能的内部状态，对每个时刻都进行这样的计算不就可以得到最有可能的内部序列的吗？还是上次我们提到的例子，我吃到了三天红烧肉，想求出隐藏的后厨序列，下面已经给出计算好的前向变量和后向变量。

那么第一天是大爷的“概率”：
$$ \alpha_{t=1}(1)*\beta_{t=1}(1) = 0.54*0.368=0.19872 $$

第一天是大叔的“概率”：
$$ \alpha_{t=1}(2)*\beta_{t=1}(2) = 0.03*0.487=0.01461 $$
而第一天不可能是大哥，所以我们可以认为第一天最有可能的是大爷。此处的概率一词加了引号是因为只是形式意义上的一种可能性的表示，但是不满足概率的归一化的要求。
以此类推，计算第二天第三天，得到最有可能的后厨序列是：大爷-大叔-大叔。这有什么问题呢？

Read More

写在前面

之前我在研一时写的隐马模型系列一因为不知道啥原因，在CSDN上的文章找不到了，然后就也一直没写，这次又看到了HMM，就争取把那篇文章复现出来。而且因为已经是第二次写了，就争取写的深入一些。

HMM有三个基本问题：

1. 评估问题：给定模型和两个状态序列，求哪个更有可能。
2. 学习问题：给定观测，求解模型参数。
3. 预测问题：给定模型和观测，求隐状态序列。

HMM被定义为一个三元组$$$(A,B,\pi)$$$，$$$\pi$$$表示初始状态概率向量，A表示状态转移概率矩阵，B表示观测概率矩阵。

两个基本假设：

1. 齐次马尔科夫性假设：只依赖于前一状态。
2. 观测独立性假设：观测只依赖于状态。

前向算法与后向算法

前向算法

刷完了DP算法之后，前向算法其实就很简单了。关键是分解该问题的最优子问题，然后写出递推式。
子问题定义为：

第 i 个隐状态是 i ，生成观测序列的前 t 个元素的概率是多少？

并将子问题定义为前向变量$$$\alpha_t(i)$$$,那么递推式就可以表示为：
$$\alpha_{t+1}(i) = \big[\sum_{j=1}^N\alpha_t^j\big]B(i,o_{t+1})$$

最后输出DP缓存中最后一列的元素和即可。这样的话，该算法的时间复杂度是$$$O(N^2T)$$$，而不是大暴搜的$$$O(N^TT)$$$。

后向算法

后向算法的子问题为：

第 i 个隐状态是 i ，生成观测序列的第 t+1 元素直到最后一个元素的概率是多少？

并将子问题定义为$$$\beta_t(i)$$$，那么递推式就可以表示为：
$$\beta_t(i)=\sum_{j=1}^NA(i,j)B(j,o_{t+1})\beta_{t+1}(j)$$

然后递推变成从前往后了，其实并没有改进算法的效率，这个算法只是为了和前向算法结合去解决后面两个问题的。具体怎么结合呢？

比如说求给定观测序列的概率可以用前向概率和后向概率结合表示：
$$P(O|\lambda) = \sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)A(i,j)B(j,o_{t+1})\beta_{t+1}(j)$$

这里引用我自己之前一篇博客中的图片解释

那么我们不使用那两个大大的求和符号的话，就可以指定某个向量两个时刻时内部两个隐状态的情况下，求输出概率了。那如果在对其中一个隐状态求边缘概率呢呢？那就可以指定某个时刻下其内部隐状态的情况下求输出概率了。那用这个概率除以$$$P(O|\lambda)$$$呢？那就会得到指定时刻内部状态是某个固定取值的后验概率。这个就是BaumWelch算法中的一个关键计算步骤。注意上面的计算都包含时间t，如果再对时间t求和，那得到的概率就不限制时刻了，而是某个内部状态或者某两个内部状态出现在隐藏序列中的概率。