原地址：http://www.cnblogs.com/freedomshe/archive/2012/04/24/2468747.html
题记：2012年4月1日回到家，南大计算机研究僧复试以后，等待着的就是独坐家中无聊的潇洒。不知哪日，无意中和未来的同学潘潘聊到了图像处理，聊到了她的论文《基于LDA的行人检测》，出于有一年半工作经验的IT男人的本能，就一起开始学习研究这篇“论文”了。众所周知，老师给学生设置论文题目的，起初都是很模糊的——自己没有思考清楚实践上的可行性和具体思路，仅从理论了解上就给学生设置一些“难以实现”的论文任务。几经修改和商讨，最后的论文实际上就是“基于SIFT+Kmeans+LDA的图片分类器的实现”了。至此，代码已经编写完毕，图片分类的效果还算满意。 ——copyright：由于是一起学习研究的结果，相关所有内容潘潘童鞋可以以第一作者身份使用！ . 实现思路
. 软件环境
. Step1——SIFT应用
. Step2——Kmeans应用
. Step3——数词频的实现
. Step4——LDA应用
. 参考 一、实现思路
分类器的功能是：输入一组图片，给定需要分类的类别数lda_k(>)；输出lda_k个文件夹，每个文件夹内的图片为一类图片。 第一步是SIFT特征提取：输入图片，输出图片的特征点集，即feature列表，每个feature代表一个图片的某个局部特征，每个feature的数据结构由一个128维浮点数组表示。至此，可以将一幅图片转换成一个feature集。 第二步是Kmeans聚类：输入是所有图片的feature集的综合，给定参数km_k代表需要聚类的类别数；输出是km_k个feature，在LDA的视角看来就是“单词表”，用“单词表”中的一个“单词”（类中的质心feature）代表kmeans聚类里面一类的所有feature。 第三步是统计词频：（对每个图片）输入是图片的feature集和“单词表”，分别计算该图片feature集中每个feature对应的“单词”，并统计每个“单词”在该feature集中出现的次数即词频；输出是词频统计数据。 最后一步是LDA训练潜在主题：输入是所有图片文件的词频统计数据，以及给定的需要训练出来的主题类别数lda_k；LDA输出参数较多，其中最有用的就是文档-主题条件概率矩阵（theta矩阵），即举证中每个元素表示P(主题k|文档m)——在文档m中，主题是k的概率——通过该概率即可判断当前文档最可能的主题，实现了将所有文档分类为lda_k个主题。 总之，理论上LDA研究的实体是一组文档，每个文档由若干单词组成，通过无监督学习，能够发现lda_k个主题，并且确定theta矩阵——文档确定的情况下生成主题k的概率，以及phi矩阵——主题确定的情况下生成单词v的概率。分类器通过SIFT算法将图片转换为若干feature，即将图片看成是“文档”feature看成是“单词”。而仅通过SIFT处理后的feature并不能直接单做“单词”作为LDA的输入，因为几乎每个feature都不一样，还需要Kmeans算法对所有图片的feature集的总和做一次聚类，得到km_k个类别的中心feature，即生成km_k个“单词”的“单词表”，并以此中心feature代替一个类别内的所有其他feature，从而将一个图片“文档”中的所有feature均在“单词表”中能够找到代表它的“单词”，这样图片就真正转换为了LDA能够处理的“文档”。 [Go Top] 二、软件环境
VS2010，MFC，C++。
安装并配置Opencv，参见VS2010+Opencv-2.4.0的配置攻略。
下载并集成SIFT源码，参见在VS2010中应用SIFT（C）源码。
下载并集成LDA源码，参见在VS2010中应用LDA（C）源码。 Kmeans为Opencv自带函数，无需应用其他源码。建立好自己的工程，集成算法源码后，工程文件夹大致结构应为下图所示： 设计好自己的例程界面，并关联好响应函数和成员变量，本例程界面如下： [Go Top] 三、Step1——SIFT应用
在该步骤内，程序依据“图片源目录”给出的图片目录路径，扫描目录内的所有图片文件，对每个执行如下操作： ...
n = _sift_features(img, &features, SIFT_INTVLS/**/, SIFT_SIGMA/*1.6*/, SIFT_CONTR_THR/*0.04*/,
SIFT_CURV_THR/**/, SIFT_IMG_DBL/**/, SIFT_DESCR_WIDTH/**/, SIFT_DESCR_HIST_BINS/**/); //SIFTfeature提取
...
export_features(out_file_name, features, n); //将features导出为文件
...
if(勾选了“保存SIFT特征图”)
{
draw_features(img, features, n); //在img图片上标记出features
cvSaveImage(out_img_name, img, NULL); //将标记后的图片保存
}
... 其中最主要的三个函数就是_sift_features(…), export_features(…), draw_features(…)均为sift源码所提供。（注：feature有两种类型——OXFD和LOWE，本程序只涉及LOWE类型，所有OXFD相关格式均自动忽略。） _sift_features(…)函数第一个参数img为传入图片的IplImage指针格式，为Opencv所定义的图片数据结构；features后面的参数均为SIFT算法的输入参数，具体含义见作者的源码注释。 需要注意和理解的是features这个参数，其指向的为一个结构体feature的数组，feature结构为： /**
Structure to represent an affine invariant image feature. The fields
x, y, a, b, c represent the affine region around the feature:
a(x-u)(x-u) + 2b(x-u)(y-v) + c(y-v)(y-v) = 1
*/
struct feature
{
double x; /**< x coord */
double y; /**< y coord */
double a; /**< Oxford-type affine region parameter */
double b; /**< Oxford-type affine region parameter */
double c; /**< Oxford-type affine region parameter */
double scl; /**< scale of a Lowe-style feature */
double ori; /**< orientation of a Lowe-style feature */
int d; /**< descriptor length */
double descr[FEATURE_MAX_D]; /**< descriptor */
int type; /**< feature type, OXFD or LOWE */
int category; /**< all-purpose feature category */
struct feature* fwd_match; /**< matching feature from forward image */
struct feature* bck_match; /**< matching feature from backmward image */
struct feature* mdl_match; /**< matching feature from model */
CvPoint2D64f img_pt; /**< location in image */
CvPoint2D64f mdl_pt; /**< location in model */
void* feature_data; /**< user-definable data */
}; 其中x,y表示feature在图片中的坐标，scl, ori表示在图中标记特征的强度和方向，descr是最重要的特征信息即128维的特征向量（FEATURE_MAX_D==）。实际上所保存的特征文件里面也只保存了feature的这些信息。features就是feature的一个数组，包含了SIFT算法所提取的图片的所有feature，返回值n表示features数组中有多少个feature元素。 export_features(…)就是将上一步中的features保存为文件，文件格式如下： 178.616459 111.621902 34.241822 1.292619 220.732993 432.886679 20.895674 -1.740551 …… 其中第一行的第一个int表示文件中feature的个数，第二个int表示feature的维数。接下来每段132个浮点数（后128个近似为整数）为一个feautre，前4个分别为x, y, scr, ori，后128个为表征feature的128维向量。 draw_features(…)就是将feature信息示例性地标注在图上了，图片如下： cvSaveImage(…)是Opencv的函数，将加过标注的图片保存为文件。 这一步实际上要做的重点就是输入为图片，输出为SIFT算法后为每张图片生成的feature文件集： 这样每张图片就相当于一个“文档”，而feature就是“文档”中的“预备单词”。 [Go Top] 四、Step2——Kmeans应用
Step1里面的feature只是“预备单词”，在成为单词之前还要通过Step2生成“单词表”和Step3将“文档”中的“预备单词”找到“单词表”中最相近的“单词”替换之（并不是真正操作上的替换，只是当成“单词表”中的“单词”统计出来而已）。 在Step2中，关键操作如下： ...
CvMat *samples=cvCreateMat(featureNum, dims, CV_32FC1); //包含所有图片的所有feature信息的矩阵，featureNum个feature，每个feature为dims（128）维向量，每一维的元素类型为32位浮点数
CvMat *clusters=cvCreateMat(featureNum, , CV_32SC1); //每个feature所在“质心”的指针（实际上本例程中没有用到该信息）
CvMat *centers=cvCreateMat(k, dims, CV_32FC1); //“质心”信息的数组，k个“质心”每个质心都是dims（128）维向量，每一维的元素类型为32位浮点数
cvSetZero(clusters); //将矩阵初始化为0
cvSetZero(centers); //将矩阵初始化为0
while(file.ReadString(strLine))
{
...
n = import_features(CIni::CStrToChar(fileName), FEATURE_LOWE, &features); //导入feature文件，n为导入的feature个数
...
//将feature文件内所有feature信息存入samples矩阵结构内
for(int i = ; i < n; i++)
{
for(int j = ; j < dims; j++)
{
samples->data.fl[temp++] = features[i].descr[j];
}
}
}
cvKMeans2(samples, k, clusters,cvTermCriteria(CV_TERMCRIT_EPS,,1.0), , (CvRNG *), KMEANS_USE_INITIAL_LABELS, centers); //Kmeans聚类
...
cvSave(CIni::CStrToChar(ini.getWordListFilePath()), centers); //保存单词表
... 其中关键函数当然是import_features(...)和cvKMeans2(...)，前者是sift源码里的方法，用来导入feature文件使之成为内存数据结构，后者是Opencv里的kmeans算法之一（cvKMeans2(...)内部调用了kmeans(...)）。 在说明他们之前首先要了解Opencv内部通用数据结构——用来存储矩阵类型的CvMat，其声明如下： typedef struct CvMat
{
int type;
int step; /* for internal use only */
int* refcount;
int hdr_refcount; union
{
uchar* ptr;
short* s;
int* i;
float* fl;
double* db;
} data; #ifdef __cplusplus
union
{
int rows;
int height;
}; union
{
int cols;
int width;
};
#else
int rows;
int cols;
#endif
} 其中最重要的就是用来存储元素信息的一维数组data，data为一个联合体，包含了各种通用类型。data虽然是一维结构，但通过rows, cols标识，可以将data看成二维的矩阵结构，type就标识了元素类型，而step就标识了一行的步长（占用字节数）。
cvCreateMat(int rows, int cols, int type);包含了三个参数，rows和cols分别表示data二维矩阵的行数和列数，type表示元素的类型。元素类型通过宏来定义，CV_32FC1中的32表示元素的位数，F表示是浮点数，Cn表示通道数，C1表示1通道，C3表示3通道（只有在类似RGB的三元表示一个像素点颜色时才用到多通道，而通道实际上就是说一个元素应该由几列来表示）。
有了上面的知识，就知道了samples矩阵是featureNum行，128列的矩阵，而矩阵元素类型为32位浮点数。featureNum的值为所有图片所有feature的总和，也就是说把所有的feature文件里面的feature信息都读入到samples矩阵内用来作为Kmeans的输入。clusters矩阵是featureNum行1列的矩阵，实际上就是featureNum个int元素的数组，每个int元素都相当于一个下标，标识对应feature（一行）属于centers中的哪个类（“质心”），而本程序里面没有使用到clusters提供的信息，featuresNum个feature分别对应于哪个“质心”是在Step3算词频里面实现的。centers就是保存“单词表”信息的矩阵了，k行128列的32位浮点矩阵，即k个feature，每个feature都是一类feature的“质心”——该类所有feature中最中间那个（通过欧式距离来计算的）feature。 关键操作来了，import_features(...)就是将feature文件读入内存，使之成为feature的内存数据结构。
int import_features(char* filename, int type, struct feature** feat)第一个参数就是feature文件的路径，第二个参数是feature的类型，本文是FEATURE_LOWE类型，第三个参数就是传出参数feature数组了，返回值为读入内存的feature个数。在接下来的操作中就需要把features里面的数据全部拷贝到samples矩阵结构内，这样所有文件的feature信息都拷贝到samples内后就完成了kmeans传入数据的工作。 cvKMeans2(...)是Opencv在kmeans(...)之上封装的一个函数，参数含义可以参照Hongquan的博客中的《OpenCV中kmean算法的实现》，本文需要说明的几个参数是：第一个参数samples——传入所有图片feature信息的总和；第二个参数k即kmeans的聚类个数（本程序通过界面的类别数设置）；第七个参数flags表示生成随机数的方式，可能每次运行程序对相同的输入输出的单词表都不同，那么就跟这个参数相关；第八个参数centers，这就是我们需要的输出了，即“单词表”信息，每一行就是一个feature也就是一个“单词”。 这样，读入Step1中输出的所有feature文件，然后综合起来并输入参数k（聚类类别数），然后就可以生成一张包含k个“单词”的“单词表”了，用Opencv的cvSave(...)函数把这张“单词表”保存为yml文件，就是本步骤的主要输出了。输出“单词表”文件wordList.yml格式类似： %YAML:1.0
wordList: !!opencv-matrix
rows:
cols:
dt: f
data: [ 1.64593754e+001, 2.11062511e+001, .,
3.48312492e+001, 2.39375000e+001, 1.28156252e+001,
1.25531254e+001, 1.16812506e+001, 4.92906265e+001,
3.61625023e+001, 2.94968758e+001, ., 1.70375004e+001,
1.12718754e+001, 1.71968746e+001, 3.13437500e+001,
3.27593765e+001, 2.25156250e+001, 2.02625008e+001,
2.26406250e+001, 2.73750000e+001, 1.90531254e+001,
2.41375008e+001, 2.85156250e+001, 1.92875004e+001,
............................................. ] [Go Top] 五、Step3——数词频的实现
前面已经说过，Step2是生成“单词表”，而Step3就是通过计算欧式距离来确定原始feature文档中每个feature对应“单词表”中的哪个“单词”，然后统计出来，生成对应的统计文件——词频文件wordF.data。主要操作如下： ...
FileStorage readfs(CIni::CStrToChar(ini.getWordListFilePath()), FileStorage::READ); //以只读的形式打开yml。
Mat wordList; //单词表矩阵
readfs[CIni::CStrToChar(CIni::removeSuffix(CIni::getFileNameFromPath(ini.getWordListFilePath())))] >> wordList; //读取单词表
...
while(“对于每个feature文档”)
{
...
int n = import_features(CIni::CStrToChar(fileName), FEATURE_LOWE, &features);
...
for(int i = ; i < n; i++)
{
...
distMin = normL2Sqr_(pa, pb, dims); //计算欧式距离
...
for(int j = ; j < wordNum; j++)
{
...
dist = normL2Sqr_(pa, pb, dims); //计算欧式距离
...
}
}
}
... 显而易见，前三个操作就是利用Opencv的函数，将Step2中生成的“单词表”读入矩阵结构（Mat）的变量wordList中。而对于每个feature文档，通过sift源码的import_features(...)函数读入以feature结构体为元素的数组features中。接下来对于每个feature计算它与wordList中哪个“单词”最接近，最后统计出该文档中包含的单词及其个数。 normL2Sqr_(pa, pb, dims)是Opencv中计算欧式距离的函数，其中pa指向一个feature，而pb指向单词表中的一个“单词”（“单词”也是feature），dims表示待计算数据的维数，在这里就是128了。 经过上面的操作，统计好每个图片文档中每个“单词”出现的次数后，统一保存为一个词频文档，即将所有图片文档的词频信息保存在一个文件wordF.data内，其文件结构为： : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 其中每一行代表一个文档（图片文档），%d:%d的结构表示单词ID:单词个数，第一个数字表示后面的元素项个数。这些数据就是Step3统计词频的输出了，也是Step4LDA运算的输入。 [Go Top] 六、Step4——LDA应用
作者对LDA的实现并不是像sift和kmeans那样由一个函数通过参数传入传出来给出，而是作者的main函数的一个实现过程。本例程LDA的主要过程为： ...
int topic_num = lda_k; //LDA分类数
struct corpus *cps;
struct est_param param;
...
cps = read_corpus(data); //读取训练集
init_param(cps,¶m,topic_num); //初始化参数
//迭代计算
while ()
{
//对每个文档使用sampling方法计算
for (int m=; m<cps->num_docs; m++)
{
...
for (int l=; l<cps->docs[m].length; l++)
{
for (int c=; c<cps->docs[m].words[l].count; c++)
{
param.z[m][word_index] = sampling(m,word_index,cps->docs[m].words[l].id,topic_num,cps,¶m,alpha,beta,p,s_talpha,vbeta); //sampling计算
...
}
}
}
if ((iter_time >= burn_in_num) && (iter_time % SAMPLE_LAG == ))
{
calcu_param(¶m, cps,topic_num,alpha,beta); //计算参数
...
}
//迭代结束条件
if (sample_time ==sample_num)
{
break;
}
}
average_param(¶m, cps,topic_num,alpha,beta,sample_num); //计算theta，phi平均值
save_model(cps,¶m,model_name,alpha,beta,topic_num,sample_num); //保存结果数据
... 对于没有接触过LDA的人要看懂函数过程还是很困难的，但对于有语言功底的程序员，如果只是应用LDA过程，只需要了解几个主要的LDA概念就行了。（其他部分请参考作者代码注释） 、LDA的输入。LDA的输入除了Step3生成的词频统计信息外，还需要一些参数。这些参数一般设置为通用参数就行了，本例程只有lda_k为界面输入，表示LDA训练的主题个数。
read_corpus(data)方法中，data为字符串，表示输入数据的路径，该方法将文件读取为struct corpus的格式。corpus结构体表示的是一个文档集，其结构可以从声明中看出： struct word
{
int id;
int count;
}; struct document
{
int id; //文档id
int num_term; //文档包含的单词个数（count的总和）
int length; //文档包含的单词类别个数（id:count结构的个数）
struct word* words;
}; struct corpus
{
struct document* docs;
int num_docs; //文档个数
int num_terms; //单词表中单词总数（实际上是所有单词中最大id+1）
}; 对比Step3中生成的词频文件，不难看懂上面这些结构体的意义。 、LDA的输出。LDA的输出包含很多数据，由save_model(...)函数输出为文件： lda.other文件保存参数alpha, beta, topic_num, sample_num的值。lda.topic_assgin文件保存z矩阵，z[m][n]==k表示文档m中的单词n所对应的主题为k。lda.theta文件保存theta矩阵，theta[m][k]表示在文档为m时，生成主题k的概率，即条件概率p(主题k|文档m)。lda.phi文件保存phi矩阵，phi[k][v]表示在主题为k时，生成单词v的概率，即条件概率p(单词v|主题k)。 其中需要关注的只有theta和phi两个矩阵，而本例程只用到了theta矩阵的信息。theta矩阵和phi矩阵的例子如下表： P(Topic_k|Doc_m) Topic1 Topic2 Doc1 0.45 0.55 Doc2 0.1 0.9 Doc3 0.6 0.4 P(Word_v|Topic_k) Word1 Word2 Word3 Word4 Topic1 0.45 0.15 0.2 0.2 Topic2 0.1 0.5 0.3 0.1 theta矩阵 phi矩阵
其中的Topic都是由LDA通过无监督学习得到的潜在主题，只需要用户告诉LDA主题数目就行了。从上面的表格中不难看出，在theta矩阵中，我们知道了每个文档生成每个主题的概率，通过比较概率大小就可以确定文档所对应的主题了。 、LDA算法及其运算所需结构体struct est_param。 struct est_param
{
int **z; //z[m][n] stands for topic assigned to nth word in mth document
double **theta; // theta[m][k] stands for the topic mixture proportion for document m
double **phi; // phi[k][v] stands for the probability of vth word in vocabulary is assigned to topic k
// count statistics
int **nd; //nd[m][k] stands for the number of words assigned to kth topic in mth document
int **nw; //nw[k][t] stands for the number of kth topic assigned to tth term
int *nd_sum; //nd_sum[m] total number of word in mth document
int *nw_sum; //nw_sum[k] total number of terms assigned to kth topic
}; 该结构体内theta和phi二维数组是我们熟悉的输出，z数组也是前面提到的输出之一，后面的四个变量是LDA运算中间过程的必备临时变量。LDA算法，如果了解了相关的变量名和LDA过程，通过下面这幅图是不难了解算法过程的： 综上，知道了每个图片文档对应哪个主题的概率最大后，就可以根据主题个数新建lda_k个文件夹，然后把分类到对应主题的图片拷贝过去，从而实现了对图片的LDA分类。分类效果如图： [Go Top] 七、参考
【Opencv】
Opencv下载：http://sourceforge.net/projects/opencvlibrary/files/
Opencv教程：http://www.opencv.org.cn/index.php 【SIFT】
SIFT源码下载：http://blogs.oregonstate.edu/hess/code/sift/
SIFT源码使用方法：http://www.open-open.com/lib/view/1325332699514 【Kmeans】Kmeans为Opencv自带函数。
OpenCV中kmean算法的实现：http://blog.hongquan.me/?p=8 【LDA】
LDA源码下载：http://code.google.com/p/lsa-lda/
gibbs的LDA实现：http://sourceforge.net/projects/gibbslda/ 【Bag-of-words模型】
SIFT算法的应用--目标识别之Bag-of-words模型：http://blog.csdn.net/v_JULY_v/article/details/6555899 [Go Top] 源码下载：http://www.pudn.com/downloads464/sourcecode/graph/texture_mapping/detail1949366.html