igllib 204 gradient

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

再看下源码:

//main.cpp
#include <igl/avg_edge_length.h>
#include <igl/barycenter.h>
#include <igl/grad.h>
#include <igl/jet.h>
#include <igl/readDMAT.h>
#include <igl/readOFF.h>
#include <igl/viewer/Viewer.h>#include <iostream>
#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 formatigl::readOFF(TUTORIAL_SHARED_PATH "/cheburashka.off", V, F);// Read scalar function values from a file, U: #V by 1VectorXd U;igl::readDMAT(TUTORIAL_SHARED_PATH "/cheburashka-scalar.dmat",U);// Compute gradient operator: #F*3 by #VSparseMatrix<double> 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 UMatrixXd GU = Map<const MatrixXd>((G*U).eval().data(),F.rows(),3);// U : 6669 * 1      U代表标量在点处的值// GU 13334 by 3 会将所有点对第i个三角形的梯度叠加起来,这样每一行就对应于一个三角形的梯度了// Compute gradient magnitudeconst VectorXd GU_mag = GU.rowwise().norm();igl::viewer::Viewer viewer;viewer.data.set_mesh(V, F);// Compute pseudocolor for original functionMatrixXd 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 barycentersMatrixXd BC;igl::barycenter(V,F,BC);// 然后算每个三角形的质心const RowVector3d black(0,0,0);viewer.data.add_edges(BC,BC+max_size*GU, black);// 在每个三角形的质心处画出梯度的方向// Hide wireframeviewer.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 <alecjacobson@gmail.com>
//
// 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 <Eigen/Geometry>
#include <vector>template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,const Eigen::PlainObjectBase<DerivedF>&F,Eigen::SparseMatrix<typename DerivedV::Scalar> &G)
{Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> eperp21(F.rows(),3), eperp13(F.rows(),3);for (int i=0;i<F.rows();++i){// renaming indices of vertices of triangles for convenienceint i1 = F(i,0);int i2 = F(i,1);int i3 = F(i,2);// #F x 3 matrices of triangle edge vectors, named after opposite verticesEigen::Matrix<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);// area of parallelogram is twice area of triangle// area of parallelogram is || v1 x v2 ||Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n  = v32.cross(v13);// This does correct l2 norm of rows, so that it contains #F list of twice// triangle areasdouble dblA = std::sqrt(n.dot(n));// now normalize normals to get unit normalsEigen::Matrix<typename DerivedV::Scalar, 1, 3> u = n / dblA;// rotate each vector 90 degrees around normaldouble 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<int> rs;rs.reserve(F.rows()*4*3);std::vector<int> cs;cs.reserve(F.rows()*4*3);std::vector<double> vs;vs.reserve(F.rows()*4*3);// row indicesfor(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 indicesfor(int r=0;r<3;r++){for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));}//列索引变成//  1 0 2 0      1 0 2 0         1 0 2 0 // valuesfor(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));// 可以发现值的前4个是x值,中间4个是y值,后面4个是z值,那么我们只需要看前面4个即可// 0 1     eperp13(i,0)//  0 0     -eperp13(i,0)// 0 2     eperp21(i,0)//  0 0     -eperp21(i,0)// 在图片就相当于在第2个点的hat function 的梯度设置为 (1/h2)e2//           在第3个点的hat function 的梯度设置为 (1/h3)e3//            在第1个点的hat function 的梯度设置为 -(1/h2)e2-(1/h3)e3 = (1/h1)e1// 证明请见:http://blog.csdn.net/seamanj/article/details/52075447//为了清晰的说明 这里加上2个三角形的情况//行索引为 0 1 0 1 0 1 0 1    2 3 2 3 2 3 2 3    4 5 4 5 4 5 4 5//列索引为 1 4 0 3 2 5 0 3    1 4 0 3 2 5 0 3    1 4 0 3 2 5 0 3 //可以看到前面8个是x部分， 然后是y, z//然后看前面8个元素 前行索引相同的分别拿出来为// 0 0 0 0  1 1 1 1// 1 0 2 0  4 3 5 3前面部分是第一个三角形的索引， 后面是另一个三角形的索引// create sparse gradient operator matrixG.resize(3*F.rows(),V.rows());std::vector<Eigen::Triplet<typename DerivedV::Scalar> > triplets;for (int i=0;i<(int)vs.size();++i){triplets.push_back(Eigen::Triplet<typename DerivedV::Scalar>(rs[i],cs[i],vs[i]));}G.setFromTriplets(triplets.begin(), triplets.end());
}#ifdef IGL_STATIC_LIBRARY
// Explicit template specialization
// template void igl::grad<double, int>(Eigen::Matrix<double, -1, -1, 0, -1,-1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&,Eigen::SparseMatrix<double, 0, int>&);
template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
//template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
template void igl::grad<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&);
#endif

it's similar for 'grad_tet' function.

igllib 204 gradient相关推荐

深度学习优化函数详解（5）-- Nesterov accelerated gradient (NAG) 优化算法
深度学习优化函数详解系列目录深度学习优化函数详解(0)– 线性回归问题深度学习优化函数详解(1)– Gradient Descent 梯度下降法深度学习优化函数详解(2)– SGD 随机梯度下降 ...
比Momentum更快：揭开Nesterov Accelerated Gradient的真面目NAG 梯度下降
d为累计梯度作为一个调参狗,每天用着深度学习框架提供的各种优化算法如Momentum.AdaDelta.Adam等,却对其中的原理不甚清楚,这样和一条咸鱼有什么分别!(误)但是我又懒得花太多时间去看 ...
PyTorch 笔记（13）— autograd（0.4 之前和之后版本差异）、Tensor（张量）、Gradient（梯度）
1. 背景简述 torch.autograd 是 PyTorch 中方便用户使用,专门开发的一套自动求导引擎,它能够根据输入和前向传播过程自动构建计算图,并执行反向传播. 计算图是现代深度学习框架 P ...
Gradient Descent和Back propagation在做什么？
Gradient Descent梯度下降实际上你要用一个Gradient Descent的方法来train一个neural network的话你应该要怎么做? 到底实际上在train neural ...
机器学习中的梯度下降（ Gradient Descent）算法
前言梯度下降(Gradient Descent,GD)算法主要分为三种:批量梯度下降(Batch Gradient Descent,BGD)算法.随机梯度下降(Stochastic Gradient ...
Gradient Descent梯度下降（透彻分析）
----------首先了解什么是梯度? 官方解释: 梯度的本意是一个向量(矢量),表示某一函数在该点处的方向导数沿着该方向取得最大值,即函数在该点处沿着该方向(此梯度的方向)变化最快,变化率最大(为 ...
python中tensor与variable_NLP实战篇之tf2中tensor、variable、gradient、ops
本文是基于tensorflow2.2.0版本,介绍了tf中变量.张量的概念,tf中梯度的计算方式和tensor相关的操作. 实战系列篇章中主要会分享,解决实际问题时的过程.遇到的问题或者使用的工具等等 ...
android 背景切换动画效果代码,关于Android shape gradient背景渐变
百度后,发现渐变色不仅可以根据xml来实现,也可以用java代码来实现,由于目前没有那么多时间,只记录xml实现的方法:以后在记录Java实现的代码. 通过Shape gradient标签来实现首先 ...
201-3-19李宏毅机器学习视频笔记七（游戏解释Gradient Descent）
视频部分: 视频7:使用帝国时代游戏解释Gradient Descent 李宏毅机器学习(2017)_哔哩哔哩 (゜-゜)つロ干杯~-bilibili https://www.bilibili.co ...

igllib 204 gradient

igllib 204 gradient相关推荐

最新文章

热门文章