From 28b40e8d44a9239efccd369c134cadbdf854fd81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Silvia=20Sell=C3=A1n?= Date: Mon, 8 Aug 2022 14:32:09 -0700 Subject: [PATCH] Docstrings (#22) --- README.md | 12 +- src/gpytoolbox/array_correspondence.py | 39 +++-- src/gpytoolbox/boundary_edges.py | 27 ++- src/gpytoolbox/boundary_loops.py | 38 +++-- src/gpytoolbox/boundary_vertices.py | 26 ++- src/gpytoolbox/cotangent_laplacian.py | 27 ++- .../cotangent_laplacian_intrinsic.py | 31 ++-- src/gpytoolbox/cotangent_weights.py | 37 ++-- src/gpytoolbox/cotangent_weights_intrinsic.py | 40 +++-- src/gpytoolbox/decimate.py | 8 +- src/gpytoolbox/doublearea.py | 12 +- src/gpytoolbox/doublearea_intrinsic.py | 32 ++-- src/gpytoolbox/edges.py | 59 ++++--- src/gpytoolbox/grad.py | 12 +- src/gpytoolbox/halfedge_edge_map.py | 67 ++++---- src/gpytoolbox/halfedge_lengths.py | 36 ++-- src/gpytoolbox/halfedge_lengths_squared.py | 44 +++-- src/gpytoolbox/halfedges.py | 35 ++-- src/gpytoolbox/in_element_aabb.py | 8 +- src/gpytoolbox/massmatrix.py | 15 +- src/gpytoolbox/massmatrix_intrinsic.py | 30 +--- src/gpytoolbox/min_quad_with_fixed.py | 161 +++++++++++------- src/gpytoolbox/normalize_points.py | 2 +- src/gpytoolbox/offset_surface.py | 8 +- src/gpytoolbox/per_face_normals.py | 10 +- src/gpytoolbox/per_vertex_normals.py | 22 +-- src/gpytoolbox/ray_mesh_intersect.py | 8 +- src/gpytoolbox/read_mesh.py | 73 +++++--- src/gpytoolbox/regular_cube_mesh.py | 21 +-- src/gpytoolbox/subdivide.py | 54 +++--- src/gpytoolbox/tip_angles.py | 38 +++-- src/gpytoolbox/tip_angles_intrinsic.py | 42 +++-- src/gpytoolbox/triangle_triangle_adjacency.py | 45 +++-- src/gpytoolbox/write_mesh.py | 63 ++++--- test/test_quadtree_laplacian.py | 113 ++++++------ 35 files changed, 761 insertions(+), 534 deletions(-) diff --git a/README.md b/README.md index 5a130538..d8f0a499 100644 --- a/README.md +++ b/README.md @@ -189,16 +189,8 @@ I = in_element_aabb(queries,V,F) # This is a C++ binding ``` --> # TO-DO -## Must do before first PyPi release -- \[Oded\] Add docstrings for `array_correspondence`, `boundary_edges`, - `boundary_loops`, `boundary_vertices`, `cotangent_laplacian_intrinsic`, - `cotangent_laplacian`, `cotangent_weights_intrinsic`, `cotangent_weights`, - `doublearea_intrinsic`, `edges`, `halfedge_edge_map`, - `halfedge_lengths_squared`, `halfedge_lengths`, `halfedges`, - `min_quad_with_fixed`, `read_mesh`, `subdivide`, `tip_angles_intrinsic`, - `tip_angles`, `triangle_triangle_adjacency`, `write_mesh`. -- \[Oded\] Add FD-type linear solve with Dirichlet conditions function. -- \[Silvia\] When Oded does this^, fix `test_quadtree_laplacian.py`. + + ## Future to-dos - Add examples to docstrings. diff --git a/src/gpytoolbox/array_correspondence.py b/src/gpytoolbox/array_correspondence.py index 1487a73c..1fbefcc4 100644 --- a/src/gpytoolbox/array_correspondence.py +++ b/src/gpytoolbox/array_correspondence.py @@ -1,21 +1,30 @@ import numpy as np def array_correspondence(A,B,axis=None): - # Computes a map from the a to the equal elements in b - # - # Inputs: - # A: #a numpy array (must be 1-dim or 2-dim) - # B: #b numpy array (must be 1-dim or 2-dim) - # Optional: - # axis: If None, will treat A and B as flat arrays. - # If a number, will check for equality of the - # entire axis, in which case the dimension of - # A and B across that axis must be equal. - # Outputs: - # f #A numpy index list mapping from the a to b, with -1 if there is - # no matching entry. - # If b contains multiple eligible entries, return an arbitrary one. - # If there are no -1s, b[f] == a + """Computes a map from A to the equal elements in B + + Parameters + ---------- + A : (a,) or (a,k) numpy array (must be 1-dim or 2-dim) + B : (b,) or (b,k) numpy array (must be 1-dim or 2-dim) + axis : int or None, optional (default None) + If None, will treat A and B as flat arrays. + If a number, will check for equality of the entire axis, in which case + the dimension of A and B across that axis must be equal. + + Returns + ------- + f : (a,) numpy int array + index list mapping from A to B, with -1 if there is no + matching entry. + If b contains multiple eligible entries, return an arbitrary one. + If there are no -1s, `b[f] == a` + + Examples + -------- + TODO + + """ if axis is None: A = A.ravel() diff --git a/src/gpytoolbox/boundary_edges.py b/src/gpytoolbox/boundary_edges.py index 3c5f4fda..8d2c5ac1 100644 --- a/src/gpytoolbox/boundary_edges.py +++ b/src/gpytoolbox/boundary_edges.py @@ -1,14 +1,25 @@ from gpytoolbox.edges import edges def boundary_edges(F): - # Given a triangle mesh with face indices F, returns all unique oriented - # boundary edges as indices into the vertex array. - # Works only on manifold meshes. - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Outputs: - # bE #bE by 2 indices of boundary edges into the vertex array + """Given a triangle mesh with face indices F, returns all unique oriented + boundary edges as indices into the vertex array. + Works only on manifold meshes. + + Parameters + ---------- + F : (m,3) numpy int array. + face index list of a triangle mesh + + Returns + ------- + bE : (be,2) numpy int array. + indices of boundary edges into the vertex array + + Examples + -------- + TODO + + """ E,b = edges(F, return_boundary_indices=True) return E[b,:] diff --git a/src/gpytoolbox/boundary_loops.py b/src/gpytoolbox/boundary_loops.py index dbee4171..231ca93d 100644 --- a/src/gpytoolbox/boundary_loops.py +++ b/src/gpytoolbox/boundary_loops.py @@ -3,21 +3,29 @@ def boundary_loops(f, allow_wrong_orientations=True): - # Computes oriented boundary loop for each boundary component of a triangle - # mesh. - # - # This function only works on manifold triangle meshes. - # TODO: assert manifoldness when running this function - # - # Input: - # F #F by 3 face index list of a triangle mesh - # Optional: - # allow_wrong_orientations whether to allow F to contain - # wrongly oriented triangles - # - # Output: - # loops list of numpy array list of boundary vertices - # in oriented loops + """Computes oriented boundary loop for each boundary component of a triangle + mesh. + This function only works on manifold triangle meshes. + + TODO: assert manifoldness when running this function + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + allow_wrong_orientations: bool, optional (default True). + whether to allow F to contain wrongly oriented triangles + + Returns + ------- + loops : list of numpy arrays that are themselves lists of boundary vertices + in oriented loops + + Examples + -------- + TODO + + """ assert f.shape[0] > 0 assert f.shape[1] == 3 diff --git a/src/gpytoolbox/boundary_vertices.py b/src/gpytoolbox/boundary_vertices.py index 451edad4..3dd0ffa1 100644 --- a/src/gpytoolbox/boundary_vertices.py +++ b/src/gpytoolbox/boundary_vertices.py @@ -2,13 +2,25 @@ from gpytoolbox.boundary_edges import boundary_edges def boundary_vertices(F): - # Given a triangle mesh with face indices F, returns the indices of all - # boundary vertices. Works only on manifold meshes. - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Outputs: - # bV #bV list of indices into F of boundary vertices + """Given a triangle mesh with face indices F, returns the indices of all + boundary vertices. + Works only on manifold meshes. + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + bV : (bv,) numpy int array + list of indices into F of boundary vertices + + Examples + -------- + TODO + + """ return np.unique(boundary_edges(F)) diff --git a/src/gpytoolbox/cotangent_laplacian.py b/src/gpytoolbox/cotangent_laplacian.py index 59719f8b..63e0e90f 100644 --- a/src/gpytoolbox/cotangent_laplacian.py +++ b/src/gpytoolbox/cotangent_laplacian.py @@ -3,14 +3,25 @@ from gpytoolbox.cotangent_laplacian_intrinsic import cotangent_laplacian_intrinsic def cotangent_laplacian(V,F): - # Builds the (pos. def.) cotangent Laplacian for a triangle mesh. - # - # Input: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indeces into V - # - # Output: - # L #V by #V scipy csr_matrix cotangent Laplacian + """Builds the (pos. def.) cotangent Laplacian for a triangle mesh. + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + L : (n,n) scipy csr_matrix + cotangent Laplacian + + Examples + -------- + TODO + + """ l_sq = halfedge_lengths_squared(V,F) return cotangent_laplacian_intrinsic(l_sq,F,n=V.shape[0]) diff --git a/src/gpytoolbox/cotangent_laplacian_intrinsic.py b/src/gpytoolbox/cotangent_laplacian_intrinsic.py index ab90bc15..3f71667a 100644 --- a/src/gpytoolbox/cotangent_laplacian_intrinsic.py +++ b/src/gpytoolbox/cotangent_laplacian_intrinsic.py @@ -3,18 +3,25 @@ from gpytoolbox.cotangent_weights_intrinsic import cotangent_weights_intrinsic def cotangent_laplacian_intrinsic(l_sq,F,n=None): - # Builds the (pos. def.) cotangent Laplacian for a triangle mesh using only - # intrinsic information (squared halfedge edge lengths). - # - # Input: - # l_sq #F by 3 numpy array of squared halfedge lengths as computed - # by halfedge_lengths_squared - # F #F by 3 int numpy array of face/edge vertex indeces into V - # Optional: - # n an integer denoting the number of vertices in the mesh - # - # Output: - # L #V by #V scipy csr_matrix cotangent Laplacian + """Builds the (pos. def.) cotangent Laplacian for a triangle mesh. + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + L : (n,n) scipy csr_matrix + cotangent Laplacian + + Examples + -------- + TODO + + """ assert F.shape[1] == 3 assert l_sq.shape == F.shape diff --git a/src/gpytoolbox/cotangent_weights.py b/src/gpytoolbox/cotangent_weights.py index b78757e0..4c4f2c1d 100644 --- a/src/gpytoolbox/cotangent_weights.py +++ b/src/gpytoolbox/cotangent_weights.py @@ -3,19 +3,30 @@ from gpytoolbox.cotangent_weights_intrinsic import cotangent_weights_intrinsic def cotangent_weights(V,F): - # Builds the cotangent weights (cotangent/2) for each halfedge. - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Input: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indeces into V - # - # Output: - # C #F by 3 numpy array of cotangent weights for each halfedge + """Builds the cotangent weights (cotangent/2) for each halfedge. + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + C : (m,3) numpy array + cotangent weights for each halfedge + + Examples + -------- + TODO + + """ l_sq = halfedge_lengths_squared(V,F) return cotangent_weights_intrinsic(l_sq,F) diff --git a/src/gpytoolbox/cotangent_weights_intrinsic.py b/src/gpytoolbox/cotangent_weights_intrinsic.py index 62a11a44..63a3cffa 100644 --- a/src/gpytoolbox/cotangent_weights_intrinsic.py +++ b/src/gpytoolbox/cotangent_weights_intrinsic.py @@ -2,21 +2,31 @@ from gpytoolbox.doublearea_intrinsic import doublearea_intrinsic def cotangent_weights_intrinsic(l_sq,F): - # Builds the cotangent weights (cotangent/2) for each halfedge using only - # intrinsic information (squared halfedge edge lengths). - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Input: - # l_sq #F by 3 numpy array of squared halfedge lengths as computed - # by halfedge_lengths_squared - # F #F by 3 int numpy array of face/edge vertex indeces into V - # - # Output: - # C #F by 3 numpy array of cotangent weights for each halfedge + """Builds the cotangent weights (cotangent/2) for each halfedge using only + intrinsic information (squared halfedge edge lengths). + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + C : (m,3) numpy array + cotangent weights for each halfedge + + Examples + -------- + TODO + + """ assert F.shape[1] == 3 assert l_sq.shape == F.shape diff --git a/src/gpytoolbox/decimate.py b/src/gpytoolbox/decimate.py index b8568f9a..028c08e5 100644 --- a/src/gpytoolbox/decimate.py +++ b/src/gpytoolbox/decimate.py @@ -7,10 +7,10 @@ def decimate(V,F,face_ratio=0.1,num_faces=None): Parameters ---------- - V : numpy double array - Matrix of vertex coordinates - F : numpy int array - Matrix of triangle indices + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh face_ratio: double, optional (default 0.1) Desired ratio of output faces to input faces num_faces : int, optional (default None) diff --git a/src/gpytoolbox/doublearea.py b/src/gpytoolbox/doublearea.py index 0ae61939..c1530c07 100644 --- a/src/gpytoolbox/doublearea.py +++ b/src/gpytoolbox/doublearea.py @@ -8,14 +8,16 @@ def doublearea(V,F=None): Parameters ---------- - V : numpy double array - Matrix of vertex coordinates - F : numpy int array, optional - Matrix of triangle indices, None if ordered polyline (by default, None) + V : (n,d) numpy array + vertex list of a polyline or triangle mesh + F : numpy int array, optional (default: None) + if None or (m,2), interpret input as ordered polyline; + if (m,3) numpy int array, interpred as face index list of a triangle + mesh Returns ------- - dblA : numpy double array + dblA : (m,) numpy double array vector of twice the (unsigned) area/length See Also diff --git a/src/gpytoolbox/doublearea_intrinsic.py b/src/gpytoolbox/doublearea_intrinsic.py index cb18e776..a6b17288 100644 --- a/src/gpytoolbox/doublearea_intrinsic.py +++ b/src/gpytoolbox/doublearea_intrinsic.py @@ -2,16 +2,28 @@ def doublearea_intrinsic(l_sq,F): - # Construct the doublearea of each element of a triangle mesh using only - # intrinsic information (squared halfedge edge lengths). - # - # Input: - # l_sq #F by 3 numpy array of squared halfedge lengths as computed - # by halfedge_lengths_squared - # F #F by 3 int numpy array of face/edge vertex indices into V - # - # Output: - # A #F vector of twice the (unsigned) area + """Construct the doublearea of each element of a line or triangle mesh. + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + dblA : (m,) numpy double array + vector of twice the (unsigned) area/length + + See Also + -------- + doublearea. + + Examples + -------- + TO-DO + """ assert F.shape == l_sq.shape assert F.shape[1]==3 diff --git a/src/gpytoolbox/edges.py b/src/gpytoolbox/edges.py index 59eb99bc..105697e6 100644 --- a/src/gpytoolbox/edges.py +++ b/src/gpytoolbox/edges.py @@ -5,33 +5,38 @@ def edges(F, return_boundary_indices=False, return_interior_indices=False, return_nonmanifold_indices=False): - # Given a triangle mesh with face indices F, returns all unique unoriented - # edges as indices into the vertex array. - # There is no particular ordering convention for edges. - # Boundary edges are guaranteed to be oriented as in F. - # - # NOTE: "nonmanifold edges" here specifically mean edges that correspond to - # 3 or more halfedges, not any complex-destroying edge in general. - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Optional: - # return_boundary_indices whether to return a list of - # indices into E denoting the - # boundary edges - # return_interior_indices whether to return a list of - # indices into E denoting the - # interior edges - # return_nonmanifold_indices whether to return a list of - # indices into E denoting the - # nonmanifold edges - # Outputs: - # E #E by 2 indices of edges into the vertex array - # Optional: - # boundary_indices list of indices into E of boundary edges - # interior_indices list of indices into E of interior edges - # nonmanifold_indices list of indices into E of nonmanifold - # edges + """Given a triangle mesh with face indices F, returns all unique unoriented + edges as indices into the vertex array. + There is no particular ordering convention for edges. + Boundary edges are guaranteed to be oriented as in F. + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + return_boundary_indices : bool, optional (default False) + whether to return a list of indices into E denoting the boundary edges + return_interior_indices : bool, optional (default False) + whether to return a list of indices into E denoting the interior edges + return_nonmanifold_indices : bool, optional (default False) + whether to return a list of indices into E denoting the nonmanifold edges + + Returns + ------- + E : (e,2) numpy int array + indices of edges into the vertex array + boundary_indices : if requested, (b,) numpy int array + list of indices into E of boundary edges + interior_indices : if requested, (i,) numpy int array + list of indices into E of interior edges + nonmanifold_indices : if requested, (nm,) numpy int array + list of indices into E of nonmanifold edges + + Examples + -------- + TODO + + """ assert F.shape[0] > 0 assert F.shape[1] == 3 diff --git a/src/gpytoolbox/grad.py b/src/gpytoolbox/grad.py index 47ca1dc1..b0c15079 100644 --- a/src/gpytoolbox/grad.py +++ b/src/gpytoolbox/grad.py @@ -12,14 +12,16 @@ def grad(V,F=None): Parameters ---------- - V : numpy double array - Matrix of vertex coordinates - F : numpy int array, optional (default None) - Matrix of triangle indices + V : (n,d) numpy array + vertex list of a polyline or triangle mesh + F : numpy int array, optional (default: None) + if None, interpret input as ordered polyline; + if (m,3) numpy int array, interpred as face index list of a triangle + mesh Returns ------- - G : scipy sparse.csr_matrix + G : (d*m,n) scipy sparse.csr_matrix Sparse FEM gradient matrix See Also diff --git a/src/gpytoolbox/halfedge_edge_map.py b/src/gpytoolbox/halfedge_edge_map.py index 2839b498..a2fd1edf 100644 --- a/src/gpytoolbox/halfedge_edge_map.py +++ b/src/gpytoolbox/halfedge_edge_map.py @@ -3,35 +3,44 @@ from gpytoolbox.array_correspondence import array_correspondence def halfedge_edge_map(F, assume_manifold=True): - # Computes unique edge indices, and a unique map from halfedges to edges, - # as well as its inverse. - # There is no particular ordering convention for edges. - # Boundary edges are guaranteed to be oriented as in F. - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Optional: - # assume_manifold if this is true, will assume that F is - # manifold, and thus every edge can be - # incident to at most two halfedges. - # The algorithm is very slow if this is - # False. - # Outputs: - # he #F by 3 by 2 halfedge list as per above conventions - # E #E by 2 indices of edges into the vertex array - # he_to_E #F by 3 index map from he to corresponding row in E - # E_to_he if assume_manifold, #E by 2 by 2 index map from e to - # corresponding row and col in he for two halfedges (or -1 - # if only one halfedge exists) - # if not assume_manifold, python list of E k by 2 index - # arrays, where i is however many halfedges are adjacent to - # each edge. - # + """Computes unique edge indices, and a unique map from halfedges to edges, + as well as its inverse. + There is no particular ordering convention for edges. + Boundary edges are guaranteed to be oriented as in F. + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + assume_manifold : bool, optional (default True) + if this is true, will assume that F is manifold, and thus every edge can + be incident to at most two halfedges. + The algorithm is very slow if this is False. + + Returns + ------- + he : (m,3,2) numpy int array + halfedge list as per above conventions + E : (e,2) numpy int array + indices of edges into the vertex array + he_to_E : (m,3) numpy int array + index map from he to corresponding row in E + E_to_he : (e,2,2) numpy int array of list of m (k,2) numpy int arrays + if assume_manifold, (e,2,2) index map from e to corresponding row and + col in he for two halfedges (or -1 if only one halfedge exists); + if not assume_manifold, python list with m entries of (k,2) arrays, + where k is however many halfedges are adjacent to each edge. + + Examples + -------- + TODO + + """ m = F.shape[0] assert m > 0 diff --git a/src/gpytoolbox/halfedge_lengths.py b/src/gpytoolbox/halfedge_lengths.py index 57125aba..705a1986 100644 --- a/src/gpytoolbox/halfedge_lengths.py +++ b/src/gpytoolbox/halfedge_lengths.py @@ -2,17 +2,29 @@ from gpytoolbox.halfedge_lengths_squared import halfedge_lengths_squared def halfedge_lengths(V,F): - # Given a triangle mesh V,F, returns the lengths of all halfedges. - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Inputs: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indices into V - # Outputs: - # l 3*#F lengths of halfedges + """Given a triangle mesh V,F, returns the lengths of all halfedges. + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + l : (m,3) numpy array + lengths of halfedges + + Examples + -------- + TODO + + """ return np.sqrt(halfedge_lengths_squared(V,F)) diff --git a/src/gpytoolbox/halfedge_lengths_squared.py b/src/gpytoolbox/halfedge_lengths_squared.py index 6e79af95..67ef42e3 100644 --- a/src/gpytoolbox/halfedge_lengths_squared.py +++ b/src/gpytoolbox/halfedge_lengths_squared.py @@ -2,22 +2,34 @@ from gpytoolbox.halfedges import halfedges def halfedge_lengths_squared(V,F): - # Given a triangle mesh V,F, returns the lengths of all halfedges, squared. - # The reason to work with squared lengths instead of just lengths is that - # lengths are computed as squared lengths, and then often used as squared - # lengths, and thus keeping track of the squared lengths circumvents a - # lossy square root followed by a square. - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Inputs: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indices into V - # Outputs: - # l_sq 3*#F squared lengths of halfedges + """Given a triangle mesh V,F, returns the lengths of all halfedges, squared. + The reason to work with squared lengths instead of just lengths is that + lengths are computed as squared lengths, and then often used as squared + lengths, and thus keeping track of the squared lengths circumvents a + lossy square root followed by a square. + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + l_sq : (m,3) numpy array + squared lengths of halfedges + + Examples + -------- + TODO + + """ assert F.shape[0] > 0 assert F.shape[1] == 3 diff --git a/src/gpytoolbox/halfedges.py b/src/gpytoolbox/halfedges.py index 1ef4371d..16fc4d40 100644 --- a/src/gpytoolbox/halfedges.py +++ b/src/gpytoolbox/halfedges.py @@ -1,18 +1,29 @@ import numpy as np def halfedges(F): - # Given a triangle mesh with face indices F, returns all oriented halfedges - # as indices into the vertex array. - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Outputs: - # he #F by 3 by 2 halfedge list as per above conventions + """Given a triangle mesh with face indices F, returns all oriented halfedges + as indices into the vertex array. + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + he : (m,3,2) numpy int array + halfedge list as per above conventions + + Examples + -------- + TODO + + """ assert F.shape[0] > 0 assert F.shape[1] == 3 diff --git a/src/gpytoolbox/in_element_aabb.py b/src/gpytoolbox/in_element_aabb.py index 2c87c215..6b848a8f 100644 --- a/src/gpytoolbox/in_element_aabb.py +++ b/src/gpytoolbox/in_element_aabb.py @@ -9,10 +9,10 @@ def in_element_aabb(queries,V,F): ---------- queries : numpy double array Matrix of query point coordinates - V : numpy double array - Matrix of vertex coordinates - F : numpy int array, optional (default None) - Matrix of triangle (2D) or tetrahedron (3D) indices + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) or (m,4) numpy int array + face/tet index list of a triangle/tet mesh Returns ------- diff --git a/src/gpytoolbox/massmatrix.py b/src/gpytoolbox/massmatrix.py index 1c6dd785..56730ffe 100755 --- a/src/gpytoolbox/massmatrix.py +++ b/src/gpytoolbox/massmatrix.py @@ -9,20 +9,23 @@ def massmatrix(V,F=None,type='voronoi'): """FEM Mass matrix - Builds the finite elements mass matrix of a triangle mesh or polyline using a piecewise linear hat function basis. + Builds the finite elements mass matrix of a triangle mesh or polyline using + a piecewise linear hat function basis. Parameters ---------- - V : numpy double array - Matrix of vertex coordinates - F : numpy int array (optional, default None) - Matrix of element indices (if None, assumed to be ordered polyline) + V : (n,d) numpy array + vertex list of a polyline or triangle mesh + F : numpy int array, optional (default: None) + if None or (m,2), interpret input as ordered polyline; + if (m,3) numpy int array, interpred as face index list of a triangle + mesh type : str, optional (default 'voronoi') Type of mass matrix computation: 'voronoi' (default), 'full' or 'barycentric' Returns ------- - M : scipy sparse.csr_matrix + M : (n,n) scipy sparse.csr_matrix Mass matrix See Also diff --git a/src/gpytoolbox/massmatrix_intrinsic.py b/src/gpytoolbox/massmatrix_intrinsic.py index 1d882311..61118050 100755 --- a/src/gpytoolbox/massmatrix_intrinsic.py +++ b/src/gpytoolbox/massmatrix_intrinsic.py @@ -8,22 +8,24 @@ def massmatrix_intrinsic(l_sq,F,n=None,type='voronoi'): """FEM intrinsic mass matrix - Builds the finite elements mass matrix for a triangle mesh using a piecewise linear hat function basis, using only intrinsic information (squared halfedge edge lengths). + Builds the finite elements mass matrix for a triangle mesh using a piecewise + linear hat function basis, using only intrinsic information (squared + halfedge edge lengths). Parameters ---------- - l_sq : numpy double array + l_sq : (m,3) numpy double array Vector of squared halfedge lengths as computed by halfedge_lengths_squared - F : numpy int array - Matrix of triangle indices into some V that is assumed to exist - n : int, optional (default None) + F : (m,3) numpy int array + face index list of a triangle mesh (into a V assumed to exist) + n : int, optional (default: None) Integer denoting the number of vertices in the mesh - type : str, optional (default 'voronoi') + type : str, optional (default: 'voronoi') Type of mass matrix computation: 'voronoi' (default), 'full' or 'barycentric' Returns ------- - M : scipy sparse.csr_matrix + M : (n,n) scipy sparse.csr_matrix Intrinsicly computed mass matrix See Also @@ -38,20 +40,6 @@ def massmatrix_intrinsic(l_sq,F,n=None,type='voronoi'): -------- TO-DO """ - # Builds the finite elements mass matrix for a triangle mesh using a - # piecewise linear hat function basis, using only intrinsic information - # (squared halfedge edge lengths). - # - # Input: - # l_sq #F by 3 numpy array of squared halfedge lengths as computed - # by halfedge_lengths_squared - # F #F by 3 int numpy array of face/edge vertex indices into V - # Optional: - # n an integer denoting the number of vertices in the mesh - # type either of 'voronoi' {default}, 'full', or 'barycentric' - # - # Output: - # M #V by #V mass matrix assert F.shape == l_sq.shape assert F.shape[1]==3 diff --git a/src/gpytoolbox/min_quad_with_fixed.py b/src/gpytoolbox/min_quad_with_fixed.py index eab7a1f7..47a41289 100644 --- a/src/gpytoolbox/min_quad_with_fixed.py +++ b/src/gpytoolbox/min_quad_with_fixed.py @@ -2,29 +2,44 @@ import scipy as sp -def min_quad_with_fixed(Q, c=None, A=None, b=None, k=None, y=None): - # Solve the following quadratic program with linear constraints: - # argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) - # A*u == b - # u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) - # - # Input: - # Q n by n symmetric sparse scipy csr_matrix. - # Will be symmetriyed if not symmetric. - # c scalar or n numpy array or (n,p) numpy array. - # Assumed to be scalar 0 if None - # A m by n sparse scipy csr_matrix. - # m=0 assumed if None. - # b scalar or m numpy array or (m,p) numpy array. - # Assumed to be scalar 0 if None. - # k o numpy int array. - # o=0 assumed if None. - # y scalar or o numpy array or (o,p) numpy array. - # Assumed to be scalar 0 if None. - # - # Output: - # u n numpy array or (n,p) numpy array. - # solution to the optimization problem. +def min_quad_with_fixed(Q, + c=None, + A=None, + b=None, + k=None, + y=None): + """Solve the following quadratic program with linear constraints: + ``` + argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) + A*u == b + u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) + ``` + + Parameters + ---------- + Q : (n,n) symmetric sparse scipy csr_matrix + This matrix will be symmetrized if not exactly symmetric. + c : None or scalar or (n,) numpy array or (n,p) numpy array + Assumed to be scalar 0 if None + A : None or (m,n) sparse scipy csr_matrix + m=0 assumed if None + b : None or scalar or (m,) numpy array or (m,p) numpy array + Assumed to be scalar 0 if None + k : None or (o,) numpy array + o=0 assumed if None + y : None or scalar or (o,) numpy array or (o,p) numpy array + Assumed to be scalar 0 if None + + Returns + ------- + u : (n,) numpy array or (n,p) numpy array + Solution to the optimization problem + + Examples + -------- + TODO + + """ return min_quad_with_fixed_precompute(Q, A, k).solve(c, b, y) @@ -32,26 +47,39 @@ def min_quad_with_fixed(Q, c=None, A=None, b=None, k=None, y=None): # This is written in snake_case on purpose, so constructing the class looks just # like calling a function. class min_quad_with_fixed_precompute: - # Prepare a precomputation object to efficiently solve the following problem: - # argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) - # A*u == b - # u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) - # - # Input: - # Q n by n symmetric sparse scipy csr_matrix. - # Will be symmetriyed if not symmetric. - # A m by n sparse scipy csr_matrix. - # m=0 assumed if None. - # k o numpy int array. - # o=0 assumed if None. - # - # TODO: Detect linearly dependent constraints in A and remove them. - # TODO: Detect constraints in A that contradict k and error. - # TODO: Allow user to specify positive definiteness. - # - # Output: - # precomputed precomputation object that canbe used to solve the problem - def __init__(self, Q, A=None, k=None): + + def __init__(self, + Q, + A=None, + k=None): + """Prepare a precomputation object to efficiently solve the following + constrained optimization problem: + ``` + argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) + A*u == b + u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) + ``` + + Parameters + ---------- + Q : (n,n) symmetric sparse scipy csr_matrix + This matrix will be symmetrized if not exactly symmetric. + A : None or (m,n) sparse scipy csr_matrix + m=0 assumed if None + k : None or (o,) numpy array + o=0 assumed if None + + Returns + ------- + precomputed : instance of class min_quad_with_fixed_precompute + precomputation object that can be used to solve the optimization problem + + Examples + -------- + TODO + + """ + self.n = Q.shape[0] assert Q.shape[1] == self.n assert self.n>0 @@ -128,24 +156,37 @@ def __init__(self, Q, A=None, k=None): self.solver = lambda x : splu.solve(x) + def solve(self, + c=None, + b=None, + y=None): + """Solve the following quadratic program with linear constraints: + ``` + argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) + A*u == b + u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) + ``` + + Parameters + ---------- + c : None or scalar or (n,) numpy array or (n,p) numpy array + Assumed to be scalar 0 if None + b : None or scalar or (m,) numpy array or (m,p) numpy array + Assumed to be scalar 0 if None + y : None or scalar or (o,) numpy array or (o,p) numpy array + Assumed to be scalar 0 if None + + Returns + ------- + u : (n,) numpy array or (n,p) numpy array + Solution to the optimization problem - # Solve the following quadratic program with linear constraints: - # argmin_u 0.5 * tr(u.transpose()*Q*u) + tr(c.transpose()*u) - # A*u == b - # u[k] == y (if y is a 1-tensor) or u[k,:] == y) (if y is a 2-tensor) - # - # Input: - # c scalar or n numpy array or (n,p) numpy array. - # Assumed to be scalar 0 if None - # b scalar or m numpy array or (m,p) numpy array. - # Assumed to be scalar 0 if None. - # y scalar or o numpy array or (o,p) numpy array. - # Assumed to be scalar 0 if None. - # - # Output: - # u n numpy array or (n,p) numpy array. - # solution to the optimization problem. - def solve(self, c=None, b=None, y=None): + Examples + -------- + TODO + + """ + def cp(x): if x is None or np.isscalar(x): return 0 diff --git a/src/gpytoolbox/normalize_points.py b/src/gpytoolbox/normalize_points.py index 7884bb11..78540fae 100644 --- a/src/gpytoolbox/normalize_points.py +++ b/src/gpytoolbox/normalize_points.py @@ -7,7 +7,7 @@ def normalize_points(v,center=None): Parameters ---------- - v : numpy double array + v : (n,d) numpy double array Matrix of point position coordinates center : numpy double array (optional, None by default) Where to center the mesh (if None, centered at zero) diff --git a/src/gpytoolbox/offset_surface.py b/src/gpytoolbox/offset_surface.py index f6b93a61..f1c192d7 100644 --- a/src/gpytoolbox/offset_surface.py +++ b/src/gpytoolbox/offset_surface.py @@ -7,10 +7,10 @@ def offset_surface(V,F,iso,grid_size=100): Parameters ---------- - V : numpy double array - Matrix of vertex coordinates - F : numpy int array - Matrix of triangle indices + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh iso : double Matrix of vertex coordinates of the second mesh grid_size : int, optional (default 100) diff --git a/src/gpytoolbox/per_face_normals.py b/src/gpytoolbox/per_face_normals.py index 7718aa34..13a019b9 100644 --- a/src/gpytoolbox/per_face_normals.py +++ b/src/gpytoolbox/per_face_normals.py @@ -9,16 +9,16 @@ def per_face_normals(V,F,unit_norm=True): Parameters ---------- - V : numpy double array - Matrix of mesh vertex coordinates - F : numpy int array - Matrix of triangle indices into V + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh unit_norm : bool, optional (default True) Whether to normalize each face's normal before outputting Returns ------- - N : numpy double array + N : (n,d) numpy double array Matrix of per-face normals See Also diff --git a/src/gpytoolbox/per_vertex_normals.py b/src/gpytoolbox/per_vertex_normals.py index f12bf85d..70188264 100644 --- a/src/gpytoolbox/per_vertex_normals.py +++ b/src/gpytoolbox/per_vertex_normals.py @@ -7,18 +7,18 @@ def per_vertex_normals(V,F): """Normal vectors to all vertices on a mesh - Computs area-weighted per-vertex unit normal vectors for a triangle mesh + Computes area-weighted per-vertex unit normal vectors for a triangle mesh Parameters ---------- - V : numpy double array - Matrix of mesh vertex coordinates - F : numpy int array - Matrix of triangle indices into V + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh Returns ------- - N : numpy double array + N : (m,d) numpy double array Matrix of per-vertex normals See Also @@ -29,16 +29,6 @@ def per_vertex_normals(V,F): -------- TODO """ - # Computs per vertex unit normal vectors for a triangle mesh - # - # Input: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indeces into V - # - # Output: - # N #F by 3 numpy array of per-vertex unit normals - # - # First compute face normals face_normals = per_face_normals(V,F,unit_norm=True) diff --git a/src/gpytoolbox/ray_mesh_intersect.py b/src/gpytoolbox/ray_mesh_intersect.py index 51a1274c..c8996785 100644 --- a/src/gpytoolbox/ray_mesh_intersect.py +++ b/src/gpytoolbox/ray_mesh_intersect.py @@ -11,10 +11,10 @@ def ray_mesh_intersect(cam_pos,cam_dir,V,F): Matrix of camera positions cam_dir : numpy double array Matrix of camera directions - V : numpy double array - Matrix of mesh vertices - F : numpy int array - Matrix of triangle indices into V + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh Returns ------- diff --git a/src/gpytoolbox/read_mesh.py b/src/gpytoolbox/read_mesh.py index 1ae626e3..f920d7fb 100644 --- a/src/gpytoolbox/read_mesh.py +++ b/src/gpytoolbox/read_mesh.py @@ -8,34 +8,51 @@ def read_mesh(file, return_UV=False, return_N=False, reader=None): - # Reads a mesh from file. - # - # If you have the approproate C++ extensions installed, this will use a fast - # C++-based reader. If you do not, this will use a slow python reader. - # - # Currently only supports triangle meshes. - # - # Input: - # file the path to the mesh - # Optional: - # fmt The file format of the mesh to open. - # If None, try to guess the format from the file - # extension. - # Supported formats: obj - # return_UV Try reading texture coordinates, if they are - # present and the file format supports it. - # return_N Try reading normal coordinates, if they are - # present and the file format supports it. - # reader Which reader engine to use. None, 'C++' or 'Python'. - # If None, will use C++ if available, and else Python. - # Output: - # V numpy array of mesh vertex positions - # F numpy array of mesh face indices into V - # UV numpy array of texture coordinates (if requested) - # Ft numpy array of mesh face indices into UV (if requested) - # N numpy array of normal coordinates (if requested) - # Fn numpy array of mesh face indices into N (if requested) - # + """Reads a mesh from file. + + If you have the approproate C++ extensions installed, this will use a fast + C++-based reader. If you do not, this will use a slow python reader. + + Currently only supports triangle meshes. + + Parameters + ------- + file : string + the path the mesh will be read from + fmt : string, optional (default: None) + The file format of the mesh to open. + If None, try to guess the format from the file extension. + Supported formats: obj + return_UV : bool, optional (default: None) + Try reading texture coordinates, if they are present and the file + format supports it. + return_N : bool, optional (default: None) + Try reading normal coordinates, if they are present and the file format + supports it. + reader : string, optional (default: None) + Which reader engine to use. None, 'C++' or 'Python'. + If None, will use C++ if available, and else Python. + + Returns + ---------- + V : (n,3) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh (into V) + UV : (n_uv,2) numpy array, if requested + vertex list for texture coordinates + Ft : (m,3) numpy int array, if requested + face index list for texture coordinates (into UV) + N : (n_n,3) numpy array, if requested + vertex list for normal coordinates + Fn : (m,3) numpy int array, if requested + face index list for normal coordinates (into N) + + Examples + -------- + TODO + + """ # Detect format if it has not been specified if fmt is None: diff --git a/src/gpytoolbox/regular_cube_mesh.py b/src/gpytoolbox/regular_cube_mesh.py index 796e8616..83665385 100644 --- a/src/gpytoolbox/regular_cube_mesh.py +++ b/src/gpytoolbox/regular_cube_mesh.py @@ -14,10 +14,10 @@ def regular_cube_mesh(gs,type='rotationally-symmetric'): Returns ------- - V : numpy double array - Matrix of tet mesh vertex coordinates - T : numpy int array - Matrix of tet vertex indices into V + V : (n,d) numpy array + vertex list of a tet mesh + T : (m,T) numpy int array + tet index list of a tet mesh See Also -------- @@ -27,19 +27,6 @@ def regular_cube_mesh(gs,type='rotationally-symmetric'): -------- TODO """ - # Generates a regular tetrahedral mesh of a one by one by one cube. - # Each grid cube is decomposed into 6 reflectionally symmetric tetrahedra - # - # Input: - # gs int number of vertices on each side - # Optional: - # type a string choosing the specific cube division scheme: 'five' for - # a division of each cube into 5 tets 'rotationally-symmetric' - # (default), 'reflectionally-symmetric' or 'hex' - # Output: - # V #V by 3 numpy array of mesh vertex positions - # T #T by 4 int numpy array of tet vertex indeces into V - # dictionary ={ 'five' : 0, diff --git a/src/gpytoolbox/subdivide.py b/src/gpytoolbox/subdivide.py index cce30787..42982999 100644 --- a/src/gpytoolbox/subdivide.py +++ b/src/gpytoolbox/subdivide.py @@ -6,27 +6,39 @@ def subdivide(V,F, method='upsample', iters=1, return_matrix=False): - # Subdivides a mesh - # - # Input: - # V #V by d numpy array of mesh vertex positions. - # Can be None, in which case the subdivision is performed only - # logically on the faces. - # F #F by k int numpy array of vertex indices into V. - # Can represent a polyline or a triangle mesh. - # Optional: - # method which method to use for subdivison. - # can be 'upsample' {default}, - # 'loop' (only triangle meshes) - # iters how many subdivision iterations to perform - # return_matrix whether to return matrix for sparse map S - # - # Output: - # Vu #Vu by d numpy array of subdivided mesh vertex positions - # Is None if V was None. - # Fu #Fu by k int numpy array of subdivided vertex indices into V. - # S If return_matrix is True, the sparse matrix such that Vu == S*V - # + """Builds the (pos. def.) cotangent Laplacian for a triangle mesh. + + Parameters + ---------- + V : (n,d) numpy array + vertex list of vertex positions + F : numpy int array + if (m,2), interpret input as ordered polyline; + if (m,3) numpy int array, interpred as face index list of a triangle + mesh + method : string, optional (default: 'upsample') + Which method to use for subdivison. + Can be 'upsample' {default}, 'loop' (only triangle meshes) + iters : int, optional (default: 1) + How many iterations of subdivision to perform. + return_matrix : bool, optional (default: False) + Whether to return the matrix for the sparse map S. + + Returns + ------- + Vu : (n_u,d) numpy array + vertex list of subdivided polyline / mesh + Fu : (m_u,2) or (m_u) numpy int array + face index list of upsampled polyline or triangle mesh + S : (n_u,n) sparse scipy csr_matrix (if requested) + sparse matrix such that `Vu == S*V`; + returned only if `return_matrix == True`. + + Examples + -------- + TODO + + """ assert iters >= 0 diff --git a/src/gpytoolbox/tip_angles.py b/src/gpytoolbox/tip_angles.py index 85b52520..5885a677 100644 --- a/src/gpytoolbox/tip_angles.py +++ b/src/gpytoolbox/tip_angles.py @@ -2,20 +2,30 @@ from gpytoolbox.halfedge_lengths_squared import halfedge_lengths_squared from gpytoolbox.tip_angles_intrinsic import tip_angles_intrinsic -def tip_angles(V, F, use_small_angle_approx=True): - # Computes the angles formed by each vertex within its respective face - # (the tip angle). - # - # Input: - # V #V by 3 numpy array of mesh vertex positions - # F #F by 3 int numpy array of face/edge vertex indeces into V - # Optional: - # use_small_angle_approx if True, uses a different, more - # more stable formula for small - # angles. - # - # Output: - # alpha #F by 3 numpy array of tip angles for each vertex +def tip_angles(V, F, + use_small_angle_approx=True): + """Computes the angles formed by each vertex within its respective face + (the tip angle). + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + use_small_angle_approx : bool, optional (default: True) + If True, uses a different, more more stable formula for small angles. + + Returns + ------- + α : (m,3) numpy array + tip angles for each vertex referenced in `F` + + Examples + -------- + TODO + + """ l_sq = halfedge_lengths_squared(V,F) return tip_angles_intrinsic(l_sq,F,use_small_angle_approx) diff --git a/src/gpytoolbox/tip_angles_intrinsic.py b/src/gpytoolbox/tip_angles_intrinsic.py index fcd19eb1..1cb1f712 100644 --- a/src/gpytoolbox/tip_angles_intrinsic.py +++ b/src/gpytoolbox/tip_angles_intrinsic.py @@ -1,21 +1,31 @@ import numpy as np -def tip_angles_intrinsic(l_sq, F, use_small_angle_approx=True): - # Computes the angles formed by each vertex within its respective face - # (the tip angle) using only intrinsic information (squared halfedge edge - # lengths). - # - # Input: - # l_sq #F by 3 numpy array of squared halfedge lengths as computed - # by halfedge_lengths_squared - # F #F by 3 int numpy array of face/edge vertex indeces into V - # Optional: - # use_small_angle_approx if True, uses a different, more - # more stable formula for small - # angles. - # - # Output: - # alpha #F by 3 numpy array of tip angles for each vertex +def tip_angles_intrinsic(l_sq, F, + use_small_angle_approx=True): + """Computes the angles formed by each vertex within its respective face + (the tip angle) using only intrinsic information (squared halfedge edge + lengths). + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + use_small_angle_approx : bool, optional (default: True) + If True, uses a different, more more stable formula for small angles. + + Returns + ------- + α : (m,3) numpy array + tip angles for each vertex referenced in `F` + + Examples + -------- + TODO + + """ + assert F.shape[1] == 3 assert l_sq.shape == F.shape diff --git a/src/gpytoolbox/triangle_triangle_adjacency.py b/src/gpytoolbox/triangle_triangle_adjacency.py index f1400976..387f25cd 100644 --- a/src/gpytoolbox/triangle_triangle_adjacency.py +++ b/src/gpytoolbox/triangle_triangle_adjacency.py @@ -3,23 +3,34 @@ from gpytoolbox.array_correspondence import array_correspondence def triangle_triangle_adjacency(F): - # Given a manifold triangle mesh with face indices F, computes adjacency - # info between triangles - # - # The ordering convention for halfedges is the following: - # [halfedge opposite vertex 0, - # halfedge opposite vertex 1, - # halfedge opposite vertex 2] - # - # Inputs: - # F #F by 3 face index list of a triangle mesh - # Outputs: - # TT #F by 3 index list specifying which face j is adjacent to which - # face i across the respective halfedge in position (i,j). - # If there is no adjacent face (boundary halfedge), the entry is - # -1 - # TTi #F by 3 index list specifying which halfedge of face j (0,1,2) - # is adjacent to i in position (i,j) + """Given a manifold triangle mesh with face indices F, computes adjacency + info between triangles + + The ordering convention for halfedges is the following: + [halfedge opposite vertex 0, + halfedge opposite vertex 1, + halfedge opposite vertex 2] + + Parameters + ---------- + F : (m,3) numpy int array + face index list of a triangle mesh + + Returns + ------- + TT : (m,3) numpy int array + Index list specifying which face j is adjacent to which face i across + the respective halfedge in position (i,j). + If there is no adjacent face (boundary halfedge), the entry is -1. + TTi : (m,3) numpy int array + Index list specifying which halfedge of face j (0,1,2) is adjacent to i + in position (i,j). + + Examples + -------- + TODO + + """ m = F.shape[0] assert m > 0 diff --git a/src/gpytoolbox/write_mesh.py b/src/gpytoolbox/write_mesh.py index a618f30a..7e3d47d1 100644 --- a/src/gpytoolbox/write_mesh.py +++ b/src/gpytoolbox/write_mesh.py @@ -11,29 +11,46 @@ def write_mesh(file, Fn=None, fmt=None, writer=None): - # Writes a mesh to a file. - # - # If you have the approproate C++ extensions installed, this will use a fast - # C++-based writer. If you do not, this will use a slow python writer. - # - # Currently only supports triangle meshes. - # - # Input: - # file the path to the mesh - # V numpy array of mesh vertex positions - # F numpy array of mesh face indices into V (0 indexing) - # Optional: - # UV numpy array of texture coordinates - # Ft numpy array of mesh face indices into UV - # N numpy array of normal coordinates - # Fn numpy array of mesh face indices into N - # fmt The file format of the mesh to open. - # If None, try to guess the format from the file - # extension. - # Supported formats: obj - # writer Which writer engine to use. None, 'C++' or 'Python'. - # If None, will use C++ if available, and else Python. - # + """Writes a mesh to a file. + + If you have the approproate C++ extensions installed, this will use a fast + C++-based writer. If you do not, this will use a slow python writer. + + Currently only supports triangle meshes. + + Parameters + ---------- + file : string + the path the mesh will be written to + V : (n,3) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh (into V) + UV : (n_uv,2) numpy array, optional (default: None) + vertex list for texture coordinates + Ft : (m,3) numpy int array, optional (default: None) + face index list for texture coordinates (into UV) + N : (n_n,3) numpy array, optional (default: None) + vertex list for normal coordinates + Fn : (m,3) numpy int array, optional (default: None) + face index list for normal coordinates (into N) + fmt : string, optional (default: None) + The file format of the mesh to write. + If None, try to guess the format from the file extension. + Supported formats: obj + writer : string, optional (default: None) + Which writer engine to use. None, 'C++' or 'Python'. + If None, will use C++ if available, and else Python. + + Returns + ------- + + + Examples + -------- + TODO + + """ # Detect format if it has not been specified if fmt is None: diff --git a/test/test_quadtree_laplacian.py b/test/test_quadtree_laplacian.py index 4972f490..a10695f9 100644 --- a/test/test_quadtree_laplacian.py +++ b/test/test_quadtree_laplacian.py @@ -8,71 +8,76 @@ class TestQuadtreeLaplacian(unittest.TestCase): # TODO WRITE WITHOUT IGL def test_laplacian_convergence(self): self.assertTrue(True) - # np.random.seed(0) - # th = 2*np.pi*np.random.rand(100,1) - # P = 0.5*np.concatenate((np.cos(th),np.sin(th)),axis=1) + np.random.seed(0) + th = 2*np.pi*np.random.rand(100,1) + P = 0.5*np.concatenate((np.cos(th),np.sin(th)),axis=1) - # C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=6,min_depth=2,vmin=np.array([-1,-1]),vmax=np.array([1,1])) - # L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) - # fun = stored_at[:,0]**2.0 - # # This will never be exactly two everywhere, but it should at least be two in most places. - # self.assertTrue(np.isclose(np.median(np.abs(L @ fun)),2.0)) + C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=6,min_depth=2,vmin=np.array([-1,-1]),vmax=np.array([1,1])) + L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) + fun = stored_at[:,0]**2.0 + # This will never be exactly two everywhere, but it should at least be two in most places. + self.assertTrue(np.isclose(np.median(np.abs(L @ fun)),2.0)) - # # This is not very satisfying so just to be safe let's solve a Poisson equation - # th = 2*np.pi*np.random.rand(500,1) - # P = 0.5*np.concatenate((np.cos(th),np.sin(th)),axis=1) - # C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=5,vmin=np.array([-1,-1]),vmax=np.array([1,1])) - # V,Q,_ = gpytoolbox.bad_quad_mesh_from_quadtree(C,W,CH) - # L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) - # gt_fun = stored_at[:,0]**2.0 - # lap_fun = 2.0 + 0.0*stored_at[:,0] - # bb = ((stored_at[:,0]>0.85) | (stored_at[:,0]<-0.85) | (stored_at[:,1]>0.85) | (stored_at[:,1]<-0.85)).nonzero()[0] - # bc = gt_fun[bb] - # Aeq = csr_matrix((0, 0), dtype=np.float64) - # Beq = np.array([], dtype=np.float64) + # This is not very satisfying so just to be safe let's solve a Poisson equation + th = 2*np.pi*np.random.rand(500,1) + P = 0.5*np.concatenate((np.cos(th),np.sin(th)),axis=1) + C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=5,vmin=np.array([-1,-1]),vmax=np.array([1,1])) + V,Q,_ = gpytoolbox.bad_quad_mesh_from_quadtree(C,W,CH) + L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) + gt_fun = stored_at[:,0]**2.0 + lap_fun = 2.0 + 0.0*stored_at[:,0] + bb = ((stored_at[:,0]>0.85) | (stored_at[:,0]<-0.85) | (stored_at[:,1]>0.85) | (stored_at[:,1]<-0.85)).nonzero()[0] + bc = gt_fun[bb] + Aeq = csr_matrix((0, 0), dtype=np.float64) + Beq = np.array([], dtype=np.float64) # u = igl.min_quad_with_fixed(L,-1.0*lap_fun,bb,bc,Aeq,Beq,False)[1] + # b = -1.0*lap_fun + # print(b.shape) + u = gpytoolbox.fixed_dof_solve(L, b=1.0*lap_fun, k=bb, y=bc) + # Plot solution: + # ps.init() + # quadtree = ps.register_surface_mesh("test quadtree",V,Q,edge_width=1) + # quadtree.add_scalar_quantity("numeric solve",u,defined_on='faces') + # quadtree.add_scalar_quantity("groundtruth",gt_fun,defined_on='faces') + # ps.register_point_cloud("boundary conditions",stored_at[bb,:]) + # ps.show() - # # Plot solution: - # # ps.init() - # # quadtree = ps.register_surface_mesh("test quadtree",V,Q,edge_width=1) - # # quadtree.add_scalar_quantity("numeric solve",u,defined_on='faces') - # # quadtree.add_scalar_quantity("groundtruth",gt_fun,defined_on='faces') - # # ps.register_point_cloud("boundary conditions",stored_at[bb,:]) - # # ps.show() - - # # Error will exist, since the lowest deapth leaf nodes are big, but it should be "reasonable" - # self.assertTrue(np.max(np.abs(u-gt_fun))<0.05) - # # It should also go down if we increase the minimum depth - # C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=6,vmin=np.array([-1,-1]),vmax=np.array([1,1])) - # L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) - # gt_fun = stored_at[:,0]**2.0 - # lap_fun = 2.0 + 0.0*stored_at[:,0] - # bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] - # bc = gt_fun[bb] + # Error will exist, since the lowest deapth leaf nodes are big, but it should be "reasonable" + self.assertTrue(np.max(np.abs(u-gt_fun))<0.05) + # It should also go down if we increase the minimum depth + C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=6,vmin=np.array([-1,-1]),vmax=np.array([1,1])) + L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) + gt_fun = stored_at[:,0]**2.0 + lap_fun = 2.0 + 0.0*stored_at[:,0] + bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] + bc = gt_fun[bb] + u = gpytoolbox.fixed_dof_solve(L, b=1.0*lap_fun, k=bb, y=bc) # u = igl.min_quad_with_fixed(L,-1.0*lap_fun,bb,bc,Aeq,Beq,False)[1] - # self.assertTrue(np.max(np.abs(u-gt_fun))<0.01) - # C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=7,vmin=np.array([-1,-1]),vmax=np.array([1,1])) - # L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) - # gt_fun = stored_at[:,0]**2.0 - # lap_fun = 2.0 + 0.0*stored_at[:,0] - # bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] - # bc = gt_fun[bb] + self.assertTrue(np.max(np.abs(u-gt_fun))<0.01) + C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=9,min_depth=7,vmin=np.array([-1,-1]),vmax=np.array([1,1])) + L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) + gt_fun = stored_at[:,0]**2.0 + lap_fun = 2.0 + 0.0*stored_at[:,0] + bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] + bc = gt_fun[bb] + u = gpytoolbox.fixed_dof_solve(L, b=1.0*lap_fun, k=bb, y=bc) # u = igl.min_quad_with_fixed(L,-1.0*lap_fun,bb,bc,Aeq,Beq,False)[1] - # self.assertTrue(np.max(np.abs(u-gt_fun))<0.005) - # # *shrug*. I don't really know how to evaluate this Laplacian, especially given some error is always expected since the stencil is not exact... + self.assertTrue(np.max(np.abs(u-gt_fun))<0.005) + # *shrug*. I don't really know how to evaluate this Laplacian, especially given some error is always expected since the stencil is not exact... - # # Does it *at least* converge if we force a regular grid? This will catch factors of two and minus sign errors: - # C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=6,min_depth=6,vmin=np.array([-1,-1]),vmax=np.array([1,1])) - # L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) - # gt_fun = stored_at[:,0]**2.0 - # lap_fun = 2.0 + 0.0*stored_at[:,0] - # bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] - # bc = gt_fun[bb] + # Does it *at least* converge if we force a regular grid? This will catch factors of two and minus sign errors: + C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=True,max_depth=6,min_depth=6,vmin=np.array([-1,-1]),vmax=np.array([1,1])) + L, stored_at = gpytoolbox.quadtree_laplacian(C,W,CH,D,A) + gt_fun = stored_at[:,0]**2.0 + lap_fun = 2.0 + 0.0*stored_at[:,0] + bb = ((stored_at[:,0]>0.8) | (stored_at[:,0]<-0.8) | (stored_at[:,1]>0.8) | (stored_at[:,1]<-0.8)).nonzero()[0] + bc = gt_fun[bb] # u = igl.min_quad_with_fixed(L,-1.0*lap_fun,bb,bc,Aeq,Beq,False)[1] - # self.assertTrue(np.isclose(np.max(np.abs(u-gt_fun)),0.0)) + u = gpytoolbox.fixed_dof_solve(L, b=1.0*lap_fun, k=bb, y=bc) + self.assertTrue(np.isclose(np.max(np.abs(u-gt_fun)),0.0)) if __name__ == '__main__': unittest.main() \ No newline at end of file