三维点云处理01-PCA代码实现

三维点云处理01-PCA代码实现

基础知识点

1.特征值与特征向量的意义
2.方阵的特征值分解（nxn）
3.一般矩阵的特征值分解（mxn）——SVD 分解
4.谱定理
5.瑞立熵

代码实现（python）

#!/usr/bin/python3

import os
import numpy as np
from pyntcloud import PyntCloud

def PCA(data,sort=True):
	'''
	PCA函数实现：
	输入：
	data：点云数据
	sort:是否排序
	输出：
	w:特征值
	v:特征向量
	'''
	
	#获取点云总数
	N = data.shape[0]
	#将点云数据转为numpy
	X = data.to_numpy()
	
	#PCA-01:将点云进行归一化处理
	mu = np.mean(X,axis=0)
	X_normalized = X - mu
	
	#PCA-02:对归一化后的点云求协方差矩阵(因为数据为Nx3，即列为属性，行为观测，所以rowvar=False, bias=True表示有偏估计）
	func = np.cov
	H = func(X_normalized, rowvar=False, bias=True)
	
	#PCA-03：对协方差矩阵进行SVD分解，获得特征值和特征向量
	eigenvalues, eigenvectors = np.linalg.eig(H)
	
	#PCA-04:获得特征值从大到小排序的索引，然后根据特征值排序索引，对特征值和特征向量进行排序
	if sort:
		sort = eigenvalues.argsort()[::-1]
		eigenvalues = eigenvalues[sort]
		eigenvectors = eigenvectors[:,sort] #列对应的是特征向量
	
	return eigenvalues,eigenvectors
	
def main(point_cloud_filename):
	'''
	主函数：
	输入：点云文件名
	'''
	#从文件中加载点云
	point_cloud_pynt = PyntCloud.from_file(point_cloud_filename)
	
	#从点云中获取点
	points = point_cloud_pynt.points
	
	#使用主成分分析获得输入点云的特征值和特征向量
	w,v = PCA(points)
	#点云主方向对应的向量
	point_cloud_vector = v[:,0]

代码实现（C++）

• pca.h

#pragma once

#include<Eigen/Dense>

#include"utils/utils.hpp"

class PCA {
public:
	enum class eigen_vector_order : uint8_t {
		ASCENDING, DESCENDING
};
//函数声明
	void input(const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic> & input_matirx);
	bool compute(eigen_vector_order order = eigen_vector_order::DESCENDING);

	const EIgen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& get_input_matrix();
	const Eigen::matrix<double, Eigen::Dynamic, Eigen::Dynamic>& get_eigen_values();
	const Eigen::matrix<double, Eigen::Dynamic, Eigen::Dynamic>& get_eigen_vector();

private:
	void sort_eigen_vector(eigen_vector_order order = eigen_vector_order::DESCENDING);

	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> _input_matrix;
	Eigen::Matrix<double,Eigen::Dynamic, Eigen::Dynamic> _centered_matrix;
	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> _covariance_matirx;

	Eigen::Matrix<double, Eigen::Dynamic, EIgen::Dynamic> _eigen_values;
	Eigen::Matrix<double, Eigen::Dynamic, EIgen::Dynamic> _eigen_vectors;
}

• pca.cpp

#include"PCA_Downsample/pca.h"

#include<iostream>
#include<vector>

void PCA::input(const Eigen::MAtrix<double,Eigen::Dynamic, Eigen::Dynamic>& input_matrix)
{
	_input_matrix = input_matrix;
}

bool PCA::compute(eigen_vector_orde)
{
	//对输入矩阵进行归一化处理
	_centered_matrix = _input_matrix.colwise() - _input_matrix.rowwise().mean();

	if(_input_matrix.cols == 1)
	{
		_covariance_matrix = (_centered_matrix * _centered_matrix.adjoint());
	}
	else
	{
		_covariance_matrix = (_centered_matrix * _centered_matrix.adjoint()) / (_input_matrix.cols() - 1);
	}
	
	//求解特征值和特征向量
	Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> eigen_solver(_covariance_matrix);
	_eigen_values = eigen_solver.eigenvalues();
	_eigen_vectors = eigen_solver.eigenvectors();

	//根据特征值大小对特征值和特征向量进行排序
	sort_eigen_vector(order);
	return true;
}

//获取输入矩阵
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& PCA::get_input_matrix()
{
	return _input_matrix;
}

//获取特征值
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& PCA::get_eigen_values()
{
	return _eigen_values;
}

//获取特征向量
const Eigen::Matrix<double, Eigen:Dynamic, Eigen::Dynamic>& PCA::get_eigen_vector()
{
	return _eigen_vectors;
}

void PCA::sort_eigen_vector(eigen_vector_order order)
{
	std::vector<std::pair<int,int>> eigen_value_index_vector;
	for(int i = 0; i <_eigen_values.size(); ++i)
	{
		eigen_value_index_vector.push_back(std::make_pair(_eigen_values(i,0),i);
	}
	if(order == eigen_vector_order::ASCENDING)
	{
		std::sort(std::begin(eigen_value_index_vector),std::end(eigen_value_index_vector),std::less<std::pair<int,int>>());
	}
	else
	{
		std::sort(std::begin(eigen_value_index_vector), std::end(eigen_value_index_vector), std::greater<std::pair<int,int>>());
	}
	auto sorted_eigen_values = Eigen::Matrix<double, Eigen::Dynamic,1>(_eigen_values.rows());
	auto sorted_eigen_vectors = Eigen::Matrix<double, Eigen::Dynamic,Eigen::Dynamic>(_eigen_vector,rows(),_eigen_vector.cols());
	for(int i = 0; i < _eigen_values.size(); ++i)
	{
		sorted_eigen_values = _eigen_values(eigen_value_index_vector[i].second,0);
		sorted_eigen_vectors.col(i) = _eigen_vectors.col(eigen_value_index_vector[i].second);
	}
	
	_eigen_values = sorted_eigen_values;
	_eigen_vectors = sorted_eigen_vectors;
}

• pca_test.cpp

#include"PCA_Downsample/pca.h"

int main(int argc, char** argv) {
	//创建输入矩阵
	//矩阵数据类型：double
	//矩阵维度：2行，3列
	Eigen::Matrix<double,2,3> data;
	//向矩阵中输入数据
	data << 1,2,3
			1,2,3;
	
	PCA pca;
	pca.input(data);
	pca.compute(PCA::eigen_vector_order::DESCENDING);

	//输出经过主成分分析后的特征值和特征向量
	std::cout<<"eigen values:\n" << pca.get_eigen_values() << std::endl;
	std::cout<<"eigen vector: \n" << pca.get_eigen_vectoe() <<std::endl;
	return 0;
}