Van Le – SGI 2023

Fellows: Matheus Araujo, Aditya Abhyankar, Sahana Kargi, Van Le, Maria Stuebner

Mentor: Martin Skrodzki

Introduction

Hyperbolic geometry has recently gained popularity as a promising model to represent non-Euclidean data, with applications ranging from natural language processing to computer vision. It also can be useful in geometry processing to verify discrete conformal equivalence of different triangulations of a surface. One important advantage of hyperbolic geometry is its capacity to represent and visualize hierarchical structures in a more natural way than Euclidean geometry. However, data visualization algorithms, such as t-distributed stochastic neighbor embedding (t-SNE), were more extensively explored in the Euclidean space than in the hyperbolic space. Specifically, spatial data structures can be leveraged to accelerate the approximation of nearest neighbors and provide significant improvements in performance. In this project, mentored by Professor Martin Skrodzi, we investigate different implementations and strategies of spatial acceleration data structures adapted for hyperbolic geometry, such as polar quadtrees, to improve the performance of nearest neighbors computation in the hyperbolic space.

Hyperbolic Geometry

Hyperbolic geometry accepts four out of five of Euclid’s postulates of geometry. An interpretation of Euclid’s fifth postulate is known as the “Playfair’s Axiom,” which states that when given a line and a point not on it, at most one line parallel to the given line can be drawn through the point. Hyperbolic geometry accepts an altered version of this axiom, which allows for an infinite number of lines parallel to the given line and through the given point. There’s no easy way to interpret this geometry in Euclidean space. Every point has constant negative Gaussian curvature equal to \(-1\). One visualization attempt is the pseudosphere, which is the shape obtained by revolving the tractrix curve about its asymptote. But this shape has singularities. Various attempts use a type of knitting technique called crocheting which results in fractal-likes shapes similar to coral reefs, but these only have finite patches with constant Gaussian curvature. Locally, we can intuit the hyperbolic plane to be a saddle surface like a Pringle, but this is just a very local approximation. There are many other models which allow us to make sense of hyperbolic geometry using projections. The one relevant to this project is the Poincaré disk model, which is obtained by perspectively projecting the hyperboloid of one sheet onto the unit disk in \(\mathbb{R}^2\).

t-SNE Algorithm

An interesting application of hyperbolic space is hierarchical data visualization, such as single-cell clustering or that of the semantic meaning of an English word vocabulary. In the Poincaré disk model, space grows exponentially as we stray further from the origin, so we can summarize a lot more data than the polynomially area-growing Euclidean space. Another intersecting property of the Poincaré disk model of hyperbolic space, as we stray further from the origin, is that the hyperbolic distance \(d(x,y)\) between points \(x,y\) approaches \(d(x,0)+d(y,0)\). This is extremely useful for visualizing tree-like data, where the closest path between two child nodes is through their parent. One general challenge in data visualization is that each point lies in very high dimensional space. A simple projection down to a lower dimension ruins clustering and spatial relationships between the data. The t-SNE algorithm is an iterative gradient-based procedure that recovers these relationships after an initial projection like that obtained via running a principal component analysis (PCA). In each iteration, the t-SNE algorithm computes an update on a per-data-point basis as a sum of “attractive forces” pulling towards spatially close points and a sum of “repulsive forces” pushing away from spatial distant points. For each data point, the algorithm computes a similarity score for every other data point in both the high and low dimensional spaces and computes the attractive and repulsive forces from these scores. If the low dimensional space is the 2D Poincaré disk, the distance used for computing the similarity scores is the hyperbolic distance:

\(d(\mathbf{u}, \mathbf{v}) = \cosh^{-1} \left(1 + 2 \cdot \frac{|| \mathbf{u} – \mathbf{v} ||^{2}}{(1 – || \mathbf{u} ||^{2})(1 – || \mathbf{v} ||^{2})} \right)\).

Quadtrees

The main bottleneck of the t-SNE algorithm happens to be the computation of repulsive forces. This normally requires an \(O(N^2)\) time complexity per iteration, which can be heavily optimized by using spatial data structures like a quadtree. Each node represents a polar rectangle on the Poincare disk (with a minimum and maximum radius and angle) and is either a leaf or has four sub-rectangles as children. Building the tree takes roughly linear time, as we add all the data points into the tree one by one, incrementally adding more and more nodes as needed, so that just one (or extremely few) points belong to each leaf node. We experimented with various splitting techniques for the generation of the four children nodes from a parent: (1) Area-based splitting, which involves computing the midpoint of the parent rectangle via

\(\text{mid}_r = \cosh^{-1}\left(\frac{\cosh(\alpha\ \text{max}_{r}) + \cosh(\alpha\ \text{min}_r)}{2}\right)\cdot\frac1\alpha\),

and

\(\text{mid}_\theta = \frac{\max_\theta + \min_\theta}{2}\),

which is the hyperbolic average of the min and max radii and a simple average of the min and max angles, and (2) Length-based splitting, in which the midpoint is computed as

\(\text{mid}_r = \frac{\min_r + \max_r}{2}\),

and

\(\text{mid}_\theta = \frac{\min_\theta + \max_\theta}{2}\),

which is just half the radius and angle of the parent. We perform the four-way split using this midpoint and the four corners of the parent rectangle. Once the tree is built, at each iterative gradient update, for each point, we approximate the repulsive forces by traversing the tree via depth-first-search, and if the Einstein midpoint,

\(m_i(\{\alpha_{ij}\}_j, \{\mathbf{v}{ij}\}_j) = \sum_j \left[\frac{\alpha{ij}\gamma(\mathbf{v}{ij})}{\sum_l \alpha{il}\gamma(\mathbf{v}{il})} \right]\mathbf{v}{ij}\),

of the points belonging to the rectangle at a given level in the search is sufficiently far away from the point, we simply use the distance to the Einstein midpoint itself in the calculation of the repulsive similarity score. Intuitively, this is a good approximation, as points far away do not have a large impact on the update anyway, so grouping together many far-away points estimate their collective impact on the update in one fell swoop.

Results

Conclusion and Future Work

Acknowledgments

We would like to express our sincere gratitude to Dr. Martin Skrodzki for his guidance and encouragement throughout this project.

Link to Repo

References

[1] Mark Gillespie, Boris Springborn, and Keenan Crane. “Discrete Conformal Equivalence of Polyhedral Surfaces”. In: ACM Trans. Graph. 40.4 (2021).

[2] Yunhui Guo et al. “Clipped hyperbolic classifiers are super-hyperbolic classifiers”. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2022, pp. 11–20.

[3] Shimizu Ryohei, Mukuta YUSUKE, and Harada Tatsuya. “Hyperbolic neural networks++”. In: International Conference on Learning Representations. Vol. 3. 2020.

Fellows: Ikaro Penha Costa, Bereket Faltamo, Sahana Kargi, Van Le

Advisor: Oded Stein, University of Southern California

Volunteers: Lucas Valenca, Shaimaa Monem

MOTIVATION

The art of assembling model kits has always been a fascinating hobby; however, the process of generating these models is often limited to surfaces that can be easily flattened. Traditionally, it takes the 3D artists weeks or even months to create a model with detailed instructions on how to assemble it. Thus, this research aims to explore the potential of developing a tool that can automate the process of capturing the curvature of more complex objects, allowing the construction of any given surface. Our ambitions are to offer professionals an unprecedented level of creative freedom by opening up the possibilities of replicating more intricate shapes and from there, pave the way for novel applications in the engineering world.

METHODOLOGY

This project employs C++ and heavily uses the geometry processing library – libigl. For this research, the example objects are an icosahedron approximating the sphere and a mesh to a boat.

Segmentation into developable shapes

The first step is to determine whether it is a developable surface or not. In mathematics, a developable surface is a smooth surface with zero Gaussian curvature. We first took the boat and the icosahedral sphere and segmented them into parts that have zero Gaussian curvature. We checked the Gaussian curvature by using the Matlab code from Day 1 of the Tutorial Week. We could also do it with C++ by using this function.

Flattening pieces via corresponding homeomorphism

Once the shape is segmented, the next step is to flatten those segmented pieces. For this, we will use the geometry-processing-parameterization library and the Least Squares conformal mapping (LSCM) method. It works by finding and preserving the local angles and lengths between pieces. It aligns the centroids within the segments and then finds the optimal rotation, scaling, and translation using a least squares approach. Here are the results of implementing it on different segments of the boat.

Before LSCM After LSCM

In addition to flattening the pieces, we also colored them according to their coordinates. We normalized the X, Y, and Z coordinates of all the vertices and colored them the coordinating RGB color. The goal is to use this colormap during reassembly so that the flattened pieces can be curved to the proper amounts.

Reassembling Pieces

In possession of the set of folded pieces, it naturally follows that each piece could be reconnected to revert this collection to the original object shape. This step is denominated reassembling instructions and essentially requires techniques of Shape Correspondence. That is, find a correspondence between a folded piece and the original shape.

Two main mathematical measures are involved in comparing the similarity of topological spaces: Hausdorff and Frechét distances. A more common implementation of bounds for Hausdorff distance is based on an energy optimization as discussed by Alec Jacobson in these tutorial notes. In this case, the Hausdorff distance is used to iteratively find the closest correspondence between a point in the source mesh and the entirety of the target mesh. This algorithm is Iterative Closest Point (ICP) and it is also available in libigl.

As seen in the image above, a segmented piece matches the correct place on the icosahedron after running the ICP algorithm. The ICP technique is, in summary, a sampling step of the mesh (VX, FX) and an optimization step in which each point in VX the objective function is to minimize the Hausdorff distance between this point and the target mesh (VY, FY). The libigl provides the point-to-normal implementation discussed in Jacobson’s notes, which is a form of Gauss-Newton optimization algorithm. This method requires an integer for the number of sample points in (VX, FX) and also an integer for the maximum number of iterations to be used in the minimization algorithm.

The number of sample points and maximum iterations is decisive in order to obtain reliable reassembling results. As a common fact in optimization, depending on the number of iterations and size of the model, the algorithm might be susceptible to being stuck at a local minimum. This is also the case with ICP, causing artifacts as the image below:

A valid approach to avoid this sinister reassembling is to generate different initial conditions so that the optimization algorithm might converge to different local minima. In this case, using the Eigen library, a set of uniformly random rigid transformations are applied to the piece before passing it to the ICP method. Since the Hausdorff distance is the objective function, the rigid transformation which causes ICP to produce the lowest Hausdorff distance is stored and then used to produce the final result.

However, research (and also life) is never linear. The approach described before would cause a few failed assertions throughout the label library. The assertions would be informed that the transformed mesh would have zero area. As the libigl area computation is robust, further investigations are necessary.

The investigations have determined the existence of overflow in the translation vector inside the ICP method. For the boat example, the entries in VX are originally in the magnitude of 10¹, and after some iterations of the ICP, the translation vector is applied to VX whose every entry would be numerically the same number in the magnitude of 10²⁸; then, causing the transformed mesh to have zero area. In this case, the proposed remediation is to ignore the transformed initial pieces which would cause large entry-wise points in VX.

Hence, after performing each ICP iteration with the extra care of managing the translation vector entries, promising results were found as shown in the image below. Observe that the piece matches a similar side of the boat but not the correct side, it is placed in the mirrored side. In this case, a possibility is that it might be necessary to add more constraints to the initial conditions.

As finding the cause of zero areas has demanded significant time, the approach of adding more constraints to the initial conditions can be left as future work. For this, one can store information on the piece boundary in the segmentation step. This information could be used to point the ICP rotation and translation in the direction of the correct final place in the target shape at each iteration step. That is, try to point ICP rigid transformation to the original place of segmentation of each piece. Hopefully, this would cause correctness of results and faster convergence.

Reassembling Instructions with Edge Coloring

When putting the pieces back together we need a way to know which parts come together. In particular, we need to know which pieces have shared edges. That is where edge coloring helps. We color-coded each cutting line with random patterns so that the pieces will take part of the pattern with them. The color coding of the icosahedron can be seen below.

Real Life Assembly by Adding Tabs and Holes

As a final step, we wanted to see how accurate our algorithm was in real life. We formatted and printed the images and assembled them. Some work is needed in ensuring that all the individual shapes are the proper size. Additionally, we added the tabs and holes to the segments so the pieces fit together.

Link to Repo

ACKNOWLEDGEMENTS

We would like to express our sincere gratitude to Professor Oded Stein for his instruction, guidance, and encouragement throughout this research.

We would also like to acknowledge the advice and recommendations of the volunteers – Lucas and Shaimaa for offering help, unique insights, and perspectives.

This research heavily used C++ and the library libigl.