梯度泊松方程网格形变

　　网格顶点坐标的3个分量当做3个独立的标量场，如此，三角面片便有3个独立的梯度场，是为网格内在属性。当顶点移动时，问题便为求解方程，

　　根据变分法得 其中Φ是为形变后的顶点坐标，W为形变后的梯度场。方程进一步用矩阵表示为 ，L为网格拉普拉斯矩阵，b为梯度场的散度.

　　1.定义梯度算子

　　设 f 为分段线性函数，网格三角面片{Xi, Xj, Xk}的顶点处有：f(Xi) = fi ; f(Xj) = fj ; f(Xk) = fk; u = (a, b)是三角面片的重心坐标，Bi, Bj, Bk是三角面片的一次拉格朗日

　　插值函数，通过线性插值，f 在三角面片上的梯度表示为。

　　函数Bi,Bj,Bk满足Bi(u) + Bj(u) + Bk(u) = 1，所以Bi(u) + Bj(u) + Bk(u)的梯度等于0，所以 f 的梯度可表示为。

　　根据重心定理有，函数B的梯度为顶点到对边的高的倒数，方向为对边指向顶点。B的梯度表示为，A表示三角形面积，⊥表示旋转90度。

　　2. 定义散度算子

　　设置向量值函数w：S——R^{3   。}

　　S表示网格，W表示三角面片上的各个顶点向量，那么W的散度为。

　　3. 拉普拉斯算子

　　

　　

　　

　　

　构建拉普拉斯矩阵

void mainNode::InitLpls()
{
    int vertexNumber = m_vecAllVertex.size();
    //构建拉普拉斯矩阵权值，方式二
    for (int i = 0; i < vertexNumber; i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        float totalWeight = 0.0;
        for (int j = 0; j < pVertex->vecNeighborEdge.size(); j++)
        {
            pEDGE pEdge = pVertex->vecNeighborEdge.at(j);
            pVERTEX pA = pEdge->pA;
            pVERTEX pB = pEdge->pB;

            pVERTEX pTarget;
            if (pA->id == pVertex->id)
                pTarget = pB;
            else
                pTarget = pA;

            std::vector<pTRIANGLE> vecTri = pEdge->vecNeighborTri;
            pVERTEX pC = NULL;
            float weight = 0.0;
            for (int k = 0; k < vecTri.size(); k++)
            {
                pTRIANGLE pTri = vecTri.at(k);
                pVERTEX p1 = pTri->pA;
                pVERTEX p2 = pTri->pB;
                pVERTEX p3 = pTri->pC;
                if ((pA->id == p1->id && pB->id == p2->id) || (pA->id == p2->id && pB->id == p1->id))
                {
                    pC = p3;
                }
                else if ((pA->id == p1->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p1->id))
                {
                    pC = p2;
                }
                else if ((pA->id == p2->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p2->id))
                {
                    pC = p1;
                }
                //开始求取权值
                float cotAngle = 0.0;
                GetCotAngle(pA->pos, pB->pos, pC->pos, cotAngle);
                weight += cotAngle;
            }
            weight /= 2.0;
            pVertex->mapWeightForOther.insert(std::pair<int, float>(pTarget->id, weight) );
            totalWeight += weight;
        }
        pVertex->fTotalWeight = totalWeight;
    }
    //计算拉普拉斯坐标。方式二
    for (int i = 0; i < vertexNumber; i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        osg::Vec3 pos = pVertex->pos;
        pos = pos * pVertex->fTotalWeight;
        osg::Vec3 otherPos = osg::Vec3(0.0, 0.0, 0.0);
        for (int j = 0; j < pVertex->vecNeighborVertex.size(); j++)
        {
            pVERTEX pNeihbor = pVertex->vecNeighborVertex.at(j);
            std::map<int, float>::iterator itor = pVertex->mapWeightForOther.find(pNeihbor->id);
            float weight = itor->second;
            otherPos += pNeihbor->pos * weight;
        }
        pVertex->lplsPos = pos - otherPos;
    }
    
    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);
        //更改后
        //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);
            //更改后
            //vectorTriplet[beginNumber[i] + j + 1] = Eigen::Triplet<float>(neighborID, i, weight);
            j++;
        }
    }

    //头一圈和顶一圈
    for (int i = 0; i < m_iAroundNumber * 2; i++)
    {
        float weight = 1.0;
        pVERTEX pVertex;
        if (i < m_iAroundNumber)
        {
            pVertex = m_vecAllVertex.at(i);
        }
        else
        {
            pVertex = m_vecAllVertex.at(i + m_iAroundNumber * 14);
        }

        int row = i + vertexNumber;
        vectorTriplet.push_back(Eigen::Triplet<float>(row, pVertex->id, weight));
    }
    
    

    m_Matrix.resize(vertexNumber + m_iAroundNumber * 2, vertexNumber);
    //m_Matrix.resize(vertexNumber, vertexNumber);
    m_Matrix.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end());
    //std::cout << m_Matrix << std::endl;
}

　计算梯度

void mainNode::CalculateGradient()
{
    std::map<int, pTRIANGLE>::iterator itorOfTri;
    for (itorOfTri = m_mapAllTri.begin(); itorOfTri != m_mapAllTri.end(); itorOfTri++)
    {
        int triID = itorOfTri->first;
        pTRIANGLE pTri = itorOfTri->second;
        osg::Vec3 A = pTri->pA->pos;
        osg::Vec3 B = pTri->pB->pos;
        osg::Vec3 C = pTri->pC->pos;

        
        osg::Vec3 AB = B - A;
        osg::Vec3 CA = A - C;
        osg::Vec3 BC = C - B;

        osg::Vec3 normal = ((B - A) ^ (C - A));
        //三角形面积
        float area = normal.length() * 0.5;
        normal.normalize();

        float times = 2.0;
        osg::Vec3 vecGradientX = ((normal ^ BC) / (times * area) * A.x()) + ((normal ^ AB) / (times * area) * C.x()) + ((normal ^ CA) / (times * area) * B.x());
        osg::Vec3 vecGradientY = ((normal ^ BC) / (times * area) * A.y()) + ((normal ^ AB) / (times * area) * C.y()) + ((normal ^ CA) / (times * area) * B.y());
        osg::Vec3 vecGradientZ = ((normal ^ BC) / (times * area) * A.z()) + ((normal ^ AB) / (times * area) * C.z()) + ((normal ^ CA) / (times * area) * B.z());
        pTri->vecGradient[0] = vecGradientX;
        pTri->vecGradient[1] = vecGradientY;
        pTri->vecGradient[2] = vecGradientZ;
    }
}

　　计算散度

void mainNode::CalculateDivergence()
{
    
    for (int i = 0; i < m_vecAllVertex.size(); i++)
    {
        pVERTEX pVertex = m_vecAllVertex.at(i);
        std::vector<pTRIANGLE> vecTri= pVertex->vecNeighborTri;

        float divergenceX, divergenceY, divergenceZ;
        divergenceX = divergenceY = divergenceZ = 0.0;

        for (int j = 0; j < vecTri.size(); j++)
        {
            pTRIANGLE pTri = vecTri.at(j);
            pVERTEX pA = pTri->pA;
            pVERTEX pB = pTri->pB;
            pVERTEX pC = pTri->pC;
            //A为主顶点
            osg::Vec3 A, B, C;
            if (pA->id == pVertex->id)
            {
                A = pA->pos;
                B = pB->pos;
                C = pC->pos;
            }
            else if (pB->id == pVertex->id)
            {
                A = pB->pos;
                B = pC->pos;
                C = pA->pos;
            }
            else if (pC->id == pVertex->id)
            {
                A = pC->pos;
                B = pA->pos;
                C = pB->pos;
            }
            //
            osg::Vec3 AB, AC;
            float cotAngleABC, cotAngleACB;
            AB = B - A;
            AC = C - A;
            GetCotAngle(A, C, B, cotAngleABC);
            GetCotAngle(A, B, C, cotAngleACB);

            for (int k = 0; k < 3; k++)
            {
                float divergenceTemp = 0 - 0.5 * (cotAngleABC * (AC * pTri->vecGradient.at(k)) + cotAngleACB * (AB * pTri->vecGradient.at(k)));
                if (k == 0)
                {
                    divergenceX += divergenceTemp;
                }
                else if (k == 1)
                {
                    divergenceY += divergenceTemp;
                }
                else if(k == 2)
                {
                    divergenceZ += divergenceTemp;
                }
            }
        }

        pVertex->vecDivergence[0] = divergenceX;
        pVertex->vecDivergence[1] = divergenceY;
        pVertex->vecDivergence[2] = divergenceZ;
    }

}

　　参考论文：

　　Mesh Editing with Poisson-Based Gradient Field Manipulation

　　三角网格曲面上的Laplace算子及其应用

　　http://blog.sina.com.cn/s/blog_73428e9a0101h4l9.html