本文共 4505 字，大约阅读时间需要 15 分钟。

证明请见:http://blog.csdn.net/seamanj/article/details/52075447

再看下源码:

//main.cpp#include 
   
    #include 
    
     #include 
     
      #include 
      
       #include 
       
        #include 
        
         #include 
         
          #include 
          
           #include "tutorial_shared_path.h"int main(int argc, char *argv[]){ using namespace Eigen; using namespace std; MatrixXd V; MatrixXi F; // Load a mesh in OFF format igl::readOFF(TUTORIAL_SHARED_PATH "/cheburashka.off", V, F); // Read scalar function values from a file, U: #V by 1 VectorXd U; igl::readDMAT(TUTORIAL_SHARED_PATH "/cheburashka-scalar.dmat",U); // Compute gradient operator: #F*3 by #V SparseMatrix
           
             G; igl::grad(V,F,G); // V : 6669 by 3 6669个点 // F : 13334 by 3 13334个面 // G : 40002 by 6669 每3行代表一个面(x,y,z), 每一列代表第j个点在第i个三角形内的梯度方向，注意x,y,z不是连在一起的， 行的排列为：x (#F行)，y (#F行), z(#F行) // Compute gradient of U MatrixXd GU = Map
            
             ((G*U).eval().data(),F.rows(),3); // U : 6669 * 1 U代表标量在点处的值 // GU 13334 by 3 会将所有点对第i个三角形的梯度叠加起来,这样每一行就对应于一个三角形的梯度了 // Compute gradient magnitude const VectorXd GU_mag = GU.rowwise().norm(); igl::viewer::Viewer viewer; viewer.data.set_mesh(V, F); // Compute pseudocolor for original function MatrixXd C; igl::jet(U,true,C); // // Or for gradient magnitude //igl::jet(GU_mag,true,C); viewer.data.set_colors(C); // Average edge length divided by average gradient (for scaling) const double max_size = igl::avg_edge_length(V,F) / GU_mag.mean(); // Draw a black segment in direction of gradient at face barycenters MatrixXd BC; igl::barycenter(V,F,BC); // 然后算每个三角形的质心 const RowVector3d black(0,0,0); viewer.data.add_edges(BC,BC+max_size*GU, black); // 在每个三角形的质心处画出梯度的方向 // Hide wireframe viewer.core.show_lines = false; viewer.launch();}
            
           
          
         
        
       
      
     
    
   

//grad.cpp// This file is part of libigl, a simple c++ geometry processing library.//// Copyright (C) 2013 Alec Jacobson 
   
    //// This Source Code Form is subject to the terms of the Mozilla Public License// v. 2.0. If a copy of the MPL was not distributed with this file, You can// obtain one at http://mozilla.org/MPL/2.0/.#include "grad.h"#include 
    
     #include 
     
      template 
      
       IGL_INLINE void igl::grad(const Eigen::PlainObjectBase
       
        &V,                     const Eigen::PlainObjectBase
        
         &F, Eigen::SparseMatrix
         
           &G){ Eigen::Matrix
          
            eperp21(F.rows(),3), eperp13(F.rows(),3); for (int i=0;i
           
             v32 = V.row(i3) - V.row(i2); Eigen::Matrix
            
              v13 = V.row(i1) - V.row(i3); Eigen::Matrix
             
               v21 = V.row(i2) - V.row(i1); // area of parallelogram is twice area of triangle // area of parallelogram is || v1 x v2 || Eigen::Matrix
              
                n = v32.cross(v13); // This does correct l2 norm of rows, so that it contains #F list of twice // triangle areas double dblA = std::sqrt(n.dot(n)); // now normalize normals to get unit normals Eigen::Matrix
               
                 u = n / dblA; // rotate each vector 90 degrees around normal double norm21 = std::sqrt(v21.dot(v21)); double norm13 = std::sqrt(v13.dot(v13)); eperp21.row(i) = u.cross(v21); eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i)));//确定方向 eperp21.row(i) *= norm21 / dblA;//确定长度 eperp13.row(i) = u.cross(v13); eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i))); eperp13.row(i) *= norm13 / dblA; //只需要算两个方向,另一个方向可以通过取负获得 //这里是第0个点,通过 -(1的梯度+2的梯度)获得 } std::vector
                
                  rs; rs.reserve(F.rows()*4*3); std::vector
                 
                   cs; cs.reserve(F.rows()*4*3); std::vector
                  
                    vs; vs.reserve(F.rows()*4*3); // row indices for(int r=0;r<3;r++)//处理x,y,z分量 { for(int j=0;j<4;j++)//每行有四个元素， 复制四次行索引 { for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i); //处理每个三角形 } } //这里的索引有点复杂,我们简化成只有一个三角形好了,那么行索引变为 // 0 0 0 0 1 1 1 1 2 2 2 2 // column indices for(int r=0;r<3;r++) { for(int i=0;i
                   
                     > triplets; for (int i=0;i<(int)vs.size();++i) { triplets.push_back(Eigen::Triplet
                    
                     (rs[i],cs[i],vs[i])); } G.setFromTriplets(triplets.begin(), triplets.end());}#ifdef IGL_STATIC_LIBRARY// Explicit template specialization// template void igl::grad
                     
                      (Eigen::Matrix
                      
                        const&, Eigen::Matrix
                       
                         const&,Eigen::SparseMatrix
                        
                         &);template void igl::grad
                         
                          
                           , Eigen::Matrix
                           
                             >(Eigen::PlainObjectBase
                            
                             
                               > const&, Eigen::PlainObjectBase
                              
                               
                                 > const&, Eigen::SparseMatrix
                                
                                 
                                  ::Scalar, 0, int>&);//template void igl::grad
                                  
                                   
                                    , Eigen::Matrix
                                    
                                      >(Eigen::PlainObjectBase
                                     
                                      
                                        > const&, Eigen::PlainObjectBase
                                       
                                        
                                          > const&, Eigen::SparseMatrix
                                         
                                          
                                           ::Scalar, 0, int>&);template void igl::grad
                                           
                                            
                                             , Eigen::Matrix
                                             
                                               >(Eigen::PlainObjectBase
                                              
                                               
                                                 > const&, Eigen::PlainObjectBase
                                                
                                                 
                                                   > const&, Eigen::SparseMatrix
                                                  
                                                   
                                                    ::Scalar, 0, int>&);#endif
                                                   
                                                  
                                                 
                                                
                                               
                                              
                                             
                                            
                                           
                                          
                                         
                                        
                                       
                                      
                                     
                                    
                                   
                                  
                                 
                                
                               
                              
                             
                            
                           
                          
                         
                        
                       
                      
                     
                    
                   
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     
    
   

it's similar for 'grad_tet' function.