GMM的EM算法实现

GMM的全名是Gaussian Mixture Model,即高斯混合模型,这个名字挺直观。

如果给了一组数据点,假设这些数据点符合iid,那么很容就能根据数据拟合出一个最佳的高斯分布。

但是如果给了一堆数据点,而他们却不符合iid,只能做到互相独立,并不能做到同分布。相反的,数据点们有可能来自不同的高斯分布,这也就是『高斯混合』模型名字的意思。这样如何能拟合出最佳的高斯分布们呢?

EM算法提供了一种方法,这种方法的思路有些像牛拉法。先是确定一个初值,然后看看当前值距离目标还有多远,然后逐渐的缩小差距(不断修正membership weight)。

具体原理见这里

这篇小短文大概记录下这次作业的内容。作业题目是生成了2个二维高斯分布混合的数据点,共2000个,然后用这两千个点拟合出2个二维高斯分布的参数。

clear;
%产生2个二维正态数据
MU1    = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2    = [-1 -1];
SIGMA2 = [1 0; 0 1];
X      = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'.');

%initial values
h_mu1 = X(randi(2000),:);
h_sigma1 = cov(X);
h_mu2 = X(randi(2000),:);
h_sigma2 = cov(X);
h_alpha1 = 0.5;
h_alpha2 = 0.5;

%simple iteration
for iter=1:100
    %E step
    w1 = mvnpdf(X,h_mu1,h_sigma1)*h_alpha1./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);
    w2 = mvnpdf(X,h_mu2,h_sigma2)*h_alpha2./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);

    %M step
    h_alpha1 = sum(w1)/2000;
    h_alpha2 = sum(w2)/2000;
    
    h_mu1 = sum([w1.*X(:,1) w1.*X(:,2)])/sum(w1);
    h_mu2 = sum([w2.*X(:,2) w2.*X(:,2)])/sum(w2);

    M1 = X-ones(2000,1)*h_mu1;
    h_sigma1 = ([sum(w1.*M1(:,1).^2) 0;0 sum(w1.*M1(:,2).^2)]+...
        sum(w1.*(M1(:,1).*M1(:,2)))*[0 1;1 0])/sum(w1);
    M2 = X-ones(2000,1)*h_mu2;
    h_sigma2 = ([sum(w2.*M2(:,1).^2) 0;0 sum(w2.*M2(:,2).^2)]+...
        sum(w2.*(M2(:,1).*M2(:,2)))*[0 1;1 0])/sum(w2);
end

figure;
[x,y] = meshgrid(linspace(-5,5,75));
p1 = mvnpdf([x(:),y(:)],h_mu1,h_sigma1);
p2 = mvnpdf([x(:),y(:)],h_mu2,h_sigma2);
surf(x,y,reshape(p1,75,75));hold on
surf(x,y,reshape(p2,75,75));

options = statset('Display','final');
obj = gmdistribution.fit(X,2,'Options',options);
figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);

绝大多数时候都没有什么问题,不过由于初始点是从X里随便选的一个数据点,极偶然会跑偏。

以下分别是数据点的分布、自行拟合出的结果、利用matlab自带拟合函数拟合出的结果。

originalgmm

根据给定的均值、协方差矩阵画出的GMM的图

xscatter

数据点的散点图,保存之后可能有点糊

fitbymyself

自己写的脚本拟合出的GMM的图

fitbymatlabtoolbox

matlab提供的拟合函数拟合出的图

迭代100次之后,高斯分布1的均值分别偏差0.2776%、4.4462%,协方差分别偏差了2.4458%、5.0146%;高斯分布2的均值偏差了2.6981%、1.2911%,协方差偏差了3.8498%、1.5517%。精度凑凑活活把...

对比之下,matlab自带的拟合函数拟合出来的高斯分布1的均值偏差2.2415%、2.2415%,协方差偏差为4.02%、1.32%;高斯分布2的均值偏差了0.4165%、0.4165%,协方差偏差为2.39%、4.62%。精度也没高到哪里去...不过均值偏差一样让我很是在意...

然而看图1、图3和图4的话,明显还是自己写的代码拟合效果更好一些啊,有些不科学。

逐帧提取《Bad Apple》至txt文件

最近实验室要搞一块液晶屏,决定也跳到俗坑里去,在液晶屏上跑一次《Bad Apple》。

不过鉴于液晶屏用哪家还没决定,所以先把前期的准备工作做好。这篇短文中将介绍下如何将《Bad Apple》中的每帧的像素点存到txt文件中。

继续按照着将大象塞到冰箱里需要三步的思路做这件事。将《Bad Apple》逐帧将像素点存入txt中也只需要三步。

1. 上网找到一个1440x1080的《Bad Apple》的视频。

当然分辨率低一些也无所谓啦,反正液晶屏那么渣。假设下载到的《Bad Apple》的视频文件名叫做“BA.mp4”。

2. 将该视频逐帧截图,并将截图按照顺序保存。

这一步可以通过两个方法实现。

2.1  第一个方法是简单的方法:用KMPlayer捕捉。

首先去下载软件。然后安装。

然后用KMPlayer打开“BA.mp4”,播放窗口右击会看到叫做“Capture”的菜单,选择"CaptureFrame Extract"(中文版本中为"捕捉画面:高级捕捉"),这个时候会打开一个窗口进行设置。

  • “Extract to”(捕获到)中设置图片输出的路径;
  • “Image format”(图像格式)中选择是输出bmp(位图)、jpeg还是png(我用的jpeg)
  • “Numbers to extract”(要捕获的数量)选择“continuously”(连续)
  • “Size to extract”(捕获尺寸)中选择“Original size”(原始尺寸)或是“Specified size”(自定义大小,注意长宽比和原视频一致。)
  • “Frames to extract”(要捕获的帧)中选择截取画面的方式,可以选择“Every frame”(所有帧),不过鉴于缩小图包体积的考虑,我选择了“In 1sec. 30 frame(s)”(在一秒内30帧)。

一切设置妥当之后一定要记得点击播放,否则是不会截图的!将视频从头到尾播放一遍,会自动的在指定的路径输出一堆图片。

比如在我的设置中就是在“E:BA”这个路径下输出了“BA.mp40000.jpg”——“BA.mp43570.jpg”,共3571张图。其中文件名中的“BA.mp4”为视频文件名,之后的“0000.jpg”——“3570.jpg”为图片的编号。

可以在cmd中用ren命令和通配符批量把文件名前缀“BA.mp4”修改掉,不过鉴于之后还要使用matlab进行下一步处理,我觉得就没有这个必要了。

2.2  接下来介绍第二步的第二种方法,这种方法比较麻烦:用opencv读取视频文件,然后逐帧输出。

鉴于opencv的开发环境配置比较繁琐,且本人在两年前使用了半年的opencv+VS2010之后就再也没有接触过这个东西,久远的记忆几乎没怎么剩下,只剩下了本人硬盘上的一堆代码。所以opencv的开发环境配置过程在此不表,出门google一下,具体过程叙述得很详细。

下面说一下opencv中截取帧的大概思路。

int pic_index = 0;  //图片编号
CvCapture* capture;
IplImage* frame;
char picname[] = "";
char filename[]="E:BA.avi";  //因为不知道mp4文件是不是能打开,记得avi是可以打开的
capture = cvCaptureFromFile(filename);
if(!capture)
	return 1;
while(1)
{
	frame = cvQueryFrame(capture);
	if(!frame)
		break;
        /*picname[] = ???; 这段懒得想了...一想类型转换就头大...*/
	cvSaveImage(picname,frame);
        pic_index++;
}
cvReleaseCapture(&capture);

最后应当也能在picname对应的路径下生成像对应的图片,不过鉴于这个方法本人并没有尝试,且opencv的开发环境配置起来实在麻烦,所以还是比较推荐用第一个比较傻瓜的方法。

3. 第三步就是将所得到的图片二值化,以在液晶中显示。

这一步按说也有两种方法,一种方法是用matlab;另一种方法还是用opencv。

3.1  先介绍用matlab的方法。

刚才提到输出了“BA.mp40000.jpg”——“BA.mp43570.jpg”多个文件。下面这段代码就是在matlab中处理这3000多张图,进行二值化的。

for i=0:3570  %图片编号
    str2 = num2str(i);  %将类型转化为str
    switch length(str2)  %判断图片编号的位数,然后再前面补0,以符合“BA.mp4xxxx.jpg”的命名规则
        case 1
            str2 = strcat('000',str2);  %一位数补仨0,以此类推
        case 2
            str2 = strcat('00',str2);
        case 3
            str2 = strcat('0',str2);
    end
I = imread(strcat('E:BABA.mp4',str2,'.jpg'));  %读取制定路径下的图片

t = graythresh(I);  %用Otsu's method(我也不知道这是啥)找到一个合适的阈值以进行二值化
I_2 = im2bw(I,t);  %进行二值化,这个时候I_2是一个存着像素点数据的矩阵,每个点0为黑,1为白。

fid = fopen(strcat('E:BA2txt',str2,'.txt'),'wt');  %将I_2的数据写到txt文件中,并且按照分辨率进行换行
for i=1:480  %此处是因为我按帧截图的分辨率是640*480,所以矩阵一共有480行
    fprintf(fid,'%d',I_2(i,:));
    if i~=480  %不加以末行判断的话最后会有一行空行
        fprintf(fid,'n');
    end
end
fclose(fid);

imwrite(I_2,strcat('E:BA2',str2,'.jpg'));  %将二值化之后的图像输出至对应路径
end

之后就会在对应路径下生成一系列记录着每个二值图片每个像素点的信息的txt文件,和一系列二值图片。

3.2  然后大概介绍下opencv中的方法。

如果第二步用了opencv,第三步还是可以继续用opencv的,这样不仅能够把开发环境配置的代价分摊到两步里,还可以使用opencv里的cvSmooth函数对图像进行平滑处理,能减少不少毛刺。不过鉴于cvSmooth函数并不属于本次的内容,所以也就在此一笔带过好了。

首先还是要读取原图片,然后大概的思路其实和matlab里的也差不多。

IplImage* blkWhtImage1;  //二值图
IplImage* grayImage = cvCreateImage(cvSize(640,480),IPL_DEPTH_8U,1);  //灰度图
cvCvtColor(srcImage,grayImage,CV_BGR2GRAY);
cvThreshold(grayImage,blkWhtImage1,0,255,CV_THRESH_BINARY|CV_THRESH_OTSU);  //转换成二值图

具体也可以参照官网文档中的示例代码进行二值化。

接下来介绍读取像素点的过程。

int color;
for(int row=0; row!=480; row++)
{
	for(int col=0; i!=640; col++)
	{
		color = cvGet2D(blkWhtImage1,row,col).val[0];  //由于blkWhtImage是一个单通道图像,所以数据存在cvScalar.val[0]中
                /*对color进行处理,按照行数和列数存到文件中即可。fopen什么的文件操作很烦,懒得看了...*/
	}
}

到此图片的像素点的信息就存到txt中去了。

结尾:当然在液晶上跑坏苹果这只是其中很很简单的一步,具体的硬件的设计以及通讯如何进行则是更麻烦的事情。

如果显示的载体是示波器则更麻烦,需要进行边界识别。在matlab中需要调用bwperim函数,在opencv中则需要调用cvFindContours函数。不过鉴于本次是在液晶上跑,就不再赘述了。