精华内容
下载资源
问答
  • kd树的构建
    2022-03-25 15:49:14

    在这里插入图片描述

    1、三维空间点代码:

    1、kdtree.h

    /**
    
      * @file   kdtree.h
      * @brief Thisis a brief description.
      * @author dongjian
      * @par   Copyright (c):
      *         All Rights Reserved
      * @date   2018:04:24 
      *  @note   mattersneeding attention
      */  
    #ifndef _5383BD42_370E_4C00_A25E_AD4403E5656A
    #define _5383BD42_370E_4C00_A25E_AD4403E5656A
    #include <vector>
    #include <queue>
    #include <iostream>
    
    namespace kdtree
    {
    	class Point //空间中的三维点
    	{
    	public:
    	    double x;
    		double y;
    		double z;
            // 末尾const的作用是 不会改变本类对象中的实例成员变量
    		double Distance(const Point& p) const;
            // 有参构造函数。
    		Point(double ox, double oy, double oz);
    		void printPoint(){
    			std::cout << "最近点:("<< this->x << "," << this->y << "," << this->z << ")"<<std::endl;
    			
    		}
    	};
    
    	class BBox
    	{
    	public:
    		double xmin, ymin, zmin;
    		double xmax, ymax, zmax;
    
    		bool Contains(const Point& p) const;
    		bool Contains(const BBox& bbox) const;
            // 相交的意思
    		bool Intersects(const BBox& bbox) const;
    
    		BBox(double x0, double x1, double y0, double y1, double z0, double z1);
    
            // 静态函数,返回值是BBox
    		static BBox UniverseBBox();
    		double ShortestDistance(const Point& p) const;
    	};
    
    	enum Dimension
    	{
    		X = 0,
    		Y = 1,
    		Z = 2
    	};
    
    	class Node
    	{
    	public:
    		Node* left;			//left child
    		Node* right;		//right child
    		int begin;			//start index [close
    		int end;				//end index  (open
    		Dimension dimension;	//cut dimension
    		double pivot;		//cut value
    
    		Node(int b, int e, Dimension dim, double piv);
    		bool IsLeaf() const;
            // 左子叶信封 什么鬼?
    		BBox LeftSubTreeEnvelope(const BBox& current_entext) const;
    		BBox RightSubTreeEnvelope(const BBox& current_entext) const;
    	};
    
    	class KDTree
    	{
    	public:
            // 设置输入,构建树
    		void SetInput(const std::vector<Point>& inputs);
    
            // 搜索
    		void RangeSearch(const BBox& bbox, std::vector<int>& indices) const ;
    
    		int NearestNeighbourSearch(const Point& searchpoint) const ;
    
    		void  NearestNeighbourSearch(const Point& searchpoint, int k, std::vector<int>& indices) const;
    
    		static int _leaf_max_size;
    	private:
    		const std::vector<Point>* pData; // 2维数组
    		std::vector<int> indices;
    		Node* pRoot;
    
    	private:
    		struct comparepair
    		{
    			bool operator()(const std::pair<int, double>&, const std::pair<int, double>&);
    		};
    
    		typedef std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double>>, comparepair> index_distance;
    
    		Node* DivideTree(const std::vector<Point>& inputs,int start,int end);
    		Dimension SelectDimension(const std::vector<Point>& inputs, int start, int end);
    		int Partition(int start, int end, Dimension dim);
    
    		double GetAt(int index, Dimension dim);
    
    		void _range_search(const BBox& search_bbox, Node* pNode, const BBox& current_extent,std::vector<int>& indices) const ;
    
    		int _nearest_neighbour(const Point& searchPoint, Node* pNode, const BBox& current_extent, double& current_shorest_dis) const ;
    
    		void _nearest_neighbour(const Point& searchPoint, int k, index_distance& ins, Node* pNode, const BBox& current_extent, double& current_shorest_dis) const ;
    	};
    
    }
    #endif //_5383BD42_370E_4C00_A25E_AD4403E5656A
    

    2、kdtree.cpp 含demo实现

    #include "kdtree.h"
    #include <algorithm>
    #include <iostream>
    #include <math.h>
    using namespace std;
    using namespace kdtree;
    
    Point::Point(double ox, double oy, double oz) :
    x(ox),
    y(oy),
    z(oz)
    {}
    
    double Point::Distance(const Point& p) const
    {
    	double dx = p.x - x;
    	double  dy = p.y - y;
    	double dz = p.z - z;
    	return sqrt(dx*dx + dy*dy + dz*dz);
    }
    
    BBox Node::LeftSubTreeEnvelope(const BBox& current_entext) const
    {
    	BBox leftRegion(current_entext);
    	switch (dimension)
    	{
    	case X:
    		leftRegion.xmax = pivot;
    		break;
    	case Y:
    		leftRegion.ymax = pivot;
    		break;
    	case Z:
    		leftRegion.zmax = pivot;
    		break;
    	}
    	return leftRegion;
    }
    
    BBox Node::RightSubTreeEnvelope(const BBox& current_entext) const
    {
    	BBox rightRegion(current_entext);
    	switch (dimension)
    	{
    	case X:
    		rightRegion.xmin = pivot;
    		break;
    	case Y:
    		rightRegion.ymin = pivot;
    		break;
    	case Z:
    		rightRegion.zmin = pivot;
    		break;
    	}
    	return rightRegion;
    }
    
    bool BBox::Contains(const BBox& bbox) const
    {
    	if (bbox.xmin < xmin)		return false;
    	if (bbox.xmax > xmax)		return false;
    	if (bbox.ymin < ymin)		return false;
    	if (bbox.ymax < ymax)		return false;
    	if (bbox.zmin < zmin)		return false;
    	if (bbox.zmax > zmax)		return false;
    	return true;
    }
    
    bool BBox::Contains(const Point& p) const
    {
    	return p.x >= xmin && p.x <= xmax && p.y >= ymin && p.y <= ymax
    		&& p.z >= zmin && p.z <= zmax;
    }
    
    bool BBox::Intersects(const BBox& bbox) const
    {
    	if (bbox.xmin > xmax || bbox.xmax < xmin) return false;
    	if (bbox.ymin > xmax || bbox.ymax < ymin) return false;
    	if (bbox.zmin > zmax || bbox.zmax < zmin) return false;
    	return true;
    }
    
    BBox BBox::UniverseBBox()
    {
    	double DOUBLE_MAX = std::numeric_limits<double>::max();
    	return BBox(-DOUBLE_MAX, DOUBLE_MAX, -DOUBLE_MAX, 
    		DOUBLE_MAX, -DOUBLE_MAX, DOUBLE_MAX);
    }
    
    BBox::BBox(double x0, double x1, double y0, double y1, double z0, double z1) :
    xmin(x0),
    xmax(x1),
    ymin(y0),
    ymax(y1),
    zmin(z0),
    zmax(z1)
    {}
    
    double BBox::ShortestDistance(const Point& p) const
    {
    	double dx = xmin - p.x > 0 ? xmin - p.x : (p.x - xmax > 0 ? p.x - xmax : 0);
    	double dy = ymin - p.y > 0 ? ymin - p.y : (p.y - ymax > 0 ? p.y - ymax : 0);
    	double dz = zmin - p.z > 0 ? zmin - p.z : (p.z - zmax > 0 ? p.z - zmax : 0);
    	return sqrt(dx*dx + dy*dy + dz*dz);
    }
    
    
    Node::Node(int b, int e, Dimension dim, double piv) :
    begin(b),
    end(e),
    dimension(dim),
    pivot(piv),
    left(nullptr),
    right(nullptr)
    {
    
    }
    
    bool Node::IsLeaf() const
    {
    	return left == nullptr && right == nullptr;
    }
    // 选择维度?
    /**
     * @brief 
     * 
     * @param inputs 输入原始数据 
     * @param start 起始索引
     * @param end 
     * @return Dimension  输入的数据中,在x y z轴上哪个跨度大,就输出哪个枚举类型
     */
    Dimension KDTree::SelectDimension(const std::vector<Point>& inputs, int start, int end)
    {
    	struct cmpobj
    	{
    		Dimension dim;
    		cmpobj(Dimension d) :
    			dim(d){}
    
    		bool operator()(const Point& p0, const Point p1) const
    		{
    			if (dim == X)
    				return p0.x < p1.x;
    			if (dim == Y)
    				return p0.y < p1.y;
    			if (dim == Z)
    				return p0.z < p1.z;
    			return false;
    		}
    	};
    	// span存的是最大-最小值
    	double span[3];
    	// 该pair 的first指向[first,last)范围内最小值元素,second指向最大值元素。
    	auto pair =  std::minmax_element(inputs.begin() + start, inputs.begin() + end,cmpobj(X));
    	span[0] = pair.second->x - pair.first->x;
    
    	pair = std::minmax_element(inputs.begin() + start, inputs.begin() + end, cmpobj(Y));
    	span[1] = pair.second->y - pair.first->y;
    
    	pair = std::minmax_element(inputs.begin() + start, inputs.begin() + end, cmpobj(Z));
    	span[2] = pair.second->z - pair.first->z;
    
    	// 获得的索引是 span里面最大元素与第一个元素之间的距离,也就是获得最大元素的索引。
    	auto index = std::distance(span, std::max_element(span, span + 3));
    	
    	return static_cast<Dimension>(index);
    }
    
    double KDTree::GetAt(int index, Dimension dim)
    {
    	auto p = (*pData)[indices[index]];
    	return  dim == X ? p.x : (dim == Y ? p.y : p.z);
    }
    
    int KDTree::Partition(int start, int end, Dimension dimension)
    {
    	int size = end - start;
    	if (size <= 0)
    	{
    		cout << "a serious error occurs " << start << "\t" << endl;
    		return -1;
    	}
    		
    	struct cmpobj
    	{
    		Dimension dim;
    		const std::vector<Point>* pData;
    		bool operator()(int i, int j) const 
    		{
    			if (X == dim)
    				return (*pData)[i].x < (*pData)[j].x;
    			if (Y == dim)
    				return (*pData)[i].y < (*pData)[j].y;
    			if (Z == dim)
    				return (*pData)[i].z < (*pData)[j].z;
    
    			return true;
    		}
    		cmpobj(Dimension dimension, const std::vector<Point>* pInputData):
    			dim(dimension),
    			pData(pInputData)
    		{}
    	};
    	
    	int median = start +size / 2;
    	std::nth_element(indices.begin() + start, indices.begin()+median , indices.begin() + end, cmpobj(dimension, pData));
    	return median;
    }
    
    int KDTree::NearestNeighbourSearch(const Point& searchpoint) const 
    {
    	double shortestDistance = std::numeric_limits<double>::max();
        return _nearest_neighbour(searchpoint, pRoot, BBox::UniverseBBox(), shortestDistance);
    }
    
    void KDTree::NearestNeighbourSearch(const Point& searchpoint, int k, std::vector<int>& indices) const 
    {
    	double shortestDistance = std::numeric_limits<double>::max();
    	indices.clear();
    	 
    	index_distance neighbours;
    	_nearest_neighbour(searchpoint, k, neighbours, pRoot, BBox::UniverseBBox(), shortestDistance);
    
    	while (!neighbours.empty())
    	{
    		indices.push_back(neighbours.top().first);
    		neighbours.pop();
    	}
    }
    
    bool KDTree::comparepair::operator()(const std::pair<int, double>& p0, const std::pair<int, double>& p1)
    {
    	return p0.second < p1.second;
    }
    
    void KDTree::_nearest_neighbour(const Point& searchPoint, int k, index_distance& ins,
    	Node* pNode, const BBox& current_extent, double& current_shorest_dis) const 
    {
    	double min_shortest_distance = current_extent.ShortestDistance(searchPoint);
    	if (min_shortest_distance >= current_shorest_dis)
    		return;
    
    	if (pNode->IsLeaf())
    	{
    		for (auto i = pNode->begin; i < pNode->end; ++i)
    		{
    			double distance = (*pData)[indices[i]].Distance(searchPoint);
    			if (ins.size() < k)
    			{
    				ins.push(pair<int, double>(indices[i], distance));	//add element
    				if (k == ins.size())
    					current_shorest_dis = ins.top().second;
    			}
    			else 
    			{
    				if (distance < current_shorest_dis)
    				{
    					ins.pop();
    					ins.push(pair<int, double>(indices[i], distance));//add element
    					current_shorest_dis = ins.top().second;
    				}
    			}
    		}
    	}
    	else
    	{
    		BBox leftRegion = pNode->LeftSubTreeEnvelope(current_extent);
    		BBox rightRegion = pNode->RightSubTreeEnvelope(current_extent);
    
    		double dis_to_left = leftRegion.ShortestDistance(searchPoint);
    		double dis_to_right = rightRegion.ShortestDistance(searchPoint);
    
    		if (dis_to_left < dis_to_right)
    		{
    			_nearest_neighbour(searchPoint,k,ins,pNode->left,leftRegion,current_shorest_dis);
    			 _nearest_neighbour(searchPoint, k,ins,pNode->right, rightRegion, current_shorest_dis);
    		}
    		else
    		{
    			_nearest_neighbour(searchPoint,k,ins, pNode->right, rightRegion, current_shorest_dis);
    			_nearest_neighbour(searchPoint,k,ins, pNode->left, leftRegion, current_shorest_dis);
    		}
    	}
    }
    
    
    int KDTree::_nearest_neighbour(const Point& searchPoint, Node* pNode, const BBox& current_extent, double& current_shorest_dis) const 
    {
    	double min_shortest_distance = current_extent.ShortestDistance(searchPoint);
    	if (min_shortest_distance >= current_shorest_dis)
    		return -1;
    
    	if (pNode->IsLeaf())
    	{
    		int shortestindex =-1;
    		for (auto i = pNode->begin; i < pNode->end;++i)
    		{
    			double distance = (*pData)[indices[i]].Distance(searchPoint);
    			if (distance < current_shorest_dis)
    			{
    				shortestindex = indices[i];
    				current_shorest_dis = distance;
    			}
    		}
    		return shortestindex;
    	}
    	else
    	{
    		BBox leftRegion = pNode->LeftSubTreeEnvelope(current_extent);
    		BBox rightRegion = pNode->RightSubTreeEnvelope(current_extent);
    
    		double dis_to_left = leftRegion.ShortestDistance(searchPoint);
    		double dis_to_right = rightRegion.ShortestDistance(searchPoint);
    
    		if (dis_to_left < dis_to_right)
    		{
    			int left = _nearest_neighbour(searchPoint, pNode->left, leftRegion, current_shorest_dis);
    			int right = _nearest_neighbour(searchPoint, pNode->right, rightRegion, current_shorest_dis);
    
    			return right == -1 ? left : right;
    		}
    		else
    		{
    			int right = _nearest_neighbour(searchPoint, pNode->right, rightRegion, current_shorest_dis);
    			int left = _nearest_neighbour(searchPoint, pNode->left, leftRegion, current_shorest_dis);
    			return left == -1 ? right : left;
    		}
    
    		return -1;
    	}
    }
    
    int KDTree::_leaf_max_size = 2; //原来是15
    
    void KDTree::SetInput(const std::vector<Point>& inputs)
    {
    	pData = &inputs;// 指针的首地址 指向了inputs的地址
    	indices.resize(inputs.size());
    	for (int i = 0, n = inputs.size(); i < n; ++i)
    		indices[i] = i;// 1 2 3 4 5 6 ....
    	// Node类型的指针
    	pRoot = DivideTree(inputs, 0, inputs.size()); //划分树
    }
    
    Node* KDTree::DivideTree(const std::vector<Point>& inputs, int start, int end)
    {
    	//cout << "build " << start << "\t" << end << endl;
    
    	int size = end - start;
    	if (size <= 0)
    		return nullptr;
    
    	Dimension dim = SelectDimension(inputs, start, end); // 选择划分维度,是按X轴划分还是y轴还是z轴
    	// 中位数
    	int median = Partition(start, end, dim); //获取中位数
    
    	Node* pNode = new Node(start, end, dim, GetAt(median, dim));
    	//递归终止条件是输入的点的个数小于叶子节点最大点个数
    	if (size > _leaf_max_size)
    	{
    		pNode->left = DivideTree(inputs, start, median);
    		pNode->right = DivideTree(inputs, median, end);
    	}
    
    	return pNode;
    }
    
    void KDTree::RangeSearch(const BBox& bbox, std::vector<int>& indices) const 
    {
    	BBox universe_bbox = BBox::UniverseBBox();
    	_range_search(bbox, pRoot, universe_bbox, indices);
    }
    
    
    void KDTree::_range_search(const BBox& search_bbox, Node* pNode, const BBox& current_extent, std::vector<int>& ins) const 
    {
    	if (nullptr == pNode)
    		return;
    
    	if (pNode->IsLeaf())
    	{
    		for (int i = pNode->begin; i < pNode->end; ++i)
    		{
    			const Point& p = (*pData)[indices[i]];
    			if (search_bbox.Contains(p))
    				ins.push_back(indices[i]);
    		}
    	}
    	else
    	{
    		//trim bounding box
    		BBox leftRegion = pNode->LeftSubTreeEnvelope(current_extent);
    
    		if (search_bbox.Contains(leftRegion))
    		{
    			for (int i = pNode->left->begin; i < pNode->left->end; ++i)
    				ins.push_back(indices[i]);
    		}
    		else if (search_bbox.Intersects(leftRegion))
    		{
    			_range_search(search_bbox, pNode->left, leftRegion, ins);
    		}
    
    		BBox rightRegion = pNode->RightSubTreeEnvelope(current_extent);
    
    		if (search_bbox.Contains(rightRegion))
    		{
    			for (int i = pNode->right->begin; i < pNode->right->end; ++i)
    				ins.push_back(indices[i]);
    		}
    		else if (search_bbox.Intersects(rightRegion))
    		{
    			_range_search(search_bbox, pNode->right, rightRegion, ins);
    		}
    	}
    }
    int main(int argc, char *argv[]){
        
        std::vector<kdtree::Point> input;
    
      
        // 1、输入数据
        for (int  i = 0; i < 10; i++)
        {
            double x,y,z;
            std::cin>>x>>y>>z;
            std::cout<<x<<y<<z<<std::endl;
            kdtree::Point p(x, y, z);
           
            input.push_back(p);
        }
    	//2、构建树
    	KDTree kdtree;
    	kdtree.SetInput(input);
    
    	//3、最近邻搜索
    	//创建寻找的点
    	int k =3;
    	Point p1(1.2 , 2.1, 3.2);
    	std::vector<int> indices;//索引
    	kdtree.NearestNeighbourSearch(p1, k, indices);
    
    	//4、打印K个最近点
    	for(int i = 0; i<k; i++ ){
    
    		input[indices[i]].printPoint();
    		std::cout << "距离是:" << p1.Distance(input[indices[i]]) << std::endl;
    	}
    
        return 0;
    }
    

    结果:

    在这里插入图片描述

    2、二维空间代码

    1、kdtree_new.h

    //
    //  KDTree.h
    //  Test
    //
    //  Created by xiuzhu on 2021/7/13.
    //
    
    #ifndef kdtree_new_h
    #define kdtree_new_h
    
    #include <vector>
    
    using namespace std;
    //template <class T>
    class KDTree {
    private:
        int key;
        vector<double> root;//树及其子树的根节点
        KDTree *parent;
        KDTree *left_child;
        KDTree *right_child;
        
    public:
        KDTree();
        ~KDTree();
        bool is_empty();//判断KD树是否为空
        bool is_leaf();//判断树是否只有一个叶子节点
        bool is_root();//判断是否是树的根节点
        bool is_left();//判断该子kd树的根结点是否是其父kd树的左结点
        bool is_right();//判断该子kd树的根结点是否是其父kd树的右结点
        vector<vector<double> > Transpose(const vector<vector<double> > &Matrix);//坐标转换
        double findMiddleValue(vector<double> vec);//查找中指
        void buildKdTree(KDTree* tree, vector<vector<double> > data, unsigned depth);//构建kd树
        void printKdTree(KDTree* tree, unsigned int depth);//逐层打印KD树
        double measureDistance(vector<double> point1, vector<double> point2, unsigned method);//计算空间中两个点的距离
        vector<double> searchNearestNeighbor(vector<double> goal, KDTree *tree);在kd树tree中搜索目标点goal的最近邻
        
    };
    
    
    #endif /* KDTree_h */
    
    

    2、kdtree_new.cpp

    //
    //  KDTree.cpp
    //  Test
    //
    //  Created by xiuzhu on 2021/7/13.
    //
    #include "kdtree_new.h"
    #include <stdio.h>
    #include <iostream>
    #include <vector>
    #include <algorithm>
    #include <string>
    #include <cmath>
    
    using namespace std;
    
    
    KDTree::KDTree():parent(nullptr),left_child(nullptr),right_child(nullptr){}
    KDTree::~KDTree(){}
    bool KDTree::is_empty()
    {
        return root.empty();
    }
    bool KDTree::is_leaf()
    {
        return (!root.empty()) && right_child == nullptr && left_child == nullptr;
    }
    bool KDTree::is_root()
    {
        return (!is_empty()) && parent == nullptr;
    }
    bool KDTree::is_left()
    {
        return parent->left_child->root == root;
    }
    bool KDTree::is_right()
    {
        return parent->right_child->root == root;
    }
    //用于转换坐标
    vector<vector<double> > KDTree::Transpose(const vector<vector<double>> &Matrix)
    {
        unsigned row = Matrix.size();
        unsigned col = Matrix[0].size();
        vector<vector<double> > Trans(col,vector<double>(row,0));
        for (unsigned i = 0; i < col; ++i)
        {
            for (unsigned j = 0; j < row; ++j)
            {
                Trans[i][j] = Matrix[j][i];
            }
        }
        return Trans;
    }
    //在不同的坐标轴上寻找中值
    
    double KDTree::findMiddleValue(vector<double> vec)
    {
        sort(vec.begin(),vec.end());
        auto pos = vec.size() / 2;
        return vec[pos];
    }
    void KDTree::buildKdTree(KDTree *tree, vector<vector<double>> data, unsigned int depth)
    {
        //样本的数量
        unsigned long samplesNum = data.size();
        //终止条件
        if (samplesNum == 0)
        {
            return;
        }
        if (samplesNum == 1)
        {
            tree->root = data[0];
            return;
        }
        //样本的维度
        unsigned long k = data[0].size();//坐标轴个数
        vector<vector<double> > transData = Transpose(data);
        //选择切分属性
        unsigned splitAttribute = depth % k;
        vector<double> splitAttributeValues = transData[splitAttribute];
        //选择切分值
        double splitValue = findMiddleValue(splitAttributeValues);
        //cout << "splitValue" << splitValue  << endl;
    
        // 根据选定的切分属性和切分值,将数据集分为两个子集
        vector<vector<double> > subset1;
        vector<vector<double> > subset2;
        for (unsigned i = 0; i < samplesNum; ++i)
        {
            if (splitAttributeValues[i] == splitValue && tree->root.empty())
                tree->root = data[i];
            else
            {
                if (splitAttributeValues[i] < splitValue)
                    subset1.push_back(data[i]);
                else
                    subset2.push_back(data[i]);
            }
        }
    
        //子集递归调用buildKdTree函数
        tree->left_child = new KDTree;
        tree->left_child->parent = tree;
        tree->right_child = new KDTree;
        tree->right_child->parent = tree;
        buildKdTree(tree->left_child, subset1, depth + 1);
        buildKdTree(tree->right_child, subset2, depth + 1);
    }
    void KDTree::printKdTree(KDTree *tree, unsigned int depth)//打印
    {
        for (unsigned i = 0; i < depth; ++i)
            cout << "\t";
    
        for (vector<double>::size_type j = 0; j < tree->root.size(); ++j)
            cout << tree->root[j] << ",";
        cout << endl;
        if (tree->left_child == nullptr && tree->right_child == nullptr )//叶子节点
            return;
        else //非叶子节点
        {
            if (tree->left_child != nullptr)
            {
                for (unsigned i = 0; i < depth + 1; ++i)
                    cout << "\t";
                cout << " left:";
                printKdTree(tree->left_child, depth + 1);
            }
    
            cout << endl;
            if (tree->right_child != nullptr)
            {
                for (unsigned i = 0; i < depth + 1; ++i)
                    cout << "\t";
                cout << "right:";
                printKdTree(tree->right_child, depth + 1);
            }
            cout << endl;
        }
    }
    //计算空间中两个点的距离
    double KDTree::measureDistance(vector<double> point1, vector<double> point2, unsigned int method){
        if (point1.size() != point2.size())
        {
            cerr << "Dimensions don't match!!" ;
            exit(1);
        }
        switch (method)
        {
            case 0://欧氏距离
                {
                    double res = 0;
                    for (vector<double>::size_type i = 0; i < point1.size(); ++i)
                    {
                        res += pow((point1[i] - point2[i]), 2);
                    }
                    return sqrt(res);
                }
            case 1://曼哈顿距离
                {
                    double res = 0;
                    for (vector<double>::size_type i = 0; i < point1.size(); ++i)
                    {
                        res += abs(point1[i] - point2[i]);
                    }
                    return res;
                }
            default:
                {
                    cerr << "Invalid method!!" << endl;
                    return -1;
                }
        }
    
    }
    //在kd树tree中搜索目标点goal的最近邻
    //输入:目标点;已构造的kd树
    //输出:目标点的最近邻
    vector<double> KDTree::searchNearestNeighbor(vector<double> goal, KDTree *tree)
    {
        /*第一步:在kd树中找出包含目标点的叶子结点:从根结点出发,
        递归的向下访问kd树,若目标点的当前维的坐标小于切分点的
        坐标,则移动到左子结点,否则移动到右子结点,直到子结点为
        叶结点为止,以此叶子结点为“当前最近点”
        */
        unsigned long k = tree->root.size();//计算出数据的维数
        unsigned d = 0;//维度初始化为0,即从第1维开始
        KDTree* currentTree = tree;
        vector<double> currentNearest = currentTree->root;
        while(!currentTree->is_leaf())
        {
            unsigned index = d % k;//计算当前维
            if (currentTree->right_child->is_empty() || goal[index] < currentNearest[index])
            {
                currentTree = currentTree->left_child;
            }
            else
            {
                currentTree = currentTree->right_child;
            }
            ++d;
        }
        currentNearest = currentTree->root;
    
        /*第二步:递归地向上回退, 在每个结点进行如下操作:
        (a)如果该结点保存的实例比当前最近点距离目标点更近,则以该例点为“当前最近点”
        (b)当前最近点一定存在于某结点一个子结点对应的区域,检查该子结点的父结点的另
        一子结点对应区域是否有更近的点(即检查另一子结点对应的区域是否与以目标点为球
        心、以目标点与“当前最近点”间的距离为半径的球体相交);如果相交,可能在另一
        个子结点对应的区域内存在距目标点更近的点,移动到另一个子结点,接着递归进行最
        近邻搜索;如果不相交,向上回退*/
    
        //当前最近邻与目标点的距离
        double currentDistance = measureDistance(goal, currentNearest, 0);
    
        //如果当前子kd树的根结点是其父结点的左孩子,则搜索其父结点的右孩子结点所代表
        //的区域,反之亦反
        KDTree* searchDistrict;
        if (currentTree->is_left())
        {
            if (currentTree->parent->right_child == nullptr)
                searchDistrict = currentTree;
            else
                searchDistrict = currentTree->parent->right_child;
        }
        else
        {
            searchDistrict = currentTree->parent->left_child;
        }
    
        //如果搜索区域对应的子kd树的根结点不是整个kd树的根结点,继续回退搜索
        while (searchDistrict->parent != nullptr)
        {
            //搜索区域与目标点的最近距离
            double districtDistance = abs(goal[(d+1)%k] - searchDistrict->parent->root[(d+1)%k]);
    
            //如果“搜索区域与目标点的最近距离”比“当前最近邻与目标点的距离”短,表明搜索
            //区域内可能存在距离目标点更近的点
            if (districtDistance < currentDistance )//&& !searchDistrict->isEmpty()
            {
    
                double parentDistance = measureDistance(goal, searchDistrict->parent->root, 0);
    
                if (parentDistance < currentDistance)
                {
                    currentDistance = parentDistance;
                    currentTree = searchDistrict->parent;
                    currentNearest = currentTree->root;
                }
                if (!searchDistrict->is_empty())
                {
                    double rootDistance = measureDistance(goal, searchDistrict->root, 0);
                    if (rootDistance < currentDistance)
                    {
                        currentDistance = rootDistance;
                        currentTree = searchDistrict;
                        currentNearest = currentTree->root;
                    }
                }
                if (searchDistrict->left_child != nullptr)
                {
                    double leftDistance = measureDistance(goal, searchDistrict->left_child->root, 0);
                    if (leftDistance < currentDistance)
                    {
                        currentDistance = leftDistance;
                        currentTree = searchDistrict;
                        currentNearest = currentTree->root;
                    }
                }
                if (searchDistrict->right_child != nullptr)
                {
                    double rightDistance = measureDistance(goal, searchDistrict->right_child->root, 0);
                    if (rightDistance < currentDistance)
                    {
                        currentDistance = rightDistance;
                        currentTree = searchDistrict;
                        currentNearest = currentTree->root;
                    }
                }
            }//end if
    
            if (searchDistrict->parent->parent != nullptr)
            {
                searchDistrict = searchDistrict->parent->is_left()?
                                searchDistrict->parent->parent->right_child:
                                searchDistrict->parent->parent->left_child;
            }
            else
            {
                searchDistrict = searchDistrict->parent;
            }
            ++d;
        }//end while
        return currentNearest;
    }
    
    const int data[6][2]={{2,3},{5,4},{9,6},{4,7},{8,1},{7,2}};
    int main()
    {
        vector<vector<double> > train(6, vector<double>(2, 0)); // 初始化为6个 值为(2,0)的坐标
        for (unsigned i = 0; i < 6; ++i)
            for (unsigned j = 0; j < 2; ++j)
                train[i][j] = data[i][j];
        KDTree* kdTree = new KDTree();
        kdTree->buildKdTree(kdTree, train, 0);
        kdTree->printKdTree(kdTree, 0);
    
        vector<double> goal;
        goal.push_back(4);
        goal.push_back(6);
        vector<double> nearestNeighbor = kdTree->searchNearestNeighbor(goal, kdTree);
        vector<double>::iterator beg = nearestNeighbor.begin();
        cout << "The nearest neighbor is: ";
        while(beg != nearestNeighbor.end()) cout << *beg++ << ",";
        cout << endl;
        return 0;
    }
    
    
    
    更多相关内容
  • 用C++实现KD树构建以及递归输出。。。
  • KD树构建

    千次阅读 2017-09-14 23:00:47
    平衡kd树构建算法流程: 输入:K维空间数据集T={x1, x2, …. xn},其中 xi={xi(1), xi(2), … xi(k)}, i=1,….N 输出:kd树1:开始:构造根结点,根结点对应于包含T的k维空间的超矩形区域。选择x(1)为坐标轴,以T中...
                              平衡kd树构建算法流程:
    

    输入:K维空间数据集T={x1, x2, …. xn},其中 xi={xi(1), xi(2), … xi(k)}, i=1,….N
    输出:kd树

    1:开始:构造根结点,根结点对应于包含T的k维空间的超矩形区域。选择x(1)为坐标轴,以T中所有实例的x(1)坐标的中位数为切分点,将根结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴x(1)垂直的超平面实现。由根结点生成深度为1的左、右子结点:左子结点对应坐标x(1)小于切分点的子区域,右子结点对应于坐标x(1)大于切分点的子区域。将落在切分超平面上的实例点保存在根结点。
    2:重复。对深度为j的结点选择x(l)为切分的坐标轴,l=j%k+1,以该结点的区域中所有实例的x(l)坐标的中位数为切分点,将该结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴x(l)垂直的超平面实现。由该结点生成深度为j+1的左、右子结点:左子结点对应坐标x(l)小于切分点的子区域,右子结点对应坐标x(l)大于切分点的子区域。将落在切分超平面上的实例点保存在该结点。

    其中k的选择有特征xi的维度或knn中k的值。

    kd树的构建!例题3.2
    给定一个二维空间的数据集:
    T = {(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)}, 构造一个平衡kd树。
    为了方便,我这里进行编号A(2,3)、B(5,4)、C(9,6)、D(4,7)、E(8,1)、F(7,2)
    初始值r=0,对应x轴。

    这里写图片描述

    首先先沿 x 坐标进行切分,我们选出 x 坐标的中位点,获取最根部节点的坐标,对数据点x坐标进行排序得:
    A(2,3)、D(4,7)、B(5,4)、F(7,2)、E(8,1)、C(9,6)

    则我们得到中位点为B或者F,我这里选择F作为我们的根结点,并作出切分(并得到左右子树),如图:
    这里写图片描述


    这里写图片描述
    这里写图片描述

    python代码:
    R = 3 # Feature dimension
    
    class Node(object):
        def __init__(self, elem, lchild=None, rchild=None):
            self.elem = elem
            self.lchild = lchild
            self.rchild = rchild
    
    def KDT(root, LR, dataSource, r):
        if dataSource==[]:
            return
        data = sorted(dataSource, key=lambda x:x[r])
        r = (r+1) % R
        length = len(data)
        node = Node(data[length/2], None, None)
        if LR==0:
            root.lchild = node
            KDT(root.lchild, 0, data[:length/2], r)
            KDT(root.lchild, 1, data[length/2 + 1:], r)
        if LR==1:
            root.rchild = node
            KDT(root.rchild, 0, data[:length/2], r)
            KDT(root.rchild, 1, data[length/2 + 1:], r)
    
    
    def InitTree(dataSource, length):
        r = 0
        if dataSource==[]:
            print "Please init dataSource."
            return None
    
        data = sorted(dataSource, key=lambda x:x[r])
        r=(r+1) % R
        root = Node(data[length/2], None, None)
    
        KDT(root, 0, data[:length/2], r)
        KDT(root, 1, data[length/2 + 1:], r)
    
        print "InitTree Done."
        return root
    
    def PreOrderTraversalTree(root):
        if root:
            print root.elem,' | ',
            PreOrderTraversalTree(root.lchild)
            PreOrderTraversalTree(root.rchild)
    
    def InOrderTraversalTree(root):
        if root:
            InOrderTraversalTree(root.lchild)
            print root.elem,' | ',
            InOrderTraversalTree(root.rchild)
    
    def PostOrderTraversalTree(root):
        if root:
            PostOrderTraversalTree(root.lchild)
            PostOrderTraversalTree(root.rchild)
            print root.elem,' | ',
    
    
    if __name__ == "__main__":
        dataSource = [(2,3, 100), (5,4,70), (9,6,55), (4,7,200), (8,1,44), (7,2,0)]
        length = len(dataSource)
        root = InitTree(dataSource, length)
    
        print "PreOrder:"
        PreOrderTraversalTree(root)
        print "\nInOrder:"
        InOrderTraversalTree(root)
        print "\nPostOrder:"
        PostOrderTraversalTree(root)

    参考资料:

    https://mshj.blog.ustc.edu.cn/?p=292
    https://www.suilengea.com/show/xazcegzcv.html
    李航 统计学习方法

    展开全文
  • 问题导入:实现k近邻算法时,主要考虑的问题是如何对训练数据进行快速k近邻搜索。这在特征空间的维数大及训练数据容量大时尤其必要。k近邻法最简单的实现是线性扫描(穷举搜索),即要计算输入实例与每一个训练...kd树

    问题导入:

    实现k近邻算法时,主要考虑的问题是如何对训练数据进行快速k近邻搜索。

    这在特征空间的维数大及训练数据容量大时尤其必要。

    k近邻法最简单的实现是线性扫描(穷举搜索),即要计算输入实例与每一个训练实例的距离。计算并存储好以后,再查找K近邻。当训练集很大时,计算非常耗时。

    为了提高kNN搜索的效率,可以考虑使用特殊的结构存储训练数据,以减小计算距离的次数。



    一、 kd树简介

    1.1 什么是kd树

    根据KNN每次需要预测一个点时,我们都需要计算训练数据集里每个点到这个点的距离,然后选出距离最近的k个点进行投票。当数据集很大时,这个计算成本非常高

    kd树:为了避免每次都重新计算一遍距离,算法会把距离信息保存在一棵树里,这样在计算之前从树里查询距离信息,尽量避免重新计算。其基本原理是,如果A和B距离很远,B和C距离很近,那么A和C的距离也很远。有了这个信息,就可以在合适的时候跳过距离远的点。

    这样优化后的算法复杂度可降低到O(DNlog(N))。感兴趣的读者可参阅论文:Bentley,J.L.,Communications of the ACM(1975)。

    1989年,另外一种称为Ball Tree的算法,在kd Tree的基础上对性能进一步进行了优化。感兴趣的读者可以搜索Five balltree construction algorithms来了解详细的算法信息。

    1.2 原理

    黄色的点作为根节点,上面的点归左子树,下面的点归右子树,接下来再不断地划分,分割的那条线叫做分割超平面(splitting hyperplane),在一维中是一个点,二维中是线,三维的是面。

    黄色节点就是Root节点,下一层是红色,再下一层是绿色,再下一层是蓝色。

    1.树的建立;

    2.最近邻域搜索(Nearest-Neighbor Lookup)

    kd树(K-dimension tree)是一种对k维空间中的实例点进行存储以便对其进行快速检索的树形数据结构。kd树是一种二叉树,表示对k维空间的一个划分,构造kd树相当于不断地用垂直于坐标轴的超平面将K维空间切分,构成一系列的K维超矩形区域。kd树的每个结点对应于一个k维超矩形区域。利用kd树可以省去对大部分数据点的搜索,从而减少搜索的计算量。

    类比“二分查找”:给出一组数据:[9 1 4 7 2 5 0 3 8],要查找8。如果挨个查找(线性扫描),那么将会把数据集都遍历一遍。而如果排一下序那数据集就变成了:[0 1 2 3 4 5 6 7 8 9],按前一种方式我们进行了很多没有必要的查找,现在如果我们以5为分界点,那么数据集就被划分为了左右两个“簇” [0 1 2 3 4]和[6 7 8 9]。

    因此,根本就没有必要进入第一个簇,可以直接进入第二个簇进行查找。把二分查找中的数据点换成k维数据点,这样的划分就变成了用超平面对k维空间的划分。空间划分就是对数据点进行分类,“挨得近”的数据点就在一个空间里面。

    二、构造方法

    (1)构造根结点,使根结点对应于K维空间中包含所有实例点的超矩形区域;

    (2)通过递归的方法,不断地对k维空间进行切分,生成子结点。在超矩形区域上选择一个坐标轴和在此坐标轴上的一个切分点,确定一个超平面,这个超平面通过选定的切分点并垂直于选定的坐标轴,将当前超矩形区域切分为左右两个子区域(子结点);这时,实例被分到两个子区域。

    (3)上述过程直到子区域内没有实例时终止(终止时的结点为叶结点)。在此过程中,将实例保存在相应的结点上。

    (4)通常,循环的选择坐标轴对空间切分,选择训练实例点在坐标轴上的中位数为切分点,这样得到的kd树是平衡的(平衡二叉树:它是一棵空树,或其左子树和右子树的深度之差的绝对值不超过1,且它的左子树和右子树都是平衡二叉树)。

    KD树中每个节点是一个向量,和二叉树按照数的大小划分不同的是,KD树每层需要选定向量中的某一维,然后根据这一维按左小右大的方式划分数据。在构建KD树时,关键需要解决2个问题:

    (1)选择向量的哪一维进行划分;

    (2)如何划分数据;

    第一个问题简单的解决方法可以是随机选择某一维或按顺序选择,但是更好的方法应该是在数据比较分散的那一维进行划分(分散的程度可以根据方差来衡量)

    第二个问题中,好的划分方法可以使构建的树比较平衡,可以每次选择中位数来进行划分。

    三、案例分析

    3.1 树的建立

    给定一个二维空间数据集:T={(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)},构造一个平衡kd树。

    (1)思路引导:

    根结点对应包含数据集T的矩形,选择x(1)轴,6个数据点的x(1)坐标中位数是6,这里选最接近的(7,2)点,以平面x(1)=7将空间分为左、右两个子矩形(子结点);接着左矩形以x(2)=4分为两个子矩形(左矩形中{(2,3),(5,4),(4,7)}点的x(2)坐标中位数正好为4),右矩形以x(2)=6分为两个子矩形,如此递归,最后得到如下图所示的特征空间划分和kd树。

    3.2 最近领域的搜索

    假设标记为星星的点是 test point, 绿色的点是找到的近似点,在回溯过程中,需要用到一个队列,存储需要回溯的点,在判断其他子节点空间中是否有可能有距离查询点更近的数据点时,做法是以查询点为圆心,以当前的最近距离为半径画圆,这个圆称为候选超球(candidate hypersphere),如果圆与回溯点的轴相交,则需要将轴另一边的节点都放到回溯队列里面来。

    样本集{(2,3),(5,4), (9,6), (4,7), (8,1), (7,2)}

    3.2.1 查找点(2.1,3.1)

    在(7,2)点测试到达(5,4),在(5,4)点测试到达(2,3),然后search_path中的结点为<(7,2),(5,4), (2,3)>,从search_path中取出(2,3)作为当前最佳结点nearest, dist为0.141;

    然后回溯至(5,4),以(2.1,3.1)为圆心,以dist=0.141为半径画一个圆,并不和超平面y=4相交,如上图,所以不必跳到结点(5,4)的右子空间去搜索,因为右子空间中不可能有更近样本点了。

    于是再回溯至(7,2),同理,以(2.1,3.1)为圆心,以dist=0.141为半径画一个圆并不和超平面x=7相交,所以也不用跳到结点(7,2)的右子空间去搜索。

    至此,search_path为空,结束整个搜索,返回nearest(2,3)作为(2.1,3.1)的最近邻点,最近距离为0.141。

    3.2.2 查找点(2,4.5)

    在(7,2)处测试到达(5,4),在(5,4)处测试到达(4,7)【优先选择在本域搜索】,然后search_path中的结点为<(7,2),(5,4), (4,7)>,从search_path中取出(4,7)作为当前最佳结点nearest, dist为3.202;

    然后回溯至(5,4),以(2,4.5)为圆心,以dist=3.202为半径画一个圆与超平面y=4相交,所以需要跳到(5,4)的左子空间去搜索。所以要将(2,3)加入到search_path中,现在search_path中的结点为<(7,2),(2, 3)>;另外,(5,4)与(2,4.5)的距离为3.04 < dist = 3.202,所以将(5,4)赋给nearest,并且dist=3.04。

    回溯至(2,3),(2,3)是叶子节点,直接平判断(2,3)是否离(2,4.5)更近,计算得到距离为1.5,所以nearest更新为(2,3),dist更新为(1.5)

    回溯至(7,2),同理,以(2,4.5)为圆心,以dist=1.5为半径画一个圆并不和超平面x=7相交, 所以不用跳到结点(7,2)的右子空间去搜索。

    至此,search_path为空,结束整个搜索,返回nearest(2,3)作为(2,4.5)的最近邻点,最近距离为1.5。

    四、总结

    • kd树的构建过程【知道】
      • 1.构造根节点
      • 2.通过递归的方法,不断地对k维空间进行切分,生成子节点
      • 3.重复第二步骤,直到子区域中没有示例时终止
      • 需要关注细节:a.选择向量的哪一维进行划分;b.如何划分数据
    • kd树的搜索过程【知道】
      • 1.二叉树搜索比较待查询节点和分裂节点的分裂维的值,(小于等于就进入左子树分支,大于就进入右子树分支直到叶子结点)
      • 2.顺着“搜索路径”找到最近邻的近似点
      • 3.回溯搜索路径,并判断搜索路径上的结点的其他子结点空间中是否可能有距离查询点更近的数据点,如果有可能,则需要跳到其他子结点空间中去搜索
      • 4.重复这个过程直到搜索路径为空

    展开全文
  • } FLANN中kd树构建(超平面的确定) 维度划分方式对查询性能的影响 网上很多其他实现会将kd树的维度划分就保持x,y,z,…,x,y,z…的顺序。实际上,并不一定要保持这样的顺序。而且,不同的划分顺序,往往能给kd树的...

    flann源码参考:flann: https://github.com/flann-lib/flannsudo apt install libflann-dev

    目录

    K-最近邻搜索(K-Nearest Neighbour,KNN)

    什么是kd树

    FLANN中kd树的构建(超平面的确定)

    维度划分方式对查询性能的影响

    大规模数据swap的方式

    FLANN中超平面确定方法

    其他超平面确定方法

    FLANN中kd树的查询

    最近距离表&搜索半径

    近似最近邻查询(Approximate Nearest Neighbor Search,ANNS)

    FLANN中kd树查询策略

    多线程并行查询


    K-最近邻搜索(K-Nearest Neighbour,KNN)

    K-最近邻(K-Nearest Neighbour, KNN)算法是一种基本分类与回归方法,属于监督学习方法,其工作机制非常简单:给定测试样本,基于某种距离度量找出训练集中与其最靠近的k个训练样本。算法的输入输出如下:

     在特征匹配的应用中,我们通常找到最接近2个最近邻点(k=2),并在这两个关键点中,若最近的距离除以次近的距离小于指定阈值ratio,则接受这一对匹配点。Low推荐ratio的取值为0.8,但Low对大量存在尺度、旋转和亮度变化的两幅图片进行匹配,实验结果表明ratio取值在0. 4~0. 6为最佳,小于0. 4会造成非常少的匹配点,大于0. 6会造成大量的错误匹配点。

    Low D G . Distinctive Image Features from Scale-Invariant Keypoints[J]. International Journal of Computer Vision, 2004.

    knn的实现:可以简单分为:索引构建和搜索策略

    比如说,kd树就是其中一个比较主流索引构建的方法,当索引构建完毕剩下的就是确定搜索策略。不同的索引构建方法会带来匹配效率和匹配效果的不同。

    常见的索引构建方法:有基于hash函数的局部敏感哈希LSH、基于kmeans的kd树(kmeans作为超平面确定方法)、层次graph等。

    【大咖讲堂】杨晓春:高维数据的近似最近邻搜索及其在跨模态检索中的应用_哔哩哔哩_bilibili

    什么是kd树

    kd树的全称是k-dimensional tree,是一种分割k维数据空间的数据结构。

    kd树是一种二叉树。在每一层上kd树沿着按照划分维度将数据分为两组,两组数据依次进行分割形成子树。分割的对象称之为超平面(hyperplane),超平面垂直于对应维度的轴。理想的超平面是对应维度的中位数(median),这样可以保证树的平衡(balance),从而降低树的深度。

    下图所示,是k=2时的一颗kd树。需要提醒的是进行划分(split)的维度的顺序可以是任意的,不一定按照x,y,z,x,y,z…的顺序进行。每一个节点都会记录划分的维度。FLANN中有划分维度选择(或者叫做,超平面确定)的算法(比如随机最大方差,或者,kmeans聚类)。 

    参考博文:KD树的主要算法以及FLANN(PCL)的实现分析 - Fun With GeometryFun With Geometry

    在FLANN中,将kd树的节点定义如下:

    ../flann/src/cpp/flann/algorithms/kdtree_index.h

    struct Node
        {
        	/**
             * Dimension used for subdivision.
             * 如果是中间节点,则为划分维度
             * 如果是叶节点的话,则为该元素在数据集中的索引
             */
            int divfeat;
            /**
             * The values used for subdivision.
             * 划分维度对应的值
             */
            DistanceType divval;
            /**
             * Point data
             * 保存节点对应的原始元素,只有叶节点才会保存
             */
            ElementType* point;
    		/**
    		* The child nodes.
    		* child1为左子树和child2为右子树
    		* 左子树节点及其子节点在划分维度的对应值均小于当前节点
    		* 右子树节点及其子节点在划分维度的对应值均大于当前节点
    		*/
    		Node* child1, *child2;
        }

    FLANN中kd树的构建(超平面的确定)

    维度划分方式对查询性能的影响

    网上很多其他实现会将kd树的维度划分就保持x,y,z,…,x,y,z…的顺序。实际上,并不一定要保持这样的顺序。而且,不同的划分顺序,往往能给kd树的效率带来影响。

    因为点数与维数是固定的,kd树的深度(depth)也是固定的,如果采用均匀划分的方式,在分散程度比较小的维度上过度划分,势必会造成分散维度较大的维度划分不充分,更容易导致“长条形”区域出现。在各维度查询概率相同的情况下,长条形肯定更容易被命中,导致被查询次数过多而效率不佳。

    以上两图为例,同样的输入数据(二维),采用y-x-y的划分方式与x-x-y的方式得到的效果,可以看到,后一种划分比较充分,从而查询涉及无关分支的情况比较小,查询效率可能会比较高(具体高不高要看查询范围是不是趋向于正方形)。

    大规模数据swap的方式

    在kd树的构建中,涉及到了频繁的依据超平面进行分割的问题。因为kd树是面向高维数据的,高维数据的点的尺寸通常比较大,如果直接针对数据本身进行swap,那么会多次调用点数据本身的拷贝赋值函数,开销比较大。有没有比较合适的的方法?

    其一就是采用指针作为数组元素,而不是对象本身,这样 指针的swap就要廉价很多了。但是这样的话,相当于对输入数据就有要求了,恐怕适用性不强。

    还有一种方法就是FLANN使用的方式,比较提倡使用。具体方法是,待构建的kd树持有原始数据points,然后我们构建一个和原始数据等长的索引数组indices,indices[i]标识数据元素在原始数据中的索引值,因此调整后的p_i=points[indices[i]],当进行分割的时候,调整的只是indices的数值,开销就要小很多了,而且原始数据的顺序并没有发生改变。

     比如,上图中将原始数据{12,45,67,89,10,72,1,-9,87}进行中值排序的时候,并不需要swap原始数据中元素本身,而是swap元素在原始数据中对应的索引。

    FLANN中超平面确定方法

    FLANN在每一次划分数据partition的时候,会随机挑选出100个元素,计算各个维度的均值和方差。找到方差最大的5个维度,并随机选出1个作为“划分维度divfeat”,该维度对应的均值作为“划分维度对应值divval”。

    利用divfeat和divval,将原始数据在划分维度divfeat上划分为,并且在索引数组ind上重新排序(flann中使用快速排序)。

         *  dataset[ind[0 , ... , lim1-1]][divfeat] < divval
         *  dataset[ind[lim1 , ... , lim2-1]][divfeat] == divval
         *  dataset[ind[lim2 , ... , count]][divfeat] > divval

    对于重新排序后的ind数组,ind[0 , ... , median]索引的原始数据在divfeat维度上的值小于等于divval;ind[median , ... , count]索引的原始数据在divfeat维度上的值大于等于divval;是否等于按照下述条件进行判断:

    1. 当 lim1 > count/2 时,median = lim1;

    2. 否则当 lim2 < count2 时,median = lim2;

    3. 否则median = count/2;

    其他超平面确定方法

    最理想的超平面是中位数,而确定中位数的算法还有其他的方法。比如说,因为而中位数并不要求数据完全排序,而是要求在中位数左侧的数据不大于中位数,中位数右侧的数据不小于中位数而已。因此还可以使用std::nth_element, 其算法复杂度为O(n),系数为2左右。有技术博客表明当输入数据分布不均匀时,采用FLANN中超平面确定方法进行划分的方式很容易导致树不均衡从而使得树的深度增加,查询效率降低。

    KD树的主要算法以及FLANN(PCL)的实现分析 - Fun With GeometryFun With Geometry

    按照其描述的方法,我对flann kd树的超平面分割方法进行改进,结果表示效率提升仅2%,且匹配准确率下降。

    其他超平面确定方法,还有k-means、Hierachical Clustering等,暂时挂起,后续会继续探索。

    FLANN中kd树的查询

    最近距离表&搜索半径

    在kd树的查询过程中,会维护一个“最近距离表”,这个表中会记录k个最近节点及其距离。并且将表中的最大距离作为搜索距离

    划分距离是指:搜索点和KD树节点在划分维度上的差的绝对值。可以设想一下,如果划分距离小于搜索半径,则存在潜在最近节点,那么后续可以在相对应的子树上继续查询。

     FLANN中采用KNNSimpleResultSet来维护“最近距离表”;addPoint的实现是一个on-line的冒泡排序。worstDist()返回表中的最远距离。

    ../flann/src/cpp/flann/util/result_set.h

    template <typename DistanceType>
    class KNNSimpleResultSet : public ResultSet<DistanceType>
    {
    public:
    	typedef DistanceIndex<DistanceType> DistIndex;
    
        /**
         * Add a point to result set
         * @param dist distance to point
         * @param index index of point
         */
        void addPoint(DistanceType dist, size_t index);
        
        DistanceType worstDist() const
        {
            return worst_distance_;
        }
    private:
        size_t capacity_;
        size_t count_;
        DistanceType worst_distance_;
        std::vector<DistIndex> dist_index_;
    };

    近似最近邻查询(Approximate Nearest Neighbor Search,ANNS

    FLANN中的ANNS实现如下:

    ../flann/src/cpp/flann/algorithms/kdtree_index.h

      /**
         *  Search starting from a given node of the tree.  Based on any mismatches at
         *  higher levels, all exemplars below this level must have a distance of
         *  at least "mindistsq".
         */
        template<bool with_removed>
        void searchLevel(ResultSet<DistanceType>& result_set, const ElementType* vec, NodePtr node, DistanceType mindist, int& checkCount, int maxCheck, float epsError, Heap<BranchSt>* heap, DynamicBitset& checked) const
        {
            if (result_set.worstDist()<mindist) {
                //			printf("Ignoring branch, too far\n");
                return;
            }
    
            /* If this is a leaf node, then do check and return. */
            if ((node->child1 == NULL)&&(node->child2 == NULL)) {
                int index = node->divfeat;
                if (with_removed) {
                	if (removed_points_.test(index)) return;
                }
                /*  Do not check same node more than once when searching multiple trees. */
                if ( checked.test(index) || ((checkCount>=maxCheck)&& result_set.full()) ) return;
    
                checked.set(index);
                checkCount++;
    
                DistanceType dist = distance_(node->point, vec, veclen_);
                result_set.addPoint(dist,index);
                return;
            }
    
            /* Which child branch should be taken first? */
            ElementType val = vec[node->divfeat];
            DistanceType diff = val - node->divval;
            NodePtr bestChild = (diff < 0) ? node->child1 : node->child2;
            NodePtr otherChild = (diff < 0) ? node->child2 : node->child1;
    
            /* Create a branch record for the branch not taken.  Add distance
                of this feature boundary (we don't attempt to correct for any
                use of this feature in a parent node, which is unlikely to
                happen and would have only a small effect).  Don't bother
                adding more branches to heap after halfway point, as cost of
                adding exceeds their value.
             */
            
            //计算划分距离
            DistanceType new_distsq = mindist + distance_.accum_dist(val, node->divval, node->divfeat);
            
            //这里进行近似判断:
            //NN过滤Node的条件是 “搜索点的划分距离d” 大于 “当前搜索半径r” ,
            //ANN通过将 r除以“大于1的系数α”,使得Node 更容易被过滤,从而加快了查询速度,
            //α = 1+eps;这里epsError就是α
            //或者说,ANN通过将d*α作为“搜索点到 Node外接矩形的距离”,这样就更容易让 “搜索点的划分距离d” 大于 “当前搜索半径r” 
            //但是得到的结果也是近似的最近点。
    
            //		if (2 * checkCount < maxCheck  ||  !result.full()) {
            if ((new_distsq*epsError < result_set.worstDist())||  !result_set.full()) {
                heap->insert( BranchSt(otherChild, new_distsq) );
            }
    
            /* Call recursively to search next level down. */
            searchLevel<with_removed>(result_set, vec, bestChild, mindist, checkCount, maxCheck, epsError, heap, checked);
        }

    FLANN中kd树查询策略

    0. 构建n个kd树:根据上文提到的随机最大方差的方法,构建n个kd树,默认为4。

    1. 对于每一个搜索点到会在n个kd树中进行查询,查询过程中维护同一个“最近距离表”

    2. 采用子树堆的方式将“搜索点到kd树节点的划分距离d” 小于 "当前搜索半径r"的对应子树分支进行缓存,在遍历过当前所有kd树之后,会对子树集进行查询

    3. 设置总的节点访问数量,用于防止过度遍历子树导致的耗时问题

    FLANN中对于kd树的数量参数设置在:../flann/src/cpp/flann/algorithms/kdtree_index.h

    struct KDTreeIndexParams : public IndexParams
    {
        KDTreeIndexParams(int trees = 4,int use_nthelement = 0)
        {
            (*this)["algorithm"] = FLANN_INDEX_KDTREE;
            (*this)["trees"] = trees;//kd树的数量
        }
    };

    多线程并行查询

    flann1.6.10没有包含多线程查询方法,flann1.8.4和flann1.9.1采用多线程,对于每个搜索点都会开启一个查询线程:

    ../flann/src/cpp/flann/algorithms/nn_index.h

        /**
         * @brief Perform k-nearest neighbor search
         * @param[in] queries The query points for which to find the nearest neighbors
         * @param[out] indices The indices of the nearest neighbors found
         * @param[out] dists Distances to the nearest neighbors found
         * @param[in] knn Number of nearest neighbors to return
         * @param[in] params Search parameters
         */
        virtual int knnSearch(const Matrix<ElementType>& queries,
        		Matrix<size_t>& indices,
        		Matrix<DistanceType>& dists,
        		size_t knn,
        		const SearchParams& params) const
        {
        	assert(queries.cols == veclen());
        	assert(indices.rows >= queries.rows);
        	assert(dists.rows >= queries.rows);
        	assert(indices.cols >= knn);
        	assert(dists.cols >= knn);
        	bool use_heap;
    
        	if (params.use_heap==FLANN_Undefined) {
        		use_heap = (knn>KNN_HEAP_THRESHOLD)?true:false;
        	}
        	else {
        		use_heap = (params.use_heap==FLANN_True)?true:false;
        	}
        	int count = 0;
    
        	if (use_heap) {
    #pragma omp parallel num_threads(params.cores)
        		{
        			KNNResultSet2<DistanceType> resultSet(knn);
    #pragma omp for schedule(static) reduction(+:count)
        			for (int i = 0; i < (int)queries.rows; i++) {
        				resultSet.clear();
        				findNeighbors(resultSet, queries[i], params);
        				size_t n = std::min(resultSet.size(), knn);
        				resultSet.copy(indices[i], dists[i], n, params.sorted);
        				indices_to_ids(indices[i], indices[i], n);
        				count += n;
        			}
        		}
        	}
        	else {
    #pragma omp parallel num_threads(params.cores)
        		{
        			KNNSimpleResultSet<DistanceType> resultSet(knn);
    #pragma omp for schedule(static) reduction(+:count)
        			for (int i = 0; i < (int)queries.rows; i++) {
        				resultSet.clear();
        				findNeighbors(resultSet, queries[i], params);
        				size_t n = std::min(resultSet.size(), knn);
        				resultSet.copy(indices[i], dists[i], n, params.sorted);
        				indices_to_ids(indices[i], indices[i], n);
        				count += n;
        			}
        		}
        	}
        	return count;
        }

    FLANN中利用SearchParams对于查询过程的参数进行管理

    通过对参数cores进行设置线程数

    ../flann/src/cpp/flann/util/params.h

    struct SearchParams
    {
        SearchParams(int checks_ = 32, float eps_ = 0.0, bool sorted_ = true ) :
        	checks(checks_), eps(eps_), sorted(sorted_)
        {
        	max_neighbors = -1;
        	use_heap = FLANN_Undefined;
        	cores = 1;
        	matrices_in_gpu_ram = false;
        }
    
        // how many leafs to visit when searching for neighbours (-1 for unlimited)
        int checks;
        // search for eps-approximate neighbours (default: 0)
        float eps;
        // only for radius search, require neighbours sorted by distance (default: true)
        bool sorted;
        // maximum number of neighbors radius search should return (-1 for unlimited)
        int max_neighbors;
        // use a heap to manage the result set (default: FLANN_Undefined)
        tri_type use_heap;
        // how many cores to assign to the search (used only if compiled with OpenMP capable compiler) (0 for auto)
        int cores;
        // for GPU search indicates if matrices are already in GPU ram
        bool matrices_in_gpu_ram;
    
    };

    小结:

    对比多种索引方法,综合考虑效率和匹配准确率,优先选择kd树。(后续可以研究一下基于“结构化graph”的索引方法,据说效率更高,但是flann没有相关实现)

    展开全文
  • 文章目录算法简介python代码实现kd树构建与搜索 算法简介 k近邻法属于监督学习,不需要训练模型(懒惰学习)。算法流程:对于测试样本,按照某种距离度量(闵可夫斯基距离、欧氏距离、曼哈顿距离、切比雪夫距离等...
  • java KD树构建

    2016-01-04 19:57:10
    k-d (k-dimensional的简称),是一种分割k维数据空间的数据结构。主要应用于多维空间关键数据的搜索(如:范围搜索和最近邻搜索)。K-D是二进制空间分割的特殊的情况。
  • K-NN之KD树构建

    2021-09-24 21:29:17
    最近在看机器学习基础,差一点就放弃写了,给我写吐了。...# 构造KD树 # 教学视频: https://www.bilibili.com/video/BV1JE41147KA?p=2&spm_id_from=pageDriver import math class Node: def
  • Kd树构造

    2022-04-04 19:39:42
    Kd树–总结 1. KdTree构造 以经典的2D样本点集合为例data = [(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)]来展示构造过程。 排序 ​ 从0维度开始,这里就是x轴的维度,以x维度对样本集合点排序。 ​ data = [(2,3)(4,7)...
  • kd树建立 MATLAB源代码

    热门讨论 2012-05-05 13:31:15
    用于建立kd树数据结构,对于做图像处理的人应该很有用
  • 通过在REFERENCE中从数据点构建kD树,并在MODEL中为每个数据点查询树,以有效的方式执行搜索。 PTS = KDRANGEQUERY( ROOT, QUERYPT, DISTLIM ) 查找存储在 kD 树 ROOT 中且在 QUERYPT 的 DISTLIM 单元内的所有点。...
  • KNN的实现——KD树

    2020-10-25 19:59:42
    KD树就是K维二叉树 我们学过的各种二叉树,使用一个键值在树中导航以执行必要的操作,二叉树中每个节点都有唯一的一个key值,通过key我们可以组织二叉查找树、平衡树、自适应树、堆等,从某种意义上讲,这是一维的...
  • kd树构建与搜索

    千次阅读 多人点赞 2016-03-15 09:13:10
    二、构建kd树之后,如今进行最近邻搜索呢? KD树的查找算法: 在k-d树中进行数据的查找也是特征匹配的重要环节,其目的是检索在k-d树中与查询点距离最近的数据点。 这里先以一个简单的实例来描述最邻近查找的...
  • 【KNN以及KD树

    2022-07-16 10:27:37
    KNN的原理以及KD树
  • kd树的构造和搜索(超详细)

    千次阅读 多人点赞 2021-03-19 20:55:29
    kd(k-dimensional)树解决...kd树类似于平衡二叉树,不过是对k维数据进行划分。 Kd树的构造: 假设特征空间T: 递归构造kd树: 输入:n个对象的特征向量组成的特征空间T={x0,x1,…,xNx_0,x_1,…,x_Nx0​,x1​,…,xN​}
  • kd树学习体会

    2019-03-30 14:42:49
    在构造1维BST时,一个1维数据依据其与的根结点和中间结点进行大小比較的结果来决定是划分到左子树还是右子。同理。我们也能够依照这种方式,将一个K维数据与Kd-tree的根结点和中间结点进行比較。仅仅只是不是...
  • 用python实现kd树构建和搜索

    千次阅读 2022-03-31 16:26:54
    前两天学习了knn算法,knn的思想很简单,不过其中提出的kd树有理解的必要。故就用python写了一个kd树代码。 个人感想是,把kd树算法实现一遍比看书看半天有用多了,而且还不会犯困(bushi 思路来自...
  • Matlab构建KD树

    2022-03-07 08:02:02
    机器学习、构建kd树
  • kd树(K-dimensional tree)

    千次阅读 2019-10-20 20:23:57
    kd-tree简介 Kd-Tree(K-dimensional tree),是一种高维索引...kd树结构类似于高维的二叉树,树中存储的是一些K维数据。在一个K维数据集合上构建一棵Kd-Tree代表了对该K维数据集合构成的K维空间的一个划分。 kd-t...
  • 5、kd树构建及查询实现 import numpy as np class TreeNode: def __init__(self, s, d): self.vec = s # 特征向量 self.Dimension = d # 即划分空间时的特征维度,这里选取方差最大的维度 self.left = None # 左子...
  • 构建kd树kd树的搜索

    2021-02-25 11:23:46
    方法2、kd树使用特殊的结构存储训练数据,以减小计算距离的次数。算法复杂度可降低到O(DNlog(N)) 基本原理 如果A和B距离很远,B和C距离很近,那么A和C的距离也很远。 构造kd树 1.构造根节点 2.通过递归的...
  • kd树java实现

    2018-09-05 10:50:20
    利用Java语言实现kd树。利用快速排序找到每一维度的中位数构建Kd树。搜索Kd树时候利用K大堆维护K个最近邻的点。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,731
精华内容 1,892
关键字:

kd树的构建