VTK学习笔记（十八）vtk GMM

1、VTK 对GMM、EM的实现：

https://github.com/daviddoria/ExpectationMaximization

因为有依赖，下载代码使用：

git clone --recursive https://github.com/daviddoria/ExpectationMaximization.git

然后使用cmake生成项目，需要指定Eigen库路径。

2、windows下编译

2.1、注意事项：需要自己实现linux才有的函数drand48(), srand48();

实现如下：

#ifndef DRAND48_H  
#define DRAND48_H  

#include <stdlib.h>  

#define MNWZ 0x100000000LL  
#define CNWZ 0xB16  
#define ANWZ 0x5DEECE66DLL  

static unsigned long long seed = 1;

double drand48(void)
{
    seed = (ANWZ * seed + CNWZ) & 0xFFFFFFFFFFFFLL;
    unsigned int x = seed >> 16;
    return  ((double)x / (double)MNWZ);

}

void srand48(unsigned int i)
{
    seed = (((long long int)i) << 16) | rand();
}

#endif  

2.2、定义宏 M_PI

#define M_PI 3.14159265358979323846

3、测试程序

#include "ExpectationMaximization.h"
#include "GaussianModel.h"

#include <iostream>
#include <random>
#include <chrono>

Eigen::MatrixXd GenerateData(const unsigned int numberOfSamplesPerGroup, const unsigned int dimensionality);
long int GetSeed(const bool random);

int main(int, char*[])
{
    int dimensionality = 1;

    // Generate some data
    Eigen::MatrixXd data = GenerateData(1000, dimensionality);

    // Initialize the model
    std::vector<Model*> models(2);

    for(unsigned int i = 0; i < models.size(); i++)
    {
      Model* model = new GaussianModel(dimensionality);
      models[i] = model;
    }

    MixtureModel mixtureModel;
    mixtureModel.SetModels(models);

    ExpectationMaximization expectationMaximization;
    expectationMaximization.SetData(data);
    expectationMaximization.SetMixtureModel(mixtureModel);
    expectationMaximization.SetMinChange(1e-5);
    expectationMaximization.SetMaxIterations(100);

    expectationMaximization.Compute();

    std::cout << "Final models:" << std::endl;
    MixtureModel finalModel = expectationMaximization.GetMixtureModel();
    for(unsigned int i = 0; i < finalModel.GetNumberOfModels(); ++i)
    {
      finalModel.GetModel(i)->Print();
    }

  return EXIT_SUCCESS;
}

Eigen::MatrixXd GenerateData(const unsigned int numberOfSamplesPerGroup, const unsigned int dimensionality)
{
  std::default_random_engine generator(GetSeed(true));

  // Mean 0, variance 1
  Eigen::MatrixXd data1(dimensionality, numberOfSamplesPerGroup);
  std::normal_distribution<double> normalDistribution1(0,sqrt(1));
  for(unsigned int i = 0; i < numberOfSamplesPerGroup; i++)
  {
    Eigen::VectorXd v(1);
    v(0) = normalDistribution1(generator);
    data1.col(i) = v;
  }

  // Mean 5, variance 2
  Eigen::MatrixXd data2(dimensionality, numberOfSamplesPerGroup);
  std::normal_distribution<double> normalDistribution2(5,sqrt(2));
  for(unsigned int i = 0; i < numberOfSamplesPerGroup; i++)
  {
    Eigen::VectorXd v(1);
    v(0) = normalDistribution2(generator);
    data2.col(i) = v;
  }

  // Concatentate the matrices horiztonally
  Eigen::MatrixXd data(dimensionality, data1.cols() + data2.cols());
  data << data1, data2;

  return data;
}


long int GetSeed(const bool random)
{
  if(random)
  {
    auto now = std::chrono::system_clock::now();
    return now.time_since_epoch().count();
  }
  else
  {
    return 0;
  }
}

测试代码的作用是生成一组测试数据，均值0，标准差1 数据1000组， 均值5，标准差2数据1000组，拼接到一起。
执行计算后输出结果如下：

ExpectationMaximization::Compute
KMeans finished in 6 iterations.
EM iteration 0...
EM iteration 1...
EM iteration 2...
EM iteration 3...
EM iteration 4...
EM iteration 5...
EM iteration 6...
EM iteration 7...
EM iteration 8...
EM iteration 9...
EM iteration 10...
EM iteration 11...
EM iteration 12...
EM iteration 13...
EM iteration 14...
EM iteration 15...
EM iteration 16...
EM iteration 17...
EM iteration 18...
EM iteration 19...
EM performed 20 iterations.
Final models:
Mean: -0.00637246
Variance: 1.05664
Mean: 5.10241
Variance: 1.87834

可见6次迭代计算出均值，20次迭代计算出方差。

输出还缺一组参数 ：权重。

修改原代码，增加属性 维度和权重系数的输出

void GaussianModel::Print() const
{
    std::cout << "Mean: " << this->Mean << std::endl;
    std::cout << "Variance: " << this->Variance << std::endl;
    std::cout << "Dimensionality: " << this->Dimensionality << std::endl;
    std::cout << "MixingCoefficient: " << this->MixingCoefficient << std::endl;
}

再次执行测试程序输出；

ExpectationMaximization::Compute
KMeans finished in 6 iterations.
EM iteration 0...
EM iteration 1...
EM iteration 2...
EM iteration 3...
EM iteration 4...
EM iteration 5...
EM iteration 6...
EM iteration 7...
EM iteration 8...
EM iteration 9...
EM iteration 10...
EM iteration 11...
EM iteration 12...
EM iteration 13...
EM iteration 14...
EM iteration 15...
EM iteration 16...
EM iteration 17...
EM iteration 18...
EM iteration 19...
EM iteration 20...
EM iteration 21...
EM iteration 22...
EM iteration 23...
EM performed 24 iterations.
Final models:
Mean: 5.00316
Variance: 2.00875
Dimensionality: 1
MixingCoefficient: 0.491851
Mean: -0.0123765
Variance: 0.995373
Dimensionality: 1
MixingCoefficient: 0.508149

本来想说是基于VTK实现的GMM算法的，测试到后来发现没有用到VTK。
如果一定要跟VTK扯上关系的话，Eigen库可以直接从VTK库中调用，而不用额外链接。

4、其他

关于drand48() and srand48(long) on windows 还找到一个实现没有测试。

参考：mortennobel/drand48.c

#define RAND48_SEED_0   (0x330e)
#define RAND48_SEED_1 (0xabcd)
#define RAND48_SEED_2 (0x1234)
#define RAND48_MULT_0 (0xe66d)
#define RAND48_MULT_1 (0xdeec)
#define RAND48_MULT_2 (0x0005)
#define RAND48_ADD (0x000b)

unsigned short _rand48_seed[3] = {
        RAND48_SEED_0,
         RAND48_SEED_1,
         RAND48_SEED_2
};
unsigned short _rand48_mult[3] = {
         RAND48_MULT_0,
         RAND48_MULT_1,
         RAND48_MULT_2
 };
unsigned short _rand48_add = RAND48_ADD;

void
 _dorand48(unsigned short xseed[3])
 {
	         unsigned long accu;
	         unsigned short temp[2];
	
	         accu = (unsigned long)_rand48_mult[0] * (unsigned long)xseed[0] +
	          (unsigned long)_rand48_add;
	         temp[0] = (unsigned short)accu;        /* lower 16 bits */
	         accu >>= sizeof(unsigned short)* 8;
	         accu += (unsigned long)_rand48_mult[0] * (unsigned long)xseed[1] +
	          (unsigned long)_rand48_mult[1] * (unsigned long)xseed[0];
	         temp[1] = (unsigned short)accu;        /* middle 16 bits */
	         accu >>= sizeof(unsigned short)* 8;
	         accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0];
	         xseed[0] = temp[0];
	         xseed[1] = temp[1];
	         xseed[2] = (unsigned short)accu;
}

double erand48(unsigned short xseed[3])
{
         _dorand48(xseed);
         return ldexp((double) xseed[0], -48) +
                ldexp((double) xseed[1], -32) +
                ldexp((double) xseed[2], -16);
}

double drand48(){
	return erand48(_rand48_seed);
}

void srand48(long seed){
	_rand48_seed[0] = RAND48_SEED_0;
	_rand48_seed[1] = (unsigned short)seed;
	_rand48_seed[2] = (unsigned short)(seed >> 16);
	_rand48_mult[0] = RAND48_MULT_0;
	_rand48_mult[1] = RAND48_MULT_1;
	_rand48_mult[2] = RAND48_MULT_2;
	_rand48_add = RAND48_ADD;
}

可能遇到的问题：

include cmath to fix ldexp declaration.

另外关于实现的理论

参考：daviddoria/ExpectationMaximization