使用while loop 构建BVH，遍历BVH （1）

> 写在前面的一些话

其实我早就想写这篇博客了，尤其是前段时间需要完成一个跟BVH打交道的算法，当时在网上搜索相关的 BVH 构建，遍历代码，全都是用 recursion syntax 写的。

的确， recursion syntax 看起来很“简洁” 和“高效”。但是，我在写程序的时候，花最多的时间是用来除错的。写起来简洁高效的代码除错的时候一点都不简洁高效。相反如果用 while loop 来写的话，我除错是容易多了。

所以在接下来的博客里我会介绍怎么用 while loop 构建 BVH 和遍历BVH。

> 下面正式开启写代码

我们先从构建BVH 开始。假设我们有一个用三角形顶点数组表示的 3D model，我们给这个 3D model 构建BVH。在我们的BVH中，每个三角形都有一个包围盒（aabb）将其包住，而 BVH 的每个叶子节点（leaf node）都指向一个aabb，且这个aabb里面只有一个三角形。

我们首先创建一个包围盒类 - Myaabb：

class Myaabb
{
public:
    int idx;  // triangle index
    vec3 nn;  // the min point of the aabb
    vec3 mm;  // the max point of the aabb
    
    Myaabb(){ idx = -1; } 

    void creat_box(vec3 vmax, vec3 vmin, int _idx)
    {
        mm = vmax;
        nn = vmin;
        idx = _idx;
    }

    // just for check the box information
    void print_box()
    {
        printf("max: %f %f %f\n", mm.x, mm.y, mm.z);
        printf("min: %f %f %f\n", nn.x, nn.y, nn.z);
        printf("tri idx: %d\n", idx);
    }

};

其中 idx 用于判断存在aabb里的三角形， nn 和 mm 分别是 aabb 最小的顶点和最大的顶点，并且是正好可以包住存在aabb 的三角形的。

函数 creat_box() 根据输入的 vmax， vmin 和 _idx 创建 aabb。

函数 print_box() 会输出aabb 的信息，主要是除错用的。

有了aabb，我们就有了构建BVH的基本元素，接下来我们创建 BVH 的node类，用于指向，如下：

class MyNode
{
public:
    MyNode *left, *right;
    Myaabb box;
    std::vector<int> idx;

    MyNode(){ memset(this, 0, sizeof(*this)); }
}

这个node 类包含了一个aabb类元素 - box，用于指明当前node 指向的 bounding box；还包含了两个 node 类型指针 - *left 和 * right，分别指向当前node 的左子树和右子树。 

如果node 为叶子节点，left 和 right 都为null， box 只包裹一个三角形。如果node 为 中间节点 或者根节点的话，left 和 right 分别指向该节点的左子节点和右子节点，box 则为包着左节点和右子节点的大包围盒。

有了aabb 和 node，我们就可以开始构建BVH了。构建过程主要包括两个步骤，首先把 3D model 的每一个三角形都创建一个包围盒，并生成一个包围盒数组，代码如下：

MyNode* build_sceneBVH( const std::vector<vec3> &vbuf, const std::vector<uint32_t> &ibuf)
{
  std::vector<Myaabb> src;  
  int i;
  int n_face = ibuf.size()/3; //the number of triangles 

  for(i=0; i<n_face; i++)
  {
    Myaabb box;
    int tidx = i;
    vec3 v0  = vbuf[ibuf[3*tidx+0]];
    vec3 v1  = vbuf[ibuf[3*tidx+1]];
    vec3 v2  = vbuf[ibuf[3*tidx+2]];

    vec3 mm = vmax(vmax(v0, v1), v2);
    vec3 nn = vmin(vmin(v0, v1), v2);

    box.creat_box(mm, nn, tidx);
    src.push_back(box);
  }
 
  return process_aabbs(src);
}

vbuf 数组存的是 3D model 三角形的顶点，ibuf 数组存的是三角形顶点对应的index。std::vector<Myaabb> src 是保存三角形包围盒的数组，我们接着就将它传进去 process_aabbs() 函数来生成BVH，代码如下：

MyNode* process_aabbs(vector<Myaabb> &src)
{
  int i, j;
  MyNode* node = new MySceneBVH; // root node of the generated BVH

  // if the 3D model only has one triangle
  if(src.size()==1) 
  {
    node->box = src[0];
    node->idx.push_back(src[0].idx);
    return node;
  }else
  {
    std::vector<Myaabb> current_src;
    MyNode *current_node = node;

    std::vector<Myaabb> mynote_srca[512];
    std::vector<Myaabb> mynote_srcb[512];
    MyNode *mynote_left[512];
    MyNode *mynote_right[512];


    int nsrca = 0, nsrcb=0;
    int nlnode = 0, nrnode=0;
    mynote_left[nlnode++] = node;
    mynote_srca[nsrca++] = src;
    
    // while loop generate BVH
    while( nsrca>0 || nsrcb>0 )
    {
      if(nsrca>0)
      {
        current_src = mynote_srca[nsrca-1];
        current_node = mynote_left[nlnode-1];
        nsrca--;
        nlnode--;
      }else if(nsrcb>0)
      {
        current_src = mynote_srcb[nsrcb-1];
        current_node = mynote_right[nrnode-1];
        nsrcb--;
        nrnode--;
      }


      vec3 mm = vec3(-FLT_MAX, -FLT_MAX, -FLT_MAX);
      vec3 nn = vec3( FLT_MAX,  FLT_MAX,  FLT_MAX);
      for(i=0;i<current_src.size(); i++)
      {
        mm = vmax(mm, current_src[i].mm);
        nn = vmin(nn, current_src[i].nn);
      }

      if(current_src.size()==1)
      {
        current_node->box = current_src[0];
        current_node->idx.push_back(current_src[0].idx);
        continue;
      }

      // divide to srca and srcb
      float b;
      int dim[3];
      vec3 a = (mm+nn)/2.f; //center
      vec3 r = (mm-nn)/2.f;
      {
        if(r.x > r.y && r.x > r.z)
        {
          dim[0] = 0;
          dim[1] = 1;
          dim[2] = 2;
          if( r.y < r.z )
            g_swap(dim[1],dim[2]);
        }else if( r.y > r.z )
        {
          dim[0] = 1;
          dim[1] = 2;
          dim[2] = 0;
          if( r.z < r.x )
            g_swap(dim[1],dim[2]);
        }else
        {
          dim[0] = 2;
          dim[1] = 0;
          dim[2] = 1;
          if( r.x < r.y )
            g_swap(dim[1],dim[2]);
        }
      }

      std::vector<Myaabb> srca, srcb;
      for( j=0; j<3; j++ )
      {
        srca.clear();
        srcb.clear();
        for(i=0; i<current_src.size(); i++)
        {
          vec3 center = (current_src[i].mm + current_src[i].nn)/2.f;
          if(  ((float*)&center)[dim[j]] < ((float*)&a)[dim[j]]  )
            srca.push_back(current_src[i]);
          else
            srcb.push_back(current_src[i]);
        }
        if( srca.size() && srcb.size() )
          break;
      }
      if(srca.size()==0 || srcb.size()==0)
      {
        srca.assign(current_src.begin(), current_src.begin()+current_src.size()/2);
        srcb.assign(current_src.begin()+current_src.size()/2, current_src.end());
      }

      current_node->box.creat_box(mm, nn, 0);

      if(srca.size()>0)
      {
        current_node->left = new MyNode;
        mynote_left[nlnode++] = current_node->left;
        mynote_srca[nsrca++] = srca;
      }
      if(srcb.size()>0)
      {
        current_node->right = new MyNode;
        mynote_srcb[nsrcb++] = srcb;
        mynote_right[nrnode++] = current_node->right;
      }
    };
    return node;
  }
}

   具体使用这个函数的代码例子如下：

    std::vector<vec3> vbuf = the 3D model vertices;
    std::vector<uint> ibuf = the indices of all vertices;

    MyNode *mysbvh;
      mysbvh = build_sceneBVH(vbuf, ibuf);

> 函数解释

process_aabbs() 这个函数传进来的是 3D model 的包围盒数组， 即src，传出去的是 BVH 的根节点， 即node。

之前就说过，我们 BVH 的叶子节点指向只包含一个三角形的aabb。为了达到这个目的，我们需要将 src 分成左包围盒数组 srca1 和右包围盒数组 srcb1，然后 srca1 又会分为左包围盒数组 srca2 和右包围盒 srcb2， srcb1 会分为左包围盒数组 srca3 和右包围盒 srcb3， 直至分到左右各只有一个包围盒。 我们一次loop 只可以分一个包围盒数组，然后产生两个待下次分割的包围盒数组。那我们下次要分哪个包围盒呢?

我们为此准备了两个工作列表，左右各一个，将每次loop 分割的包围盒数组都分别放到各自的工作列表中，然后我们选择一直分左边的包围盒数组，即下次loop分的包围盒数组是左工作列表里最后一个包围盒数组，直至左工作列表里没有了，我们再分记录在右工作列表的最后一个包围盒数组，这里要注意的是，这个数组可能又会分成左包围盒数组和右包围盒数组，那么下次loop开始我们又会继续拿左工作列表里的数组进行分割，而不是取右工作列表的数组进行分割。当左工作列表和右工作列表都没有东西了，就证明工作完成了，while loop 结束了，BVH 生成了，收工！ 

具体到函数执行呢，在第一次 loop 时，我们先把 3D model 所有包围盒生成一个最大包围盒，然后我们比较每一个小包围盒与当前最大包围盒的中心，将小包围盒记录到左节点的包围盒数组（srca）或者右结点包围盒数组（srcb)， 在 loop 最后，我们让 node 指向最大包围盒，给node 生成左子节点和右子节点，并记录到对应的子节点数组中（mynote_left 和 mynote_right）, srca 和 srcb 也会记录到对应的包围盒数组的数组（mynote_srca和 mynote_right， 也就是我们的左右工作列表）。在下一次loop开始，我们会从子节点数组中拿出一个 node （current_node），再从对应的包围盒数组的数组中拿出一个包围盒数组(current_src)，重复上面的操作，while loop 结束后，我们就得到了想要的BVH。

> 参数解释

current_src 是我们每一次 loop 处理的包围盒数组，在 loop 最开始的时候它是 3D model 所有包围盒的数组；在之后的 loop 里，它被设定为之前 loop 划分的左包围盒数组的数组（mynote_srca，如果里面有value的话）里的最后一个数组；如果 mynote_srca size 为0， current_src 就会被设定为右包围盒数组的数组（mynote_srcb）里的最后一个数组。

current_node 是我们每一次 loop 处理的节点，在loop最开始的时候它被设置为root node；在之后的 loop 中，它会先从mynote_left 中查找是否有可以处理的node （即继续处理左子树的节点），不然的话就处理 右子树 mynote_right 中的node。

mynote_srca 和 mynote_srcb 分别是记录 current_node  左子树和右子树下需要处理的aabb 数组, 我们每次开始都先从处理 mynote_srca 下的 aabb 开启，在函数最开始的时候，mynote_srca 存的是所有包围盒的数组，即 src；在每次while loop 结束后，mynote_srca 和 mynote_srcb 会分别记录下此次分出来的左子树 和右子树的包围盒数组 （srca 和 srcb），供之后的while loop 处理。

nsrca 和 nsrcb 分别记录mynote_srca 和 mynote_srcb 里还有多少项要处理，也是我们终止while loop 的判断条件。 nlnode 和 nrnode 功能类似于 nsrca 和 nsrcb。

mm 和 nn 是current_node 所指向的 aabb 的最大和最小顶点。
dim 这个数组是我把单个包围盒分左右的一个参考依据，但具体为什么这么做我忘了，之后补上。。。

> 总结

写到这里，用 while loop 构建 BVH 的内容基本上就完结了。总结一下呢，我这个构建方法是自上而下构建的，效率可能也没有彻底优化，而且代码可能看起来觉得好长；但是用上 while loop，有错的时候我们可以很快速找到是哪里出了错， 在我看来，啰嗦但能让我快速除错的代码比让我看起来简洁却不好除错的代码更有用。

这一篇文章我写的好啰嗦，但有这些啰嗦的内容，下一篇我们介绍如何用 while loop 来遍历BVH 就会简单很多，下篇见！

> vmax, vmin 和 g_swap function

vec3 vmax(const vec3 &a, const vec3 &b)
{
  vec3 c;
    c.x = (a.x > b.x) ? a.x : b.x;
    c.y = (a.y > b.y) ? a.y : b.y;
    c.z = (a.z > b.z) ? a.z : b.z;
  return c;
}

vec3 vmin(const vec3 &a, const vec3 &b)
{
  vec3 c;
    c.x = (a.x < b.x) ? a.x : b.x;
    c.y = (a.y < b.y) ? a.y : b.y;
    c.z = (a.z < b.z) ? a.z : b.z;
  return c;
}

void g_swap( int &a, int &b )
{ 
  int c;  
    c = a;  
    a = b;  
    b = c; 
}