第五节《VTK 点与网格模型求交处理技巧》

        使用VTK做项目常常需要用到单点或点集与网格模型求交计算，常用手段是使用vtkOBBTree进行求交计算。我基于VTK提供的算法总结了3种点与网格求交计算的方法，第一种vtkOBBTree直线与网格求交计算；第二种局部网格直线求交计算；第三种基于深度缓存求交计算；

        第一种vtkOBBTree直线与网格求交计算，vtkOBBTree是用于生成模型OBB树的对象数据结构，它基于模型创建一个有向包围盒，对包围盒空间划分了更多层次的空间，直线与模型相交能够快速锁定某个区域，进行求交计算，用法如下：

// polyData 待求交计算的网格模型

vtkSmartPointer<vtkOBBTree> obbTree = vtkSmartPointer<vtkOBBTree>::New();
obbTree->SetDataSet(polyData);
obbTree->BuildLocator();
double p1[3];
double p2[3]; 
vtkSmartPointer<vtkPoints> intersectPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkIdList> intersectCells = vtkSmartPointer<vtkIdList>::New();
obbTree->IntersectWithLine(p1, p2, intersectPoints, intersectCells);
参数p1 线段起点
参数p2 线段终点
参数intersectPoints 输出交点
参数intersectCells 输出交到的三角形

vtk提供的接口是一条线段与模型求交，输出结果为线段内所有的交点和所有的网格单元，实际应用要比这个复杂一些，大部分应用我们希望输入是一个点p和一个向量n来表示一条直线，输出有多种情况，1.输出距离p点最近的交点；2.输出投影点前方并且距离p点最近的交点；3.输出投影点前方并且距离p点最近，并且交到的网格单元正面；4.输出投影点前方并且距离p点最近，并且交到的网格单元背面。下面我展示一个函数来完成上述的过程：此函数所在位置是FreePolyData类里面，使用前需要执行buildOBBTree（）；

    enum MAP_TYPE
    {
        DIST_MIN,                       // 距离最近的交点
        DIST_MIN_LINEFRONT,             // 投影点前方，距离最近的点
        DIST_MIN_LINEFRONT_CELLFRONT,   // 投影点前方，距离最近的点，交到三角形正面
        DIST_MIN_LINEFRONT_CELLBACK     // 投影点前方，距离最近的点，交到三角形背面
    };
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, MAP_TYPE mapType, FreeVector &outPut)
{
    FreeVector p1 = p;
    FreeVector p2 = p;
    p1.move(n, -100000);
    p2.move(n, 100000);
    vtkSmartPointer<vtkPoints> intersectPoints = vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkIdList> intersectCells = vtkSmartPointer<vtkIdList>::New();
    if(obbTree == NULL)
    {
        buildOBBTree();
    }
    obbTree->IntersectWithLine(p1.GetData(), p2.GetData(), intersectPoints, intersectCells);

    int outPutID = -1;
    double dist = 10000000;
    for(int i = 0; i < intersectPoints->GetNumberOfPoints(); i++)
    {
        double resP[3];
        intersectPoints->GetPoint(i, resP);
        int cellID = intersectCells->GetId(i);
        FreeVector theP(resP);

        if(mapType == DIST_MIN) // 所有交点
        {
            double d = theP.point3Distance(p);
            if(d < dist)
            {
                dist = d;
                outPut = theP;
                outPutID = cellID;
            }
        } else if(mapType == DIST_MIN_LINEFRONT) { // 投影点前方的
            FreeVector n1 = p2 - p;
            FreeVector n2 = theP - p;
            if(n1*n2 > 0)
            {
                double d = theP.point3Distance(p);
                if(d < dist)
                {
                    dist = d;
                    outPut = theP;
                    outPutID = cellID;
                }
            }
        } else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT){
            FreeVector n1 = p2 - p;
            FreeVector n2 = theP - p;
            FreeVector cellNormal = getCellNormal(cellID);
            if(n1*n2 > 0 && cellNormal*n2 < 0)
            {
                double d = theP.point3Distance(p);
                if(d < dist)
                {
                    dist = d;
                    outPut = theP;
                    outPutID = cellID;
                }
            }
        } else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT) {
            FreeVector n1 = p2 - p;
            FreeVector n2 = theP - p;
            FreeVector cellNormal = getCellNormal(cellID);
            if(n1*n2 > 0 && cellNormal*n2 > 0)
            {
                double d = theP.point3Distance(p);
                if(d < dist)
                {
                    dist = d;
                    outPut = theP;
                    outPutID = cellID;
                }
            }
        }
    }
    return outPutID;
}

        第二种局部网格直线求交计算：上述采用vtkOBBTree进行求交计算，vtkOBBTree构建了层级结构，加快了求交计算方法，但有些时候我们的数据量过大时，采用vtkOBBTree计算依然不能满足计算时间，并且我们待计算的点可能在网格附近，此时采用第二种求交计算策略恰当好处。局部网格求交计算讲述的是，网格模型周围的点或点集向网格模型上的投影；

        适用场景，例如我们想在模型表面显示一条参数化曲线，曲线控制点在模型表面，当进行拟合计算后，曲线的拟合点会脱离模型表面，此时需要将曲线的拟合点向模型表面进行投影，这样才能保证显示的曲线在模型表面，下面请看局部点或点集投影算法：

/**
 * @brief FreePolyData::pointMappingToPolyData 点按照n方向与模型的投影
 * @param p 点
 * @param n 投影方向
 * @param range 投影范围
 * @param outPut 输出交点
 * @return 输出投影到的Cell ID
 */
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, double range, FreeVector &outPut)
{
    if(kdTree == NULL)
    {
        buildKdtree();
    }
    FreeVector p1 = p;
    FreeVector p2 = p;
    p1.move(n, -100);
    p2.move(n, 100);
    vtkSmartPointer<vtkIdList> result = vtkSmartPointer<vtkIdList>::New();
    kdTree->FindPointsWithinRadius(range, p.GetData(), result);
    std::vector<bool> pointMark(polyData->GetNumberOfPoints(), false);
    for(int i = 0; i < result->GetNumberOfIds(); i++)
    {
        pointMark[result->GetId(i)] = true;
    }

    for(int i = 0; i < polyData->GetNumberOfCells(); i++)
    {
        vtkCell* cell = polyData->GetCell(i);
        int v1 = cell->GetPointId(0);
        int v2 = cell->GetPointId(1);
        int v3 = cell->GetPointId(2);
        if(pointMark[v1]||pointMark[v2]||pointMark[v3])
        {

            double t;
            double intersectionCoordinates[3];
            double parametricCoordinates[3];
            int subId;
            int res = cell->IntersectWithLine(p1.GetData(), p2.GetData(), 0.0001, t, intersectionCoordinates, parametricCoordinates, subId);
            if(res == 1)
            {
                outPut = FreeVector(intersectionCoordinates[0], intersectionCoordinates[1], intersectionCoordinates[2]);
                return i;
            }

        }

    }
    return -1;
}

        需要使用vtkKdtree寻找附近点集，然后与找到的点集所在的Cell进行求交计算，范围越大计算越慢，范围越小计算越快。此函数所在位置是FreePolyData类里面，使用前需要执行buildKdtree（）；

        第三种基于深度缓存求交计算：上述两种方法求交计算还是比较耗时的，不适用大量点云计算的，如果有大量点云需要求交计算应该怎么处理呢，我的想法是想利用vtk中的相机投影矩阵进行计算，将相机视口2D坐标转换到3D坐标，获取深度信息，即可完成投影计算，条件是点集所有点的投影方向是一致的。

        1.构建局部坐标系：将投影方向当做Z轴构建坐标系

        2.将点和待求交的模型变换到局部坐标系内

        3.计算模型包围盒，计算窗口尺寸

        4.创建相机视口，并将模型数据天骄到视口中

        5.创建基于深度缓存拾取器，计算每个点的交点

        6.将结果点集变回原来坐标系，没投上的点至0,0,0

这里函数输出点集与输出点集个数相等，顺序相同，没投上的点是（0,0,0）；

/**
 * @brief pointsDepthMapping 点集与模型求交
 * @param inputPoints 输入点集
 * @param mapNormal 投影方向
 * @param polyData 模型
 * @param outPoints 输出点集，个数与输入对应，没投上是0，0，0
 */
void pointsDepthMapping(std::vector<FreeVector> &inputPoints, FreeVector mapNormal, vtkPolyData *polyData, std::vector<FreeVector> &outPoints)
{
    vtkSmartPointer<vtkPolyData> tPolyData = vtkSmartPointer<vtkPolyData>::New();
    tPolyData->DeepCopy(polyData);
    // 1.构建局部坐标系：将投影方向当做Z轴构建坐标系
    FreeVector zDir, xDir, yDir;
    zDir = -mapNormal;
    zDir.Normalize();
    constructCoordinateByZAxle(zDir, xDir, yDir); // 利用Z轴构建一个坐标系
    vtkSmartPointer<vtkMatrix4x4> m = vtkSmartPointer<vtkMatrix4x4>::New();
    m->SetElement(0, 0, xDir[0]);
    m->SetElement(1, 0, xDir[1]);
    m->SetElement(2, 0, xDir[2]);
    m->SetElement(0, 1, yDir[0]);
    m->SetElement(1, 1, yDir[1]);
    m->SetElement(2, 1, yDir[2]);
    m->SetElement(0, 2, zDir[0]);
    m->SetElement(1, 2, zDir[1]);
    m->SetElement(2, 2, zDir[2]);
    vtkSmartPointer<vtkMatrix4x4> m2 = vtkSmartPointer<vtkMatrix4x4>::New();
    m2->DeepCopy(m);
    m2->Invert();
    // 2.将点和待求交的模型变换到局部坐标系内
    transformationCurveData(m2, inputPoints);
    transformationObjectData(m2, tPolyData);
    // 3.计算模型包围盒,计算窗口尺寸
    tPolyData->ComputeBounds();
    double* bound = tPolyData->GetBounds();
    double modelW = bound[1] - bound[0];
    double modelH = bound[3] - bound[2];
    int xres = 0;
    int yres = 0;
    double pix = 2000;
    if(modelW > modelH)
    {
        xres = pix;
        yres = pix*(modelH/modelW);
    } else {
        xres = pix*(modelW/modelH);
        yres = pix;
    }
    // 4.创建相机视口
    vtkSmartPointer<vtkRenderWindow> render_win = vtkSmartPointer<vtkRenderWindow>::New ();
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New ();
    render_win->SetPosition(100000, 100000); // 起始位置移除电脑屏幕外，防止界面看到
    render_win->AddRenderer (renderer);
    render_win->SetSize (xres, yres);
    renderer->SetBackground (1.0, 0, 0);
    //render view
    vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New ();
    mapper->SetInputData(tPolyData);
    mapper->Update ();
    vtkSmartPointer<vtkActor> actor_view = vtkSmartPointer<vtkActor>::New ();
    actor_view->SetMapper (mapper);
//    renderer->SetActiveCamera (cam);
    renderer->AddActor (actor_view);
    renderer->ResetCamera();
    renderer->ResetCameraClippingRange();
    renderer->GetActiveCamera()->ParallelProjectionOn();
    render_win->Render ();
    // 5.创建基于深度拾取器
    vtkSmartPointer<vtkWorldPointPicker> worldPicker = vtkSmartPointer<vtkWorldPointPicker>::New ();
    std::vector<bool> mark(inputPoints.size(), false);
    for(int i = 0; i < inputPoints.size(); i++)
    {
        FreeVector p = inputPoints[i];
        double worldCoord[3] = {p[0], p[1], p[2]};
        vtkSmartPointer<vtkCoordinate> pCoorPress = vtkSmartPointer<vtkCoordinate>::New();
        pCoorPress->SetCoordinateSystemToWorld();
        pCoorPress->SetValue(worldCoord);
        int *dispCoord = pCoorPress->GetComputedDisplayValue(renderer);

        int x = dispCoord[0];
        int y = dispCoord[1];
        float value = render_win->GetZbufferDataAtPoint(x, y);
        if(value == 1)
        {
            outPoints.push_back(FreeVector(0, 0, 0));
            mark[i] = true;
            continue;
        }
        worldPicker->Pick (x, y, value, renderer); // 计算交点
        double coords[3] = {0};
        worldPicker->GetPickPosition (coords);
        p[2] = coords[2];
        outPoints.push_back(p);
    }
    // 开始将点变换回去
    // 旋转回去
    transformationCurveData(m, outPoints);
    for(int i = 0; i < outPoints.size(); i++) // 没投上的点置零
    {
        FreeVector& p = outPoints[i];
        if(mark[i])
        {
            p = FreeVector(0, 0, 0);
        }
    }
}

        需要注意的是，此方法计算出来的结果精度低于前两种方法，但是此方法最大有点就是计算速度极快，其中创建投影窗口会浪费几十毫秒时间。此方法函数在FreeWorld库中FreeMath中。

        FreeWorld源码安装配置下载链接请参考本专题第六节。