flann算法中的kd树构建代码解析

构建的总体函数
    void buildIndexImpl()
    {
        // Create a permutable array of indices to the input vectors.
        vind_.resize(size_);  //size是点的个数参数   vind_为存放所有点索引的std::vector<int>变量 (kdtree_single_index.h 中)
        for (size_t i = 0; i < size_; i++) {
            vind_[i] = i;  //vind_向量初始化
        }
        computeBoundingBox(root_bbox_);  //根部包围盒root_bbox_ 是一个vector vector中 每个元素是struct struct中的两个参数为 T类型的 low 和high 两个变量
        root_node_ = divideTree(0, size_, root_bbox_ );   // construct the tree
        if (reorder_) {
            data_ = flann::Matrix<ElementType>(new ElementType[size_*veclen_],  size_, veclen_);
            for (size_t i=0; i<size_; ++i) {
                std::copy(points_[vind_[i]], points_[vind_[i]]+veclen_, data_[i]);
            }
        }
    }
重点函数computeBoundingBox(root_bbox_);分析:
    该函数是计算根部包围盒的
       我的理解 veclen_应该是要构建成kd树的数据的维数  比如三维点云数据  该参数应该是3  下面以点云数据为例 分析代码
    points_存储所有点云的数据   例如points_[0][0]表示点云数据中第一个点的x坐标数据,points_[0][1]表示点云数据中第一个点的y坐标数据,points_[0][2]表示点云数据中第一个点的z坐标数据
    bbox[i].low和bbox[i].high 表示点云最外层包围盒的第i维的最小值和最大值
    void computeBoundingBox(BoundingBox& bbox)
    {
        bbox.resize(veclen_);
        //首先对vector<struct>bbox  进行初始化
        for (size_t i=0; i<veclen_; ++i) {
            bbox[i].low = (DistanceType)points_[0][i];
            bbox[i].high = (DistanceType)points_[0][i];
        }
        //循环整个点云数据找到 能够完整包围整个点云的最小的外层包围长方体,每一个维度包含一个最大值一个最小值
        for (size_t k=1; k<size_; ++k) {
            for (size_t i=0; i<veclen_; ++i) {
                if (points_[k][i]<bbox[i].low) bbox[i].low =  (DistanceType)points_[k][i];  //k为点云中点的索引值  i为该点坐标维度索引
                if (points_[k][i]>bbox[i].high) bbox[i].high =  (DistanceType)points_[k][i];
            }
        }
    }
root_bbox_参数分析:
kdtree_single_index.h 中
678行
    /**
     * Root bounding box根部包围盒
     */
    BoundingBox root_bbox_;

334行
typedef std::vector<Interval> BoundingBox; 

320行
struct Interval
    {
        DistanceType low, high;
        
    private:
        template <typename Archive>
        void serialize(Archive& ar)
        {
            ar & low;
            ar & high;
        }
        friend struct serialization::access;
    };

flann.hpp
80行
    typedef typename Distance::ResultType DistanceType;

dist.h
79行
    typedef typename Accumulator<T>::Type ResultType;
50-63行
template<typename T>
struct Accumulator { typedef T Type; };
template<>
struct Accumulator<unsigned char>  { typedef float Type; };
template<>
struct Accumulator<unsigned short> { typedef float Type; };
template<>
struct Accumulator<unsigned int> { typedef float Type; };
template<>
struct Accumulator<char>   { typedef float Type; };
template<>
struct Accumulator<short>  { typedef float Type; };
template<>
struct Accumulator<int> { typedef float Type; };

重点函数root_node_ = divideTree(0, size_, root_bbox_ );   // construct the tree分析:
kdtree_single_index.h 中
673行: 
   /**
     * Array of k-d trees used to find neighbours.
     */
    NodePtr root_node_;
317行:
    typedef Node* NodePtr;
262行:
    struct Node
    {
       /**
        * Indices of points in leaf node  叶子节点中点的索引
        */
       int left, right;
       /**
        * Dimension used for subdivision.  用于划分的尺寸
        */
       int divfeat;
       /**
        * The values used for subdivision. 用于划分的值
        */
       DistanceType divlow, divhigh;
        /**
         * The child nodes.
         */
        Node* child1, * child2;
        
        ~Node()
        {
            if (child1) child1->~Node();
            if (child2) child2->~Node();
        }
    private:
        template<typename Archive>
        void serialize(Archive& ar)
        {
            typedef KDTreeSingleIndex<Distance> Index;
            Index* obj = static_cast<Index*>(ar.getObject());
            ar & left;
            ar & right;
            ar & divfeat;
            ar & divlow;
            ar & divhigh;
            bool leaf_node = false;
            if (Archive::is_saving::value) {
                leaf_node = ((child1==NULL) && (child2==NULL));
            }
            ar & leaf_node;
            if (!leaf_node) {
                if (Archive::is_loading::value) {
                    child1 = new(obj->pool_) Node();
                    child2 = new(obj->pool_) Node();
                }
                ar & *child1;
                ar & *child2;
            }
        }
        friend struct serialization::access;
    };
    /**创建一个树状节点,将vecs的列表从vind[first]细分到vind[last]。 该例程在每个子列表上都被递归调用。在pTree位置放置一个指向这个新树节点的指针。
     * Create a tree node that subdivides the list of vecs from vind[first] 
     * to vind[last].  The routine is called recursively on each sublist.
     * Place a pointer to this new tree node in the location pTree.
     *
     * Params: pTree = the new node to create
     *                  first = index of the first vector
     *                  last = index of the last vector
     */
    NodePtr divideTree(int left, int right, BoundingBox& bbox)
    {
        NodePtr node = new (pool_) Node(); // allocate memory
        /* If too few exemplars remain, then make this a leaf node. */
                //当前待分配的区域的点数的大小小于叶节点可容纳的最大点数时,认为当前节点为叶节点,
                //注意此处的left是和right 是索引vector(vind_)中点云索引的索引值(这个索引vector中放置的是原始点云中点云数据在点云中的索引值)points_为原始点云数据  
                // 则 points_[vind_[left]][i] 表示points_中第n个点的第i维的值   这个点是索引值n  这个索引值为vind_向量中第left+1个值的内容
        if ( (right-left) <= leaf_max_size_) {
            //因为是个叶节点 ,所以该叶节点的左孩子和有孩子的树 都是空的  并且这个叶节点的取值范围是left~right
            //然后这个叶节点及所有的他的子节点 所在区域的包围盒 是两个for循环计算得到的。 
            node->child1 = node->child2 = NULL;    /* Mark as leaf node. */
            node->left = left;
            node->right = right;
            // compute bounding-box of leaf points
            for (size_t i=0; i<veclen_; ++i) {
                bbox[i].low = (DistanceType)points_[vind_[left]][i];
                bbox[i].high = (DistanceType)points_[vind_[left]][i];
            }
            for (int k=left+1; k<right; ++k) {
                for (size_t i=0; i<veclen_; ++i) {
                    if (bbox[i].low>points_[vind_[k]][i])  bbox[i].low=(DistanceType)points_[vind_[k]][i];
                    if (bbox[i].high<points_[vind_[k]][i])  bbox[i].high=(DistanceType)points_[vind_[k]][i];
                }
            }
        }
        else {
            int idx;
            int cutfeat;
            DistanceType cutval;
            middleSplit(&vind_[0]+left, right-left, idx, cutfeat, cutval, bbox);
            node->divfeat = cutfeat;
            BoundingBox left_bbox(bbox);
            left_bbox[cutfeat].high = cutval;
            node->child1 = divideTree(left, left+idx, left_bbox);
            BoundingBox right_bbox(bbox);
            right_bbox[cutfeat].low = cutval;
            node->child2 = divideTree(left+idx, right, right_bbox);
            node->divlow = left_bbox[cutfeat].high;
            node->divhigh = right_bbox[cutfeat].low;
            for (size_t i=0; i<veclen_; ++i) {
               bbox[i].low = std::min(left_bbox[i].low, right_bbox[i].low);
               bbox[i].high = std::max(left_bbox[i].high, right_bbox[i].high);
            }
        }
        return node;
    }
pool_
kdtree_single_index.h 中
690行:
/** 池化内存分配器。当有大量小的内存分配时,使用池式内存分配器比直接分配内存更有效率。
     * Pooled memory allocator.
     *
     * Using a pooled memory allocator is more efficient
     * than allocating memory directly when there is a large
     * number small of memory allocations.
     */
    PooledAllocator pool_;
PooledAllocator类的具体实现与解释在 allocator.h中
NodePtr node = new (pool_) Node(); // allocate memory 
上一句是定位new重载的一个用法 (X)
不是定位new 是operator new 的一种带参重载 带的参数是PooledAllocator类对象而已
PooledAllocator类实现的功能有些类似定位new
其实new的全写是:new (expression-list)_optional Type (constructor arguments)_optional
pMemoryBlock = new (m_uiGrowNum, m_uiUnitSize) MemoryBlock(m_uiGrowNum, m_uiUnitSize);
就相当于调用 operator new(sizeof(tagMemoryBlock), m_uiGrowNum, m_uiUnitSize)
然后再调用constructor tagMemoryBlock(m_uiGrowNum, m_uiUnitSize)
具体功能及实现:
在allocator.h中的 197行
inline void* operator new (std::size_t size, flann::PooledAllocator& allocator)
{
    return allocator.allocateMemory(size) ;
}
void middleSplit(int* ind, int count, int& index, int& cutfeat, DistanceType&  cutval, const BoundingBox& bbox)
    {
        // find the largest span from the approximate bounding box
                //1.从近似的边界盒中找出最大的跨度。 
        ElementType max_span = bbox[0].high-bbox[0].low;
        cutfeat = 0;
        cutval = (bbox[0].high+bbox[0].low)/2;
        for (size_t i=1; i<veclen_; ++i) {
            ElementType span = bbox[i].high-bbox[i].low;
            if (span>max_span) {
                max_span = span;
                cutfeat = i;
                cutval = (bbox[i].high+bbox[i].low)/2;
            }
        }
        // compute exact span on the found dimension
        //2.在找到的维度上计算精确的跨度
        ElementType min_elem, max_elem;
        computeMinMax(ind, count, cutfeat, min_elem, max_elem);
        cutval = (min_elem+max_elem)/2;
        max_span = max_elem - min_elem;
        // check if a dimension of a largest span exists
            //3.检查是否存在一个最大跨度的尺寸
        size_t k = cutfeat;
        for (size_t i=0; i<veclen_; ++i) {
            if (i==k) continue;
            ElementType span = bbox[i].high-bbox[i].low;
            if (span>max_span) {
                computeMinMax(ind, count, i, min_elem, max_elem);
                span = max_elem - min_elem;
                if (span>max_span) {
                    max_span = span;
                    cutfeat = i;
                    cutval = (min_elem+max_elem)/2;
                }
            }
        }
        //以上三个步骤不断检测是否找到的维度是否正确,以及精确的跨度。
                //上述步骤解算出来的值为cutfeat为找到的哪个维度所包含的点在该维度上跨度最大。
              //min_elem表示该维度上的最小值
                //max_elem表示该维度上的最大值
        //cutval表示该维度上的最大值和最小值的平均值
        int lim1, lim2;
        planeSplit(ind, count, cutfeat, cutval, lim1, lim2);
        if (lim1>count/2) index = lim1;
        else if (lim2<count/2) index = lim2;
        else index = count/2;
        
        assert(index > 0 && index < count);
    }
//计算一堆点的数据,在dim这个维度上的最大值和最小值
//这一堆点数据 数量为count个,在原始点云中的索引存在了ind所指向空间向后数count个空间内(类型为int)
void computeMinMax(int* ind, int count, int dim, ElementType& min_elem,  ElementType& max_elem)
    {
        min_elem = points_[ind[0]][dim];
        max_elem = points_[ind[0]][dim];
        for (int i=1; i<count; ++i) {
            ElementType val = points_[ind[i]][dim];
            if (val<min_elem) min_elem = val;
            if (val>max_elem) max_elem = val;
        }
    }
/**在'cutval'位置用垂直于'cutfeat'尺寸的轴线对点列表进行细分。
//列表细分 调整索引位置的主要调整函数 其实就是找中位数的函数
重要  要仔细的读清楚
     *  Subdivide the list of points by a plane perpendicular on axe corresponding
     *  to the 'cutfeat' dimension at 'cutval' position.
     *  
     *  On return:
     *  dataset[ind[0..lim1-1]][cutfeat]<cutval
     *  dataset[ind[lim1..lim2-1]][cutfeat]==cutval
     *  dataset[ind[lim2..count]][cutfeat]>cutval
     */函数就是将索引列表中该包围盒中的索引调整至以cutval为分界线 左侧索引的点在该维的数据小于cutval 右侧索引的点在该维的数据大于cutval    该函数功能做了一个循环调整位置 
    void planeSplit(int* ind, int count, int cutfeat, DistanceType cutval, int&  lim1, int& lim2)
    {
        int left = 0;
        int right = count-1;
        for (;; ) {
            while (left<=right && points_[ind[left]][cutfeat]<cutval) ++left;
            while (left<=right && points_[ind[right]][cutfeat]>=cutval) --right;
            if (left>right) break;
            std::swap(ind[left], ind[right]); ++left; --right;
        }
        lim1 = left;
        right = count-1;
        for (;; ) {
            while (left<=right && points_[ind[left]][cutfeat]<=cutval) ++left;
            while (left<=right && points_[ind[right]][cutfeat]>cutval) --right;
            if (left>right) break;
            std::swap(ind[left], ind[right]); ++left; --right;
        }
        lim2 = left;
    }
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值