statismo代码分析--statismo-build-gp-model | HyperPlane

statismo代码分析--statismo-build-gp-model

发表于 | 更新于 | 阅读次数：

使用高斯过程描述形变场，输入的数据只需要参考帧就行。需要设定kernel，由参考帧和kernel生成model

命令行参数

  • -t, --type TYPE : 指定模型是shape还是deformation
  • -d, --dimensionality : 输入数据的维度，当TYPE为deformation时有效
  • -k, --kernel KERNEL : 指定高斯过程的kernel，默认只支持gaussian
  • -p, --parameters KERNEL_PARAMETERS : kernel的参数，gaussian只有一个sigma参数
  • -s, --scale SCALE : 作用在kernel上的scaling参数
  • -n, --numberofbasisfunctions NUMBER_OF_BASIS_FUNCTIONS : 这个model需要的参数个数
  • -r, --reference REFERENCE_FILE : reference文件路径
  • -o, --output-file OUTPUT_FILE : 输出文件路径
  • -m, --input-model MODEL_FILE : 已有model的文件路径，用于扩展当前已有的model

kernel的代码

在文件 utils/statismo-build-gp-model-kernels.h 中添加kernel，在这里kernel是对应的一对点计算得到的一个值，并非对应点的每个元素得到的每个一个对应向量

截取guassian kernel代码

在对应的头文件中定义类似的模板类，内联重载运算符()设定具体的kernel函数

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
// 用法，opt.vKernelParameters指定kernel类型
typedef typename DataType::PointType PointType;
typedef boost::scoped_ptr<const statismo::ScalarValuedKernel>
    MatrixPointerType;
  MatrixPointerType pKernel;
if (isShapeModel == true) {
  pKernel.reset((statismo::ScalarValuedKernel *)
                      it->second.createKernelShape(opt.vKernelParameters));
} else {
  if (Dimenstionality == Dimensionality2D) {
    pKernel.reset((statismo::ScalarValuedKernel *)it->second
                      .createKernel2DDeformation(opt.vKernelParameters));
  } else {
    pKernel.reset((statismo::ScalarValuedKernel *)it->second
                      .createKernel3DDeformation(opt.vKernelParameters));
  }
}
//////////////////////////////////////////////////////////////////////////
template <class TPoint>
class GaussianKernel : public statismo::ScalarValuedKernel {
public:
  // 这两个必须要有，并且不能边
  typedef typename TPoint::CoordRepType CoordRepType;
  typedef vnl_vector VectorType;

  // 只有一个参数sigma
  GaussianKernel(double sigma)
      : m_sigma(sigma), m_sigma2(-1.0 / (sigma * sigma)) {}

  // 内联重载操作符()，定义kernel计算
  // kernel是一个二元函数，返回kernel函数值
  inline double operator()(const TPoint &x, const TPoint &y) const {
    VectorType xv = x.GetVnlVector();
    VectorType yv = y.GetVnlVector();

    VectorType r = yv - xv;

    return exp((double)dot_product(r, r) * m_sigma2);
  }

  std::string GetKernelInfo() const {
    std::ostringstream os;
    os << "GaussianKernel(" << m_sigma << ")";
    return os.str();
  }

private:
  double m_sigma;
  double m_sigma2;
};

融合kernel

代码在文件KernelCombinators.h中，详细内容请在文件中具体查阅，这里简单说明一下
如果读取已有model，则需要读取model的kernel，并与当前设置的kernel融合，这个程序就是简单的求和
kernel的求和

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
// 用法
pModelBuildingKernel.reset(new statismo::SumKernel(
          pStatModelKernel.get(), pScaledKernel.get()));
/////////////////////////////////////////////////////////////////////////////////
/**
 * A (matrix valued) kernel, which represents the sum of two matrix valued
 * kernels.
 * 这个代码就是简单的将kernel求和，也是这个程序使用的
 */
template <class TPoint> class SumKernel : public MatrixValuedKernel {
public:
  typedef MatrixValuedKernel MatrixValuedKernelType;

  SumKernel(const MatrixValuedKernelType *lhs,
            const MatrixValuedKernelType *rhs)
      : MatrixValuedKernelType(lhs->GetDimension()), m_lhs(lhs), m_rhs(rhs) {
    if (lhs->GetDimension() != rhs->GetDimension()) {
      throw StatisticalModelException(
          "Kernels in SumKernel must have the same dimensionality");
    }
  }

  // 这里求和
  MatrixType operator()(const TPoint &x, const TPoint &y) const {
    return (*m_lhs)(x, y) + (*m_rhs)(x, y);
  }

  std::string GetKernelInfo() const {
    std::ostringstream os;
    os << m_lhs->GetKernelInfo() << " + " << m_rhs->GetKernelInfo();
    return os.str();
  }

private:
  const MatrixValuedKernelType *m_lhs;
  const MatrixValuedKernelType *m_rhs;
};

kernel相乘

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
/**
 * A (matrix valued) kernel, which represents the product of two matrix valued
 * kernels.
 */

template <class TPoint>
class ProductKernel : public MatrixValuedKernel {

public:
  typedef MatrixValuedKernel MatrixValuedKernelType;

  ProductKernel(const MatrixValuedKernelType *lhs,
                const MatrixValuedKernelType *rhs)
      : MatrixValuedKernelType(lhs->GetDimension()), m_lhs(lhs), m_rhs(rhs) {
    if (lhs->GetDimension() != rhs->GetDimension()) {
      throw StatisticalModelException(
          "Kernels in SumKernel must have the same dimensionality");
    }
  }

  // 将kernel逐元素相乘
  MatrixType operator()(const TPoint &x, const TPoint &y) const {
    return (*m_lhs)(x, y) * (*m_rhs)(x, y);
  }

  std::string GetKernelInfo() const {
    std::ostringstream os;
    os << m_lhs->GetKernelInfo() << " * " << m_rhs->GetKernelInfo();
    return os.str();
  }

private:
  const MatrixValuedKernelType *m_lhs;
  const MatrixValuedKernelType *m_rhs;
};

设定kernel

kernel可以单独设置也可以将多个kernel融合设置，上面贴出了单个kernel设定的模板和kernel融合的模板，这里贴一下其他的设定

根据维度设定kernel的每个分块
因为kernel计算出来的值是一个常数，所以需要修正其维度

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
 /**
 * Takes a scalar valued kernel and creates a matrix valued kernel of the given
 * dimension. The new kernel models the output components as independent, i.e.
 * if K(x,y) is a scalar valued Kernel, the matrix valued kernel becomes
 * Id*K(x,y), where Id is an identity matrix of dimensionality d.
 */
template <class TPoint>
class UncorrelatedMatrixValuedKernel : public MatrixValuedKernel {
public:
  typedef MatrixValuedKernel MatrixValuedKernelType;

  UncorrelatedMatrixValuedKernel(const ScalarValuedKernel *scalarKernel,
                                 unsigned dimension)
      : MatrixValuedKernelType(dimension), m_kernel(scalarKernel),
        m_ident(MatrixType::Identity(dimension, dimension)) {}

  // m_ident是根据维度生成的单位矩阵
  MatrixType operator()(const TPoint &x, const TPoint &y) const {

    return m_ident * (*m_kernel)(x, y);
  }

  virtual ~UncorrelatedMatrixValuedKernel() {}

  std::string GetKernelInfo() const {
    std::ostringstream os;
    os << "UncorrelatedMatrixValuedKernel(" << (*m_kernel).GetKernelInfo()
       << ", " << this->m_dimension << ")";
    return os.str();
  }

private:
  const ScalarValuedKernel *m_kernel;
  MatrixType m_ident;
};

设定kernel的尺度

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
/**
 * A (matrix valued) kernel, which represents a scalar multiple of a matrix
 * valued kernel.
 */

template <class TPoint> class ScaledKernel : public MatrixValuedKernel {
public:
  typedef MatrixValuedKernel MatrixValuedKernelType;

  ScaledKernel(const MatrixValuedKernelType *kernel, double scalingFactor)
      : MatrixValuedKernelType(kernel->GetDimension()), m_kernel(kernel),
        m_scalingFactor(scalingFactor) {}

  // 设定kernel的洗漱
  MatrixType operator()(const TPoint &x, const TPoint &y) const {
    return (*m_kernel)(x, y) * m_scalingFactor;
  }

  std::string GetKernelInfo() const {
    std::ostringstream os;
    os << (*m_kernel).GetKernelInfo() << " * " << m_scalingFactor;
    return os.str();
  }

private:
  const MatrixValuedKernelType *m_kernel;
  double m_scalingFactor;
};

使用LowRankGP构建model

这个程序只需要reference和kernel就可以了，PCA的方法需要读取很多数据

根据kernel生成model，平均值为根据reference维度得到的0向量

1
2
3
4
5
6
7
8
typedef itk::LowRankGPModelBuilder ModelBuilderType;
typename ModelBuilderType::Pointer gpModelBuilder = ModelBuilderType::New();
gpModelBuilder->SetRepresenter(pRepresenter);

typedef itk::StatisticalModel StatisticalModelType;
typename StatisticalModelType::Pointer pModel;
pModel = gpModelBuilder->BuildNewModel(pMean, *pModelBuildingKernel.get(),
                                      opt.iNrOfBasisFunctions);

相关的数学推导有时间再进行