基于nanoflann 的欧氏聚类

1、原理

         三维空间欧氏距离聚类算法中，涉及到唯一距离参数d。当点之间距离小于距离d 时，表明两点局部相连，属于同一簇点集。其中点之间相连满足如下性质：若点A 与点B 相连，点B 与点C 相连，则点A 与点C 也相连。下图展示了基于欧氏距离聚类的原理，当距离阈值设置合理时，使用空间欧氏距离聚类算法对点进行聚类，所有的点被划分成左右两簇点集，并且两簇点云中任意两点相连。即对于一簇点中任意点，总能从该簇点中找到离该点距离小于d的另外一点。对于左边点集中C点，右边点集中距离C点最近的D点，两点之间距离大于距离阈值，经过聚类后C点与D点归属于两簇不同的点集。

2、源代码展示

直接将nanoflann.hpp放入到工程下，直接运行

#include<vector>
#include <fstream>
#include<iostream>
using namespace std;
#include"nanoflann.hpp"
#include"IO.h"
#include<sstream>
using namespace nanoflann;

typedef KDTreeSingleIndexAdaptor < L2_Simple_Adaptor<double, PointCloud2<double>>,
	PointCloud2<double>, 3 > kdTree;

struct Point3D
{
	double x;
	double y;
	double z;
};

//基于欧氏聚类的点云聚类
vector<vector<Point3D>> EucliClusters(vector<Point3D> points, double dis)
{
	vector<vector<Point3D>> result;

	vector<int> flag(points.size(), 0);
	PointCloud2<double> kdPts;
	kdPts.pts.resize(points.size());

	for (int i = 0; i < points.size(); i++)
	{
		kdPts.pts[i].x = points[i].x;
		kdPts.pts[i].y = points[i].y;
		kdPts.pts[i].z = points[i].z;
	}

	kdTree* currKDTree = new kdTree(3, kdPts, KDTreeSingleIndexAdaptorParams(10));

	currKDTree->buildIndex();
	std::vector<std::pair<size_t, double>> pairs;
	SearchParams params;
	double quertPt[3];//查询点

	for (int i = 0; i < points.size(); i++)
	{
		if (flag[i]==0)//当前点没有被搜索过
		{
			quertPt[0] = points[i].x;
			quertPt[1] = points[i].y;
			quertPt[2] = points[i].z;

			const size_t numPairs = currKDTree->radiusSearch(&quertPt[0], dis*dis, pairs, params);

			vector<Point3D> onecluster;//当前簇点集
			vector<Point3D> seed;//种子点集
			for (size_t j = 0; j < numPairs; ++j)
			{
				int id = pairs[j].first;
				if (flag[id] == 0)//未增长的点
				{
					flag[id] = 1;//标示标记过

					Point3D onepoint = points[pairs[j].first];
					seed.push_back(onepoint);
					onecluster.push_back(onepoint);
				}
			}

			//while循环直至所有点
			while (seed.size() >= 1)
			{
				Point3D onepoint = seed.back();//取最后一个元素
				seed.pop_back();//将最后一个元素删除

				quertPt[0] = onepoint.x;
				quertPt[1] = onepoint.y;
				quertPt[2] = onepoint.z;

				const size_t numPairs = currKDTree->radiusSearch(&quertPt[0], dis*dis, pairs, params);

				for (int j = 0; j < numPairs; j++)
				{
					int id = pairs[j].first;

					if (flag[id] == 0)//未增长的点
					{
						flag[id] = 1;//标示标记过

						Point3D onepoint = points[pairs[j].first];
						seed.push_back(onepoint);
						onecluster.push_back(onepoint);
					}
				}
			}
			result.push_back(onecluster);
		}
	}
	return result;
}

void main()
{
	//ifstream infile("..\\测试数据\\half落叶松.txt", ios::in);
	ifstream infile("..\\测试数据\\trees.xyz", ios::in);
	PointCloud2<double> kdPts;
	char line[128];


	vector<Point3D> points;
	Point3D pt;
	while (infile.getline(line, sizeof(line)))
	{
		std::stringstream word(line);
		word >> pt.x;
		word >> pt.y;
		word >> pt.z;

		points.push_back(pt);
	}
	double dis = 0.3;
	vector<vector<Point3D>> clusters = EucliClusters(points, dis);

	ofstream outfile("D:\\separation_two_trees.txt", ios::out);
	srand((int)time(NULL));
	for (int i = 0; i < clusters.size();i++)
	{
		double r = rand() % 255;
		double g = rand() % 255;
		double b = rand() % 255;
		for (int j = 0; j < clusters[i].size();j++)
		{
			outfile << fixed << setprecision(3) << clusters[i][j].x << " " << clusters[i][j].y << " " << clusters[i][j].z << " " <<
				fixed << setprecision(0) << r << " " << g << " " << b << endl;
		}
	}
	outfile.close();
	system("pause");
}