Lin-Canny算法

求解凸包间的最近点对是计算几何中一个非常有用的算法，经常被用在诸如碰撞检测、物理引擎等图形学相关的领域，而且该算法的效率对于最终整个系统的效能有着相当关键的制约。常规的对于几何体之间的最近点对求解一般即是暴力的遍历算法，这种效率太过于低下，不具有普遍的适用意义，尤其是在一个较大的碰撞系统里边。因而，Lin以及Canny等人提出了渐进式求解的Lin-Canny算法。

        Lin-Canny（简称LC）算法的主要用途即是求解两个凸包几何体间的最近距离（或最近点对等）。它的实现基础首先是将几何体元进行分类，即：点、边（三维线段）、面（有限面片），进而将两个凸包中的几何体元组成基本的几何体元对：点-点、点-边、点-面、边-边、边-面、面-面。如此一来，对于两个几何体之间的最近点对可能出现的位置即归结到了上述六种体元对之间。

        LC的高效实现中，另外一个关键的概念即是三维Voronoi区域（3D Voronoi Area， VA），它是二维Voronoi区域的升维扩展。对于一种由基于几何体元组成的体元对A和B来说，A 的Voronoi区域即是指，在此区域内定可以获得B上到A的最近点对，而在此区域则不一定能得到这样的最近点对。Voronoi区域是LC算法相对于暴力遍历算法进行效率提升的一个重要基础。一般来说，对于凸包上的点所对应的VA就是由它所关联到的若干条边所对应的在该点上的垂面所包围而成，VA的方向为凸包在该点上由内向外的方向。对于每条边所对应的VA则是由四个平面包围而成，其中两个平面为该边所关联的两个面片在该条边处所对应的垂面，另外两个则是由该边的过其两个端点的垂面。对于组成凸包的面片所对应的VA则是由垂直于该面片且过组成该面片的几条边的平面包围而成。对于三维空间中的三种基本体元：点、边、面所对应的各自的Voronoi区域可见下图所示。

                                      点和边的VA区域示意

面片的VA区域示意

    得到了凸包中的每个基本几何体元所对应的VA之后，即可进行基于VA的LC算法。首先，来了解一下基于VA的基本体元对之间的最近点对的求解方法。

1.       点-点

对于来自不同几何体的两个顶点，判断其是否是最近点对。如果这两个点组成了最近点对，那么它们之间必定满足: A处于B所对应的VA域内；B处于A所对应的VA域内。

2.       点-边

对于由点和边组成的体元对，首先求解出点到该边（即线段）间的最近点对，接着判断该最近点对是否是原始凸包的最近点对。假设点A到边E上的最近点为P，则判断标准为：A处于E所对应的VA域内；P处于A所对应的VA域内。如果条件满足，则可以得到由边点计算而来的最近点对（A,P）。

3.       点-面

同样，此种情况需要先求解出点到面片上的最近点对，接着进行VA检验。假设A到面片F上的最近点为P：A处于F所对应的VA域内；P处于A所对应的VA域内。

有了上述三种基于体元对间的VA判断以后，基本上相当于已经实现了LC算法，对于其它三种基本体元对的判断，均是做相应的处理之后转为上述三种情况中的一种，并接着进行算法的判断操作，直至得到最终的解。

4.       边-边

首先，求解出边-边（线段-线段）间的最近点对，这里要保证得到的最近点对均处于线段之上。进而分别以各个点与另外一条边为体元对，做2情况的判断。

5.       边-面

同样，求解出边-面间的最近点对，类同于情况4归结为前3种情况中的一种后继续进行判断。

6.       面-面

最后，对于面-面间的情况。一般情况下来说，对于两个凸包间的最近点对分别出现于个集合的面片上的情况较少，因而在此处的判断方法同样是做类型变换。首先，需要对特殊情况进行一些处理，判断两个面片是否平行，如果平行，则判断其间中否有重叠的情况，这里的重叠是指在平面的法向量上观察两面片是否有相交，若有，则重叠部分中任取一对点即是相应的最近点对。若两个面片不平行或是平行而无重叠，则求出两个面片之间的最近边对来组成新的几何体元对接着进行之后的检测。

        LC算法的实现主要以上述几种基本体元之间的情形做为基本情形进而进行渐进式的近点对求解，而且每两次基本体元对之间的转换必定使得求得到有效距离变得更小，直到求出最终的最近点对时得到最小的凸包间距离。这种实现方法也正是LC算法相对于传统暴力算法的效率提高优势所在。

        实现主要用Lin-Canny算法求解两个凸包几何体之间的最近点对，并不能直接求解非凸包体间的最近体元。但是可以对于非凸包的几何体做一些预处理，比如凸部件的分解等，进行将其组织成一棵空间划分树，再使用该算法来求解其间的距离。另外，改进的LC算法也可以用于求解两个几何体元间的最近距离，包括发生碰撞检测时的干涉距离等，这也就为碰撞检测以及物理引擎系统中的某些关键计算提出了一种很好的解决方案。

C语言代码片段

POLYHEDRON *createPolyhedron(libName, name)
     char *libName, *name;
{
  POLYHEDRON *newP;
  int i;

  for (i = 0; i < polyhedronLibraryCount && 
              strcmp(polyhedronLibrary[i].name, libName) ; i++);
  if (i == polyhedronLibraryCount) {
    printf("***** error:  can't find %s in library\n", libName);
    return NULL;
  }

  newP = allocPolyhedron;
  strcpy(newP->name, name);
  mat4copy(mat4I, newP->pose);
  newP->verts = polyhedronLibrary[i].verts;
  newP->edges = polyhedronLibrary[i].edges;
  newP->faces = polyhedronLibrary[i].faces;

  return newP;
}

/*
  Compute the voronoi region of a vertex.  This creates a list of plane
  nodes (one for each incident edge) and points the vertex's cone field
  to this list.
*/

void computeVertCone(v)
     VERTEX *v;
{
  PNODE *pn;
  FNODE *fn;
  EDGE *e;
  VECT3 u, tmpV;

  for (fn = v->edges; fn; fn = fn->next) {
    e = fn->f.e;
    /* compute vector about which to rotate vertex planes REWRITE
     this vector points toward the right of edge */
    vectAdd(e->fl->plane, e->fr->plane, tmpV);
    vectXprod(e->u, tmpV, u);
    vectNormalize(u, u);  /* bug fix:  Moritz Breipohl, 21.10.96 */

    /* construct vertex plane */
    pn = allocPnode;
    if (e->v1 == v) {
      tweakPlaneNorm(u, VERT_FLARE, e->u, pn->plane);
      vectNeg(pn->plane, pn->plane);
    }
    else tweakPlaneNorm(u, -VERT_FLARE, e->u, pn->plane);
    pn->plane[3] = - vectDotProd(v->coords, pn->plane);
    pn->nbr = e;
    addPlane(pn, &v->cone);

    /* construct edge plane */
    pn = allocPnode;
    if (e->v1 == v) tweakPlaneNorm(u, EDGE_V_FLARE, e->u, pn->plane);
    else {
      tweakPlaneNorm(u, -EDGE_V_FLARE, e->u, pn->plane);
      vectNeg(pn->plane, pn->plane);
    }
    pn->plane[3] = - vectDotProd(v->coords, pn->plane);
    pn->nbr = v;
    addPlane(pn, &e->cone);
  }
}


/*
  Compute the voronoi region of a face.  This creates a list of plane
  nodes (one for each edge of the face) and points the faces's cone field
  to this list.  The final plane of the voronoi region is the plane of
  the face itself.  This plane is stored explicitly when the face is
  created.
*/

void computeFaceCone(f)
     FACE *f;
{
  PNODE *pn;
  EDGE *e;
  VERTEX *v;
  FNODE *fnE, *fnV;
  VECT3 norm;
  
  /* we only do side planes; 
     base plane is already stored explicitly in f->plane */
  for (fnE = f->edges, fnV = f->verts; fnE; fnE = fnE->next, fnV = fnV->next) {
    e = fnE->f.e;
    v = fnV->f.v;
    /*
      (e->v1 = v) <==> e goes CCW around f; f is left face of e
      (e->v2 = v) <==> e goes CW around f; f is right face of e
    */

    /* construct face plane */
    vectXprod(f->plane, e->u, norm);
    pn = allocPnode;
    if (e->v1 == v) tweakPlaneNorm(e->u, FACE_FLARE, norm, pn->plane);
    else {
      tweakPlaneNorm(e->u, -FACE_FLARE, norm, pn->plane);
      vectNeg(pn->plane, pn->plane);
    }
    pn->plane[3] = - vectDotProd(v->coords, pn->plane);
    pn->nbr = e;
    addPlane(pn, &f->cone);

    /* construct edge plane */
    pn = allocPnode;
    if (e->v1 == v) {
      tweakPlaneNorm(e->u, EDGE_F_FLARE, norm, pn->plane);
      vectNeg(pn->plane, pn->plane);
    }
    else tweakPlaneNorm(e->u, -EDGE_F_FLARE, norm, pn->plane);
    pn->plane[3] = - vectDotProd(v->coords, pn->plane);
    pn->nbr = f;
    addPlane(pn, &e->cone);
  }
}


/*
  build the voronoi region structure for an entire polyhedron
*/

void buildCones(p)
     POLYHEDRON *p;
{
  FNODE *fn;

  for (fn = p->verts; fn; fn = fn->next) computeVertCone(fn->f.v);
  for (fn = p->faces; fn; fn = fn->next) computeFaceCone(fn->f.f);
}