diff --git a/.gitignore b/.gitignore index f0fb7da9c..3fa0d594f 100644 --- a/.gitignore +++ b/.gitignore @@ -98,3 +98,5 @@ python/builddebug tutorial/.idea python/buildstatic tutorial/cmake-build-debug +.vscode/ +.idea/ diff --git a/coding-guidelines.html b/coding-guidelines.html index 5352d4c52..6f7ec37c6 100644 --- a/coding-guidelines.html +++ b/coding-guidelines.html @@ -1,14 +1,11 @@ - + libigl - - - - + @@ -60,3 +57,4 @@ + diff --git a/include/igl/exact_geodesic.cpp b/include/igl/exact_geodesic.cpp new file mode 100644 index 000000000..0ffa8b039 --- /dev/null +++ b/include/igl/exact_geodesic.cpp @@ -0,0 +1,3224 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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 "exact_geodesic.h" + +//Copyright (C) 2008 Danil Kirsanov, MIT License +//Code from https://code.google.com/archive/p/geodesic/ +// Compiled into a single file by Zhongshi Jiang + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace igl{ +namespace geodesic{ + +//#include "geodesic_constants_and_simple_functions.h" + +//double const GEODESIC_INF = std::numeric_limits::max(); +double const GEODESIC_INF = 1e100; + +//in order to avoid numerical problems with "infinitely small" intervals, +//we drop all the intervals smaller than SMALLEST_INTERVAL_RATIO*edge_length +double const SMALLEST_INTERVAL_RATIO = 1e-6; +//double const SMALL_EPSILON = 1e-10; + + +inline double cos_from_edges(double const a, //compute the cosine of the angle given the lengths of the edges + double const b, + double const c) +{ + assert(a>1e-50); + assert(b>1e-50); + assert(c>1e-50); + + double result = (b*b + c*c - a*a)/(2.0*b*c); + result = std::max(result, -1.0); + return std::min(result, 1.0); +} + +inline double angle_from_edges(double const a, //compute the cosine of the angle given the lengths of the edges + double const b, + double const c) +{ + return acos(cos_from_edges(a,b,c)); +} + +template +inline bool read_mesh_from_file(char* filename, + Points& points, + Faces& faces) +{ + std::ifstream file(filename); + assert(file.is_open()); + if(!file.is_open()) return false; + + unsigned num_points; + file >> num_points; + assert(num_points>=3); + + unsigned num_faces; + file >> num_faces; + + points.resize(num_points*3); + for(typename Points::iterator i=points.begin(); i!=points.end(); ++i) + { + file >> *i; + } + + faces.resize(num_faces*3); + for(typename Faces::iterator i=faces.begin(); i!=faces.end(); ++i) + { + file >> *i; + } + file.close(); + + return true; +} + +// #include "geodesic_memory" +template //quickly allocates multiple elements of a given type; no deallocation +class SimlpeMemoryAllocator +{ +public: + typedef T* pointer; + + SimlpeMemoryAllocator(unsigned block_size = 0, + unsigned max_number_of_blocks = 0) + { + reset(block_size, + max_number_of_blocks); + }; + + ~SimlpeMemoryAllocator(){}; + + void reset(unsigned block_size, + unsigned max_number_of_blocks) + { + m_block_size = block_size; + m_max_number_of_blocks = max_number_of_blocks; + + + m_current_position = 0; + + m_storage.reserve(max_number_of_blocks); + m_storage.resize(1); + m_storage[0].resize(block_size); + }; + + pointer allocate(unsigned const n) //allocate n units + { + assert(n < m_block_size); + + if(m_current_position + n >= m_block_size) + { + m_storage.push_back( std::vector() ); + m_storage.back().resize(m_block_size); + m_current_position = 0; + } + pointer result = & m_storage.back()[m_current_position]; + m_current_position += n; + + return result; + }; +private: + std::vector > m_storage; + unsigned m_block_size; //size of a single block + unsigned m_max_number_of_blocks; //maximum allowed number of blocks + unsigned m_current_position; //first unused element inside the current block +}; + + +template //quickly allocates and deallocates single elements of a given type +class MemoryAllocator +{ +public: + typedef T* pointer; + + MemoryAllocator(unsigned block_size = 1024, + unsigned max_number_of_blocks = 1024) + { + reset(block_size, + max_number_of_blocks); + }; + + ~MemoryAllocator(){}; + + void clear() + { + reset(m_block_size, + m_max_number_of_blocks); + } + + void reset(unsigned block_size, + unsigned max_number_of_blocks) + { + m_block_size = block_size; + m_max_number_of_blocks = max_number_of_blocks; + + assert(m_block_size > 0); + assert(m_max_number_of_blocks > 0); + + m_current_position = 0; + + m_storage.reserve(max_number_of_blocks); + m_storage.resize(1); + m_storage[0].resize(block_size); + + m_deleted.clear(); + m_deleted.reserve(2*block_size); + }; + + pointer allocate() //allocates single unit of memory + { + pointer result; + if(m_deleted.empty()) + { + if(m_current_position + 1 >= m_block_size) + { + m_storage.push_back( std::vector() ); + m_storage.back().resize(m_block_size); + m_current_position = 0; + } + result = & m_storage.back()[m_current_position]; + ++m_current_position; + } + else + { + result = m_deleted.back(); + m_deleted.pop_back(); + } + + return result; + }; + + void deallocate(pointer p) //allocate n units + { + if(m_deleted.size() < m_deleted.capacity()) + { + m_deleted.push_back(p); + } + }; + +private: + std::vector > m_storage; + unsigned m_block_size; //size of a single block + unsigned m_max_number_of_blocks; //maximum allowed number of blocks + unsigned m_current_position; //first unused element inside the current block + + std::vector m_deleted; //pointers to deleted elemets +}; + + +class OutputBuffer +{ +public: + OutputBuffer(): + m_num_bytes(0) + {} + + void clear() + { + m_num_bytes = 0; + m_buffer = std::shared_ptr(); + } + + template + T* allocate(unsigned n) + { + double wanted = n*sizeof(T); + if(wanted > m_num_bytes) + { + unsigned new_size = (unsigned) ceil(wanted / (double)sizeof(double)); + m_buffer = std::shared_ptr(new double[new_size]); + m_num_bytes = new_size*sizeof(double); + } + + return (T*)m_buffer.get(); + } + + template + T* get() + { + return (T*)m_buffer.get(); + } + + template + unsigned capacity() + { + return (unsigned)floor((double)m_num_bytes/(double)sizeof(T)); + }; + +private: + + std::shared_ptr m_buffer; + unsigned m_num_bytes; +}; + + + + +class Vertex; +class Edge; +class Face; +class Mesh; +class MeshElementBase; + +typedef Vertex* vertex_pointer; +typedef Edge* edge_pointer; +typedef Face* face_pointer; +typedef Mesh* mesh_pointer; +typedef MeshElementBase* base_pointer; + +template //simple vector that stores info about mesh references +class SimpleVector //for efficiency, it uses an outside memory allocator +{ +public: + SimpleVector(): + m_size(0), + m_begin(NULL) + {}; + + typedef Data* iterator; + + unsigned size(){return m_size;}; + iterator begin(){return m_begin;}; + iterator end(){return m_begin + m_size;}; + + template + void set_allocation(DataPointer begin, unsigned size) + { + assert(begin != NULL || size == 0); + m_size = size; + m_begin = (iterator)begin; + } + + Data& operator[](unsigned i) + { + assert(i < m_size); + return *(m_begin + i); + } + + void clear() + { + m_size = 0; + m_begin = NULL; + } + +private: + unsigned m_size; + Data* m_begin; +}; + +enum PointType +{ + VERTEX, + EDGE, + FACE, + UNDEFINED_POINT +}; + +class MeshElementBase //prototype of vertices, edges and faces +{ +public: + typedef SimpleVector vertex_pointer_vector; + typedef SimpleVector edge_pointer_vector; + typedef SimpleVector face_pointer_vector; + + MeshElementBase(): + m_id(0), + m_type(UNDEFINED_POINT) + {}; + + vertex_pointer_vector& adjacent_vertices(){return m_adjacent_vertices;}; + edge_pointer_vector& adjacent_edges(){return m_adjacent_edges;}; + face_pointer_vector& adjacent_faces(){return m_adjacent_faces;}; + + unsigned& id(){return m_id;}; + PointType type(){return m_type;}; + +protected: + vertex_pointer_vector m_adjacent_vertices; //list of the adjacent vertices + edge_pointer_vector m_adjacent_edges; //list of the adjacent edges + face_pointer_vector m_adjacent_faces; //list of the adjacent faces + + unsigned m_id; //unique id + PointType m_type; //vertex, edge or face +}; + +class Point3D //point in 3D and corresponding operations +{ +public: + Point3D(){}; + Point3D(Point3D* p) + { + x() = p->x(); + y() = p->y(); + z() = p->z(); + }; + + double* xyz(){return m_coordinates;}; + double& x(){return *m_coordinates;}; + double& y(){return *(m_coordinates+1);}; + double& z(){return *(m_coordinates+2);}; + + void set(double new_x, double new_y, double new_z) + { + x() = new_x; + y() = new_y; + z() = new_z; + } + + void set(double* data) + { + x() = *data; + y() = *(data+1); + z() = *(data+2); + } + + double distance(double* v) + { + double dx = m_coordinates[0] - v[0]; + double dy = m_coordinates[1] - v[1]; + double dz = m_coordinates[2] - v[2]; + + return sqrt(dx*dx + dy*dy + dz*dz); + }; + + double distance(Point3D* v) + { + return distance(v->xyz()); + }; + + void add(Point3D* v) + { + x() += v->x(); + y() += v->y(); + z() += v->z(); + }; + + void multiply(double v) + { + x() *= v; + y() *= v; + z() *= v; + }; + +private: + double m_coordinates[3]; //xyz +}; + +class Vertex: public MeshElementBase, public Point3D +{ +public: + Vertex() + { + m_type = VERTEX; + }; + + ~Vertex(){}; + + bool& saddle_or_boundary(){return m_saddle_or_boundary;}; +private: + //this flag speeds up exact geodesic algorithm + bool m_saddle_or_boundary; //it is true if total adjacent angle is larger than 2*PI or this vertex belongs to the mesh boundary +}; + + +class Face: public MeshElementBase +{ +public: + Face() + { + m_type = FACE; + }; + + ~Face(){}; + + edge_pointer opposite_edge(vertex_pointer v); + vertex_pointer opposite_vertex(edge_pointer e); + edge_pointer next_edge(edge_pointer e, vertex_pointer v); + + double vertex_angle(vertex_pointer v) + { + for(unsigned i=0; i<3; ++i) + { + if(adjacent_vertices()[i]->id() == v->id()) + { + return m_corner_angles[i]; + } + } + assert(0); + return 0; + } + + double* corner_angles(){return m_corner_angles;}; + +private: + double m_corner_angles[3]; //triangle angles in radians; angles correspond to vertices in m_adjacent_vertices +}; + +class Edge: public MeshElementBase +{ +public: + Edge() + { + m_type = EDGE; + }; + + ~Edge(){}; + + double& length(){return m_length;}; + + face_pointer opposite_face(face_pointer f) + { + if(adjacent_faces().size() == 1) + { + assert(adjacent_faces()[0]->id() == f->id()); + return NULL; + } + + assert(adjacent_faces()[0]->id() == f->id() || + adjacent_faces()[1]->id() == f->id()); + + return adjacent_faces()[0]->id() == f->id() ? + adjacent_faces()[1] : adjacent_faces()[0]; + }; + + vertex_pointer opposite_vertex(vertex_pointer v) + { + assert(belongs(v)); + + return adjacent_vertices()[0]->id() == v->id() ? + adjacent_vertices()[1] : adjacent_vertices()[0]; + }; + + bool belongs(vertex_pointer v) + { + return adjacent_vertices()[0]->id() == v->id() || + adjacent_vertices()[1]->id() == v->id(); + } + + bool is_boundary(){return adjacent_faces().size() == 1;}; + + vertex_pointer v0(){return adjacent_vertices()[0];}; + vertex_pointer v1(){return adjacent_vertices()[1];}; + + void local_coordinates(Point3D* point, + double& x, + double& y) + { + double d0 = point->distance(v0()); + if(d0 < 1e-50) + { + x = 0.0; + y = 0.0; + return; + } + + double d1 = point->distance(v1()); + if(d1 < 1e-50) + { + x = m_length; + y = 0.0; + return; + } + + x = m_length/2.0 + (d0*d0 - d1*d1)/(2.0*m_length); + y = sqrt(std::max(0.0, d0*d0 - x*x)); + return; + } + +private: + double m_length; //length of the edge +}; + +class SurfacePoint:public Point3D //point on the surface of the mesh +{ +public: + SurfacePoint(): + m_p(NULL) + {}; + + SurfacePoint(vertex_pointer v): //set the surface point in the vertex + SurfacePoint::Point3D(v), + m_p(v) + {}; + + SurfacePoint(face_pointer f): //set the surface point in the center of the face + m_p(f) + { + set(0,0,0); + add(f->adjacent_vertices()[0]); + add(f->adjacent_vertices()[1]); + add(f->adjacent_vertices()[2]); + multiply(1./3.); + }; + + SurfacePoint(edge_pointer e, //set the surface point in the middle of the edge + double a = 0.5): + m_p(e) + { + double b = 1 - a; + + vertex_pointer v0 = e->adjacent_vertices()[0]; + vertex_pointer v1 = e->adjacent_vertices()[1]; + + x() = b*v0->x() + a*v1->x(); + y() = b*v0->y() + a*v1->y(); + z() = b*v0->z() + a*v1->z(); + }; + + SurfacePoint(base_pointer g, + double x, + double y, + double z, + PointType t = UNDEFINED_POINT): + m_p(g) + { + set(x,y,z); + }; + + void initialize(SurfacePoint const& p) + { + *this = p; + } + + ~SurfacePoint(){}; + + PointType type(){return m_p ? m_p->type() : UNDEFINED_POINT;}; + base_pointer& base_element(){return m_p;}; +protected: + base_pointer m_p; //could be face, vertex or edge pointer +}; + +inline edge_pointer Face::opposite_edge(vertex_pointer v) +{ + for(unsigned i=0; i<3; ++i) + { + edge_pointer e = adjacent_edges()[i]; + if(!e->belongs(v)) + { + return e; + } + } + assert(0); + return NULL; +} + +inline vertex_pointer Face::opposite_vertex(edge_pointer e) +{ + for(unsigned i=0; i<3; ++i) + { + vertex_pointer v = adjacent_vertices()[i]; + if(!e->belongs(v)) + { + return v; + } + } + assert(0); + return NULL; +} + +inline edge_pointer Face::next_edge(edge_pointer e, vertex_pointer v) +{ + assert(e->belongs(v)); + + for(unsigned i=0; i<3; ++i) + { + edge_pointer next = adjacent_edges()[i]; + if(e->id() != next->id() && next->belongs(v)) + { + return next; + } + } + assert(0); + return NULL; +} + +struct HalfEdge //prototype of the edge; used for mesh construction +{ + unsigned face_id; + unsigned vertex_0; //adjacent vertices sorted by id value + unsigned vertex_1; //they are sorted, vertex_0 < vertex_1 +}; + +inline bool operator < (const HalfEdge &x, const HalfEdge &y) +{ + if(x.vertex_0 == y.vertex_0) + { + return x.vertex_1 < y.vertex_1; + } + else + { + return x.vertex_0 < y.vertex_0; + } +} + +inline bool operator != (const HalfEdge &x, const HalfEdge &y) +{ + return x.vertex_0 != y.vertex_0 || x.vertex_1 != y.vertex_1; +} + +inline bool operator == (const HalfEdge &x, const HalfEdge &y) +{ + return x.vertex_0 == y.vertex_0 && x.vertex_1 == y.vertex_1; +} + +struct edge_visible_from_source +{ + unsigned source; + edge_pointer edge; +}; + +class Mesh +{ +public: + Mesh() + {}; + + ~Mesh(){}; + + template + void initialize_mesh_data(unsigned num_vertices, + Points& p, + unsigned num_faces, + Faces& tri); //build mesh from regular point-triangle representation + + template + void initialize_mesh_data(Points& p, Faces& tri); //build mesh from regular point-triangle representation + + std::vector& vertices(){return m_vertices;}; + std::vector& edges(){return m_edges;}; + std::vector& faces(){return m_faces;}; + + unsigned closest_vertices(SurfacePoint* p, + std::vector* storage = NULL); //list vertices closest to the point + +private: + + void build_adjacencies(); //build internal structure of the mesh + bool verify(); //verifies connectivity of the mesh and prints some debug info + + typedef void* void_pointer; + void_pointer allocate_pointers(unsigned n) + { + return m_pointer_allocator.allocate(n); + } + + std::vector m_vertices; + std::vector m_edges; + std::vector m_faces; + + SimlpeMemoryAllocator m_pointer_allocator; //fast memory allocating for Face/Vertex/Edge cross-references +}; + +inline unsigned Mesh::closest_vertices(SurfacePoint* p, + std::vector* storage) +{ + assert(p->type() != UNDEFINED_POINT); + + if(p->type() == VERTEX) + { + if(storage) + { + storage->push_back(static_cast(p->base_element())); + } + return 1; + } + else if(p->type() == FACE) + { + if(storage) + { + vertex_pointer* vp= p->base_element()->adjacent_vertices().begin(); + storage->push_back(*vp); + storage->push_back(*(vp+1)); + storage->push_back(*(vp+2)); + } + return 2; + } + else if(p->type() == EDGE) //for edge include all 4 adjacent vertices + { + edge_pointer edge = static_cast(p->base_element()); + + if(storage) + { + storage->push_back(edge->adjacent_vertices()[0]); + storage->push_back(edge->adjacent_vertices()[1]); + + for(unsigned i = 0; i < edge->adjacent_faces().size(); ++i) + { + face_pointer face = edge->adjacent_faces()[i]; + storage->push_back(face->opposite_vertex(edge)); + } + } + return 2 + edge->adjacent_faces().size(); + } + + assert(0); + return 0; +} + +template +void Mesh::initialize_mesh_data(Points& p, Faces& tri) //build mesh from regular point-triangle representation +{ + assert(p.size() % 3 == 0); + unsigned const num_vertices = p.size() / 3; + assert(tri.size() % 3 == 0); + unsigned const num_faces = tri.size() / 3; + + initialize_mesh_data(num_vertices, p, num_faces, tri); +} + +template +void Mesh::initialize_mesh_data(unsigned num_vertices, + Points& p, + unsigned num_faces, + Faces& tri) +{ + unsigned const approximate_number_of_internal_pointers = (num_vertices + num_faces)*4; + unsigned const max_number_of_pointer_blocks = 100; + m_pointer_allocator.reset(approximate_number_of_internal_pointers, + max_number_of_pointer_blocks); + + m_vertices.resize(num_vertices); + for(unsigned i=0; iadjacent Faces + std::vector count(m_vertices.size()); //count adjacent vertices + for(unsigned i=0; iid(); + assert(vertex_id < m_vertices.size()); + count[vertex_id]++; + } + } + + for(unsigned i=0; iadjacent_faces()[count[v->id()]++] = &f; + } + } + + //find all edges + //i.e. find all half-edges, sort and combine them into edges + std::vector half_edges(m_faces.size()*3); + unsigned k = 0; + for(unsigned i=0; iid(); + unsigned vertex_id_2 = f.adjacent_vertices()[(j+1) % 3]->id(); + half_edges[k].vertex_0 = std::min(vertex_id_1, vertex_id_2); + half_edges[k].vertex_1 = std::max(vertex_id_1, vertex_id_2); + + k++; + } + } + std::sort(half_edges.begin(), half_edges.end()); + + unsigned number_of_edges = 1; + for(unsigned i=1; iadjacent Vertices and Faces + m_edges.resize(number_of_edges); + unsigned edge_id = 0; + for(unsigned i=0; idistance(e.adjacent_vertices()[1]); + assert(e.length() > 1e-100); //algorithm works well with non-degenerate meshes only + + if(i != half_edges.size()-1 && half_edges[i] == half_edges[i+1]) //double edge + { + e.adjacent_faces().set_allocation(allocate_pointers(2),2); + e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id]; + e.adjacent_faces()[1] = &m_faces[half_edges[i+1].face_id]; + i += 2; + } + else //single edge + { + e.adjacent_faces().set_allocation(allocate_pointers(1),1); //one adjucent faces + e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id]; + i += 1; + } + } + + // Vertices->adjacent Edges + std::fill(count.begin(), count.end(), 0); + for(unsigned i=0; iid()]++; + count[e.adjacent_vertices()[1]->id()]++; + } + for(unsigned i=0; iadjacent_edges()[count[v->id()]++] = &e; + } + } + + // Faces->adjacent Edges + for(unsigned i=0; iid()]<3); + f->adjacent_edges()[count[f->id()]++] = &e; + } + } + + //compute angles for the faces + for(unsigned i=0; ilength(); + } + + double angle = angle_from_edges(abc[0], abc[1], abc[2]); + assert(angle>1e-5); //algorithm works well with non-degenerate meshes only + + f.corner_angles()[j] = angle; + sum += angle; + } + assert(std::abs(sum - igl::PI) < 1e-5); //algorithm works well with non-degenerate meshes only + } + + //define m_turn_around_flag for vertices + std::vector total_vertex_angle(m_vertices.size()); + for(unsigned i=0; iid()] += f.corner_angles()[j]; + } + } + + for(unsigned i=0; i 2.0*igl::PI - 1e-5); + } + + for(unsigned i=0; isaddle_or_boundary() = true; + e.adjacent_vertices()[1]->saddle_or_boundary() = true; + } + } + + assert(verify()); +} + +inline bool Mesh::verify() //verifies connectivity of the mesh and prints some debug info +{ + std::cout << std::endl; + // make sure that all vertices are mentioned at least once. + // though the loose vertex is not a bug, it most likely indicates that something is wrong with the mesh + std::vector map(m_vertices.size(), false); + for(unsigned i=0; iadjacent_vertices()[0]->id()] = true; + map[e->adjacent_vertices()[1]->id()] = true; + } + assert(std::find(map.begin(), map.end(), false) == map.end()); + + //make sure that the mesh is connected trough its edges + //if mesh has more than one connected component, it is most likely a bug + std::vector stack(1,&m_faces[0]); + stack.reserve(m_faces.size()); + + map.resize(m_faces.size()); + std::fill(map.begin(), map.end(), false); + map[0] = true; + + while(!stack.empty()) + { + face_pointer f = stack.back(); + stack.pop_back(); + + for(unsigned i=0; i<3; ++i) + { + edge_pointer e = f->adjacent_edges()[i]; + face_pointer f_adjacent = e->opposite_face(f); + if(f_adjacent && !map[f_adjacent->id()]) + { + map[f_adjacent->id()] = true; + stack.push_back(f_adjacent); + } + } + } + assert(std::find(map.begin(), map.end(), false) == map.end()); + + //print some mesh statistics that can be useful in debugging + // std::cout << "mesh has " << m_vertices.size() + // << " vertices, " << m_faces.size() + // << " faces, " << m_edges.size() + // << " edges\n"; + + unsigned total_boundary_edges = 0; + double longest_edge = 0; + double shortest_edge = 1e100; + for(unsigned i=0; iset(data); + unsigned type = (unsigned) data[3]; + unsigned id = (unsigned) data[4]; + + + if(type == 0) //vertex + { + point->base_element() = &mesh->vertices()[id]; + } + else if(type == 1) //edge + { + point->base_element() = &mesh->edges()[id]; + } + else //face + { + point->base_element() = &mesh->faces()[id]; + } +} + +inline void fill_surface_point_double(geodesic::SurfacePoint* point, + double* data, + long mesh_id) +{ + data[0] = point->x(); + data[1] = point->y(); + data[2] = point->z(); + data[4] = point->base_element()->id(); + + if(point->type() == VERTEX) //vertex + { + data[3] = 0; + } + else if(point->type() == EDGE) //edge + { + data[3] = 1; + } + else //face + { + data[3] = 2; + } +} + +class Interval; +class IntervalList; +typedef Interval* interval_pointer; +typedef IntervalList* list_pointer; + +class Interval //interval of the edge +{ +public: + + Interval(){}; + ~Interval(){}; + + enum DirectionType + { + FROM_FACE_0, + FROM_FACE_1, + FROM_SOURCE, + UNDEFINED_DIRECTION + }; + + double signal(double x) //geodesic distance function at point x + { + assert(x>=0.0 && x <= m_edge->length()); + + if(m_d == GEODESIC_INF) + { + return GEODESIC_INF; + } + else + { + double dx = x - m_pseudo_x; + if(m_pseudo_y == 0.0) + { + return m_d + std::abs(dx); + } + else + { + return m_d + sqrt(dx*dx + m_pseudo_y*m_pseudo_y); + } + } + } + + double max_distance(double end) + { + if(m_d == GEODESIC_INF) + { + return GEODESIC_INF; + } + else + { + double a = std::abs(m_start - m_pseudo_x); + double b = std::abs(end - m_pseudo_x); + + return a > b ? m_d + sqrt(a*a + m_pseudo_y*m_pseudo_y): + m_d + sqrt(b*b + m_pseudo_y*m_pseudo_y); + } + } + + void compute_min_distance(double stop) //compute min, given c,d theta, start, end. + { + assert(stop > m_start); + + if(m_d == GEODESIC_INF) + { + m_min = GEODESIC_INF; + } + else if(m_start > m_pseudo_x) + { + m_min = signal(m_start); + } + else if(stop < m_pseudo_x) + { + m_min = signal(stop); + } + else + { + assert(m_pseudo_y<=0); + m_min = m_d - m_pseudo_y; + } + } + //compare two intervals in the queue + bool operator()(interval_pointer const x, interval_pointer const y) const + { + if(x->min() != y->min()) + { + return x->min() < y->min(); + } + else if(x->start() != y->start()) + { + return x->start() < y->start(); + } + else + { + return x->edge()->id() < y->edge()->id(); + } + } + + double stop() //return the endpoint of the interval + { + return m_next ? m_next->start() : m_edge->length(); + } + + double hypotenuse(double a, double b) + { + return sqrt(a*a + b*b); + } + + void find_closest_point(double const x, + double const y, + double& offset, + double& distance); //find the point on the interval that is closest to the point (alpha, s) + + double& start(){return m_start;}; + double& d(){return m_d;}; + double& pseudo_x(){return m_pseudo_x;}; + double& pseudo_y(){return m_pseudo_y;}; + double& min(){return m_min;}; + interval_pointer& next(){return m_next;}; + edge_pointer& edge(){return m_edge;}; + DirectionType& direction(){return m_direction;}; + bool visible_from_source(){return m_direction == FROM_SOURCE;}; + unsigned& source_index(){return m_source_index;}; + + void initialize(edge_pointer edge, + SurfacePoint* point = NULL, + unsigned source_index = 0); + +protected: + double m_start; //initial point of the interval on the edge + double m_d; //distance from the source to the pseudo-source + double m_pseudo_x; //coordinates of the pseudo-source in the local coordinate system + double m_pseudo_y; //y-coordinate should be always negative + double m_min; //minimum distance on the interval + + interval_pointer m_next; //pointer to the next interval in the list + edge_pointer m_edge; //edge that the interval belongs to + unsigned m_source_index; //the source it belongs to + DirectionType m_direction; //where the interval is coming from +}; + +struct IntervalWithStop : public Interval +{ +public: + double& stop(){return m_stop;}; +protected: + double m_stop; +}; + +class IntervalList //list of the of intervals of the given edge +{ +public: + IntervalList(){m_first = NULL;}; + ~IntervalList(){}; + + void clear() + { + m_first = NULL; + }; + + void initialize(edge_pointer e) + { + m_edge = e; + m_first = NULL; + }; + + interval_pointer covering_interval(double offset) //returns the interval that covers the offset + { + assert(offset >= 0.0 && offset <= m_edge->length()); + + interval_pointer p = m_first; + while(p && p->stop() < offset) + { + p = p->next(); + } + + return p;// && p->start() <= offset ? p : NULL; + }; + + void find_closest_point(SurfacePoint* point, + double& offset, + double& distance, + interval_pointer& interval) + { + interval_pointer p = m_first; + distance = GEODESIC_INF; + interval = NULL; + + double x,y; + m_edge->local_coordinates(point, x, y); + + while(p) + { + if(p->min()find_closest_point(x, y, o, d); + if(d < distance) + { + distance = d; + offset = o; + interval = p; + } + } + p = p->next(); + } + }; + + unsigned number_of_intervals() + { + interval_pointer p = m_first; + unsigned count = 0; + while(p) + { + ++count; + p = p->next(); + } + return count; + } + + interval_pointer last() + { + interval_pointer p = m_first; + if(p) + { + while(p->next()) + { + p = p->next(); + } + } + return p; + } + + double signal(double x) + { + interval_pointer interval = covering_interval(x); + + return interval ? interval->signal(x) : GEODESIC_INF; + } + + interval_pointer& first(){return m_first;}; + edge_pointer& edge(){return m_edge;}; +private: + interval_pointer m_first; //pointer to the first member of the list + edge_pointer m_edge; //edge that owns this list +}; + +class SurfacePointWithIndex : public SurfacePoint +{ +public: + unsigned index(){return m_index;}; + + void initialize(SurfacePoint& p, unsigned index) + { + SurfacePoint::initialize(p); + m_index = index; + } + + bool operator()(SurfacePointWithIndex* x, SurfacePointWithIndex* y) const //used for sorting + { + assert(x->type() != UNDEFINED_POINT && y->type() !=UNDEFINED_POINT); + + if(x->type() != y->type()) + { + return x->type() < y->type(); + } + else + { + return x->base_element()->id() < y->base_element()->id(); + } + } + +private: + unsigned m_index; +}; + +class SortedSources : public std::vector +{ +private: + typedef std::vector sorted_vector_type; +public: + typedef sorted_vector_type::iterator sorted_iterator; + typedef std::pair sorted_iterator_pair; + + sorted_iterator_pair sources(base_pointer mesh_element) + { + m_search_dummy.base_element() = mesh_element; + + return equal_range(m_sorted.begin(), + m_sorted.end(), + &m_search_dummy, + m_compare_less); + } + + void initialize(std::vector& sources) //we initialize the sources by copie + { + resize(sources.size()); + m_sorted.resize(sources.size()); + for(unsigned i=0; ilength(); + if(std::abs(hs+hc) < local_epsilon) + { + if(rs<=m_start) + { + r = m_start; + d_out = signal(m_start) + std::abs(rs - m_start); + } + else if(rs>=end) + { + r = end; + d_out = signal(end) + fabs(end - rs); + } + else + { + r = rs; + d_out = signal(rs); + } + } + else + { + double ri = (rs*hc + hs*rc)/(hs+hc); + + if(riend) + { + r = end; + d_out = signal(end) + hypotenuse(end - rs, hs); + } + else + { + r = ri; + d_out = m_d + hypotenuse(rc - rs, hc + hs); + } + } + } + + +inline void Interval::initialize(edge_pointer edge, + SurfacePoint* source, + unsigned source_index) +{ + m_next = NULL; + //m_geodesic_previous = NULL; + m_direction = UNDEFINED_DIRECTION; + m_edge = edge; + m_source_index = source_index; + + m_start = 0.0; + //m_stop = edge->length(); + if(!source) + { + m_d = GEODESIC_INF; + m_min = GEODESIC_INF; + return; + } + m_d = 0; + + if(source->base_element()->type() == VERTEX) + { + if(source->base_element()->id() == edge->v0()->id()) + { + m_pseudo_x = 0.0; + m_pseudo_y = 0.0; + m_min = 0.0; + return; + } + else if(source->base_element()->id() == edge->v1()->id()) + { + m_pseudo_x = stop(); + m_pseudo_y = 0.0; + m_min = 0.0; + return; + } + } + + edge->local_coordinates(source, m_pseudo_x, m_pseudo_y); + m_pseudo_y = -m_pseudo_y; + + compute_min_distance(stop()); +} + + + +// #include "geodesic_algorithm_base.h" +class GeodesicAlgorithmBase +{ +public: + enum AlgorithmType + { + EXACT, + DIJKSTRA, + SUBDIVISION, + UNDEFINED_ALGORITHM + }; + + GeodesicAlgorithmBase(geodesic::Mesh* mesh): + m_type(UNDEFINED_ALGORITHM), + m_max_propagation_distance(1e100), + m_mesh(mesh) + {}; + + virtual ~GeodesicAlgorithmBase(){}; + + virtual void propagate(std::vector& sources, + double max_propagation_distance = GEODESIC_INF, //propagation algorithm stops after reaching the certain distance from the source + std::vector* stop_points = NULL) = 0; //or after ensuring that all the stop_points are covered + + virtual void trace_back(SurfacePoint& destination, //trace back piecewise-linear path + std::vector& path) = 0; + + void geodesic(SurfacePoint& source, + SurfacePoint& destination, + std::vector& path); //lazy people can find geodesic path with one function call + + void geodesic(std::vector& sources, + std::vector& destinations, + std::vector >& paths); //lazy people can find geodesic paths with one function call + + virtual unsigned best_source(SurfacePoint& point, //after propagation step is done, quickly find what source this point belongs to and what is the distance to this source + double& best_source_distance) = 0; + + virtual void print_statistics() //print info about timing and memory usage in the propagation step of the algorithm + { + std::cout << "propagation step took " << m_time_consumed << " seconds " << std::endl; + }; + + AlgorithmType type(){return m_type;}; + + virtual std::string name(); + + geodesic::Mesh* mesh(){return m_mesh;}; +protected: + + void set_stop_conditions(std::vector* stop_points, + double stop_distance); + double stop_distance() + { + return m_max_propagation_distance; + } + + AlgorithmType m_type; // type of the algorithm + + typedef std::pair stop_vertex_with_distace_type; + std::vector m_stop_vertices; // algorithm stops propagation after covering certain vertices + double m_max_propagation_distance; // or reaching the certain distance + + geodesic::Mesh* m_mesh; + + double m_time_consumed; //how much time does the propagation step takes + double m_propagation_distance_stopped; //at what distance (if any) the propagation algorithm stopped +}; + +inline double length(std::vector& path) +{ + double length = 0; + if(!path.empty()) + { + for(unsigned i=0; i& path) +{ + std::cout << "number of the points in the path = " << path.size() + << ", length of the path = " << length(path) + << std::endl; +} + +inline std::string GeodesicAlgorithmBase::name() +{ + switch(m_type) + { + case EXACT: + return "exact"; + case DIJKSTRA: + return "dijkstra"; + case SUBDIVISION: + return "subdivision"; + default: + case UNDEFINED_ALGORITHM: + return "undefined"; + } +} + +inline void GeodesicAlgorithmBase::geodesic(SurfacePoint& source, + SurfacePoint& destination, + std::vector& path) //lazy people can find geodesic path with one function call +{ + std::vector sources(1, source); + std::vector stop_points(1, destination); + double const max_propagation_distance = GEODESIC_INF; + + propagate(sources, + max_propagation_distance, + &stop_points); + + trace_back(destination, path); +} + +inline void GeodesicAlgorithmBase::geodesic(std::vector& sources, + std::vector& destinations, + std::vector >& paths) //lazy people can find geodesic paths with one function call +{ + double const max_propagation_distance = GEODESIC_INF; + + propagate(sources, + max_propagation_distance, + &destinations); //we use desinations as stop points + + paths.resize(destinations.size()); + + for(unsigned i=0; i* stop_points, + double stop_distance) +{ + m_max_propagation_distance = stop_distance; + + if(!stop_points) + { + m_stop_vertices.clear(); + return; + } + + m_stop_vertices.resize(stop_points->size()); + + std::vector possible_vertices; + for(unsigned i = 0; i < stop_points->size(); ++i) + { + SurfacePoint* point = &(*stop_points)[i]; + + possible_vertices.clear(); + m_mesh->closest_vertices(point, &possible_vertices); + + vertex_pointer closest_vertex = NULL; + double min_distance = 1e100; + for(unsigned j = 0; j < possible_vertices.size(); ++j) + { + double distance = point->distance(possible_vertices[j]); + if(distance < min_distance) + { + min_distance = distance; + closest_vertex = possible_vertices[j]; + } + } + assert(closest_vertex); + + m_stop_vertices[i].first = closest_vertex; + m_stop_vertices[i].second = min_distance; + } +} + + + +class GeodesicAlgorithmExact : public GeodesicAlgorithmBase +{ +public: + GeodesicAlgorithmExact(geodesic::Mesh* mesh): + GeodesicAlgorithmBase(mesh), + m_memory_allocator(mesh->edges().size(), mesh->edges().size()), + m_edge_interval_lists(mesh->edges().size()) + { + m_type = EXACT; + + for(unsigned i=0; iedges()[i]); + } + }; + + ~GeodesicAlgorithmExact(){}; + + void propagate(std::vector& sources, + double max_propagation_distance = GEODESIC_INF, //propagation algorithm stops after reaching the certain distance from the source + std::vector* stop_points = NULL); //or after ensuring that all the stop_points are covered + + void trace_back(SurfacePoint& destination, //trace back piecewise-linear path + std::vector& path); + + unsigned best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source + double& best_source_distance); + + void print_statistics(); + +private: + typedef std::set IntervalQueue; + + void update_list_and_queue(list_pointer list, + IntervalWithStop* candidates, //up to two candidates + unsigned num_candidates); + + unsigned compute_propagated_parameters(double pseudo_x, + double pseudo_y, + double d, //parameters of the interval + double start, + double end, //start/end of the interval + double alpha, //corner angle + double L, //length of the new edge + bool first_interval, //if it is the first interval on the edge + bool last_interval, + bool turn_left, + bool turn_right, + IntervalWithStop* candidates); //if it is the last interval on the edge + + void construct_propagated_intervals(bool invert, + edge_pointer edge, + face_pointer face, //constructs iNew from the rest of the data + IntervalWithStop* candidates, + unsigned& num_candidates, + interval_pointer source_interval); + + double compute_positive_intersection(double start, + double pseudo_x, + double pseudo_y, + double sin_alpha, + double cos_alpha); //used in construct_propagated_intervals + + unsigned intersect_intervals(interval_pointer zero, + IntervalWithStop* one); //intersecting two intervals with up to three intervals in the end + + interval_pointer best_first_interval(SurfacePoint& point, + double& best_total_distance, + double& best_interval_position, + unsigned& best_source_index); + + bool check_stop_conditions(unsigned& index); + + void clear() + { + m_memory_allocator.clear(); + m_queue.clear(); + for(unsigned i=0; iid()]; + }; + + void set_sources(std::vector& sources) + { + m_sources.initialize(sources); + } + + void initialize_propagation_data(); + + void list_edges_visible_from_source(MeshElementBase* p, + std::vector& storage); //used in initialization + + long visible_from_source(SurfacePoint& point); //used in backtracing + + void best_point_on_the_edge_set(SurfacePoint& point, + std::vector const& storage, + interval_pointer& best_interval, + double& best_total_distance, + double& best_interval_position); + + void possible_traceback_edges(SurfacePoint& point, + std::vector& storage); + + bool erase_from_queue(interval_pointer p); + + IntervalQueue m_queue; //interval queue + + MemoryAllocator m_memory_allocator; //quickly allocate and deallocate intervals + std::vector m_edge_interval_lists; //every edge has its interval data + + enum MapType {OLD, NEW}; //used for interval intersection + MapType map[5]; + double start[6]; + interval_pointer i_new[5]; + + unsigned m_queue_max_size; //used for statistics + unsigned m_iterations; //used for statistics + + SortedSources m_sources; +}; + +inline void GeodesicAlgorithmExact::best_point_on_the_edge_set(SurfacePoint& point, + std::vector const& storage, + interval_pointer& best_interval, + double& best_total_distance, + double& best_interval_position) +{ + best_total_distance = 1e100; + for(unsigned i=0; ifind_closest_point(&point, + offset, + distance, + interval); + + if(distance < best_total_distance) + { + best_interval = interval; + best_total_distance = distance; + best_interval_position = offset; + } + } +} + +inline void GeodesicAlgorithmExact::possible_traceback_edges(SurfacePoint& point, + std::vector& storage) +{ + storage.clear(); + + if(point.type() == VERTEX) + { + vertex_pointer v = static_cast(point.base_element()); + for(unsigned i=0; iadjacent_faces().size(); ++i) + { + face_pointer f = v->adjacent_faces()[i]; + storage.push_back(f->opposite_edge(v)); + } + } + else if(point.type() == EDGE) + { + edge_pointer e = static_cast(point.base_element()); + for(unsigned i=0; iadjacent_faces().size(); ++i) + { + face_pointer f = e->adjacent_faces()[i]; + + storage.push_back(f->next_edge(e,e->v0())); + storage.push_back(f->next_edge(e,e->v1())); + } + } + else + { + face_pointer f = static_cast(point.base_element()); + storage.push_back(f->adjacent_edges()[0]); + storage.push_back(f->adjacent_edges()[1]); + storage.push_back(f->adjacent_edges()[2]); + } +} + + +inline long GeodesicAlgorithmExact::visible_from_source(SurfacePoint& point) //negative if not visible +{ + assert(point.type() != UNDEFINED_POINT); + + if(point.type() == EDGE) + { + edge_pointer e = static_cast(point.base_element()); + list_pointer list = interval_list(e); + double position = std::min(point.distance(e->v0()), e->length()); + interval_pointer interval = list->covering_interval(position); + //assert(interval); + if(interval && interval->visible_from_source()) + { + return (long)interval->source_index(); + } + else + { + return -1; + } + } + else if(point.type() == FACE) + { + return -1; + } + else if(point.type() == VERTEX) + { + vertex_pointer v = static_cast(point.base_element()); + for(unsigned i=0; iadjacent_edges().size(); ++i) + { + edge_pointer e = v->adjacent_edges()[i]; + list_pointer list = interval_list(e); + + double position = e->v0()->id() == v->id() ? 0.0 : e->length(); + interval_pointer interval = list->covering_interval(position); + if(interval && interval->visible_from_source()) + { + return (long)interval->source_index(); + } + } + + return -1; + } + + assert(0); + return 0; +} + +inline double GeodesicAlgorithmExact::compute_positive_intersection(double start, + double pseudo_x, + double pseudo_y, + double sin_alpha, + double cos_alpha) +{ + assert(pseudo_y < 0); + + double denominator = sin_alpha*(pseudo_x - start) - cos_alpha*pseudo_y; + if(denominator<0.0) + { + return -1.0; + } + + double numerator = -pseudo_y*start; + + if(numerator < 1e-30) + { + return 0.0; + } + + if(denominator < 1e-30) + { + return -1.0; + } + + return numerator/denominator; +} + +inline void GeodesicAlgorithmExact::list_edges_visible_from_source(MeshElementBase* p, + std::vector& storage) +{ + assert(p->type() != UNDEFINED_POINT); + + if(p->type() == FACE) + { + face_pointer f = static_cast(p); + for(unsigned i=0; i<3; ++i) + { + storage.push_back(f->adjacent_edges()[i]); + } + } + else if(p->type() == EDGE) + { + edge_pointer e = static_cast(p); + storage.push_back(e); + } + else //VERTEX + { + vertex_pointer v = static_cast(p); + for(unsigned i=0; iadjacent_edges().size(); ++i) + { + storage.push_back(v->adjacent_edges()[i]); + } + + } +} + +inline bool GeodesicAlgorithmExact::erase_from_queue(interval_pointer p) +{ + if(p->min() < GEODESIC_INF/10.0)// && p->min >= queue->begin()->first) + { + assert(m_queue.count(p)<=1); //the set is unique + + IntervalQueue::iterator it = m_queue.find(p); + + if(it != m_queue.end()) + { + m_queue.erase(it); + return true; + } + } + + return false; +} + +inline unsigned GeodesicAlgorithmExact::intersect_intervals(interval_pointer zero, + IntervalWithStop* one) //intersecting two intervals with up to three intervals in the end +{ + assert(zero->edge()->id() == one->edge()->id()); + assert(zero->stop() > one->start() && zero->start() < one->stop()); + assert(one->min() < GEODESIC_INF/10.0); + + double const local_epsilon = SMALLEST_INTERVAL_RATIO*one->edge()->length(); + + unsigned N=0; + if(zero->min() > GEODESIC_INF/10.0) + { + start[0] = zero->start(); + if(zero->start() < one->start() - local_epsilon) + { + map[0] = OLD; + start[1] = one->start(); + map[1] = NEW; + N = 2; + } + else + { + map[0] = NEW; + N = 1; + } + + if(zero->stop() > one->stop() + local_epsilon) + { + map[N] = OLD; //"zero" interval + start[N++] = one->stop(); + } + + start[N+1] = zero->stop(); + return N; + } + + double const local_small_epsilon = 1e-8*one->edge()->length(); + + double D = zero->d() - one->d(); + double x0 = zero->pseudo_x(); + double x1 = one->pseudo_x(); + double R0 = x0*x0 + zero->pseudo_y()*zero->pseudo_y(); + double R1 = x1*x1 + one->pseudo_y()*one->pseudo_y(); + + double inter[2]; //points of intersection + char Ninter=0; //number of the points of the intersection + + if(std::abs(D)local_small_epsilon) + { + inter[0] = (R1 - R0)/(2.*denom); //one solution + Ninter = 1; + } + } + else + { + double D2 = D*D; + double Q = 0.5*(R1-R0-D2); + double X = x0 - x1; + + double A = X*X - D2; + double B = Q*X + D2*x0; + double C = Q*Q - D2*R0; + + if (std::abs(A)local_small_epsilon) + { + inter[0] = -C/B; //one solution + Ninter = 1; + } + } + else + { + double det = B*B-A*C; + if(det>local_small_epsilon*local_small_epsilon) //two roots + { + det = sqrt(det); + if(A>0.0) //make sure that the roots are ordered + { + inter[0] = (-B - det)/A; + inter[1] = (-B + det)/A; + } + else + { + inter[0] = (-B + det)/A; + inter[1] = (-B - det)/A; + } + + if(inter[1] - inter[0] > local_small_epsilon) + { + Ninter = 2; + } + else + { + Ninter = 1; + } + } + else if(det>=0.0) //single root + { + inter[0] = -B/A; + Ninter = 1; + } + } + } + //---------------------------find possible intervals--------------------------------------- + double left = std::max(zero->start(), one->start()); //define left and right boundaries of the intersection of the intervals + double right = std::min(zero->stop(), one->stop()); + + double good_start[4]; //points of intersection within the (left, right) limits +"left" + "right" + good_start[0] = left; + char Ngood_start=1; //number of the points of the intersection + + for(char i=0; i left + local_epsilon && x < right - local_epsilon) + { + good_start[Ngood_start++] = x; + } + } + good_start[Ngood_start++] = right; + + MapType mid_map[3]; + for(char i=0; isignal(mid) <= one->signal(mid) ? OLD : NEW; + } + + //-----------------------------------output---------------------------------- + N = 0; + if(zero->start() < left - local_epsilon) //additional "zero" interval + { + if(mid_map[0] == OLD) //first interval in the map is already the old one + { + good_start[0] = zero->start(); + } + else + { + map[N] = OLD; //"zero" interval + start[N++] = zero->start(); + } + } + + for(long i=0;istop() > one->stop() + local_epsilon) + { + if(N==0 || map[N-1] == NEW) + { + map[N] = OLD; //"zero" interval + start[N++] = one->stop(); + } + } + + start[0] = zero->start(); // just to make sure that epsilons do not damage anything + //start[N] = zero->stop(); + + return N; +} + +inline void GeodesicAlgorithmExact::initialize_propagation_data() +{ + clear(); + + IntervalWithStop candidate; + std::vector edges_visible_from_source; + for(unsigned i=0; ibase_element(), + edges_visible_from_source); + + for(unsigned j=0; jlength(); + candidate.compute_min_distance(candidate.stop()); + candidate.direction() = Interval::FROM_SOURCE; + + update_list_and_queue(interval_list(e), &candidate, 1); + } + } +} + +inline void GeodesicAlgorithmExact::propagate(std::vector& sources, + double max_propagation_distance, //propagation algorithm stops after reaching the certain distance from the source + std::vector* stop_points) +{ + set_stop_conditions(stop_points, max_propagation_distance); + set_sources(sources); + initialize_propagation_data(); + + clock_t start = clock(); + + unsigned satisfied_index = 0; + + m_iterations = 0; //for statistics + m_queue_max_size = 0; + + IntervalWithStop candidates[2]; + + while(!m_queue.empty()) + { + m_queue_max_size = std::max(static_cast(m_queue.size()), m_queue_max_size); + + unsigned const check_period = 10; + if(++m_iterations % check_period == 0) //check if we covered all required vertices + { + if (check_stop_conditions(satisfied_index)) + { + break; + } + } + + interval_pointer min_interval = *m_queue.begin(); + m_queue.erase(m_queue.begin()); + edge_pointer edge = min_interval->edge(); + list_pointer list = interval_list(edge); + + assert(min_interval->d() < GEODESIC_INF); + + bool const first_interval = min_interval->start() == 0.0; + //bool const last_interval = min_interval->stop() == edge->length(); + bool const last_interval = min_interval->next() == NULL; + + bool const turn_left = edge->v0()->saddle_or_boundary(); + bool const turn_right = edge->v1()->saddle_or_boundary(); + + for(unsigned i=0; iadjacent_faces().size(); ++i) //two possible faces to propagate + { + if(!edge->is_boundary()) //just in case, always propagate boundary edges + { + if((i == 0 && min_interval->direction() == Interval::FROM_FACE_0) || + (i == 1 && min_interval->direction() == Interval::FROM_FACE_1)) + { + continue; + } + } + + face_pointer face = edge->adjacent_faces()[i]; //if we come from 1, go to 2 + edge_pointer next_edge = face->next_edge(edge,edge->v0()); + + unsigned num_propagated = compute_propagated_parameters(min_interval->pseudo_x(), + min_interval->pseudo_y(), + min_interval->d(), //parameters of the interval + min_interval->start(), + min_interval->stop(), //start/end of the interval + face->vertex_angle(edge->v0()), //corner angle + next_edge->length(), //length of the new edge + first_interval, //if it is the first interval on the edge + last_interval, + turn_left, + turn_right, + candidates); //if it is the last interval on the edge + bool propagate_to_right = true; + + if(num_propagated) + { + if(candidates[num_propagated-1].stop() != next_edge->length()) + { + propagate_to_right = false; + } + + bool const invert = next_edge->v0()->id() != edge->v0()->id(); //if the origins coinside, do not invert intervals + + construct_propagated_intervals(invert, //do not inverse + next_edge, + face, + candidates, + num_propagated, + min_interval); + + update_list_and_queue(interval_list(next_edge), + candidates, + num_propagated); + } + + if(propagate_to_right) + { + //propogation to the right edge + double length = edge->length(); + next_edge = face->next_edge(edge,edge->v1()); + + num_propagated = compute_propagated_parameters(length - min_interval->pseudo_x(), + min_interval->pseudo_y(), + min_interval->d(), //parameters of the interval + length - min_interval->stop(), + length - min_interval->start(), //start/end of the interval + face->vertex_angle(edge->v1()), //corner angle + next_edge->length(), //length of the new edge + last_interval, //if it is the first interval on the edge + first_interval, + turn_right, + turn_left, + candidates); //if it is the last interval on the edge + + if(num_propagated) + { + bool const invert = next_edge->v0()->id() != edge->v1()->id(); //if the origins coinside, do not invert intervals + + construct_propagated_intervals(invert, //do not inverse + next_edge, + face, + candidates, + num_propagated, + min_interval); + + update_list_and_queue(interval_list(next_edge), + candidates, + num_propagated); + } + } + } + } + + m_propagation_distance_stopped = m_queue.empty() ? GEODESIC_INF : (*m_queue.begin())->min(); + clock_t stop = clock(); + m_time_consumed = (static_cast(stop)-static_cast(start))/CLOCKS_PER_SEC; + +/* for(unsigned i=0; ifirst(); + assert(p->start() == 0.0); + while(p->next()) + { + assert(p->stop() == p->next()->start()); + assert(p->d() < GEODESIC_INF); + p = p->next(); + } + }*/ +} + + +inline bool GeodesicAlgorithmExact::check_stop_conditions(unsigned& index) +{ + double queue_distance = (*m_queue.begin())->min(); + if(queue_distance < stop_distance()) + { + return false; + } + + while(index < m_stop_vertices.size()) + { + vertex_pointer v = m_stop_vertices[index].first; + edge_pointer edge = v->adjacent_edges()[0]; //take any edge + + double distance = edge->v0()->id() == v->id() ? + interval_list(edge)->signal(0.0) : + interval_list(edge)->signal(edge->length()); + + if(queue_distance < distance + m_stop_vertices[index].second) + { + return false; + } + + ++index; + } + return true; +} + + +inline void GeodesicAlgorithmExact::update_list_and_queue(list_pointer list, + IntervalWithStop* candidates, //up to two candidates + unsigned num_candidates) +{ + assert(num_candidates <= 2); + //assert(list->first() != NULL); + edge_pointer edge = list->edge(); + double const local_epsilon = SMALLEST_INTERVAL_RATIO * edge->length(); + + if(list->first() == NULL) + { + interval_pointer* p = &list->first(); + IntervalWithStop* first; + IntervalWithStop* second; + + if(num_candidates == 1) + { + first = candidates; + second = candidates; + first->compute_min_distance(first->stop()); + } + else + { + if(candidates->start() <= (candidates+1)->start()) + { + first = candidates; + second = candidates+1; + } + else + { + first = candidates+1; + second = candidates; + } + assert(first->stop() == second->start()); + + first->compute_min_distance(first->stop()); + second->compute_min_distance(second->stop()); + } + + if(first->start() > 0.0) + { + *p = m_memory_allocator.allocate(); + (*p)->initialize(edge); + p = &(*p)->next(); + } + + *p = m_memory_allocator.allocate(); + memcpy(*p,first,sizeof(Interval)); + m_queue.insert(*p); + + if(num_candidates == 2) + { + p = &(*p)->next(); + *p = m_memory_allocator.allocate(); + memcpy(*p,second,sizeof(Interval)); + m_queue.insert(*p); + } + + if(second->stop() < edge->length()) + { + p = &(*p)->next(); + *p = m_memory_allocator.allocate(); + (*p)->initialize(edge); + (*p)->start() = second->stop(); + } + else + { + (*p)->next() = NULL; + } + return; + } + + bool propagate_flag; + + for(unsigned i=0; ifirst(); + assert(p->start() == 0.0); + + while(p != NULL && p->stop() - local_epsilon < q->start()) + { + p = p->next(); + } + + while(p != NULL && p->start() < q->stop() - local_epsilon) //go through all old intervals + { + unsigned const N = intersect_intervals(p, q); //interset two intervals + + if(N == 1) + { + if(map[0]==OLD) //if "p" is always better, we do not need to update anything) + { + if(previous) //close previous interval and put in into the queue + { + previous->next() = p; + previous->compute_min_distance(p->start()); + m_queue.insert(previous); + previous = NULL; + } + + p = p->next(); + + } + else if(previous) //extend previous interval to cover everything; remove p + { + previous->next() = p->next(); + erase_from_queue(p); + m_memory_allocator.deallocate(p); + + p = previous->next(); + } + else //p becomes "previous" + { + previous = p; + interval_pointer next = p->next(); + erase_from_queue(p); + + memcpy(previous,q,sizeof(Interval)); + + previous->start() = start[0]; + previous->next() = next; + + p = next; + } + continue; + } + + //update_flag = true; + + Interval swap(*p); //used for swapping information + propagate_flag = erase_from_queue(p); + + for(unsigned j=1; jnext() = p; + previous->compute_min_distance(previous->stop()); + m_queue.insert(previous); + previous = NULL; + } + i_new[0] = p; + p->next() = i_new[1]; + p->start() = start[0]; + } + else if(previous) //extend previous interval to cover everything; remove p + { + i_new[0] = previous; + previous->next() = i_new[1]; + m_memory_allocator.deallocate(p); + previous = NULL; + } + else //p becomes "previous" + { + i_new[0] = p; + memcpy(p,q,sizeof(Interval)); + + p->next() = i_new[1]; + p->start() = start[0]; + } + + assert(!previous); + + for(unsigned j=1; jnext() = swap.next(); + } + else + { + current_interval->next() = i_new[j+1]; + } + + current_interval->start() = start[j]; + } + + for(unsigned j=0; jcompute_min_distance(current_interval->stop()); //compute minimal distance + + if(map[j]==NEW || (map[j]==OLD && propagate_flag)) + { + m_queue.insert(current_interval); + } + } + } + + p = swap.next(); + } + + if(previous) //close previous interval and put in into the queue + { + previous->compute_min_distance(previous->stop()); + m_queue.insert(previous); + previous = NULL; + } + } +} + +inline unsigned GeodesicAlgorithmExact::compute_propagated_parameters(double pseudo_x, + double pseudo_y, + double d, //parameters of the interval + double begin, + double end, //start/end of the interval + double alpha, //corner angle + double L, //length of the new edge + bool first_interval, //if it is the first interval on the edge + bool last_interval, + bool turn_left, + bool turn_right, + IntervalWithStop* candidates) //if it is the last interval on the edge +{ + assert(pseudo_y<=0.0); + assert(dstart() = 0.0; + p->stop() = L; + p->d() = d - pseudo_x; + p->pseudo_x() = 0.0; + p->pseudo_y() = 0.0; + return 1; + } + else if(last_interval && pseudo_x >= end) + { + p->start() = 0.0; + p->stop() = L; + p->d() = d + pseudo_x-end; + p->pseudo_x() = end*cos(alpha); + p->pseudo_y() = -end*sin(alpha); + return 1; + } + else if(pseudo_x >= begin && pseudo_x <= end) + { + p->start() = 0.0; + p->stop() = L; + p->d() = d; + p->pseudo_x() = pseudo_x*cos(alpha); + p->pseudo_y() = -pseudo_x*sin(alpha); + return 1; + } + else + { + return 0; + } + } + + double sin_alpha = sin(alpha); + double cos_alpha = cos(alpha); + + //important: for the first_interval, this function returns zero only if the new edge is "visible" from the source + //if the new edge can be covered only after turn_over, the value is negative (-1.0) + double L1 = compute_positive_intersection(begin, + pseudo_x, + pseudo_y, + sin_alpha, + cos_alpha); + + if(L1 < 0 || L1 >= L) + { + if(first_interval && turn_left) + { + p->start() = 0.0; + p->stop() = L; + p->d() = d + sqrt(pseudo_x*pseudo_x + pseudo_y*pseudo_y); + p->pseudo_y() = 0.0; + p->pseudo_x() = 0.0; + return 1; + } + else + { + return 0; + } + } + + double L2 = compute_positive_intersection(end, + pseudo_x, + pseudo_y, + sin_alpha, + cos_alpha); + + if(L2 < 0 || L2 >= L) + { + p->start() = L1; + p->stop() = L; + p->d() = d; + p->pseudo_x() = cos_alpha*pseudo_x + sin_alpha*pseudo_y; + p->pseudo_y() = -sin_alpha*pseudo_x + cos_alpha*pseudo_y; + + return 1; + } + + p->start() = L1; + p->stop() = L2; + p->d() = d; + p->pseudo_x() = cos_alpha*pseudo_x + sin_alpha*pseudo_y; + p->pseudo_y() = -sin_alpha*pseudo_x + cos_alpha*pseudo_y; + assert(p->pseudo_y() <= 0.0); + + if(!(last_interval && turn_right)) + { + return 1; + } + else + { + p = candidates + 1; + + p->start() = L2; + p->stop() = L; + double dx = pseudo_x - end; + p->d() = d + sqrt(dx*dx + pseudo_y*pseudo_y); + p->pseudo_x() = end*cos_alpha; + p->pseudo_y() = -end*sin_alpha; + + return 2; + } +} + +inline void GeodesicAlgorithmExact::construct_propagated_intervals(bool invert, + edge_pointer edge, + face_pointer face, //constructs iNew from the rest of the data + IntervalWithStop* candidates, + unsigned& num_candidates, + interval_pointer source_interval) //up to two candidates +{ + double edge_length = edge->length(); + double local_epsilon = SMALLEST_INTERVAL_RATIO * edge_length; + + //kill very small intervals in order to avoid precision problems + if(num_candidates == 2) + { + double start = std::min(candidates->start(), (candidates+1)->start()); + double stop = std::max(candidates->stop(), (candidates+1)->stop()); + if(candidates->stop()-candidates->start() < local_epsilon) // kill interval 0 + { + *candidates = *(candidates+1); + num_candidates = 1; + candidates->start() = start; + candidates->stop() = stop; + } + else if ((candidates+1)->stop() - (candidates+1)->start() < local_epsilon) + { + num_candidates = 1; + candidates->start() = start; + candidates->stop() = stop; + } + } + + IntervalWithStop* first; + IntervalWithStop* second; + if(num_candidates == 1) + { + first = candidates; + second = candidates; + } + else + { + if(candidates->start() <= (candidates+1)->start()) + { + first = candidates; + second = candidates+1; + } + else + { + first = candidates+1; + second = candidates; + } + assert(first->stop() == second->start()); + } + + if(first->start() < local_epsilon) + { + first->start() = 0.0; + } + if(edge_length - second->stop() < local_epsilon) + { + second->stop() = edge_length; + } + + //invert intervals if necessary; fill missing data and set pointers correctly + Interval::DirectionType direction = edge->adjacent_faces()[0]->id() == face->id() ? + Interval::FROM_FACE_0 : + Interval::FROM_FACE_1; + + if(!invert) //in this case everything is straighforward, we do not have to invert the intervals + { + for(unsigned i=0; inext() = (i == num_candidates - 1) ? NULL : candidates + i + 1; + p->edge() = edge; + p->direction() = direction; + p->source_index() = source_interval->source_index(); + + p->min() = 0.0; //it will be changed later on + + assert(p->start() < p->stop()); + } + } + else //now we have to invert the intervals + { + for(unsigned i=0; inext() = (i == 0) ? NULL : candidates + i - 1; + p->edge() = edge; + p->direction() = direction; + p->source_index() = source_interval->source_index(); + + double length = edge_length; + p->pseudo_x() = length - p->pseudo_x(); + + double start = length - p->stop(); + p->stop() = length - p->start(); + p->start() = start; + + p->min() = 0; + + assert(p->start() < p->stop()); + assert(p->start() >= 0.0); + assert(p->stop() <= edge->length()); + } + } +} + + +inline unsigned GeodesicAlgorithmExact::best_source(SurfacePoint& point, //quickly find what source this point belongs to and what is the distance to this source + double& best_source_distance) +{ + double best_interval_position; + unsigned best_source_index; + + best_first_interval(point, + best_source_distance, + best_interval_position, + best_source_index); + + return best_source_index; +} + +inline interval_pointer GeodesicAlgorithmExact::best_first_interval(SurfacePoint& point, + double& best_total_distance, + double& best_interval_position, + unsigned& best_source_index) +{ + assert(point.type() != UNDEFINED_POINT); + + interval_pointer best_interval = NULL; + best_total_distance = GEODESIC_INF; + + if(point.type() == EDGE) + { + edge_pointer e = static_cast(point.base_element()); + list_pointer list = interval_list(e); + + best_interval_position = point.distance(e->v0()); + best_interval = list->covering_interval(best_interval_position); + if(best_interval) + { + //assert(best_interval && best_interval->d() < GEODESIC_INF); + best_total_distance = best_interval->signal(best_interval_position); + best_source_index = best_interval->source_index(); + } + } + else if(point.type() == FACE) + { + face_pointer f = static_cast(point.base_element()); + for(unsigned i=0; i<3; ++i) + { + edge_pointer e = f->adjacent_edges()[i]; + list_pointer list = interval_list(e); + + double offset; + double distance; + interval_pointer interval; + + list->find_closest_point(&point, + offset, + distance, + interval); + + if(interval && distance < best_total_distance) + { + best_interval = interval; + best_total_distance = distance; + best_interval_position = offset; + best_source_index = interval->source_index(); + } + } + + //check for all sources that might be located inside this face + SortedSources::sorted_iterator_pair local_sources = m_sources.sources(f); + for(SortedSources::sorted_iterator it=local_sources.first; it != local_sources.second; ++it) + { + SurfacePointWithIndex* source = *it; + double distance = point.distance(source); + if(distance < best_total_distance) + { + best_interval = NULL; + best_total_distance = distance; + best_interval_position = 0.0; + best_source_index = source->index(); + } + } + } + else if(point.type() == VERTEX) + { + vertex_pointer v = static_cast(point.base_element()); + for(unsigned i=0; iadjacent_edges().size(); ++i) + { + edge_pointer e = v->adjacent_edges()[i]; + list_pointer list = interval_list(e); + + double position = e->v0()->id() == v->id() ? 0.0 : e->length(); + interval_pointer interval = list->covering_interval(position); + if(interval) + { + double distance = interval->signal(position); + + if(distance < best_total_distance) + { + best_interval = interval; + best_total_distance = distance; + best_interval_position = position; + best_source_index = interval->source_index(); + } + } + } + } + + if(best_total_distance > m_propagation_distance_stopped) //result is unreliable + { + best_total_distance = GEODESIC_INF; + return NULL; + } + else + { + return best_interval; + } +} + +inline void GeodesicAlgorithmExact::trace_back(SurfacePoint& destination, //trace back piecewise-linear path + std::vector& path) +{ + path.clear(); + double best_total_distance; + double best_interval_position; + unsigned source_index = std::numeric_limits::max(); + interval_pointer best_interval = best_first_interval(destination, + best_total_distance, + best_interval_position, + source_index); + + if(best_total_distance >= GEODESIC_INF/2.0) //unable to find the right path + { + return; + } + + path.push_back(destination); + + if(best_interval) //if we did not hit the face source immediately + { + std::vector possible_edges; + possible_edges.reserve(10); + + while(visible_from_source(path.back()) < 0) //while this point is not in the direct visibility of some source (if we are inside the FACE, we obviously hit the source) + { + SurfacePoint& q = path.back(); + + possible_traceback_edges(q, possible_edges); + + interval_pointer interval; + double total_distance; + double position; + + best_point_on_the_edge_set(q, + possible_edges, + interval, + total_distance, + position); + + //std::cout << total_distance + length(path) << std::endl; + assert(total_distancesource_index(); + + edge_pointer e = interval->edge(); + double local_epsilon = SMALLEST_INTERVAL_RATIO*e->length(); + if(position < local_epsilon) + { + path.push_back(SurfacePoint(e->v0())); + } + else if(position > e->length()-local_epsilon) + { + path.push_back(SurfacePoint(e->v1())); + } + else + { + double normalized_position = position/e->length(); + path.push_back(SurfacePoint(e, normalized_position)); + } + } + } + + SurfacePoint& source = static_cast(m_sources[source_index]); + if(path.back().distance(&source) > 0) + { + path.push_back(source); + } +} + +inline void GeodesicAlgorithmExact::print_statistics() +{ + GeodesicAlgorithmBase::print_statistics(); + + unsigned interval_counter = 0; + for(unsigned i=0; i +IGL_INLINE void igl::exact_geodesic( + const Eigen::MatrixBase &V, + const Eigen::MatrixBase &F, + const Eigen::MatrixBase &VS, + const Eigen::MatrixBase &FS, + const Eigen::MatrixBase &VT, + const Eigen::MatrixBase &FT, + Eigen::PlainObjectBase &D) +{ + assert(V.cols() == 3 && F.cols() == 3 && "Only support 3D triangle mesh"); + assert(VS.cols() ==1 && FS.cols() == 1 && VT.cols() == 1 && FT.cols() ==1 && "Only support one dimensional inputs"); + std::vector points(V.rows() * V.cols()); + std::vector faces(F.rows() * F.cols()); + for (int i = 0; i < points.size(); i++) + { + points[i] = V(i / 3, i % 3); + } + for (int i = 0; i < faces.size(); i++) + { + faces[i] = F(i / 3, i % 3); + } + + igl::geodesic::Mesh mesh; + mesh.initialize_mesh_data(points, faces); + igl::geodesic::GeodesicAlgorithmExact exact_algorithm(&mesh); + + std::vector source(VS.rows() + FS.rows()); + std::vector target(VT.rows() + FT.rows()); + for (int i = 0; i < VS.rows(); i++) + { + source[i] = (igl::geodesic::SurfacePoint(&mesh.vertices()[VS(i)])); + } + for (int i = 0; i < FS.rows(); i++) + { + source[i] = (igl::geodesic::SurfacePoint(&mesh.faces()[FS(i)])); + } + + for (int i = 0; i < VT.rows(); i++) + { + target[i] = (igl::geodesic::SurfacePoint(&mesh.vertices()[VT(i)])); + } + for (int i = 0; i < FT.rows(); i++) + { + target[i] = (igl::geodesic::SurfacePoint(&mesh.faces()[FT(i)])); + } + + exact_algorithm.propagate(source); + std::vector path; + D.resize(target.size(), 1); + for (int i = 0; i < target.size(); i++) + { + exact_algorithm.trace_back(target[i], path); + D(i) = igl::geodesic::length(path); + } +} + +#ifdef IGL_STATIC_LIBRARY +template void igl::exact_geodesic, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); +#endif diff --git a/include/igl/exact_geodesic.h b/include/igl/exact_geodesic.h new file mode 100644 index 000000000..cec59de69 --- /dev/null +++ b/include/igl/exact_geodesic.h @@ -0,0 +1,55 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2017 Zhongshi Jiang +// +// 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/. + +#ifndef IGL_EXACT_GEODESIC_H +#define IGL_EXACT_GEODESIC_H + +#include "igl_inline.h" +#include + +namespace igl +{ + // Exact geodesic algorithm for triangular mesh with the implementation from https://code.google.com/archive/p/geodesic/, + // and the algorithm first described by Mitchell, Mount and Papadimitriou in 1987 + // + // Inputs: + // V #V by 3 list of 3D vertex positions + // F #F by 3 list of mesh faces + // VS #VS by 1 vector specifying indices of source vertices + // FS #FS by 1 vector specifying indices of source faces + // VT #VT by 1 vector specifying indices of target vertices + // FT #FT by 1 vector specifying indices of target faces + // Output: + // D #VT+#FT by 1 vector of geodesic distances of each target w.r.t. the nearest one in the source set + // + // Note: + // Specifying a face as target/source means its center. + // + template < + typename DerivedV, + typename DerivedF, + typename DerivedVS, + typename DerivedFS, + typename DerivedVT, + typename DerivedFT, + typename DerivedD> + IGL_INLINE void exact_geodesic( + const Eigen::MatrixBase &V, + const Eigen::MatrixBase &F, + const Eigen::MatrixBase &VS, + const Eigen::MatrixBase &FS, + const Eigen::MatrixBase &VT, + const Eigen::MatrixBase &FT, + Eigen::PlainObjectBase &D); +} + +#ifndef IGL_STATIC_LIBRARY +# include "exact_geodesic.cpp" +#endif + +#endif \ No newline at end of file diff --git a/index.html b/index.html index a6d95df4f..e9c4fc5da 100644 --- a/index.html +++ b/index.html @@ -1,14 +1,11 @@ - + libigl - - - - + @@ -22,10 +19,10 @@

Get started with:

+
git clone --recursive https://github.com/libigl/libigl.git
-

libigl is a simple C++ geometry processing library. We have a wide functionality including construction of sparse discrete differential geometry @@ -39,7 +36,7 @@ like MATLAB.

just include igl headers (e.g. #include <igl/cotmatrix.h>) and run. Each header file contains a single function (e.g. igl/cotmatrix.h contains igl::cotmatrix()). Most are tailored to operate on a generic triangle mesh -stored in an n-by–3 matrix of vertex positions V and an m-by–3 matrix of +stored in an n-by-3 matrix of vertex positions V and an m-by-3 matrix of triangle indices F.

Optionally the library may also be pre-compiled into a statically @@ -121,14 +118,14 @@ libigl depends only on the Eigen librar

Libigl compartmentalizes its optional dependences via its directory organization in the include/ folder. All header files located directly in the include/igl/ folder have only stl and Eigen as dependencies. For example, -all of the headers that depend on CGAL are located in include/igl/cgal. For a -full list of optional dependencies check optional/CMakeLists.txt.

+all of the headers that depend on CGAL are located in include/igl/copyleft/cgal. +For a full list of optional dependencies check optional/CMakeLists.txt.

GCC and the Optional CGAL Dependency

-

The include/igl/cgal/*.h headers depend on CGAL. It has come to our attention -that CGAL does not work properly with GCC 4.8. To the best of our knowledge, -GCC 4.7 and clang will work correctly.

+

The include/igl/copyleft/cgal/*.h headers depend on CGAL. It has come to +our attention that CGAL does not work properly with GCC 4.8. To the best of +our knowledge, GCC 4.7 and clang will work correctly.

OpenMP and Windows

@@ -210,7 +207,7 @@ BibTeX entry:

title = {{libigl}: A simple {C++} geometry processing library}, author = {Alec Jacobson and Daniele Panozzo and others}, note = {http://libigl.github.io/libigl/}, - year = {2016}, + year = {2017}, } @@ -221,56 +218,59 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a few labs/companies/institutions using libigl:

Contact

Libigl is a group endeavor led by Alec Jacobson and Daniele -Panozzo. Please contact -us if you have +Panozzo. Please contact +us if you have questions or comments. For troubleshooting, please post an issue on github.

-

If you’re using libigl in your projects, quickly drop us a -note. Tell us who you +

If you’re using libigl in your projects, quickly drop us a +note. Tell us who you are and what you’re using it for. This helps us apply for funding and justify spending time maintaining this.

@@ -279,11 +279,10 @@ page.

-

2017 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan -Zhou, Sebastian Koch, Amir Vaxman, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De -Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.

+

2017 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan Zhou, Sebastian Koch, Jeremie Dumas, Amir Vaxman, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.

Please see individual files for appropriate copyright notices.

+ diff --git a/optional/index.html b/optional/index.html index fe9fa820a..51b7168e1 100644 --- a/optional/index.html +++ b/optional/index.html @@ -1,14 +1,11 @@ - + libigl - - - - + @@ -24,7 +21,7 @@ improve your compilation times.

Libigl is developed most often on Mac OS X, though has current users in Linux and Windows.

-

Linux/Mac OS X/Cygwin

+

Linux/Mac OS X/Cygwin

Libigl may also be compiled to a static library. This is advantageous when building a project with libigl, since when used as an header-only library can @@ -40,7 +37,7 @@ cmake -DCMAKE_BUILD_TYPE=Release ../optional make -

Warnings

+

Warnings

You should expect to see a few linker warnings of the form:

@@ -50,7 +47,7 @@ make

These are (admittedly unpopular) functions that have never been used by us statically so we haven’t explicit instantiations (yet).

-

External

+

External

Finally there are a number of external libraries that we include in ./external/ because they are either difficult to obtain or they have been @@ -58,7 +55,7 @@ patched for easier use with libigl. Please see the respective readmes in those directories or build the tutorial using cmake, which will recursively build all dependencies.

-
Installing Embree 2.0
+
Installing Embree 2.0

To build the embree library and executables on Mac OS X issue:

@@ -74,72 +71,72 @@ make #sudo make install -

Extras

+

Extras

-

bbw

+

bbw

This library extra contains functions for computing Bounded Biharmonic Weights, can be used with and without the mosek extra via the IGL_NO_MOSEK macro.

-

boolean

+

boolean

This library extra contains functions for computing mesh-mesh booleans, depending on CGAL and optionally Cork.

-

cgal

+

cgal

This library extra utilizes CGAL’s efficient and exact intersection and proximity queries.

-

embree

+

embree

This library extra utilizes embree’s efficient ray tracing queries.

-

matlab

+

matlab

This library extra provides support for reading and writing .mat workspace files, interfacing with Matlab at run time and compiling mex functions.

-

mosek

+

mosek

This library extra utilizes mosek’s efficient interior-point solver for quadratic programs.

-

png

+

png

This library extra uses libpng and YImage to read and write .png files.

-

tetgen

+

tetgen

This library extra provides a simplified wrapper to the tetgen 3d tetrahedral meshing library.

-

Triangle

+

Triangle

This library extra provides a simplified wrapper to the triangle 2d triangle meshing library.

-

viewer

+

viewer

This library extra utilizes glfw and glew to open an opengl context and launch a simple mesh viewer.

-

xml

+

xml

This library extra utilizes tinyxml2 to read and write serialized classes containing Eigen matrices and other standard simple data-structures.

-

Development

+

Development

Further documentation for developers is listed in style_guidelines.html.

-

License

+

License

See LICENSE.txt

-

Zipping

+

Zipping

Zip this directory without .git litter and binaries using:

@@ -275,28 +272,28 @@ exposing templated functions.
  • Easy incorporation: This can be easily incorporated - into external projects.
    • +into external projects.

Drawbacks of compressed .h/.cpp pair

  • Hard to debug/edit: The compressed files are - automatically generated. They’re huge and should not be edited. Thus - debugging and editing are near impossible.

    • +automatically generated. They’re huge and should not be edited. Thus +debugging and editing are near impossible.

  • Compounded dependencies: - An immediate disadvantage of this - seems to be that even to use a single function (e.g. - cotmatrix), compiling and linking against - igl.cpp will require linking to all of libigl’s - dependencies (OpenGL, GLUT, - AntTweakBar, BLAS). However, because all - dependencies other than Eigen should be encapsulated between - #ifndef guards (e.g. #ifndef IGL_NO_OPENGL, it - is possible to ignore certain functions that have such dependencies.

    • +An immediate disadvantage of this +seems to be that even to use a single function (e.g. +cotmatrix), compiling and linking against +igl.cpp will require linking to all of libigl’s +dependencies (OpenGL, GLUT, +AntTweakBar, BLAS). However, because all +dependencies other than Eigen should be encapsulated between +#ifndef guards (e.g. #ifndef IGL_NO_OPENGL, it +is possible to ignore certain functions that have such dependencies.

  • Long compile: - Compiling igl.cpp takes a long time and isn’t easily parallelized (no make - -j12 equivalent).

    • +Compiling igl.cpp takes a long time and isn’t easily parallelized (no make +-j12 equivalent).

Here’s a tiny test example using igl.h and igl.cpp. Save the following in test.cpp:

@@ -334,7 +331,7 @@ return (argc>=2 && igl::read_triangle_mesh(argv[1],V,F)?0:1); 
grep OpenGL `grep -L IGL_NO_OPENGL include/igl/*`
-

Optional

+

Optional

  • OpenGL (disable with IGL_NO_OPENGL) @@ -343,7 +340,7 @@ return (argc>=2 && igl::read_triangle_mesh(argv[1],V,F)?0:1);
  • OpenGL >= 4 (enable with IGL_OPENGL_4)
  • AntTweakBar (disable with IGL_NO_ANTTWEAKBAR) Last tested 1.16 (see - libigl/external/AntTweakBar)
    • +libigl/external/AntTweakBar)
  • GLEW Windows and Linux
  • OpenMP
  • libpng libiglpng extra only
    • @@ -361,7 +358,7 @@ return (argc>=2 && igl::read_triangle_mesh(argv[1],V,F)?0:1);
  • CoMiSo libcomiso extra only
    • -

    Optional (included in external/)

    +

    Optional (included in external/)

    • TetGen libigltetgen extra only
      • @@ -373,3 +370,4 @@ return (argc>=2 && igl::read_triangle_mesh(argv[1],V,F)?0:1); + diff --git a/python/iglhelpers.py b/python/iglhelpers.py index bf8d9c1b4..0695bb533 100644 --- a/python/iglhelpers.py +++ b/python/iglhelpers.py @@ -11,15 +11,15 @@ import pyigl as igl def p2e(m): if isinstance(m, np.ndarray): - if not m.flags['C_CONTIGUOUS']: - raise TypeError("p2e only support C-contiguous order") - if m.dtype.type == np.int32: - return igl.eigen.MatrixXi(m) - elif m.dtype.type == np.float64: - return igl.eigen.MatrixXd(m) + if not (m.flags['C_CONTIGUOUS'] or m.flags['F_CONTIGUOUS']): + raise TypeError('p2e support either c-order or f-order') + if m.dtype.type in [np.int32, np.int64]: + return igl.eigen.MatrixXi(m.astype(np.int32)) + elif m.dtype.type in [np.float64, np.float32]: + return igl.eigen.MatrixXd(m.astype(np.float64)) elif m.dtype.type == np.bool: return igl.eigen.MatrixXb(m) - raise TypeError("p2e only support dtype float64, int32 and bool") + raise TypeError("p2e only support dtype float64/32, int64/32 and bool") if sparse.issparse(m): # convert in a dense matrix with triples coo = m.tocoo() diff --git a/python/py_doc.cpp b/python/py_doc.cpp index c5f255baf..7174e51e1 100644 --- a/python/py_doc.cpp +++ b/python/py_doc.cpp @@ -598,6 +598,23 @@ const char *__doc_igl_embree_reorient_facets_raycast = R"igl_Qu8mg5v7(// Orient // Outputs: // I #F list of whether face has been flipped // C #F list of patch ID (output of bfs_orient > manifold patches))igl_Qu8mg5v7"; +const char *__doc_igl_exact_geodesic = R"igl_Qu8mg5v7( + // Exact geodesic algorithm for triangular mesh with the implementation from https://code.google.com/archive/p/geodesic/, + // and the algorithm first described by Mitchell, Mount and Papadimitriou in 1987 + // + // Inputs: + // V #V by 3 list of 3D vertex positions + // F #F by 3 list of mesh faces + // VS #VS by 1 vector specifying indices of source vertices + // FS #FS by 1 vector specifying indices of source faces + // VT #VT by 1 vector specifying indices of target vertices + // FT #FT by 1 vector specifying indices of target faces + // Output: + // D #VT+#FT by 1 vector of geodesic distances of each target w.r.t. the nearest one in the source set + // + // Note: + // Specifying a face as target/source means its center. + //)igl_Qu8mg5v7"; const char *__doc_igl_find_cross_field_singularities = R"igl_Qu8mg5v7(// Inputs: // V #V by 3 eigen Matrix of mesh vertex 3D positions // F #F by 3 eigen Matrix of face (quad) indices @@ -1453,3 +1470,33 @@ const char *__doc_igl_readPLY= R"igl_Qu8mg5v7(// Read a mesh from an ascii ply f // N double matrix of corner normals #N by 3 // UV #V by 2 texture coordinates // Returns true on success, false on errors)igl_Qu8mg5v7"; +const char *__doc_igl_seam_edges=R"igl_Qu8mg5v7(// Finds all UV-space boundaries of a mesh. + // + // Inputs: + // V #V by dim list of positions of the input mesh. + // TC #TC by 2 list of 2D texture coordinates of the input mesh + // F #F by 3 list of triange indices into V representing a + // manifold-with-boundary triangle mesh + // FTC #F by 3 list of indices into TC for each corner + // Outputs: + // seams Edges where the forwards and backwards directions have different + // texture coordinates, as a #seams-by-4 matrix of indices. Each row is + // organized as [ forward_face_index, forward_face_vertex_index, + // backwards_face_index, backwards_face_vertex_index ] such that one side + // of the seam is the edge: + // F[ seams( i, 0 ), seams( i, 1 ) ], F[ seams( i, 0 ), (seams( i, 1 ) + 1) % 3 ] + // and the other side is the edge: + // F[ seams( i, 2 ), seams( i, 3 ) ], F[ seams( i, 2 ), (seams( i, 3 ) + 1) % 3 ] + // boundaries Edges with only one incident triangle, as a #boundaries-by-2 + // matrix of indices. Each row is organized as + // [ face_index, face_vertex_index ] + // such that the edge is: + // F[ boundaries( i, 0 ), boundaries( i, 1 ) ], F[ boundaries( i, 0 ), (boundaries( i, 1 ) + 1) % 3 ] + // foldovers Edges where the two incident triangles fold over each other + // in UV-space, as a #foldovers-by-4 matrix of indices. + // Each row is organized as [ forward_face_index, forward_face_vertex_index, + // backwards_face_index, backwards_face_vertex_index ] + // such that one side of the foldover is the edge: + // F[ foldovers( i, 0 ), foldovers( i, 1 ) ], F[ foldovers( i, 0 ), (foldovers( i, 1 ) + 1) % 3 ] + // and the other side is the edge: + // F[ foldovers( i, 2 ), foldovers( i, 3 ) ], F[ foldovers( i, 2 ), (foldovers( i, 3 ) + 1) % 3 ])igl_Qu8mg5v7"; diff --git a/python/py_doc.h b/python/py_doc.h index 3f7678dc0..2e76267f0 100644 --- a/python/py_doc.h +++ b/python/py_doc.h @@ -48,6 +48,7 @@ extern const char *__doc_igl_eigs; extern const char *__doc_igl_embree_ambient_occlusion; extern const char *__doc_igl_embree_line_mesh_intersection; extern const char *__doc_igl_embree_reorient_facets_raycast; +extern const char *__doc_igl_exact_geodesic; extern const char *__doc_igl_find_cross_field_singularities; extern const char *__doc_igl_fit_rotations; extern const char *__doc_igl_fit_rotations_planar; @@ -99,6 +100,7 @@ extern const char *__doc_igl_readTGF; extern const char *__doc_igl_read_triangle_mesh; extern const char *__doc_igl_remove_duplicate_vertices; extern const char *__doc_igl_rotate_vectors; +extern const char *__doc_igl_seam_edges; extern const char *__doc_igl_setdiff; extern const char *__doc_igl_signed_distance; extern const char *__doc_igl_signed_distance_pseudonormal; @@ -125,4 +127,4 @@ extern const char *__doc_igl_winding_number_2; extern const char *__doc_igl_writeMESH; extern const char *__doc_igl_writeOBJ; extern const char *__doc_igl_writePLY; -extern const char *__doc_igl_readPLY; +extern const char *__doc_igl_readPLY; \ No newline at end of file diff --git a/python/py_igl.cpp b/python/py_igl.cpp index 5eaf9a977..5ec827282 100644 --- a/python/py_igl.cpp +++ b/python/py_igl.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -102,6 +103,7 @@ #include #include #include +#include void python_export_igl(py::module &m) { @@ -141,6 +143,7 @@ void python_export_igl(py::module &m) #include "py_igl/py_edge_lengths.cpp" #include "py_igl/py_edge_topology.cpp" #include "py_igl/py_eigs.cpp" +#include "py_igl/py_exact_geodesic.cpp" #include "py_igl/py_find_cross_field_singularities.cpp" #include "py_igl/py_fit_rotations.cpp" #include "py_igl/py_floor.cpp" @@ -199,4 +202,5 @@ void python_export_igl(py::module &m) #include "py_igl/py_writeOBJ.cpp" #include "py_igl/py_writePLY.cpp" #include "py_igl/py_readPLY.cpp" +#include "py_igl/py_seam_edges.cpp" } diff --git a/python/py_igl/py_exact_geodesic.cpp b/python/py_igl/py_exact_geodesic.cpp new file mode 100644 index 000000000..b05ac744b --- /dev/null +++ b/python/py_igl/py_exact_geodesic.cpp @@ -0,0 +1,24 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. + + +m.def("exact_geodesic", [] +( + const Eigen::MatrixXd &V, + const Eigen::MatrixXi &F, + const Eigen::MatrixXi &VS, + const Eigen::MatrixXi &FS, + const Eigen::MatrixXi &VT, + const Eigen::MatrixXi &FT, + Eigen::MatrixXd &D +) +{ + return igl::exact_geodesic(V, F, VS,FS,VT,FT, D); +}, __doc_igl_exact_geodesic, +py::arg("V"), py::arg("F"), py::arg("VS"), py::arg("FS"), py::arg("VT"), py::arg("FT"), py::arg("D")); + diff --git a/python/py_igl/py_seam_edges.cpp b/python/py_igl/py_seam_edges.cpp new file mode 100755 index 000000000..c3464012b --- /dev/null +++ b/python/py_igl/py_seam_edges.cpp @@ -0,0 +1,21 @@ +m.def("seam_edges", [] +( + const Eigen::MatrixXd& V, + const Eigen::MatrixXd& TC, + const Eigen::MatrixXi& F, + const Eigen::MatrixXi& FTC, + Eigen::MatrixXi& seams, + Eigen::MatrixXi& boundaries, + Eigen::MatrixXi& foldovers +) +{ +return igl::seam_edges( V, TC, F, FTC, seams, boundaries, foldovers); +}, __doc_igl_seam_edges, +py::arg("V"), +py::arg("TC"), +py::arg("F"), +py::arg("FTC"), +py::arg("seams"), +py::arg("boundaries"), +py::arg("foldovers")); + diff --git a/style-guidelines.html b/style-guidelines.html index 316e7588f..e81d7af75 100644 --- a/style-guidelines.html +++ b/style-guidelines.html @@ -1,14 +1,11 @@ - + libigl - - - - + @@ -18,7 +15,7 @@ style guidelines for developers of the library, but also acts as best-practices for users.

      -

      One function, one .h/.cpp pair

      +

      One function, one .h/.cpp pair

      The structure of libigl is very flat and function-based. For every function/sub-routine, create a single .h and .cpp file. For example, if you have @@ -61,7 +58,7 @@ namespace igl } #ifndef IGL_STATIC_LIBRARY -# include "example_fun.cpp" +#include "example_fun.cpp" #endif #endif @@ -185,21 +182,21 @@ Eigen::SparseMatrix<Atype> adjacency_matrix(const ... & F);

      Templating with Eigen

      Functions taking Eigen dense matrices/arrays as inputs and outputs (but not -return arguments), should template on top of Eigen::PlainObjectBase. Each +return arguments), should template on top of Eigen::MatrixBase. Each parameter should be derived using its own template.

      For example,

      template <typename DerivedV, typename DerivedF, typename DerivedBC>
 void barycenter(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
-  const Eigen::PlainObjectBase<DerivedBC> & BC);
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedBC> & BC);

      The Derived* template encodes the scalar type (e.g. double, int), the number of rows and cols at compile time, and the data storage (Row-major vs. -column-major).

      +column-major).

      Returning Eigen types is discouraged. In cases where the size and scalar type are a fixed and matching function of an input Derived* template, then @@ -219,7 +216,7 @@ output-argument version and call that. So a full implementation looks like:

      template <typename DerivedV, typename DerivedW>
 void fit_to_unit_cube(
-  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedV> & V,
   Eigen::PlainObjectBase<DerivedW> & W);
 template <typename DerivedV>
 void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase<DerivedV> & V);
@@ -229,7 +226,7 @@ void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase<DerivedV> &
 
 
      template <typename DerivedV, typename DerivedW>
 void fit_to_unit_cube(
-  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedV> & V,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
   W = (V.rowwise()-V.colwise().minCoeff()).array() /
@@ -237,7 +234,7 @@ void fit_to_unit_cube(
 }
 
 template <typename DerivedV>
-void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase<DerivedV> & V)
+void DerivedV fit_to_unit_cube(const Eigen::MatrixBase<DerivedV> & V)
 {
   DerivedV W;
   fit_to_unit_cube(V,W);
@@ -274,7 +271,7 @@ to mean a list of faces/triangles.
      
 
      Classes should be avoided. When naming a class use CamelCase (e.g.
 SortableRow.h).
      
 
-
      Enum naming conversion
      
+
      Enum naming conversion
      
 
 
      Enums types should be placed in the appropriate igl:: namespace and should be
 named in CamelCase (e.g. igl::SolverStatus) and instances should be named in
@@ -392,3 +389,4 @@ implementation so we’re keeping it as long as possible and profitable.
      
 
 
 
+
diff --git a/tutorial/206_GeodesicDistance/CMakeLists.txt b/tutorial/206_GeodesicDistance/CMakeLists.txt
new file mode 100644
index 000000000..97a9d6786
--- /dev/null
+++ b/tutorial/206_GeodesicDistance/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(206_GeodesicDistance)
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)
diff --git a/tutorial/206_GeodesicDistance/main.cpp b/tutorial/206_GeodesicDistance/main.cpp
new file mode 100755
index 000000000..ead43f0f5
--- /dev/null
+++ b/tutorial/206_GeodesicDistance/main.cpp
@@ -0,0 +1,73 @@
+#include 
+#include 
+#include 
+#include 
+#include 
+
+#include 
+#include "tutorial_shared_path.h"
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+void plotMeshDistance(igl::opengl::glfw::Viewer& viewer, const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, const Eigen::VectorXd& d, const double strip_size )
+{
+    // Rescale the function depending on the strip size
+    Eigen::VectorXd f = (d/strip_size);
+
+    // The function should be 1 on each integer coordinate
+    f = (f*M_PI).array().sin().abs();
+
+    // Compute per-vertex colors
+    Eigen::MatrixXd C;
+    igl::colormap(igl::COLOR_MAP_TYPE_INFERNO,f,false,C);
+
+    // Plot the mesh
+    viewer.data().set_mesh(V, F);
+    viewer.data().set_colors(C);
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // Load a mesh in OFF format
+  igl::readOBJ(TUTORIAL_SHARED_PATH "/armadillo.obj", V, F);
+
+
+  igl::opengl::glfw::Viewer viewer;
+  // Plot a distance when a vertex is picked
+  viewer.callback_mouse_down =
+  [&](igl::opengl::glfw::Viewer& viewer, int, int)->bool
+  {
+      int fid;
+      Eigen::Vector3f bc;
+      // Cast a ray in the view direction starting from the mouse position
+      double x = viewer.current_mouse_x;
+      double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+      if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core.view * viewer.core.model,
+                                  viewer.core.proj, viewer.core.viewport, V, F, fid, bc))
+      {
+          int max;
+          bc.maxCoeff(&max);
+          int vid = F(fid,max);
+          Eigen::VectorXi VS,FS,VT,FT;
+          // The selected vertex is the source
+          VS.resize(1);
+          VS << vid;
+          // All vertices are the targets
+          VT.setLinSpaced(V.rows(),0,V.rows()-1);
+          Eigen::VectorXd d;
+          igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
+
+          plotMeshDistance(viewer,V,F,d,0.05);
+      }
+      return false;
+  };
+  viewer.data().set_mesh(V,F);
+
+  cout << "Press [space] to smooth." << endl;;
+  cout << "Press [r] to reset." << endl;;
+  return viewer.launch();
+}
diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt
index 3ca47bb50..9c993a873 100644
--- a/tutorial/CMakeLists.txt
+++ b/tutorial/CMakeLists.txt
@@ -64,6 +64,7 @@ if(TUTORIALS_CHAPTER2)
   add_subdirectory("203_CurvatureDirections")
   add_subdirectory("204_Gradient")
   add_subdirectory("205_Laplacian")
+  add_subdirectory("206_GeodesicDistance")
 endif()
 
 # Chapter 3
diff --git a/tutorial/images/geodesicdistance.jpg b/tutorial/images/geodesicdistance.jpg
new file mode 100644
index 000000000..eefc4b9c2
Binary files /dev/null and b/tutorial/images/geodesicdistance.jpg differ
diff --git a/tutorial/tutorial.html b/tutorial/tutorial.html
index 1dcb725e3..8f9d830ca 100644
--- a/tutorial/tutorial.html
+++ b/tutorial/tutorial.html
@@ -1,15 +1,12 @@
 
-
+
 
 	
 	libigl Tutorial
 	
 	
 	
-
-
-
-
+   
 
 
 
@@ -19,6 +16,7 @@
 
 
      
 
+
      
 
      
 
 
      Libigl is an open source C++ library for geometry processing research and
@@ -48,7 +46,7 @@ lecture notes links to a cross-platform example application.
      
 
    • 107 Multiple Meshes
      •
  • Chapter 2: Discrete Geometric Quantities and - Operators +Operators
  • Chapter 3: Matrices and Linear Algebra @@ -96,7 +95,7 @@ lecture notes links to a cross-platform example application.

  • 404 Dual Quaternion Skinning
  • 405 As-rigid-as-possible
  • 406 Fast automatic skinning - transformations +transformations -

    Chapter 1

    +

    Chapter 1

    We introduce libigl with a series of self-contained examples. The purpose of each example is to showcase a feature of libigl while applying to a practical @@ -166,24 +165,24 @@ concepts of libigl and introduce a simple mesh viewer that allows to visualize a surface mesh and its attributes. All the tutorial examples are cross-platform and can be compiled on MacOSX, Linux and Windows.

    -

    libigl design principles

    +

    libigl design principles

    Before getting into the examples, we summarize the main design principles in libigl:

    1. No complex data types. We mostly use matrices and vectors. This greatly - favors code reusability and forces the function authors to expose all the - parameters used by the algorithm.

      2. +favors code reusability and forces the function authors to expose all the +parameters used by the algorithm.

    2. Minimal dependencies. We use external libraries only when necessary and - we wrap them in a small set of functions.

      3. +we wrap them in a small set of functions.

    3. Header-only. It is straight forward to use our library since it is only - one additional include directory in your project. (if you are worried about - compilation speed, it is also possible to build the library as a static - library)

      4. +one additional include directory in your project. (if you are worried about +compilation speed, it is also possible to build the library as a static +library)

    4. Function encapsulation. Every function (including its full - implementation) is contained in a pair of .h/.cpp files with the same name of - the function.

      5. +implementation) is contained in a pair of .h/.cpp files with the same name of +the function.

    Downloading libigl

    @@ -249,7 +248,7 @@ in 32 bit on windows.

    We provide a blank project example showing how to use libigl and cmake. Feel free and encouraged to copy or fork this project as a way of starting a new personal project using libigl.

    -

    Mesh representation

    +

    Mesh representation

    libigl uses the Eigen library to encode vector and matrices. We suggest that you keep the @@ -303,7 +302,7 @@ Similarly, a mesh can be written in an OBJ file using:

    Example 101 contains a simple mesh converter from OFF to OBJ format.

    -

    Visualizing surfaces

    +

    Visualizing surfaces

    Libigl provides an glfw-based OpenGL 3.2 viewer to visualize surfaces, their properties and additional debugging information.

    @@ -339,13 +338,13 @@ Please see the documentation in Viewer.h for more details.

    -(Example 102) loads and draws a +<img src=
    (Example 102) loads and draws a mesh.
    -

    Interaction with keyboard and mouse

    +

    Interaction with keyboard and mouse

    Keyboard and mouse events triggers callbacks that can be registered in the viewer. The viewer supports the following callbacks:

    @@ -398,7 +397,7 @@ control the camera directly in your code.

    the viewer’s callbacks. See the Viewer_plugin for more details.

    -

    Scalar field visualization

    +

    Scalar field visualization

    Colors and normals can be associated to faces or vertices using the set_colors function:

    @@ -427,7 +426,7 @@ igl::jet(Z,true,C); vertex) and the second calls a libigl functions that converts a scalar field to colors. The second parameter of jet normalizes the scalar field to lie between 0 and 1 before applying the transfer function.

    -(Example 104) igl::jet converts a scalar field to a +<img src=
    (Example 104) igl::jet converts a scalar field to a color field.
    @@ -437,7 +436,7 @@ color field. types and can be easily reused for many different tasks. Not committing to heavy data structures types favors simplicity, ease of use and reusability.

    -

    Overlays

    +

    Overlays

    In addition to plotting the surface, the viewer supports the visualization of points, lines and text labels: these overlays can be very helpful while developing geometric processing algorithms to plot debug information.

    @@ -467,13 +466,13 @@ Eigen::Vector3d M = V.colwise().maxCoeff();
    -(Example 105) The bounding box of a mesh is shown +<img src=
    (Example 105) The bounding box of a mesh is shown using overlays.
    -

    Viewer Menu

    +

    Viewer Menu

    As of version 1.2 the viewer uses a new menu and completely replaces AntTweakBar. It is based on the @@ -528,22 +527,24 @@ viewer.ngui->addVariable<bool>("bool",[&](bool val) {

    -(Example 106) The UI of the viewer can be easily customized. -
    (Example 106) The UI of the viewer can be easily customized.
    +([Example 106](106_ViewerMenu/main.cpp)) The UI of the viewer can be easily +customized. +
    (Example 106) The UI of the viewer can be easily +customized.
    -

    Multiple Meshes

    +

    Multiple Meshes

    Libigl’s igl::opengl::glfw::Viewer provides basic support for rendering multiple meshes.

    -

    Which mesh is selected is controled via the viewer.selected_data_index +

    Which mesh is selected is controlled via the viewer.selected_data_index field. By default it his is set to 0, so in the typical case of a single mesh viewer.data() returns the igl::ViewerData corresponding to the one and only mesh.

    -(Example 107) The igl::opengl::glfw::Viewer +<img src=
    (Example 107) The igl::opengl::glfw::Viewer @@ -555,12 +556,15 @@ colors.

    This chapter illustrates a few discrete quantities that libigl can compute on a mesh and the libigl functions that construct popular discrete differential -geometry operators. It also provides an introduction to basic drawing and coloring routines of our viewer.

    +geometry operators. It also provides an introduction to basic drawing and +coloring routines of our viewer.

    Normals

    Surface normals are a basic quantity necessary for rendering a surface. There -are a variety of ways to compute and store normals on a triangle mesh. Example 201 demonstrates how to compute and visualize normals with libigl.

    +are a variety of ways to compute and store normals on a triangle mesh. Example +201 demonstrates how to compute and visualize normals +with libigl.

    Per-face

    @@ -613,7 +617,7 @@ implements a simple scheme which computes corner normals as averages of normals of faces incident on the corresponding vertex which do not deviate by more than a specified dihedral angle (e.g. 20°).

    -The Normals example computes per-face (left), per-vertex (middle) and +<img src=
    The Normals example computes per-face (left), per-vertex (middle) and per-corner (right) normals
    @@ -630,8 +634,8 @@ principal curvatures:

    not the surface’s embedding.

    Intuitively, Gaussian curvature tells how locally spherical or elliptic the -surface is ( \(k_G>0\) ), how locally saddle-shaped or hyperbolic the surface -is ( \(k_G<0\) ), or how locally cylindrical or parabolic ( \(k_G=0\) ) the +surface is ( \(k_G>0\) ), how locally saddle-shaped or hyperbolic the surface +is ( \(k_G<0\) ), or how locally cylindrical or parabolic ( \(k_G=0\) ) the surface is.

    In the discrete setting, one definition for a “discrete Gaussian curvature” @@ -640,13 +644,13 @@ on a triangle mesh is via a vertex’s angular deficit:

    \(k_G(v_i) = 2π - \sum\limits_{j\in N(i)}θ_{ij},\)

    where \(N(i)\) are the triangles incident on vertex \(i\) and \(θ_{ij}\) is the angle -at vertex \(i\) in triangle \(j\) [1].

    +at vertex \(i\) in triangle \(j\) (1).

    Just like the continuous analog, our discrete Gaussian curvature reveals elliptic, hyperbolic and parabolic vertices on the domain, as demonstrated in Example 202.

    -The GaussianCurvature example computes discrete Gaussian curvature and +<img src=
    The GaussianCurvature example computes discrete Gaussian curvature and visualizes it in pseudocolor.
    @@ -670,7 +674,7 @@ normal:

    \(-\Delta \mathbf{x} = H \mathbf{n}.\)

    It is easy to compute this on a discrete triangle mesh in libigl using the -cotangent Laplace-Beltrami operator [1].

    +cotangent Laplace-Beltrami operator (1).

    #include <igl/cotmatrix.h>
 #include <igl/massmatrix.h>
@@ -687,16 +691,16 @@ H = HN.rowwise().norm(); //up to sign
 
 
    Combined with the angle defect definition of discrete Gaussian curvature, one
 can define principal curvatures and use least squares fitting to find
-directions [1].
    
+directions (1).
    
 
 
    Alternatively, a robust method for determining principal curvatures is via
-quadric fitting [2]. In the neighborhood around every vertex, a
+quadric fitting (2). In the neighborhood around every vertex, a
 best-fit quadric is found and principal curvature values and directions are
 analytically computed on this quadric (Example
 203).
    
 
 
    
-The CurvatureDirections example computes principal curvatures via quadric
+<img src=
 
    The CurvatureDirections example computes principal curvatures via quadric
@@ -726,8 +730,8 @@ linear on incident triangles.
    
 of the hat functions:
    
 
 
    \(\nabla f(\mathbf{x}) \approx
- \nabla \sum\limits_{i=1}^n \phi_i(\mathbf{x})\, f_i =
- \sum\limits_{i=1}^n \nabla \phi_i(\mathbf{x})\, f_i.\)
    
+\nabla \sum\limits_{i=1}^n \phi_i(\mathbf{x})\, f_i =
+\sum\limits_{i=1}^n \nabla \phi_i(\mathbf{x})\, f_i.\)
    
 
 
    This reveals that the gradient is a linear function of the vector of \(f_i\)
 values. Because the \(\phi_i\) are linear in each triangle, their gradients are
@@ -738,12 +742,12 @@ as a matrix multiplication taking vertex values to triangle values:
    
 
 
    where \(\mathbf{f}\) is \(n\times 1\) and \(\mathbf{G}\) is an \(md\times n\) sparse
 matrix. This matrix \(\mathbf{G}\) can be derived geometrically, e.g.
-[ch. 2, 3].
+(ch. 2, 3).
 Libigl’s grad function computes \(\mathbf{G}\) for
 triangle and tetrahedral meshes (Example 204):
    
 
 
    
-The Gradient example computes gradients of an input function on a mesh and
+<img src=
 
    The Gradient example computes gradients of an input function on a mesh and
 visualizes the vector field.
    
@@ -759,9 +763,9 @@ gradient (or equivalently the Laplacian of a function is the trace of its
 Hessian):
    
 
 
    \(\Delta f =
- \frac{\partial^2 f}{\partial x^2} +
- \frac{\partial^2 f}{\partial y^2} +
- \frac{\partial^2 f}{\partial z^2}.\)
    
+\frac{\partial^2 f}{\partial x^2} +
+\frac{\partial^2 f}{\partial y^2} +
+\frac{\partial^2 f}{\partial z^2}.\)
    
 
 
    The Laplace-Beltrami operator generalizes this to surfaces.
    
 
@@ -772,9 +776,9 @@ simultaneously from FEM, DEC and applying divergence theorem to vertex
 one-rings. As a linear operator taking vertex values to vertex values, the
 Laplacian \(\mathbf{L}\) is a \(n\times n\) matrix with elements:
    
 
-
    \(L_{ij} = \begin{cases}j \in N(i) &\cot \alpha_{ij} + \cot \beta_{ij},\\
-j \notin N(i) & 0,\\
-i = j & -\sum\limits_{k\neq i} L_{ik},
+
    \(L_{ij} = \begin{cases}j \in N(i) &\cot \alpha_{ij} + \cot \beta_{ij},\\
+j \notin N(i) & 0,\\
+i = j & -\sum\limits_{k\neq i} L_{ik},
 \end{cases}\)
    
 
 
    where \(N(i)\) are the vertices adjacent to (neighboring) vertex \(i\), and
@@ -813,16 +817,16 @@ of the Dirichlet energy (sum of squared gradients):
    
 
    Libigl implements discrete “cotangent” Laplacians for triangles meshes and
 tetrahedral meshes, building both with fast geometric rules rather than “by the
 book” FEM construction which involves many (small) matrix inversions, cf.
-[4].
    
+(4)[].
    
 
 
    The operator applied to mesh vertex positions amounts to smoothing by flowing
 the surface along the mean curvature normal direction (Example 205). Note that this is equivalent to minimizing surface area.
    
 
 
    
-The Laplacian example computes conformalized mean curvature flow using the
-cotangent Laplacian .
+The `Laplacian` example computes conformalized mean curvature flow using the
+cotangent Laplacian [#kazhdan_2012][].
 
    The Laplacian example computes conformalized mean curvature flow using the
-cotangent Laplacian [5].
    
+cotangent Laplacian (5)[].
 
    
 
 
    Mass matrix
    
@@ -831,7 +835,7 @@ cotangent Laplacian [5
 values to vertex values. From an FEM point of view, it is a discretization of
 the inner-product: it accounts for the area around each vertex. Consequently,
 \(\mathbf{M}\) is often a diagonal matrix, such that \(M_{ii}\) is the barycentric
-or voronoi area around vertex \(i\) in the mesh [1]. The inverse of
+or voronoi area around vertex \(i\) in the mesh (1)[]. The inverse of
 this matrix is also very useful as it transforms integrated quantities into
 point-wise quantities, e.g.:
    
 
@@ -859,8 +863,8 @@ Green’s identity (ignoring boundary conditions for the moment):
    
 
    Or in matrix form which is immediately translatable to code:
    
 
 
    \(\mathbf{f}^T \mathbf{G}^T \mathbf{T} \mathbf{G} \mathbf{f} =
-  \mathbf{f}^T \mathbf{M} \mathbf{M}^{-1} \mathbf{L} \mathbf{f} =
-  \mathbf{f}^T \mathbf{L} \mathbf{f}.\)
    
+\mathbf{f}^T \mathbf{M} \mathbf{M}^{-1} \mathbf{L} \mathbf{f} =
+\mathbf{f}^T \mathbf{L} \mathbf{f}.\)
    
 
 
    So we have that \(\mathbf{L} = \mathbf{G}^T \mathbf{T} \mathbf{G}\). This also
 hints that we may consider \(\mathbf{G}^T\) as a discrete divergence operator,
@@ -868,6 +872,34 @@ since the Laplacian is the divergence of the gradient. Naturally, \(n \times md\) sparse matrix which takes vector values stored at triangle faces
 to scalar divergence values at vertices.
    
 
+
    Geodesic
    
+
+
    The discrete geodesic distance between two points is the length of the shortest path between then restricted to the surface. For triangle meshes, such a path is made of a set of segments which can be either edges of the mesh or crossing a triangle.
    
+
+
    Libigl includes a wrapper for the exact geodesic algorithm (6) developed by Danil Kirsanov (https://code.google.com/archive/p/geodesic/), exposing it through an Eigen-based API. The function
    
+
+
    igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
+
    
+
+
    computes the closest geodesic distances of each vertex in VT or face in FT, from the source vertices VS or faces FS of the input mesh V,F. The output is writted in the vector d, which lists first the distances for the vertices in VT, and then for the faces in FT. For example, if you want to compute the distance from the vertex with id vid, to all vertices of F you can use:
    
+
+
    Eigen::VectorXi VS,FS,VT,FT;
+// The selected vertex is the source
+VS.resize(1);
+VS << vid;
+// All vertices are the targets
+VT.setLinSpaced(V.rows(),0,V.rows()-1);
+Eigen::VectorXd d;
+igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
+
    
+
+
    
+[Example 206](206_GeodesicDistance/main.cpp) allows to interactively pick the source vertex and displays the distance using a periodic color pattern.
+
+
    Example 206 allows to interactively pick the source vertex and displays the distance using a periodic color pattern.
+
    
+
    
+
 
    Chapter 3: Matrices and linear algebra
    
 
 
    Libigl relies heavily on the Eigen library for dense and sparse linear algebra
@@ -911,7 +943,7 @@ functionality is provided in libigl using slice_into:
    -The example Slice shows how to use igl::slice to change the colors for +<img src=
    The example Slice shows how to use igl::slice to change the colors for triangles on a mesh.
    @@ -953,11 +985,11 @@ X(I(i,j),j);. That is, I reveals how X is sorte

    Analogous functions are available in libigl for: max, min, and unique.

    -The example Sort shows how to use igl::sortrows to +<img src= +order ([Example 302](302_Sort/main.cpp))." />
    The example Sort shows how to use igl::sortrows to -pseudocolor triangles according to their barycenters’ sorted +pseudocolor triangles according to their barycenters' sorted order (Example 302).
    @@ -974,187 +1006,187 @@ functionality as common Matlab functions.

    - Name - Description + Name + Description - igl::all - Whether all elements are non-zero (true) + igl::all + Whether all elements are non-zero (true) - igl::any - Whether any elements are non-zero (true) + igl::any + Whether any elements are non-zero (true) - igl::cat - Concatenate two matrices (especially useful for dealing with Eigen sparse matrices) + igl::cat + Concatenate two matrices (especially useful for dealing with Eigen sparse matrices) - igl::ceil - Round entries up to nearest integer + igl::ceil + Round entries up to nearest integer - igl::cumsum - Cumulative sum of matrix elements + igl::cumsum + Cumulative sum of matrix elements - igl::colon - Act like Matlab’s :, similar to Eigen’s LinSpaced + igl::colon + Act like Matlab’s :, similar to Eigen’s LinSpaced - igl::components - Connected components of graph (cf. Matlab’s graphconncomp) + igl::components + Connected components of graph (cf. Matlab’s graphconncomp) - igl::count - Count non-zeros in rows or columns + igl::count + Count non-zeros in rows or columns - igl::cross - Cross product per-row + igl::cross + Cross product per-row - igl::cumsum - Cumulative summation + igl::cumsum + Cumulative summation - igl::dot - dot product per-row + igl::dot + dot product per-row - igl::eigs - Solve sparse eigen value problem + igl::eigs + Solve sparse eigen value problem - igl::find - Find subscripts of non-zero entries + igl::find + Find subscripts of non-zero entries - igl::floor - Round entries down to nearest integer + igl::floor + Round entries down to nearest integer - igl::histc - Counting occurrences for building a histogram + igl::histc + Counting occurrences for building a histogram - igl::hsv_to_rgb - Convert HSV colors to RGB (cf. Matlab’s hsv2rgb) + igl::hsv_to_rgb + Convert HSV colors to RGB (cf. Matlab’s hsv2rgb) - igl::intersect - Set intersection of matrix elements. + igl::intersect + Set intersection of matrix elements. - igl::isdiag - Determine whether matrix is diagonal + igl::isdiag + Determine whether matrix is diagonal - igl::ismember - Determine whether elements in A occur in B + igl::ismember + Determine whether elements in A occur in B - igl::jet - Quantized colors along the rainbow. + igl::jet + Quantized colors along the rainbow. - igl::max - Compute maximum entry per row or column + igl::max + Compute maximum entry per row or column - igl::median - Compute the median per column + igl::median + Compute the median per column - igl::min - Compute minimum entry per row or column + igl::min + Compute minimum entry per row or column - igl::mod - Compute per element modulo + igl::mod + Compute per element modulo - igl::mode - Compute the mode per column + igl::mode + Compute the mode per column - igl::null - Compute the null space basis of a matrix + igl::null + Compute the null space basis of a matrix - igl::nchoosek - Compute all k-size combinations of n-long vector + igl::nchoosek + Compute all k-size combinations of n-long vector - igl::orth - Orthogonalization of a basis + igl::orth + Orthogonalization of a basis - igl::parula - Generate a quantized colormap from blue to yellow + igl::parula + Generate a quantized colormap from blue to yellow - igl::pinv - Compute Moore-Penrose pseudoinverse + igl::pinv + Compute Moore-Penrose pseudoinverse - igl::randperm - Generate a random permutation of [0,…,n–1] + igl::randperm + Generate a random permutation of [0,…,n-1] - igl::rgb_to_hsv - Convert RGB colors to HSV (cf. Matlab’s rgb2hsv) + igl::rgb_to_hsv + Convert RGB colors to HSV (cf. Matlab’s rgb2hsv) - igl::repmat - Repeat a matrix along columns and rows + igl::repmat + Repeat a matrix along columns and rows - igl::round - Per-element round to whole number + igl::round + Per-element round to whole number - igl::setdiff - Set difference of matrix elements + igl::setdiff + Set difference of matrix elements - igl::setunion - Set union of matrix elements + igl::setunion + Set union of matrix elements - igl::setxor - Set exclusive “or” of matrix elements + igl::setxor + Set exclusive “or” of matrix elements - igl::slice - Slice parts of matrix using index lists: (cf. Matlab’s B = A(I,J)) + igl::slice + Slice parts of matrix using index lists: (cf. Matlab’s B = A(I,J)) - igl::slice_mask - Slice parts of matrix using boolean masks: (cf. Matlab’s B = A(M,N)) + igl::slice_mask + Slice parts of matrix using boolean masks: (cf. Matlab’s B = A(M,N)) - igl::slice_into - Slice left-hand side of matrix assignment using index lists (cf. Matlab’s B(I,J) = A) + igl::slice_into + Slice left-hand side of matrix assignment using index lists (cf. Matlab’s B(I,J) = A) - igl::sort - Sort elements or rows of matrix + igl::sort + Sort elements or rows of matrix - igl::speye - Identity as sparse matrix + igl::speye + Identity as sparse matrix - igl::sum - Sum along columns or rows (of sparse matrix) + igl::sum + Sum along columns or rows (of sparse matrix) - igl::unique - Extract unique elements or rows of matrix + igl::unique + Extract unique elements or rows of matrix @@ -1184,31 +1216,31 @@ the right-hand side.

    vertices come first and then boundary vertices:

    \[\left(\begin{array}{cc} - \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\\ - \mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{array}\right) - \left(\begin{array}{c} - \mathbf{z}_{in}\\ - \mathbf{z}_{b}\end{array}\right) = - \left(\begin{array}{c} - \mathbf{0}_{in}\\ - \mathbf{z}_{bc}\end{array}\right)\]

    +\mathbf{L}_{in,in} & \mathbf{L}_{in,b}\\ +\mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{array}\right) +\left(\begin{array}{c} +\mathbf{z}_{in}\\ +\mathbf{z}_{b}\end{array}\right) = +\left(\begin{array}{c} +\mathbf{0}_{in}\\ +\mathbf{z}_{bc}\end{array}\right)\]

    The bottom block of equations is no longer meaningful so we’ll only consider the top block:

    \[\left(\begin{array}{cc} - \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\end{array}\right) - \left(\begin{array}{c} - \mathbf{z}_{in}\\ - \mathbf{z}_{b}\end{array}\right) = - \mathbf{0}_{in}\]

    +\mathbf{L}_{in,in} & \mathbf{L}_{in,b}\end{array}\right) +\left(\begin{array}{c} +\mathbf{z}_{in}\\ +\mathbf{z}_{b}\end{array}\right) = +\mathbf{0}_{in}\]

    We can move the known values to the right-hand side:

    \[\mathbf{L}_{in,in} - \mathbf{z}_{in} = - - \mathbf{L}_{in,b} - \mathbf{z}_{b}\]

    +\mathbf{z}_{in} = - +\mathbf{L}_{in,b} +\mathbf{z}_{b}\]

    Finally we can solve this equation for the unknown values at interior vertices \(\mathbf{z}_{in}\).

    @@ -1223,7 +1255,7 @@ the linear algebra above directly. Then we can slice the solution into rows of Z corresponding to the interior vertices (Example 303).

    -The LaplaceEquation example solves a Laplace equation with Dirichlet +<img src=
    The LaplaceEquation example solves a Laplace equation with Dirichlet boundary conditions.
    @@ -1239,7 +1271,7 @@ energy subject to the same boundary conditions:

    On our discrete mesh, recall that this becomes

    \(\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{G}^T \mathbf{D} - \mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}\)

    +\mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}\)

    The general problem of minimizing some energy over a mesh subject to fixed value boundary conditions is so wide spread that libigl has a dedicated api for @@ -1249,26 +1281,28 @@ solving such systems.

    common constraints:

    \[\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} + - \mathbf{z}^T \mathbf{B} + \text{constant},\]

    +\mathbf{z}^T \mathbf{B} + \text{constant},\]

    subject to

    \[\mathbf{z}_b = \mathbf{z}_{bc} \text{ and } \mathbf{A}_{eq} \mathbf{z} = - \mathbf{B}_{eq},\]

    +\mathbf{B}_{eq},\]

    where

      -
    • \(\mathbf{Q}\) is a (usually sparse) \(n \times n\) positive semi-definite +
    • \(\mathbf{Q}\) is a (usually sparse) \(n \times n\) positive semi-definite matrix of quadratic coefficients (Hessian),
      • -
    • \(\mathbf{B}\) is a \(n \times 1\) vector of linear coefficients,
      • -
    • \(\mathbf{z}_b\) is a \(|b| \times 1\) portion of +
    • \(\mathbf{B}\) is a \(n \times 1\) vector of linear coefficients,
      • +
    • \(\mathbf{z}_b\) is a \(|b| \times 1\) portion of + \(\mathbf{z}\) corresponding to boundary or fixed vertices,
      • -
    • \(\mathbf{z}_{bc}\) is a \(|b| \times 1\) vector of known values corresponding to +
    • \(\mathbf{z}_{bc}\) is a \(|b| \times 1\) vector of known values corresponding to + \(\mathbf{z}_b\),
      • -
    • \(\mathbf{A}_{eq}\) is a (usually sparse) \(m \times n\) matrix of linear +
    • \(\mathbf{A}_{eq}\) is a (usually sparse) \(m \times n\) matrix of linear equality constraint coefficients (one row per constraint), and
      • -
    • \(\mathbf{B}_{eq}\) is a \(m \times 1\) vector of linear equality constraint +
    • \(\mathbf{B}_{eq}\) is a \(m \times 1\) vector of linear equality constraint right-hand side values.
    @@ -1303,7 +1337,7 @@ with active linear equality constraints. Specifically let’s solve the bi-Laplace equation or equivalently minimize the Laplace energy:

    \[\Delta^2 z = 0 \leftrightarrow \mathop{\text{minimize }}\limits_z \frac{1}{2} - \int\limits_S (\Delta z)^2 dA\]

    +\int\limits_S (\Delta z)^2 dA\]

    subject to fixed value constraints and a linear equality constraint:

    @@ -1323,41 +1357,44 @@ general a \(m \times 1\) vector of variables

    \[\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} + - \mathbf{z}^T \mathbf{B} + \text{constant} + \lambda^T\left(\mathbf{A}_{eq} - \mathbf{z} - \mathbf{B}_{eq}\right)\]

    +\mathbf{z}^T \mathbf{B} + \text{constant} + \lambda^T\left(\mathbf{A}_{eq} +\mathbf{z} - \mathbf{B}_{eq}\right)\]

    This can be rewritten in a more familiar form by stacking \(\mathbf{z}\) and \(\lambda\) into one \((m+n) \times 1\) vector of unknowns:

    -

    \[\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, - \frac{1}{2} - \left( - \mathbf{z}^T - \lambda^T - \right) - \left( - \begin{array}{cc} - \mathbf{Q} & \mathbf{A}_{eq}^T\\ - \mathbf{A}_{eq} & 0 - \end{array} - \right) - \left( - \begin{array}{c} - \mathbf{z}\\ - \lambda - \end{array} - \right) + - \left( - \mathbf{z}^T - \lambda^T - \right) - \left( - \begin{array}{c} - \mathbf{B}\\ - -\mathbf{B}_{eq} - \end{array} - \right) - + \text{constant}\]

    +

    $$\mathop{\text{find saddle }}_{\mathbf{z},\lambda}, +\frac{1}{2} +\left( +\mathbf{z}T +\lambdaT +\right) +\left( +\begin{array}{cc} +\mathbf{Q} & \mathbf{A}{eq}T\ +\mathbf{A}{eq} & 0 +\end{array} +\right) +\left( +\begin{array}{c} +\mathbf{z}\ +\lambda +\end{array} +\right) + +\left( +\mathbf{z}T +\lambdaT +\right) +\left( +\begin{array}{c} +\mathbf{B}\ +-\mathbf{B}_{eq} +\end{array} +\right)

    + +
      +
    • \text{constant}$$
      • +

    Differentiating with respect to \(\left( \mathbf{z}^T \lambda^T \right)\) reveals a linear system and we can solve for \(\mathbf{z}\) and \(\lambda\). The only @@ -1368,12 +1405,12 @@ different factorization technique (LDLT rather than LLT): libigl’s the presence of linear equality constraints (Example 304).

    -The example LinearEqualityConstraints first solves with just fixed value +<img src=
    The example LinearEqualityConstraints first solves with just fixed value -constraints (left: 1 and –1 on the left hand and foot respectively), then +constraints (left: 1 and -1 on the left hand and foot respectively), then solves with an additional linear equality constraint (right: points on right hand and foot constrained to be equal).
    @@ -1420,11 +1457,11 @@ igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
    - Example 305 uses an active set solver to optimize -discrete biharmonic kernels at multiple scales +<img src=
    Example 305 uses an active set solver to optimize -discrete biharmonic kernels [6] at multiple scales +discrete biharmonic kernels (7)[] at multiple scales .
    @@ -1437,7 +1474,7 @@ eigen value problem:

    where \(A\) is a sparse symmetric matrix and \(B\) is a sparse positive definite matrix. Most commonly in geometry processing, we let \(A=L\) the cotangent -Laplacian and \(B=M\) the per-vertex mass matrix (e.g. [7]). +Laplacian and \(B=M\) the per-vertex mass matrix (e.g. (8)[]). Typically applications will make use of the low frequency eigen modes. Analogous to the Fourier decomposition, a function \(f\) on a surface can be represented via its spectral decomposition of the eigen modes of the @@ -1457,7 +1494,7 @@ eigen vector satisfying:

    \(\mathbf{L} \phi_i = \lambda_i \mathbf{M} \phi_i\).

    -

    Note that Vallet & Levy [7] propose solving a symmetrized +

    Note that Vallet & Levy (8)[] propose solving a symmetrized standard eigen problem \(\mathbf{M}^{-1/2}\mathbf{L}\mathbf{M}^{-1/2} \phi_i = \lambda_i \phi_i\). Libigl implements a generalized eigen problem solver so this unnecessary symmetrization can be avoided.

    @@ -1465,9 +1502,9 @@ this unnecessary symmetrization can be avoided.

    Often the sum above is truncated to the first \(k\) eigen vectors. If the low frequency modes are chosen, i.e. those corresponding to small \(\lambda_i\) values, then this truncation effectively regularizes \(\mathbf{f}\) to smooth, -slowly changing functions over the mesh (e.g. [8]). Modal +slowly changing functions over the mesh (e.g. (9)[]). Modal analysis and model subspaces have been used frequently in real-time deformation -(e.g. [9]).

    +(e.g. (10)[]).

    In Example 306), the first 5 eigen vectors of the discrete Laplace-Beltrami operator are computed and displayed in @@ -1484,7 +1521,7 @@ igl::eigs(L,M,5,igl::EIGS_TYPE_SM,U,S);

    -(Example 306) Low frequency eigen vectors +<img src=
    (Example 306) Low frequency eigen vectors @@ -1511,9 +1548,9 @@ partial differential equation.

    There are many flavors of these techniques, but a prototypical subset are those that consider solutions to the bi-Laplace equation, that is a biharmonic -function [10]. This fourth-order PDE provides sufficient +function (11)[]. This fourth-order PDE provides sufficient flexibility in boundary conditions to ensure \(C^1\) continuity at handle -constraints (in the limit under refinement) [11].

    +constraints (in the limit under refinement) (12)[].

    Biharmonic surfaces

    @@ -1585,7 +1622,7 @@ U = V+D;
    -The BiharmonicDeformation example deforms a statues head as a biharmonic +<img src=
    The BiharmonicDeformation example deforms a statue’s head as a biharmonic @@ -1609,10 +1646,10 @@ terms of the original positions \(\mathbf{x}\) and the \(\mathbf{x}' = \mathbf{x} - \mathbf{d}\):

    \(\int\limits_S \|\Delta (\mathbf{x}' - \mathbf{x})\|^2 dA = \int\limits_S - \|\Delta \mathbf{x}' - \Delta \mathbf{x})\|^2 dA.\)

    +\|\Delta \mathbf{x}' - \Delta \mathbf{x})\|^2 dA.\)

    In the early work of Sorkine et al., the quantities \(\Delta \mathbf{x}'\) and -\(\Delta \mathbf{x}\) were dubbed “differential coordinates” [12]. +\(\Delta \mathbf{x}\) were dubbed “differential coordinates” (13)[]. Their deformations (without linearized rotations) is thus equivalent to biharmonic deformation fields.

    @@ -1632,7 +1669,7 @@ igl::harmonic(V,F,b,bc,k,Z);
    -The PolyharmonicDeformation example deforms a flat domain (left) into a bump as a +<img src=
    The PolyharmonicDeformation example deforms a flat domain (left) into a bump as a solution to various \(k\)-harmonic PDEs.
    @@ -1650,7 +1687,7 @@ rotations.

    computes its new location as a linear combination of bone transformations:

    \(\mathbf{x}' = \sum\limits_{i = 1}^m w_i(\mathbf{x}) \mathbf{T}_i - \left(\begin{array}{c}\mathbf{x}_i\\1\end{array}\right),\)

    +\left(\begin{array}{c}\mathbf{x}_i\\1\end{array}\right),\)

    where \(w_i(\mathbf{x})\) is the scalar weight function of the ith bone evaluated at \(\mathbf{x}\) and \(\mathbf{T}_i\) is the bone transformation as a \(4 \times 3\) @@ -1680,15 +1717,15 @@ any handle structure such as a cage, collection of points, selected regions, etc.).

    Bounded biharmonic weights are one such technique that casts weight computation -as a constrained optimization problem [13]. The weights enforce +as a constrained optimization problem (14)[]. The weights enforce smoothness by minimizing the familiar Laplacian energy:

    \(\sum\limits_{i = 1}^m \int_S (\Delta w_i)^2 dA\)

    subject to constraints which enforce interpolation of handle constraints:

    -

    \(w_i(\mathbf{x}) = \begin{cases} 1 & \text{ if } \mathbf{x} \in H_i\\ 0 & - \text{ otherwise } \end{cases},\)

    +

    \(w_i(\mathbf{x}) = \begin{cases} 1 & \text{ if } \mathbf{x} \in H_i\\ 0 & +\text{ otherwise } \end{cases},\)

    where \(H_i\) is the ith handle, and constraints which enforce non-negativity, parition of unity and encourage sparsity:

    @@ -1699,7 +1736,7 @@ parition of unity and encourage sparsity:

    set solver or by calling out to Mosek.

    -The example BoundedBiharmonicWeights computes weights for a tetrahedral +<img src=
    The example BoundedBiharmonicWeights computes weights for a tetrahedral @@ -1717,22 +1754,22 @@ about the \(z\)-axis. Intuitively one might expect to but instead the blend is a degenerate matrix scaling the \(x\) and \(y\) coordinates by zero:

    -

    \(0.5\left(\begin{array}{ccc}0&-1&0\\1&0&0\\0&0&1\end{array}\right)+ - 0.5\left(\begin{array}{ccc}0&1&0\\-1&0&0\\0&0&1\end{array}\right)= - \left(\begin{array}{ccc}0&0&0\\0&0&0\\0&0&1\end{array}\right)\)

    +

    \(0.5\left(\begin{array}{ccc}0&-1&0\\1&0&0\\0&0&1\end{array}\right)+ +0.5\left(\begin{array}{ccc}0&1&0\\-1&0&0\\0&0&1\end{array}\right)= +\left(\begin{array}{ccc}0&0&0\\0&0&0\\0&0&1\end{array}\right)\)

    In practice, this means the shape shrinks and collapses in regions where bone weights overlap: near joints.

    -

    Dual quaternion skinning presents a solution [14]. This method +

    Dual quaternion skinning presents a solution (15). This method represents rigid transformations as a pair of unit quaternions, \(\hat{\mathbf{q}}\). The linear blend skinning formula is replaced with a linear blend of dual quaternions:

    \(\mathbf{x}' = - \cfrac{\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}} - {\left\|\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}\right\|} - \mathbf{x},\)

    +\cfrac{\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}} +{\left\|\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}\right\|} +\mathbf{x},\)

    where \(\hat{\mathbf{q}_i}\) is the dual quaternion representation of the rigid transformation of bone \(i\). The normalization forces the result of the linear @@ -1752,8 +1789,8 @@ igl::dqs(V,W,vQ,vT,U);

    -The example DualQuaternionSkinning compares linear blend skinning (top) to dual -quaternion skinning (bottom), highlighting LBSs candy wrapper effect (middle) +<img src=
    The example DualQuaternionSkinning compares linear blend skinning (top) to dual quaternion skinning (bottom), highlighting LBS’s candy wrapper effect (middle) @@ -1781,15 +1818,15 @@ They work by comparing the deformation of a mesh vertex to its rest position rotated to a new coordinate frame which best matches the deformation. The non-linearity stems from the mutual dependence of the deformation and the best-fit rotation. These techniques are often labeled -“as-rigid-as-possible” as they penalize the sum of all local deformations’ +“as-rigid-as-possible” as they penalize the sum of all local deformations' deviations from rotations.

    To arrive at such an energy, let’s consider a simple per-triangle energy:

    \(E_\text{linear}(\mathbf{X}') = \sum\limits_{t \in T} a_t \sum\limits_{\{i,j\} - \in t} w_{ij} \left\| - \left(\mathbf{x}'_i - \mathbf{x}'_j\right) - - \left(\mathbf{x}_i - \mathbf{x}_j\right)\right\|^2\)

    +\in t} w_{ij} \left\| +\left(\mathbf{x}'_i - \mathbf{x}'_j\right) - +\left(\mathbf{x}_i - \mathbf{x}_j\right)\right\|^2\)

    where \(\mathbf{X}'\) are the mesh’s unknown deformed vertex positions, \(t\) is a triangle in a list of triangles \(T\), \(a_t\) is the area of triangle \(t\) and @@ -1808,10 +1845,11 @@ for each triangle \(t\) which are constrained to be ro rewritten, this time comparing deformed edge vectors to their rotated rest counterparts:

    -

    \(E_\text{arap}(\mathbf{X}',\{\mathbf{R}_1,\dots,\mathbf{R}_{|T|}\}) = \sum\limits_{t \in T} a_t \sum\limits_{\{i,j\} - \in t} w_{ij} \left\| - \left(\mathbf{x}'_i - \mathbf{x}'_j\right)- - \mathbf{R}_t\left(\mathbf{x}_i - \mathbf{x}_j\right)\right\|^2.\)

    +

    $E_\text{arap}(\mathbf{X}',{\mathbf{R}1,\dots,\mathbf{R}{|T|}}) = \sum\limits_{t \in T} a_t \sum\limits_{{i,j}

    + +

    \in t} w_{ij} \left| +\left(\mathbf{x}‘_i - \mathbf{x}’_j\right)- +\mathbf{R}_t\left(\mathbf{x}_i - \mathbf{x}_j\right)\right|2.$

    The separation into the primary vertex position variables \(\mathbf{X}'\) and the rotations \(\{\mathbf{R}_1,\dots,\mathbf{R}_{|T|}\}\) lead to strategy for @@ -1825,11 +1863,11 @@ the energy, thus we may safely iterate them until convergence.

    The different flavors of “as-rigid-as-possible” depend on the dimension and codimension of the domain and the edge-sets \(T\). The proposed surface -manipulation technique by Sorkine and Alexa [15], considers \(T\) to +manipulation technique by Sorkine and Alexa (16)[], considers \(T\) to be the set of sets of edges emanating from each vertex (spokes). Later, Chao et al. derived the relationship between “as-rigid-as-possible” mesh energies and co-rotational elasticity considering 0-codimension elements as edge-sets: -triangles in 2D and tetrahedra in 3D [16]. They also showed how +triangles in 2D and tetrahedra in 3D (17)[]. They also showed how Sorkine and Alexa’s edge-sets are not a discretization of a continuous energy, proposing instead edge-sets for surfaces containing all edges of elements incident on a vertex (spokes and rims). They show that this amounts to @@ -1856,10 +1894,10 @@ certain constraints on the positions of vertices in b, we may call:

    Libigl’s implementation of as-rigid-as-possible deformation takes advantage of the highly optimized singular value decomposition code from McAdams et al. -[17] which leverages SSE intrinsics.

    +(18)[] which leverages SSE intrinsics.

    -The example AsRigidAsPossible deforms a surface as if it were made of an +<img src=
    The example AsRigidAsPossible deforms a surface as if it were made of an elastic material
    @@ -1875,7 +1913,7 @@ the case of the as-rigid-as-possible optimization, the bottleneck is typically the large number of polar decompositions necessary to recover best fit rotations for each edge-set (i.e. for each triangle, tetrahedron, or vertex cell). Even if this code is optimized, the number of primary degrees of freedom -is tied to the discretization level, despite the deformations’ low frequency +is tied to the discretization level, despite the deformations' low frequency behavior.

    This invites two routes toward fast non-linear optimization. First, is it @@ -1895,7 +1933,7 @@ linear blend skinning in matrix form:

    replaced by a linear combination of a small number of degrees of freedom in the \((3+1)m \times 3\) stack of transposed “handle” transformations. Swapping in \(\mathbf{M}\mathbf{T}\) for \(\mathbf{X}'\) in the ARAP energies above immediately -sees performance gains during the global solve step as \(m << n\).

    +sees performance gains during the global solve step as \(m << n\).

    The complexity of the local step—fitting rotations—is still bound to the original mesh discretization. However, if the skinning is well behaved, @@ -1908,7 +1946,7 @@ clustered edge-sets show diminishing returns on the deformation quality so we may choose a small number of clusters, proportional to the number of skinning weight functions (rather than the number of discrete mesh vertices).

    -

    This proposed deformation model [18], can simultaneously be seen as a +

    This proposed deformation model (19)[], can simultaneously be seen as a fast, subspace optimization for ARAP and as an automatic method for finding the best skinning transformation degrees of freedom.

    @@ -1917,11 +1955,11 @@ the skinning transformations associated with handles. To fix a transformation entirely we simply add the constraint:

    \(\left(\begin{array}{cccc} - 1 & 0 & 0 & 0\\ - 0 & 1 & 0 & 0\\ - 0 & 0 & 1 & 0\\ - 0 & 0 & 0 & 1\end{array}\right) - \mathbf{T}_i^T = \hat{\mathbf{T}}_i^T,\)

    +1 & 0 & 0 & 0\\ +0 & 1 & 0 & 0\\ +0 & 0 & 1 & 0\\ +0 & 0 & 0 & 1\end{array}\right) +\mathbf{T}_i^T = \hat{\mathbf{T}}_i^T,\)

    where \(\hat{\mathbf{T}}_i^T\) is the \((3+1) \times 3\) transposed fixed transformation for handle \(i\).

    @@ -1940,10 +1978,10 @@ user.

    freeing the translation component (producing a “chickenhead” effect):

    \(\left(\begin{array}{cccc} - 1&0&0&0\\ - 0&1&0&0\\ - 0&0&1&0\end{array}\right) - \mathbf{T}_i^T = \hat{\mathbf{L}}_i^T,\)

    +1&0&0&0\\ +0&1&0&0\\ +0&0&1&0\end{array}\right) +\mathbf{T}_i^T = \hat{\mathbf{L}}_i^T,\)

    where \(\hat{\mathbf{L}}_i^T\) is the fixed \(3 \times 3\) linear part of the transformation at handle \(i\).

    @@ -1971,7 +2009,7 @@ biharmonic distance embedding.

    redundant) clustering of the per-triangle edge-sets.

    -The example FastAutomaticSkinningTransformations compares a full (slow) +<img src= @@ -1987,7 +2025,7 @@ rotation edge sets (right of middle), to the very fast subpsace method propagating full affine transformations at handles (bones, points, regions, etc.) to the rest of the shape via weights. Another deformation framework, called “generalized barycentric coordinates”, is a special case of linear blend -skinning [19]: transformations are restricted to +skinning (20)[]: transformations are restricted to pure translations and weights are required to retain affine precision. This latter requirement means that we can write the rest-position of any vertex in the mesh as the weighted combination of the control handle locations:

    @@ -2003,19 +2041,19 @@ position of each point of the shape to be the weighted combination of the

    There are many different flavors of “generalized barycentric coordinates” (see table in “Automatic Methods” section, -[19]). The vague goal of “generalized barycentric +(20)[]). The vague goal of “generalized barycentric coordinates” is to capture as many properties of simplicial barycentric coordinates (e.g. for triangles in 2D and tetrahedral in 3D) for larger sets of points or polyhedra. Some generalized barycentric coordinates can be computed in closed form; others require optimization-based precomputation. Nearly all flavors require connectivity information describing how the control points form a external polyhedron around the input shape: a cage. However, a recent -techinique does not require a cage [20]. This method ensures +techinique does not require a cage (21)[]. This method ensures affine precision during optimization over weights of a smoothness energy with affine functions in its kernel:

    \(\mathop{\text{min}}_\mathbf{W}\,\, \text{trace}(\frac{1}{2}\mathbf{W}^T \mathbf{A} - \mathbf{W}), \text{subject to: } \mathbf{C} = \mathbf{W}\mathbf{C}\)

    +\mathbf{W}), \text{subject to: } \mathbf{C} = \mathbf{W}\mathbf{C}\)

    subject to interpolation constraints at selected vertices. If \(\mathbf{A}\) has affine functions in its kernel—that is, if \(\mathbf{A}\mathbf{V} = 0\)—then @@ -2040,12 +2078,12 @@ the integral-average of \(\mathbf{K}\) applied to a fu the mesh:

    \(\mathbf{A} = (\mathbf{M}^{-1}\mathbf{K})^2_\mathbf{M} = \mathbf{K}^T \mathbf{M}^{-1} - \mathbf{K}.\)

    +\mathbf{K}.\)

    Since the Laplacian \(\mathbf{K}\) is a second-order derivative it measures zero on affine functions, thus \(\mathbf{A}\) has affine functions in its null space. A short derivation proves that this implies \(\mathbf{W}\) will be affine precise (see -[20]).

    +(21)[]).

    Minimizers of this “squared Laplacian” energy are in some sense discrete biharmonic functions. Thus they’re dubbed “biharmonic coordinates” (not the @@ -2060,7 +2098,7 @@ handles):

    -(Example 407) shows a physics +<img src= @@ -2070,7 +2108,7 @@ points for a biharmonic coordinates deformation of the blue high-resolution mesh.
    -

    Chapter 5: Parametrization

    +

    Chapter 5: Parametrization

    In computer graphics, we denote as surface parametrization a map from the surface to \(\mathbf{R}^2\). It is usually encoded by a new set of 2D @@ -2096,9 +2134,9 @@ genus. They initially cut the mesh in multiple patches that can be separately pa

  • Global seamless parametrization: these are global parametrization algorithm that hides the seams, making the parametrization “continuous”, under specific assumptions that we will discuss later.

    • -

    Harmonic parametrization

    +

    Harmonic parametrization

    -

    Harmonic parametrization [21] is a single patch, fixed boundary parametrization +

    Harmonic parametrization (22)[] is a single patch, fixed boundary parametrization algorithm that computes the 2D coordinates of the flattened mesh as two harmonic functions.

    @@ -2133,7 +2171,7 @@ functions is designed to be reusable in other parametrization algorithms.

    mesh (Example 501).

    -(Example 501) Harmonic parametrization. (left) +<img src=
    (Example 501) Harmonic parametrization. (left) @@ -2141,9 +2179,9 @@ mesh with texture, (right) UV parametrization with texture
    -

    Least squares conformal maps

    +

    Least squares conformal maps

    -

    Least squares conformal maps parametrization [22] minimizes the +

    Least squares conformal maps parametrization (23)[] minimizes the conformal (angular) distortion of the parametrization. Differently from harmonic parametrization, it does not need to have a fixed boundary.

    @@ -2151,7 +2189,7 @@ harmonic parametrization, it does not need to have a fixed boundary.

    \[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \int_X \frac{1}{2}| \nabla \mathbf{u}^{\perp} - \nabla \mathbf{v} |^2 dA \]

    -

    which can be rewritten in matrix form as [23]:

    +

    which can be rewritten in matrix form as (24)[]:

    \[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \frac{1}{2} [\mathbf{u},\mathbf{v}]^t (L_c - 2A) [\mathbf{u},\mathbf{v}] \]

    @@ -2185,15 +2223,15 @@ case we do not need to fix the boundary. To remove the null space of the energy vertices to two arbitrary positions. The full source code is provided in Example 502.

    -(Example 502) LSCM parametrization. (left) mesh +<img src=
    (Example 502) LSCM parametrization. (left) mesh with texture, (right) UV parametrization
    -

    As-rigid-as-possible parametrization

    +

    As-rigid-as-possible parametrization

    -

    As-rigid-as-possible parametrization [24] is a powerful single-patch, +

    As-rigid-as-possible parametrization (25)[] is a powerful single-patch, non-linear algorithm to compute a parametrization that strives to preserve distances (and thus angles). The idea is very similar to ARAP surface deformation: each triangle is mapped to the plane trying to preserve its @@ -2210,7 +2248,7 @@ parametrization. Similarly to LSCM, the boundary is free to deform to minimize the distortion.

    -(Example 503) As-Rigid-As-Possible parametrization. +<img src=
    (Example 503) As-Rigid-As-Possible parametrization. @@ -2218,7 +2256,7 @@ texture" /> texture
    -

    N-rotationally symmetric tangent fields

    +

    N-rotationally symmetric tangent fields

    The design of tangent fields is a basic tool used to design guidance fields for uniform quadrilateral and hexahedral remeshing. Libigl contains an @@ -2245,7 +2283,7 @@ the triangle mesh (output_field), plus the singularities of the field

    The singularities are vertices where the field vanishes (highlighted in red in -the figure above). igl::nrosy can also generate N-RoSy fields [25], +the figure above). igl::nrosy can also generate N-RoSy fields (26)[], which are a generalization of vector fields where in every face the vector is defined up to a constant rotation of \(2\pi / N\). As can be observed in the following figure, the singularities of the fields generated with different @@ -2259,11 +2297,11 @@ N are of different types and they appear in different positions.

    We demonstrate how to call and plot N-RoSy fields in Example 504, where the degree of the field can be change pressing the number keys. igl::nrosy implements the algorithm proposed in -[26]. N-RoSy fields can also be interpolated with the algorithm -proposed in [27], see Section npolyvectorfields for more details +(27)[]. N-RoSy fields can also be interpolated with the algorithm +proposed in (28)[], see Section npolyvectorfields for more details (igl::n_polyvector).

    -

    Global, seamless integer-grid parametrization

    +

    Global, seamless integer-grid parametrization

    The previous parametrization methods were focusing on creating parametrizations of surface patches aimed at texture mapping or baking of other surface @@ -2271,7 +2309,7 @@ properties such as normals and high-frequency details. Global, seamless parametrization aims at parametrizing complex shapes with a parametrization that is aligned with a given set of directions for the purpose of surface remeshing. In libigl, we provide a reference implementation of the pipeline -proposed in the mixed integer quadrangulation paper [26].

    +proposed in the mixed integer quadrangulation paper (27)[].

    The first step involves the design of a 4-RoSy field (sometimes called cross field) that describes the alignment of the edges of the desired quadrilateral @@ -2350,7 +2388,7 @@ input cross field.

    We hide the seams by adding integer constraints to the Poisson problem -that align the isolines on both sides of each seam [26].

    +that align the isolines on both sides of each seam (27).

    Seamless Poisson parametrization. @@ -2369,12 +2407,12 @@ it contains many overlaps.

    libQEx (not included in libigl). The full pipeline is implemented in Example 505.

    -

    Anisotropic remeshing

    +

    Anisotropic remeshing

    Anisotropic and non-uniform quad remeshing is important to concentrate the elements in the regions with more details. It is possible to extend the MIQ quad meshing framework to generate anisotropic quad meshes using a mesh -deformation approach [28].

    +deformation approach (29)[].

    The input of the anisotropic remeshing algorithm is a sparse set of constraints that define the shape and scale of the desired quads. This can be encoded as a @@ -2430,11 +2468,11 @@ generate the UV parametrization, but other algorithms could be applied: the only desiderata is that the generated quad mesh should be as isotropic as possible.

    -

    N-PolyVector fields

    +

    N-PolyVector fields

    N-RoSy vector fields can be further generalized to represent arbitrary vector-sets, with arbitrary angles between them and with arbitrary lengths -[29]. This generalization is called N-PolyVector field, and +(30)[]. This generalization is called N-PolyVector field, and libigl provides the function igl::n_polyvector to design them starting from a sparse set of constraints (Example 507).

    @@ -2447,12 +2485,12 @@ sparse set of constraints (Example 507 -

    Globally optimal direction fields [27] are a special case of +

    Globally optimal direction fields (28)[] are a special case of PolyVector fields. If the constraints are taken from an N-RoSy field, igl::n_polyvector generates a field that is equivalent, after normalization, to a globally optimal direction field.

    -

    Conjugate vector fields

    +

    Conjugate vector fields

    Two tangent vectors lying on a face of a triangle mesh are conjugate if

    @@ -2461,13 +2499,13 @@ to a globally optimal direction field.

    This condition is very important in architectural geometry: The faces of an infinitely dense quad mesh whose edges are aligned with a conjugate field are planar. Thus, a quad mesh whose edges follow a conjugate field are easier to -planarize [30].

    +planarize (31).

    Finding a conjugate vector field that satisfies given directional constraints is a standard problem in architectural geometry, which can be tackled by deforming a Poly-Vector field to the closest conjugate field.

    -

    This algorithm [29] alternates a global step, which enforces +

    This algorithm (30) alternates a global step, which enforces smoothness, with a local step, that projects the field on every face to the closest conjugate field (Example 508).

    @@ -2478,10 +2516,10 @@ closest conjugate field (Example 508). (right).
    -

    Planarization

    +

    Planarization

    A quad mesh can be transformed in a planar quad mesh with Shape-Up -[31], a local/global approach that uses the global step to enforce +(32), a local/global approach that uses the global step to enforce surface continuity and the local step to enforce planarity.

    Example 509 planarizes a quad mesh until it @@ -2496,7 +2534,7 @@ igl::palanarize (right). The colors represent the planarity of the quads.

    -

    Integrable PolyVector Fields

    +

    Integrable PolyVector Fields

    Vector-field guided surface parameterization is based on the idea of designing the gradients of the parameterization functions (which are tangent vector fields @@ -2505,8 +2543,8 @@ on the surface) instead of the functions themselves. Thus, vector-set fields parameterization (and subsequent remeshing) need to be integrable: it must be possible to break them down into individual vector fields that are gradients of scalar functions. Fields obtained by most smoothness-based design methods (eg. -[25], [27], [29], [26], -[28]) do not have this property. In [32], a method +(26)[], (28)[], (30)[], (27)[], +(29)[]) do not have this property. In (33)[], a method for creating integrable polyvector fields was introduced. This method takes as input a given field and improves its integrability by removing the vector field curl, thus turning it into a gradient of a function (Example @@ -2526,11 +2564,11 @@ continuous variables only. This is done using coefficients of appropriately defined polynomials. The parameterizations generated by the resulting fields are exactly aligned to the field directions and contain no inverted triangles.

    -

    General N-PolyVector fields

    +

    General N-PolyVector fields

    While mostly applicable for the design of symmetric fields (i.e. fields that comprise of vector sets with symmetries between them at each point, e.g. N-RoSy -or frame-fields), the framework presented in [29] can be used to +or frame-fields), the framework presented in (30)[] can be used to design completely general fields, with possibly no such symmetries. For example, one can design fields that at each point comprise of an arbitrary number of vectors, not required to be collinear - as opposed e.g. to the case of the 4 @@ -2542,10 +2580,12 @@ function igl::n_polyvector_general, and is illustrated in the example ( +mesh faces, but is only shown on a subset for clarity. +" />

    Interpolation of a general field with 3 (left) and 9 vectors per point field from a sparse set of random constraints (in red). The field is defined on all -mesh faces, but is only shown on a subset for clarity.
    +mesh faces, but is only shown on a subset for clarity. +

    The design of these general directional fields (also called vector-set fields) @@ -2562,14 +2602,14 @@ to the particular nature of the polynomial that applies in that case (two coefficients are 0).

    For a complete categorization of fields used in various applications (including -these general ones) see Vaxman et al. 2016 [33].

    +these general ones) see Vaxman et al. 2016 (34).

    -

    Chapter 6: External libraries

    +

    Chapter 6: External libraries

    An additional positive side effect of using matrices as basic types is that it is easy to exchange data between libigl and other software and libraries.

    -

    State serialization

    +

    State serialization

    Geometry processing applications often require a considerable amount of computational time and/or manual input. Serializing the state of the application @@ -2714,7 +2754,7 @@ common to have to do small changes to figures, and being able to serialize the entire state just before you take screenshots will save you many painful hours before a submission deadline.

    -

    Mixing Matlab code

    +

    Mixing Matlab code

    Libigl can be interfaced with Matlab to offload numerically heavy computation to a Matlab script. The major advantage of this approach is that you will be @@ -2840,7 +2880,7 @@ L = sparse(LIJV(:,1),LIJV(:,2),LIJV(:,3));

    which is easily copied and pasted into Matlab for debugging, etc.

    -

    Calling libigl functions from Matlab

    +

    Calling libigl functions from Matlab

    It is also possible to call libigl functions from matlab, compiling them as MEX functions. This can be used to offload to C++ code the computationally @@ -2851,7 +2891,7 @@ We plan to provide wrappers for all our functions in the future, if you are interested in this feature (or if you want to help implementing it) please let us know.

    -

    Triangulation of closed polygons

    +

    Triangulation of closed polygons

    The generation of high-quality triangle and tetrahedral meshes is a very common task in geometry processing. We provide wrappers in libigl to @@ -2876,7 +2916,7 @@ in its interior) is triangulated.

    Triangulation of the interior of a polygon.
    -

    Tetrahedralization of closed surfaces

    +

    Tetrahedralization of closed surfaces

    Similarly, the interior of a closed manifold surface can be tetrahedralized using the function igl::tetrahedralize which wraps the Tetgen library (Example @@ -2890,7 +2930,7 @@ using the function igl::tetrahedralize which wraps the Tetgen libra

    Tetrahedralization of the interior of a surface mesh.
    -

    Baking ambient occlusion

    +

    Baking ambient occlusion

    Ambient occlusion is a rendering technique used to calculate the exposure of each point in a surface @@ -2927,7 +2967,7 @@ occlusion." /> occlusion.

    -

    Screen Capture

    +

    Screen Capture

    Libigl supports read and writing to .png files via the stb image code.

    @@ -2951,14 +2991,14 @@ igl::png::writePNG(R,G,B,A,"out.png");

    In Example 607 a scene is rendered in a temporary png and used to texture a quadrilateral.

    -

    Locally Injective Maps

    +

    Locally Injective Maps

    Extreme deformations or parametrizations with high-distortion might flip elements. This is undesirable in many applications, and it is possible to avoid it by introducing a non-linear constraints that guarantees that the area of every element remain positive.

    -

    Libigl can be used to compute Locally Injective Maps [34] using a variety of +

    Libigl can be used to compute Locally Injective Maps (35)[] using a variety of deformation energies. A simple deformation of a 2D grid is computed in Example 608.

    @@ -2969,7 +3009,7 @@ editing plus the anti-flipping constraints (right)." /> editing plus the anti-flipping constraints (right).
    -

    Boolean operations on meshes

    +

    Boolean operations on meshes

    Constructive solid geometry (CSG) is a technique to define a complex surface as the result of a number of set operations on solid regions of space: union, @@ -2983,23 +3023,27 @@ function \(a(\mathbf{x})\) “returns true”. operations are straightforward. For example, the union of solids \(A\) and \(B\) is simply

    -

    \(A \cup B = \{\mathbf{x} \left.\right| - a(\mathbf{x}) \text{ or } b(\mathbf{x})\},\)

    +

    $A \cup B = {\mathbf{x} \left.\right|

    + +

    a(\mathbf{x}) \text{ or } b(\mathbf{x})},$

    the intersection is

    -

    \(A \cap B = \{\mathbf{x} \left.\right| - a(\mathbf{x}) \text{ and } b(\mathbf{x})\},\)

    +

    $A \cap B = {\mathbf{x} \left.\right|

    + +

    a(\mathbf{x}) \text{ and } b(\mathbf{x})},$

    the difference \(A\) minus \(B\) is

    -

    \(A \setminus B = \{\mathbf{x} \left.\right| - a(\mathbf{x}) \text{ and _not_ } b(\mathbf{x})\},\)

    +

    $A \setminus B = {\mathbf{x} \left.\right|

    + +

    a(\mathbf{x}) \text{ and not } b(\mathbf{x})},$

    and the symmetric difference (XOR) is

    -

    \(A \triangle B = \{\mathbf{x} \left.\right| - \text{either } a(\mathbf{x}) \text{ or } b(\mathbf{x}) \text{ but not both }\}.\)

    +

    $A \triangle B = {\mathbf{x} \left.\right|

    + +

    \text{either } a(\mathbf{x}) \text{ or } b(\mathbf{x}) \text{ but not both }}.$

    Stringing together many of these operations, one can design quite complex shapes. A typical CSG library might only keep explicit base-case @@ -3012,7 +3056,7 @@ compute robustly with boundary representations, but are nonetheless useful.

    To compute a boolean operation on a triangle mesh with vertices VA and triangles FA and another mesh VB and FB, libigl first computes a unified -“mesh arrangement” (see [35]) with vertices V and triangles F where all triangle-triangle +“mesh arrangement” (see (36)[]) with vertices V and triangles F where all triangle-triangle intersections have been “resolved”. That is, edges and vertices are added exactly at the intersection lines, so the resulting non-manifold mesh (V,F) has no self-intersections.

    @@ -3031,10 +3075,10 @@ intersection) the boundary of the corresponding cells are extracted.

    The following figure shows each boolean operation on two meshes.

    -The example Boolean conducts +<img src=
    The example Boolean conducts boolean operations on the Cheburashka (red) and Knight (green). From left @@ -3056,7 +3100,7 @@ together coincident vertices, maintaining original triangle orientations.

    cork, which is typically faster, but is not always robust.

    -

    CSG Tree

    +

    CSG Tree

    The previous section discusses using igl::copyleft::cgal::mesh_boolean to compute the result of a single boolean @@ -3092,13 +3136,13 @@ operations

    then the final composite.

    -Example 610 computes complex CSG Tree operation on 5 +<img src=
    Example 610 computes complex CSG Tree operation on 5 input meshes.
    -

    Mesh Statistics

    +

    Mesh Statistics

    Libigl contains various mesh statistics, including face angles, face areas and the detection of singular vertices, which are vertices with more or less than 6 @@ -3133,7 +3177,7 @@ the angles are to 60 degrees the more stable will the optimization be. In this case, it is clear that the mesh is of bad quality and it will probably result in artifacts if used for solving PDEs.

    -

    Generalized Winding Number

    +

    Generalized Winding Number

    The problem of tetrahedralizing the interior of closed watertight surface mesh is a difficult, but well-posed problem (see our Tetgen wrappers). But @@ -3155,7 +3199,7 @@ mesh and which are outside. That is, which should be kept and which should be removed.

    The “Generalized Winding Number” is a robust method for determined -inside and outside for troublesome meshes [36]. The generalized +inside and outside for troublesome meshes (37)[]. The generalized winding number with respect to (V,F) at some point \(\mathbf{p} \in \mathcal{R}^3\) is defined as scalar function:

    @@ -3174,7 +3218,7 @@ oriented), then \(w(\mathbf{p})\) tends smoothly towar more inside (V,F), and toward 0 as \(\mathbf{p}\) is more outside.

    -Example 702 computes the +<img src=Mesh Decimation +

    Mesh Decimation

    The study of mesh simplification or decimation is nearly as old as meshes themselves. Given a high resolution mesh with too many triangles, find a “well @@ -3198,14 +3242,14 @@ methods are fairly advanced.

    One family of mesh decimation methods operates by successively remove elements from the mesh. In particular, Hoppe advocates for successively remove or rather -collapsing edges [37]. The generic form of this technique is to +collapsing edges (38)[]. The generic form of this technique is to construct a sequence of n meshes from the initial high-resolution mesh \(M_0\) to the lowest resolution mesh \(M_n\) by collapsing a single edge:

    \(M_0 \mathop{\longrightarrow}_\text{edge collapse} - M_1 \mathop{\longrightarrow}_\text{edge collapse} - \dots \mathop{\longrightarrow}_\text{edge collapse} - M_{n-1} \mathop{\longrightarrow}_\text{edge collapse} M_n.\)

    +M_1 \mathop{\longrightarrow}_\text{edge collapse} +\dots \mathop{\longrightarrow}_\text{edge collapse} +M_{n-1} \mathop{\longrightarrow}_\text{edge collapse} M_n.\)

    Hoppe’s original method and subsequent follow-up works propose various ways to choose the next edge to collapse in this sequence. Using a cost-based paradigm, @@ -3265,13 +3309,13 @@ drops below m=1000.

    One can also scratch deeper inside the decimation loop and call igl::collapse_edge directly. In order to operate efficiently, this routine needs more than the usual (V,F) mesh representation. We need E a list of -edge indices, where E.row(i) --> [s,d]; we need EMAP which maps the +edge indices, where E.row(i) --> [s,d]; we need EMAP which maps the “half”-edges of each triangle in F to its corresponding edge in E so that -E.row(EMAP(f+i*F.rows)) --> [s,d] if the edge across from the ith corner of the +E.row(EMAP(f+i*F.rows)) --> [s,d] if the edge across from the ith corner of the fth face is [s,d] (up to orientation); we need EF and EI which keep track of the faces incident on each edge and across from which corner of those faces the edges appears, so that EF(e,o) = f and EI(e,o) = i means that the edge -E.row(e) --> [s,d] appears in the fth face across from its ith corner (for +E.row(e) --> [s,d] appears in the fth face across from its ith corner (for o=0 the edge orientations should match, for o=1 the orientations are opposite).

    @@ -3324,7 +3368,7 @@ model. queue based approach with the simple shortest-edge-midpoint cost/placement strategy discussed above.

    -

    Signed Distances

    +

    Signed Distances

    In the Generalized Winding Number section, we examined a robust method for determining whether points lie inside or outside @@ -3399,15 +3443,15 @@ tree.squared_distance(V,F,P,sqrD,I,C);

    Finally, from the closest point or the winding number it’s possible to sign this distance. In igl::signed_distance we provide two methods for signing: -the so-called “pseudo-normal test” [38] and the generalized -winding number [36].

    +the so-called “pseudo-normal test” (39)[] and the generalized +winding number (37)[].

    The pseudo-normal test (see also igl::pseudonormal_test) assumes the input mesh is a watertight (closed, non-self-intersecting, manifold) mesh. Then given a query point \(\mathbf{q}\) and its closest point \(\mathbf{c} \in (V,F)\), it carefully chooses an outward normal \(\mathbf{n}\) at \(\mathbf{c}\) so that \(\text{sign}(\mathbf{q}-\mathbf{c})\cdot \mathbf{n}\) reveals whether -\(\mathbf{q}\) is inside \((V,F)\): –1, or outside: +1. This is a fast \(O(1)\) test +\(\mathbf{q}\) is inside \((V,F)\): -1, or outside: +1. This is a fast \(O(1)\) test once \(\mathbf{c}\) is located, but may fail if V,F is not watertight.

    An alternative is to use the generalized winding @@ -3428,13 +3472,13 @@ contains signed (unsquared) distances and the extra output N (only with the pseudo-normal test.

    -Example 704 computes signed distance on +<img src=
    Example 704 computes signed distance on slices through the bunny.
    -

    Marching Cubes

    +

    Marching Cubes

    Often 3D data is captured as scalar field defined over space \(f(\mathbf{x}) : \mathcal{R}^3 \rightarrow \mathcal{R}\). Lurking within this field, @@ -3444,7 +3488,7 @@ iso-surface at value \(v\) is composed of all points < processing is to extract an iso-surface as a triangle mesh for further mesh-based processing or visualization. This is referred to as iso-contouring.

    -

    “Marching Cubes” [39] is a famous +

    “Marching Cubes” (40) is a famous method for iso-contouring tri-linear functions \(f\) on a regular lattice (aka grid). The core idea of this method is to contour the iso-surface passing through each cell (if it does at @@ -3461,7 +3505,7 @@ input scalar field S sampled at vertex locations GV of

    -(Example 705) samples signed distance to the +<img src=Facet Orientation +

    Facet Orientation

    Models from the web occasionally arrive unorientated in the sense that the orderings of each triangles vertices do not consistently agree. Determining @@ -3495,7 +3539,7 @@ enforce a consistent facet orientation in the output faces FF.

    For (closed or nearly closed) surfaces representing the boundary of a solid object, libigl provides a routine to reorient faces so that the vertex ordering corresponds to a counter-clockwise ordering of the vertices with a -right-hand-rule normal pointing outward. This method [40] assumes +right-hand-rule normal pointing outward. This method (41)[] assumes that most of the universe is empty. That is, most points in space are outside of the solid object than inside. @@ -3510,10 +3554,10 @@ occluded (lighter, i.e., facing more void space).

    The boolean vector I reveals which rows of F have been flipped in FF.

    -(Example 706) loads a truck model with +<img src=
    (Example 706) loads a truck model with inconsistent orientations (back facing triangles shown darker). Orientable @@ -3522,7 +3566,7 @@ Alternatively, each individual triangle is considered a “patch” (mid and oriented outward independently.
    -

    Swept Volume

    +

    Swept Volume

    The swept volume \(S\) of a moving solid object \(A\) can be defined as any point in space such that at one moment in time the point lies inside the solid. In other @@ -3547,20 +3591,21 @@ error.

    volume is by exploiting an alternative definition of the swept volume based on signed distances:

    -

    \(S = \left\{ \mathbf{p}\ \middle| \ d(\mathbf{p},\partial S) < 0 \right\} = \left\{ \mathbf{p}\ -\middle|\ -\min\limits_{t \in [0,1]} d(\mathbf{p},f(t)\ \partial A) < 0 \right\}\)

    +

    $S = \left{ \mathbf{p} \middle|  d(\mathbf{p},\partial S) < 0 \right} = \left{ \mathbf{p}
    +\middle|

    + +

    \min\limits_{t \in [0,1]} d(\mathbf{p},f(t) \partial A) < 0 \right}$

    If \(\partial A\) is a triangle mesh, then we can approximate this by 1) discretizing time at a finite step of steps \([0,\Delta t,2\Delta t, \dots, 1]\) and by 2) discretizing space with a regular grid and representing the distance field using trilinear interpolation of grid values. Finally the output mesh, \(\partial S\) is approximated by contouring using Marching Cubes -[39].

    +(40).

    This method is similar to one described by Schroeder et al. in 1994 -[41], and the one used in conjunction with boolean operations by -Garg et al. 2016 [42].

    +(42), and the one used in conjunction with boolean operations by +Garg et al. 2016 (43).

    In libigl, if your input solid’s surface is represented by (V,F) then the output surface mesh will be (SV,SF) after calling:

    @@ -3573,7 +3618,7 @@ volume, greater than zero to approximate a positive offset of the swept volume or less than zero to approximate a negative offset.

    -(Example 707) computes +<img src=
    (Example 707) computes @@ -3581,7 +3626,7 @@ the surface of the swept volume (silver) of the bunny model undergoing a rigid motion (gold).
    -

    Picking

    +

    Picking

    Picking vertices and faces using the mouse is very common in geometry processing applications. While this might seem a simple operation, its @@ -3610,15 +3655,15 @@ Hierarchy constructed by Embree, and fid and vid are the picked face and vertex, respectively.

    -(Example 708) Picking via ray casting. The selected +<img src=
    (Example 708) Picking via ray casting. The selected vertices are colored in red.
    -

    Vector Field Visualization

    +

    Vector Field Visualization

    -

    Vector fields on surfaces are commonly visualized by tracing streamlines. Libigl +

    Vector fields on surfaces are commonly visualized by tracing [streamlines] (https://en.wikipedia.org/wiki/Streamlines,_streaklines,_and_pathlines). Libigl supports the seeding and tracing of streamlines, for both simple vector fields and for N-rosy fields. The seeds for the streamlines are initialized using streamlines_init, and the lines are traced using streamlines_next. Each call to streamlines_next extends @@ -3626,13 +3671,13 @@ each line by one triangle, allowing interactive rendering of the traced lines, a in Example 709.

    -(Example 709) Interactive streamlines tracing. +([Example 709](709_VectorFieldVisualizer/main.cpp)) Interactive streamlines tracing.
    (Example 709) Interactive streamlines tracing.
    -

    Scalable Locally Injective Maps

    +

    Scalable Locally Injective Maps

    -

    The Scalable Locally Injective Maps [43] algorithm allows to +

    The Scalable Locally Injective Maps (44) algorithm allows to compute locally injective maps on massive datasets. The algorithm shares many similarities with ARAP, but uses a reweighting scheme to minimize arbitrary distortion energies, including those that prevent the introduction of flips.

    @@ -3652,7 +3697,7 @@ using the SLIM algorithm in 10 iterations." /> using the SLIM algorithm in 10 iterations.
    -

    Subdivision surfaces

    +

    Subdivision surfaces

    Given a coarse mesh (aka cage) with vertices V and faces F, one can createa higher-resolution mesh with more vertices and faces by subdividing every @@ -3676,7 +3721,7 @@ a finer and finer mesh.

    The subdivision method of igl::loop is not in plane. The vertices of the refined mesh are moved to weight combinations of their neighbors: the mesh is -smoothed as it is refined [44]. This and other smooth subdivision +smoothed as it is refined (45). This and other smooth subdivision methods can be understood as generalizations of spline curves to surfaces. In particular the Loop subdivision method will converge to a \(C^1\) surface as we consider the limit of recursive applications of subdivision. Away from @@ -3695,14 +3740,14 @@ the carrier surfaces with extreme bias.

    The original coarse mesh and three different subdivision methods: -igl::upsample, igl::loop and -igl::false_barycentric_subdivision. +`igl::upsample`, `igl::loop` and +`igl::false_barycentric_subdivision`." />
    The original coarse mesh and three different subdivision methods: igl::upsample, igl::loop and igl::false_barycentric_subdivision.
    -

    Data smoothing

    +

    Data smoothing

    A noisy function \(f\) defined on a surface \(\Omega\) can be smoothed using an energy minimization that balances a smoothing term \(E_S\) with a quadratic @@ -3733,7 +3778,7 @@ of the function to be perpendicular to the boundary, the Hessian energy gives an unbiased result.

    -(Example 712) From left to right: a function +<img src= @@ -3743,31 +3788,31 @@ with the Laplacian energy and zero Neumann boundary conditions, and the result of smoothing with the Hessian energy.
    -

    Miscellaneous

    +

    Miscellaneous

    Libigl contains a wide variety of geometry processing tools and functions for dealing with meshes and the linear algebra related to them: far too many to discuss in this introductory tutorial. We’ve pulled out a couple of the interesting functions in this chapter to highlight.

    -

    Outlook for continuing development

    +

    Outlook for continuing development

    Libigl is in active development, and we plan to focus on the following features in the next months:

    • A better and more consistent documentation, plus extending this tutorial - to cover more libigl features.

      • +to cover more libigl features.

    • Implement a mixed-integer solver which only uses Eigen to remove the - dependency on CoMiSo.

      • +dependency on CoMiSo.

    • Improve the robustness and performance of the active set QP solver. In - particular, handle linearly dependent constraints.

      • +particular, handle linearly dependent constraints.

    • Implement more mesh analysis functions, including structural analysis for - masonry and 3D-printability analysis.

      • +masonry and 3D-printability analysis.

    • Increase support for point clouds and general polygonal meshes.

      • -

    • What would you like to see in libigl? Contact - us! or post a feature - request.

      • +

    • What would you like to see in libigl? Contact +us! or post a feature +request.

    We encourage you to contribute to the library and to report problems and bugs. @@ -3776,262 +3821,310 @@ repository and to open a our github repository.

    -
    +
      -

    1. Mark Meyer, Mathieu Desbrun, Peter Schröder and Alan H. Barr, - Discrete Differential-Geometry Operators for Triangulated - 2-Manifolds, - 2003.

      +
    2. +

      Mark Meyer, Mathieu Desbrun, Peter Schröder and Alan H. Barr, +Discrete Differential–Geometry Operators for Triangulated +2–Manifolds, +2003.  ↩

      3. -

    3. Daniele Panozzo, Enrico Puppo, Luigi Rocca, Efficient - Multi-scale Curvature and Crease - Estimation, - 2010.

      +
    4. +

      Daniele Panozzo, Enrico Puppo, Luigi Rocca, Efficient +Multi–scale Curvature and Crease +Estimation, +2010.  ↩

      5. -

    5. Alec Jacobson, - Algorithms and Interfaces for Real-Time Deformation of 2D and 3D - Shapes, - 2013.

      +
    6. +

      Alec Jacobson, +Algorithms and Interfaces for Real–Time Deformation of 2D and 3D +Shapes, +2013.  ↩

      7. -

    7. Andrei Sharf, Thomas Lewiner, Gil Shklarski, Sivan Toledo, and - Daniel Cohen-Or. Interactive topology-aware surface - reconstruction, - 2007.

      +
    8. +

      Andrei Sharf, Thomas Lewiner, Gil Shklarski, Sivan Toledo, and +Daniel Cohen–Or. Interactive topology–aware surface +reconstruction, +2007.  ↩

      9. -

    9. Michael Kazhdan, Jake Solomon, Mirela Ben-Chen, - Can Mean-Curvature Flow Be Made - Non-Singular, - 2012.

      +
    10. +

      Michael Kazhdan, Jake Solomon, Mirela Ben–Chen, +Can Mean–Curvature Flow Be Made +Non–Singular, +2012.  ↩

      11. -

    11. Raid M. Rustamov, Multiscale Biharmonic - Kernels, 2011.

      +
    12. +

      Joseph S. B. Mitchell, David M. Mount, Christos H. Papadimitriou. The Discrete Geodesic Problem, 1987  ↩

      13. -

    13. Bruno Vallet and Bruno Lévy. Spectral Geometry Processing with - Manifold - Harmonics, - 2008.

      +
    14. +

      Raid M. Rustamov, Multiscale Biharmonic +Kernels, 2011.  ↩

      15. -

    15. Klaus Hildebrandt, Christian Schulz, Christoph von - Tycowicz, and Konrad Polthier. Interactive Surface Modeling using Modal - Analysis, - 2011.

      +
    16. +

      Bruno Vallet and Bruno Lévy. Spectral Geometry Processing with +Manifold +Harmonics, +2008.  ↩

      17. -

    17. Jernej Barbic and Doug James. Real-Time Subspace Integration - for St.Venant-Kirchhoff Deformable - Models, - 2005.

      +
    18. +

      Klaus Hildebrandt, Christian Schulz, Christoph von +Tycowicz, and Konrad Polthier. Interactive Surface Modeling using Modal +Analysis, +2011.  ↩

      19. -

    19. Matrio Botsch and Leif Kobbelt. - An Intuitive Framework for Real-Time Freeform - Modeling, - 2004.

      +
    20. +

      Jernej Barbic and Doug James. Real–Time Subspace Integration +for St.Venant–Kirchhoff Deformable +Models, +2005.  ↩

      21. -

    21. Alec Jacobson, Elif Tosun, Olga Sorkine, and Denis - Zorin. Mixed Finite Elements for Variational Surface - Modeling, - 2010.

      +
    22. +

      Matrio Botsch and Leif Kobbelt. +An Intuitive Framework for Real–Time Freeform +Modeling, +2004.  ↩

      23. -

    23. Olga Sorkine, Yaron Lipman, Daniel Cohen-Or, Marc Alexa, - Christian Rössl and Hans-Peter Seidel. Laplacian Surface - Editing, 2004.

      +
    24. +

      Alec Jacobson, Elif Tosun, Olga Sorkine, and Denis +Zorin. Mixed Finite Elements for Variational Surface +Modeling, +2010.  ↩

      25. -

    25. Alec Jacobson, Ilya Baran, Jovan Popović, and Olga Sorkine. - Bounded Biharmonic Weights for Real-Time - Deformation, - 2011.

      +
    26. +

      Olga Sorkine, Yaron Lipman, Daniel Cohen–Or, Marc Alexa, +Christian Rössl and Hans–Peter Seidel. Laplacian Surface +Editing, 2004.  ↩

      27. -

    27. Ladislav Kavan, Steven Collins, Jiri Zara, and Carol O’Sullivan. - Geometric Skinning with Approximate Dual Quaternion - Blending, - 2008.

      +
    28. +

      Alec Jacobson, Ilya Baran, Jovan Popović, and Olga Sorkine. +Bounded Biharmonic Weights for Real–Time +Deformation, +2011.  ↩

      29. -

    29. Olga Sorkine and Marc Alexa. As-rigid-as-possible Surface - Modeling, 2007.

      +
    30. +

      Ladislav Kavan, Steven Collins, Jiri Zara, and Carol O'Sullivan. +Geometric Skinning with Approximate Dual Quaternion +Blending, +2008.  ↩

      31. -

    31. Isaac Chao, Ulrich Pinkall, Patrick Sanan, Peter Schröder. - A Simple Geometric Model for Elastic - Deformations, - 2010.

      +
    32. +

      Olga Sorkine and Marc Alexa. As–rigid–as–possible Surface +Modeling, 2007.  ↩

      33. -

    33. Alexa McAdams, Andrew Selle, Rasmus Tamstorf, Joseph Teran, - Eftychios Sifakis. Computing the Singular Value Decomposition of 3x3 - matrices with minimal branching and elementary floating point - operations, - 2011.

      +
    34. +

      Isaac Chao, Ulrich Pinkall, Patrick Sanan, Peter Schröder. +A Simple Geometric Model for Elastic +Deformations, +2010.  ↩

      35. -

    35. Alec Jacobson, Ilya Baran, Ladislav Kavan, Jovan Popović, and - Olga Sorkine. Fast Automatic Skinning - Transformations, - 2012.

      +
    36. +

      Alexa McAdams, Andrew Selle, Rasmus Tamstorf, Joseph Teran, +Eftychios Sifakis. Computing the Singular Value Decomposition of 3x3 +matrices with minimal branching and elementary floating point +operations, +2011.  ↩

      37. -

    37. Alec Jacobson, Zhigang Deng, Ladislav Kavan, - J.P. Lewis. Skinning: Real-Time Shape - Deformation, - 2014.

      +
    38. +

      Alec Jacobson, Ilya Baran, Ladislav Kavan, Jovan Popović, and +Olga Sorkine. Fast Automatic Skinning +Transformations, +2012.  ↩

      39. -

    39. Yu Wang, Alec Jacobson, Jernej Barbic, Ladislav Kavan. Linear - Subspace Design for Real-Time Shape - Deformation, - 2015

      +
    40. +

      Alec Jacobson, Zhigang Deng, Ladislav Kavan, +J.P. Lewis. Skinning: Real–Time Shape +Deformation, +2014.  ↩

      41. -

    41. Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, Werner - Stuetzle. Multiresolution Analysis of Arbitrary - Meshes, 2005.

      +
    42. +

      Yu Wang, Alec Jacobson, Jernej Barbic, Ladislav Kavan. Linear +Subspace Design for Real–Time Shape +Deformation, +2015  ↩

      43. -

    43. Bruno Lévy, Sylvain Petitjean, Nicolas Ray, Jérome Maillot. - Least Squares Conformal Maps, for Automatic Texture Atlas - Generation,, 2002.

      +
    44. +

      Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, Werner +Stuetzle. Multiresolution Analysis of Arbitrary +Meshes, 2005.  ↩

      45. -

    45. Patrick Mullen, Yiying Tong, Pierre Alliez, Mathieu Desbrun. - Spectral Conformal - Parameterization, 2008.

      +
    46. +

      Bruno Lévy, Sylvain Petitjean, Nicolas Ray, Jérome Maillot. +Least Squares Conformal Maps, for Automatic Texture Atlas +Generation,, 2002.  ↩

      47. -

    47. Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, Steven J. Gortler. - A Local/Global Approach to Mesh - Parameterization, 2008.

      +
    48. +

      Patrick Mullen, Yiying Tong, Pierre Alliez, Mathieu Desbrun. +Spectral Conformal +Parameterization, 2008.  ↩

      49. -

    49. Nicolas Ray, Bruno Vallet, Wan Chiu Li, Bruno Lévy. - N-Symmetry Direction Field - Design, - 2008.

      +
    50. +

      Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, Steven J. Gortler. +A Local/Global Approach to Mesh +Parameterization, 2008.  ↩

      51. -

    51. David Bommes, Henrik Zimmer, Leif Kobbelt. - Mixed-integer - quadrangulation, - 2009.

      +
    52. +

      Nicolas Ray, Bruno Vallet, Wan Chiu Li, Bruno Lévy. +N–Symmetry Direction Field +Design, +2008.  ↩

      53. -

    53. Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter - Schröder. Globally Optimal Direction - Fields, - 2013.

      +
    54. +

      David Bommes, Henrik Zimmer, Leif Kobbelt. +Mixed–integer +quadrangulation, +2009.  ↩

      55. -

    55. Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga - Sorkine-Hornung. Frame Fields: Anisotropic and Non-Orthogonal Cross - Fields, - 2014.

      +
    56. +

      Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter +Schröder. Globally Optimal Direction +Fields, +2013.  ↩

      57. -

    57. Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga - Sorkine-Hornung. Designing N-PolyVector Fields with Complex - Polynomials, 2014

      +
    58. +

      Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga +Sorkine–Hornung. Frame Fields: Anisotropic and Non–Orthogonal Cross +Fields, +2014.  ↩

      59. -

    59. Yang Liu, Weiwei Xu, Jun Wang, Lifeng Zhu, Baining Guo, Falai Chen, Guoping - Wang. General Planar Quadrilateral Mesh Design Using Conjugate Direction - Field, - 2008.

      +
    60. +

      Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga +Sorkine–Hornung. Designing N–PolyVector Fields with Complex +Polynomials, 2014  ↩

      61. -

    61. Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly - Shape-Up: Shaping Discrete Geometry with - Projections, 2012

      +
    62. +

      Yang Liu, Weiwei Xu, Jun Wang, Lifeng Zhu, Baining Guo, Falai Chen, Guoping +Wang. General Planar Quadrilateral Mesh Design Using Conjugate Direction +Field, +2008.  ↩

      63. -

    63. Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga - Sorkine-Hornung. Integrable PolyVector Fields, 2015

      +
    64. +

      Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly +Shape–Up: Shaping Discrete Geometry with +Projections, 2012  ↩

      65. -

    65. Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, - David Bommes, Klaus Hildebrandt, Mirela Ben-Chen. Directional Field - Synthesis, Design, and - Processing, - 2016

      +
    66. +

      Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga +Sorkine–Hornung. Integrable PolyVector Fields, 2015  ↩

      67. -

    67. Christian Schüller, Ladislav Kavan, Daniele Panozzo, Olga - Sorkine-Hornung. Locally Injective - Mappings, 2013.

      +
    68. +

      Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, +David Bommes, Klaus Hildebrandt, Mirela Ben–Chen. Directional Field +Synthesis, Design, and +Processing, +2016  ↩

      69. -

    69. Qingnan Zhou, Eitan Grinspun, Denis Zorin. Mesh Arrangements for - Solid - Geometry, - 2016

      +
    70. +

      Christian Schüller, Ladislav Kavan, Daniele Panozzo, Olga +Sorkine–Hornung. Locally Injective +Mappings, 2013.  ↩

      71. -

    71. Alec Jacobson, Ladislav Kavan, and Olga Sorkine. - Robust Inside-Outside Segmentation using Generalized Winding - Numbers, - 2013.

      +
    72. +

      Qingnan Zhou, Eitan Grinspun, Denis Zorin. Mesh Arrangements for +Solid +Geometry, +2016  ↩

      73. -

    73. Hugues Hoppe. Progressive - Meshes, 1996

      +
    74. +

      Alec Jacobson, Ladislav Kavan, and Olga Sorkine. +Robust Inside–Outside Segmentation using Generalized Winding +Numbers, +2013.  ↩

      75. -

    75. J Andreas Baerentzen and Henrik Aanaes. +

    76. +

      Hugues Hoppe. Progressive +Meshes, 1996  ↩

      +
      77. + +
    77. +

      J Andreas Baerentzen and Henrik Aanaes. Signed distance computation using the angle weighted pseudonormal, - 2005.

      +2005.  ↩

      78. -

    78. W.E. Lorensen and Harvey E. Cline. Marching cubes: A high - resolution 3d surface construction - algorithm, - 1987.

      +
    79. +

      W.E. Lorensen and Harvey E. Cline. Marching cubes: A high +resolution 3d surface construction +algorithm, +1987.  ↩

      80. -

    80. Kenshi Takayama, Alec Jacobson, Ladislav Kavan, Olga - Sorkine-Hornung. A Simple Method for Correcting Facet Orientations in - Polygon Meshes Based on Ray - Casting, - 2014.

      +
    81. +

      Kenshi Takayama, Alec Jacobson, Ladislav Kavan, Olga +Sorkine–Hornung. A Simple Method for Correcting Facet Orientations in +Polygon Meshes Based on Ray +Casting, +2014.  ↩

      82. -

    82. William J. Schroeder, William E. Lorensen, and Steve - Linthicum. Implicit Modeling of Swept Surfaces and - Volumes, - 1994.

      +
    83. +

      William J. Schroeder, William E. Lorensen, and Steve +Linthicum. Implicit Modeling of Swept Surfaces and +Volumes, +1994.  ↩

      84. -

    84. Akash Garg, Alec Jacobson, Eitan Grinspun. Computational Design - of - Reconfigurables, - 2016

      +
    85. +

      Akash Garg, Alec Jacobson, Eitan Grinspun. Computational Design +of +Reconfigurables, +2016  ↩

      86. -

    86. Michael Rabinovich, Roi Poranne, Daniele Panozzo, Olga - Sorkine-Hornung. Scalable Locally Injective - Mappings, 2016.

      +
    87. +

      Michael Rabinovich, Roi Poranne, Daniele Panozzo, Olga +Sorkine–Hornung. Scalable Locally Injective +Mappings, 2016.  ↩

      88. -

    88. Charles Loop. Smooth Subdivision Surfaces Based on - Triangles, - 1987.

      +
    89. +

      Charles Loop. Smooth Subdivision Surfaces Based on +Triangles, +1987.  ↩

    - + diff --git a/tutorial/tutorial.md b/tutorial/tutorial.md index 96b0c1c12..41fe9810e 100644 --- a/tutorial/tutorial.md +++ b/tutorial/tutorial.md @@ -44,10 +44,11 @@ lecture notes links to a cross-platform example application. * [202 Gaussian Curvature](#gaussiancurvature) * [203 Curvature Directions](#curvaturedirections) * [204 Gradient](#gradient) - * [204 Laplacian](#laplacian) + * [205 Laplacian](#laplacian) * [Mass matrix](#massmatrix) * [Alternative construction of Laplacian](#alternativeconstructionoflaplacian) + * [206 Geodesic Distance](#geodesic) * [Chapter 3: Matrices and Linear Algebra](#chapter3:matricesandlinearalgebra) * [301 Slice](#slice) * [302 Sort](#sort) @@ -790,6 +791,28 @@ since the Laplacian is the divergence of the gradient. Naturally, $\mathbf{G}^T$ a $n \times md$ sparse matrix which takes vector values stored at triangle faces to scalar divergence values at vertices. +## Geodesic + +The discrete geodesic distance between two points is the length of the shortest path between then restricted to the surface. For triangle meshes, such a path is made of a set of segments which can be either edges of the mesh or crossing a triangle. + +Libigl includes a wrapper for the exact geodesic algorithm [#mitchell_1987] developed by Danil Kirsanov (https://code.google.com/archive/p/geodesic/), exposing it through an Eigen-based API. The function +```cpp +igl::exact_geodesic(V,F,VS,FS,VT,FT,d); +``` +computes the closest geodesic distances of each vertex in VT or face in FT, from the source vertices VS or faces FS of the input mesh V,F. The output is writted in the vector d, which lists first the distances for the vertices in VT, and then for the faces in FT. For example, if you want to compute the distance from the vertex with id ```vid```, to all vertices of F you can use: +```cpp +Eigen::VectorXi VS,FS,VT,FT; +// The selected vertex is the source +VS.resize(1); +VS << vid; +// All vertices are the targets +VT.setLinSpaced(V.rows(),0,V.rows()-1); +Eigen::VectorXd d; +igl::exact_geodesic(V,F,VS,FS,VT,FT,d); +``` +![[Example 206](206_GeodesicDistance/main.cpp) allows to interactively pick the source vertex and displays the distance using a periodic color pattern. +](images/geodesicdistance.jpg) + # Chapter 3: Matrices and linear algebra Libigl relies heavily on the Eigen library for dense and sparse linear algebra routines. Besides geometry processing routines, libigl has linear algebra @@ -3524,3 +3547,4 @@ pseudonormal](https://www.google.com/search?q=Signed+distance+computation+using+ Solid Geometry](https://www.google.com/search?q=Mesh+Arrangements+for+Solid+Geometry), 2016 +[#mitchell_1987]: Joseph S. B. Mitchell, David M. Mount, Christos H. Papadimitriou. [The Discrete Geodesic Problem](https://www.google.com/search?q=The+Discrete+Geodesic+Problem), 1987