拉普拉斯算子网格形变

　　网格顶点的拉普拉斯坐标定义为 ，公式中di代表顶点Vi的一环邻域顶点数量

　　网格的拉普拉斯坐标用矩阵表示， ，  然后通过网格的拉普拉斯坐标求解为稀疏线性方程组便可以得到形变后的网格顶点

　　

　　初始化拉普拉斯矩阵

void mainNode::InitLpls()
{
    int vertexNumber = m_vecAllVertex.size();
    //计算拉普拉斯坐标。方式一
    for (int i = 0; i < vertexNumber; i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        osg::Vec3 neighborPos(0.0, 0.0, 0.0);
        for (int j = 0; j < pVertex->vecNeighborVertex.size(); j++)
        {
            neighborPos += pVertex->vecNeighborVertex.at(j)->pos;
        }
        pVertex->lplsPos = pVertex->pos - neighborPos / pVertex->vecNeighborVertex.size();
    }

    //计算顶点一阶临接顶点的权值
    for (int i = 0; i < vertexNumber; i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        int neighborNumber = pVertex->vecNeighborVertex.size();
        for (int j = 0; j < neighborNumber; j++)
        {
            pVERTEX pNeighbor = pVertex->vecNeighborVertex.at(j);
            float weight = float(1) / float(neighborNumber);
            pVertex->mapWeightForOther.insert(std::pair<int, float>(pNeighbor->id, weight) );
            pVertex->fTotalWeight = 1.0;
        }
    }

    //构建拉普拉斯矩阵
    int count = 0;
    std::vector<int> beginNumber(vertexNumber);
    for (int i = 0; i < vertexNumber; i++)
    {
        beginNumber[i] = count;
        count += m_vecAllVertex[i]->vecNeighborVertex.size() + 1;
    }

    std::vector<Eigen::Triplet<float> > vectorTriplet(count);
    for (int i = 0; i < vertexNumber; i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        vectorTriplet[beginNumber[i]] = Eigen::Triplet<float>(i, i, pVertex->fTotalWeight);

        int j = 0;
        std::map<int, float>::iterator itor;
        for(itor = pVertex->mapWeightForOther.begin(); itor != pVertex->mapWeightForOther.end(); itor++)
        {
            int neighborID = itor->first;
            float weight = itor->second;
            vectorTriplet[beginNumber[i] + j + 1] = Eigen::Triplet<float>(i, neighborID, -weight);
            j++;
        }
    }

    //加入移动点以及不动点，这里便是调参的地方
    for (int i = 0; i < 2; i++)
    {
        pVERTEX pVertex;
        if (i == 0)
        {
            pVertex = m_vecAllVertex.at(i);
        }
        else
        {
            pVertex = m_vecAllVertex.at(vertexNumber - 1);
        }

        int row = i + vertexNumber;
        vectorTriplet.push_back(Eigen::Triplet<float>(row, pVertex->id, pVertex->fTotalWeight));
    }
    
    m_vectorTriplet = vectorTriplet;
    m_Matrix.resize(vertexNumber + 2, vertexNumber);
    m_Matrix.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end());
}

　　求解稀疏线性方程组

void mainNode::CalculateMaitrxLpls()
{
    Eigen::SparseMatrix<float> Matrix = m_Matrix;
    Eigen::SparseMatrix<float> MatrixTranspose = Matrix.transpose();
    Eigen::SparseMatrix<float> MatrixLpls = MatrixTranspose * Matrix;

    Eigen::VectorXf targetX, targetY, targetZ;
    int vertexNumber = m_vecAllVertex.size();
    
    targetX.resize(vertexNumber + 2);
    targetY.resize(vertexNumber + 2);
    targetZ.resize(vertexNumber + 2);
    
    for (int i = 0; i < vertexNumber + 2; i++)
    {
        if (i < vertexNumber)
        {
            targetX[i] = m_vecAllVertex.at(i)->lplsPos.x();
            targetY[i] = m_vecAllVertex.at(i)->lplsPos.y();
            targetZ[i] = m_vecAllVertex.at(i)->lplsPos.z();
        }
        else if (i < vertexNumber + 1)
        {
            //第一个的顶点坐标分量
            targetX[i] = m_vecAllVertex.at(i - vertexNumber)->pos.x();
            targetY[i] = m_vecAllVertex.at(i - vertexNumber)->pos.y();
            targetZ[i] = m_vecAllVertex.at(i - vertexNumber)->pos.z();
        }
        else
        {
            //最后一个的顶点坐标分量
            targetX[i] = m_vecAllVertex.at(vertexNumber - 1)->pos.x();
            targetY[i] = m_vecAllVertex.at(vertexNumber - 1)->pos.y();
            targetZ[i] = m_vecAllVertex.at(vertexNumber - 1)->pos.z();
        }
    }
    
    Eigen::SimplicialCholesky<Eigen::SparseMatrix<float> > MatrixCholesky(MatrixLpls);
    Eigen::VectorXf resX, resY, resZ;
    resX = MatrixTranspose * targetX;
    resY = MatrixTranspose * targetY;
    resZ = MatrixTranspose * targetZ;

    Eigen::VectorXf X, Y, Z;
    X = MatrixCholesky.solve(resX);
    //std::cout << X << std::endl;
    Y = MatrixCholesky.solve(resY);
    //std::cout << Y << std::endl;
    Z = MatrixCholesky.solve(resZ);
    //std::cout << Z << std::endl;
    for (int i = 0; i < m_vecAllVertex.size(); i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        float x, y, z;
        x = X[i];
        y = Y[i];
        z = Z[i];

        pVertex->pos = osg::Vec3(X[i], Y[i], Z[i]);
    }
}

　　参考论文：

　　Differential Coordinates for Interactive Mesh Editing

　　基于微分坐标的网格morphing