Note: Although this is an ongoing work, this report documents our progress between the official 2 weeks of the project. (August 2, 2021 – August 13, 2021)
The past 2 weeks at SGI, we have been working with David Palmer on investigating a novel Bayesian approach towards the angular synchronization problem. This blog post is written to document our work and share a sneak peek into our research.
Introduction
Consider a set of unknown absolute orientations \(\{q_1, q_2, \ldots, q_n\}\) with respect to some fixed basis. The problem of angular synchronization deals with the accurate estimation of these orientations from noisy observations of their relative offsets \(O_{i, j}\), up to a constant additive phase. We are particularly interested in estimating these “true” orientations in the presence of many outlier measurements. Our interest in this topic stems from the fact that the angular synchronization problem arises in various avenues of science, including reconstruction problems in computer vision, ranking problems, sensor network localization, structural biology, etc. In our work, we study this problem from a Bayesian perspective, by modelling the observed data as a mixture between noisy observations and outliers. We also investigate the problem of continuous label switching, a global ambiguity that arises from the lack of knowledge about the basis of the absolute orientations \(q_i\). Finally, we experiment on a novel Riemannian gradient descent method for alleviating this continuous label switching problem and provide our observations herein.
Brief Primer on Bayesian Inference
Before going deeper, we’ll briefly discuss Bayesian inference. At the heart of Bayesian inference lies the celebrated Bayes’ rule (read \(a|b\) as “a given b”):
In our problem, \(q\) and \(O\) denote the true orientations that we are estimating and the noisy observations with outliers respectively. We are interested in finding the posterior distribution (or at least samples from it) over the ground truth \(q\) given the noisy observations \(O\).
The denominator (i.e., the evidence or partition function) is an integral over all \(q\)s. Exactly evaluating this integral is often intractable if \(q\) lies on some continuous manifold, as in our problem. This makes drawing samples from the posterior becomes hard.
One way to avoid computing the partition function is a sampling method called Markov Chain Monte Carlo (MCMC). Intuitively, the posterior is approximated by a markov chain whose transitions can be computed using a simpler distribution called the proposal distribution. Successive samples are then accepted or rejected based on an acceptance probability designed to ensure convergence to the posterior distribution in the limit of infinite samples. Simply put, after enough samples are drawn using MCMC, they will look like the samples from the posterior\(P(q|O) \propto P(O|q) \cdot P(q)\) without requiring us to calculate the intractable normalization \(P(O) = \int\limits_q P(O|q)\cdot P(q)\). In our work, we use Hamiltonian Monte Carlo (HMC), an efficient variant of MCMC, which uses Hamiltonian Dynamics to propose the next sample. From an implementation perspective, we use the built-in HMC sampler in Stan for drawing samples.
Mixture Model
As mentioned before, we model the noisy observation as a mixture model of true distribution and outliers. This is denoted by (Equation 1):
where \(\eta_{i j}\) is the additive noise to our true observation, \(q_i q_j^T\) is the relative orientation between the \(i^\textrm{th}\) and \(j^\textrm{th}\) objects, and \(\textrm{Uniform}(\textrm{SO}(D))\) is the uniform distribution over the rotation group \(\textrm{SO}(D)\), representing our outlier distribution. \(\textrm{SO}(D)\) is the space where every element represents a D-dimensional rotation. This mixture model serves as the likelihood \(P(O|q)\) for our Bayesian framework.
Sampled Result
Ground truth samples of \(q_i\)
\(\hat{q}_i \sim p(q|O)\) (estimated \(q_i\)) sampled from our posterior \(p(q|O)\)
The orientations \(\hat{q}_i\) sampled from the posterior look significantly rotated with respect to the original samples. Notice this is a global rotation since all the samples are rotated equally. This problem of global ambiguity of absolute orientations \(q_i\) arises from the fact that the relative orientations \(q_i q_j^T\) and \(\tilde{q}_i \tilde{q_j}^T\) of two different set of vectors can be the same even if the absolute orientations are different. The following section goes over this and provides a sneak peek into our solution for alleviating this problem.
Continuous Label Switching
A careful observation of our problem formulation (Equation 1) would reveal the problem is invariant to transformation of the absolute orientations as long as the relative orientations are preserved. Consider the 2 pairs of observations in the figure below. (Blue and Red; Yellow and Green) Let the absolute orientations be \(q_1\), \(q_2\), \(\tilde{q}_1\) and \(\tilde{q}_2\) and relative orientations between the pairs be \(R_{12}\) and \(\tilde{R}_{12}\). As the absolute orientations \(q_1\) and \(q_2\) are equally rotated by a rotation matrix \(A\), the relative orientation between them \(R_{12}\) is preserved.
More formally, Let \(A \in \textrm{SO}(D)\) be a random orientation matrix in D-dimensions. The following equation demonstrates how rotating two absolute orientations \(q_i\) and \(q_j\) by a rotation matrix \(A\) preserves the relative orientations — which in turn gives rise to a global ambiguity in our framework.
\[
R_{ij} = q_i q_j = q_i A A^T q_j = (q_i A) (A^T q_j) = \tilde{q_i} \tilde{q_j}
\]
Since our inputs to the model are relative orientations, this ambiguity (known as label switching) causes our Bayesian estimates to come randomly rotated by some rotation \(A\).
Proposed Solution
Based on Monteiller et al., in this project we explored a novel solution for alleviating this problem. The intuition is that we believe the unknown ground truth is close to the posterior samples up to a global rotation. Hence we try to approximate the ground truth by starting out with a random guess and optimizing for the alignment map between the estimate and the ground truth. Using this alignment map, and the posterior samples, we iteratively update the guess, using a custom Riemannian Stochastic Gradient Descent over \(\textrm{SO}(D)\).
Algorithm
Start with a random guess \(\mu \sim \mathrm{Uniform}(\mathrm{SO}(D))\).
Sample \(\hat{q} \sim P(q | O)\), where \(P(q | O)\) is the posterior.
Find the global ambiguity \(R\), between \(\hat{q}\) and \(\mu\). This can be obtained by solving for \(\mathrm{argmin}_R \, \| \mu – R \hat{q}\|_\mathrm{F}\).
Move \(\mu\) along the shortest geodesic toward \(\hat{q}\).
Repeat Steps 2 – 4 until convergence. Convergence is detected by a threshold on geodesic distance.
Results
We use this method to estimate the mean of the posterior over \(\mathrm{SO}(2)\) and plot the results (i.e. 2D orientations) in the complex plane as shown below.
Original Sample
Sampled Posterior
Optimized Mean Posterior
The proposed optimization proceedure is able to successfully re-align the posterior samples by alleviating the continuous label switching problem.
Conclusion
In conclusion, we study the rotation synchronisation problem from a Bayesian perspective. We explore a custom Riemannian Gradient Descent procedure and perform experiments in the \(\mathrm{SO}(2)\) case. The current method is tested on a simple toy dataset. As future work we are interested in improving our Bayesian model and benchmarking it against the current state-of-the-art. There are certain performance bottlenecks in our current architecture, which constrain us to test only on \(\mathrm{SO}(2)\). In the future, we are also interested in carrying out experiments more thoroughly in various dimensions. While the current MCMC procedure we are using does not account for the non-Euclidean geometry of the space of orientations, \(\mathrm{SO}(D)\), we are looking into replacing it with Riemannian versions of MCMC.
By Zhecheng Wang, Zeltzyn Guadalupe Montes Rosales, and Olga Guțan
Note: This post describes work that has occurred between August 9 and August 20. The project will continue for a third week; more details to come.
I. Introduction
For the past two weeks we had the pleasure of working with Nicholas Sharp, Etienne Vouga, Josh Vekhter, and Erik Amezquita. We learned about a special type of minimal surfaces: triply-periodicminimal surfaces. Their name stems from their repeating pattern. Broadly speaking, a minimal surface minimizes its surface area. This is equivalent to having zero mean curvature. A triply-periodic minimal surface (TPMS) is a surface in \(\mathbb{R}^{3}\) that is invariant under a rank-3 lattice of translations.
Figure 1. (Left) A Minimal Surface [source] and (right) a TPMS [source].
Let’s talk about nonmanifold surfaces. “Manifold” is a geometry term that means: every local region of the surface looks like the plane (more formally — it is homeomorphic to a subset of Euclidean space). Non-manifold then allows for parts of the surface that do not look like the plane, such as T-junctions. Within the context of triangle meshes, a nonmanifold surface is a surface where more than 3 faces share an edge.
II. What We Did
First, we read and studied the 1993 paper by Pinkhall and Polthier that describes the algorithm for generating minimal surfaces. Our next goal was to generate minimal surfaces. Initially, we used pinned (Dirichlet) boundary conditions and regular manifold shapes. After ensuring that our code worked on manifold surfaces, we tested it on non-manifold input.
Additionally, our team members learned how to use Blender. It has been a very enjoyable process, and the work was deeply satisfying, because of the embedded mathematical ideas intertwined with the artistic components.
III. Reading the Pinkall Paper
“Computing Discrete Minimal Surfaces and Their Conjugates,” by Ulrich Pinkall and Konrad Polthier, is the classic paper on this subject; it introduces the iterative scheme we used to find minimal surfaces. Reading this paper was our first step in this project. The algorithm that finds a discrete locally area-minimizing surface is as follows:
Take the initial surface \(M_0\) with boundary \(∂M_0\) as first approximation of M. Set i to 0.
Compute the next surface \(M_i\) by solving the linear Dirichlet problem $$ \min_{M} \frac{1}{2}\int_{M_{i}}|\nabla (f:M_{i} \to M)|^{2}$$
Set i to i+1 and continue with step 2.
The stopping condition is \(|\text{area}(M_i)-\text{area}(M_{i+1})|<\epsilon\). In our case, we used a maximum number of iterations, set by the user, as a stopping condition.
There are additional subtleties that must be considered (such as “what to do with degenerate triangles?”), but since we did not implement them — their discussion is beyond the scope of this post.
IV. Adding Periodic Boundary Conditions
This is, at its core, an optimization problem. To ensure that the optimization works, the boundary conditions have to be periodic instead of fixed in space. This is because we are enforcing a set of boundary conditions on periodic shapes — that is, tiling in a 3D space.
IV(a). Matching Vertices
First, we need to check every two pairs of vertices in the mesh. We are looking to see if they have identical coordinates in two dimensions, but are separated by exactly two units in the third dimension. When we find such pairs of vertices, we classify them to \(G_x\), \(G_y\), or \(G_z\). Note that we only store unique pairs.
IV(b). Laplacian Smoothness at the Boundary Vertices
Instead of using the discrete Laplacian, now we introduce a sparse matrix K to adjust our smoothness term based on the new boundary:
Next, we construct the matrix K, which is a sparse square matrix of dimension #vertices by #vertices. To do so, we set \(K(i,i) = 1\), \(K(i,j) = 1\), and set the \(j\)th row entries to 0 for every pair of unique matched boundary vertices. For every interior vertex \(i\), we set \(K(i,i) = 1\).
IV(c). Aligning Boundary Vertices
Now we no longer want to pin boundary vertices to their original location in space. Instead, we want to allow our vertices to move, while the opposite sides of the boundary still match. To do that, we need to adjust the existing constraint term and to include additional linear constraints \(Ax=b\).
Therefore, we add two sets of linear constraints to our linear system:
For any pair of boundary vertices distanced by 2 units in one direction, the new coordinates should differ by 2 units.
For any pair of boundary vertices matched in the two other directions, the new coordinates should differ by 0.
We construct a selection matrix \(A\), which is a #pairs of boundary vertices in \(G\) by #vertices sparse matrix, to get the distance between any pair of boundary vertices. For every \(r\)th row, \(A(r,i)=1, A(r,j)=-1.\)$
Then we need to construct 3 \(b\) vectors, each of which is a sparse square matrix of the size [# of vector pairs of boundary vertices in G for a 3D coordinate (x,y,z)]. Based on whether at one given moment we are working with \(G_x\), \(G_y\), or \(G_z\), we enter \(2\) for those selected pairs, and \(0\) for the rest.
V. Correct Outcomes
Below, we can see the algorithm being correctly implemented. Each video represents a different mesh.
VI. Aesthetically Pleasing Bugs
Nothing is perfect, and coding in Matlab is no exception. We went through many iterations of our code before we got a functional version. Below are some examples of cool-looking bugs we encountered along the way, while testing the code on (what has become) one of our favorite shapes. Each video represents a different bug applied onto the same mesh.
VII. Conclusion and Future Work
Further work may include studying the physical properties of nonmanifold TPMS. It may also include additional basic structural simulations for physical properties, and establishing a comparison between the results for nonmanifold surfaces and the existing results for manifold surfaces.
Additional goals may be of computational or algebraic nature. For example, one can write scripts to generate many possible initial conditions, then use code to convert the surfaces with each of these conditions into minimal surfaces. An algebraic goal may be to enumerate all possible possibly-nonmanifold structures and perhaps categorize them based on their properties. This is, in fact, an open problem.
The possibilities are truly endless, and potential directions depend on the interests of the group of researchers undertaking this project further.
Most branches of geometry are viewed as either (1) subsets of pure mathematics or (2) disciplines derived from problems in engineering, architecture, and computational mathematics. But what is geometry itself?
Our work in the Summer Geometry Institute has been in the modern field of geometry processing. In this post, we will consider the philosophical aspects of geometry, and what they tell us about geometry processing specifically. The discussion has an inherent Eurocentric bias, and is, unfortunately, limited to Caucasian male philosophers.
The Roots of Geometry
In the eighteenth century, mathematics was primarily used to understand the natural sciences. For example, it was perceived as a language to better articulate physics, rather than as having intrinsic value. At that time, geometry was perceived as a form of applied mathematics. But where does the philosophy of geometry have its roots?
First we have Kant (1724-1804), who believed that our grasp of euclidean geometry must be instinctive. Opinions about this vary, based on each individual’s definitions of a priori knowledge and empiricism. Helmholtz, for example, agreed with Kant.
However, geometry processing incorporates not only euclidean geometry, but also calculus, linear algebra, differential equations, and topology. In addition, geometry processing includes a framework of algorithmic and optimization techniques. The knowledge of these topics varies from person to person, and therefore geometry processing can not be known a priori.
Back to the 18th century: as logical empiricism became popular, the progress of philosophy of geometry stalled. Geometry was reduced to a system of definitions and conventions; logic and mathematics were given priority. As a result, the last big wave of interest in the philosophy of geometry was in the 1920s — 100 years ago, and 116 years after Kant!
Space as a Mathematical Concept
Mathematical conclusions depend on the structure we confer to a space. The structure determines the permissible elements, as well as the operations that can be applied to the elements. Without imposing structure, it is difficult to impose “rules,” and the potency of math is, thus, diluted. Most spaces in euclidean geometry exist; in fact, they are embedded in the reality we live in. With other kinds of geometry, we must define the properties of the space first. Then, we can ask: what are the conditions for creating a space? What kinds of spaces, and therefore geometries, can we have? (Our personal favorite is the geometry of chance.)
Helmholtz (1821-1894) said the only allowable spaces are those that maintain constant curvature. This was eventually challenged by Einstein’s (1879-1955) Theory of Relativity. Weyl (1885-1955) came up with intuitive spaces, where it is possible to compare lengths only if two edges share the same spatio-temporal point. Similarly to the earlier-mentioned Kant (1724-1804), Carnap (1891-1970) claimed that intuition of spaces requires an n-dimensional topological space. As the philosophy of geometry evolved, so did our understanding of space.
“Waiting for Godot” in a Hyperbolic Space?
(The authors would like to gratefully acknowledge the primary source for the following sections.)
So what is geometry? Many people study it or use it… but what exactly are we interacting with?
As an exercise of imagination, consider being a two-dimensional (2D) entity in a 2D universe. In fact, you don’t even have to imagine, somebody already wrote about this. The book is, unfortunately, at times misogynistic, but the geometry is intriguing to consider. Back to the exercise of imagination: since everything in the 2D universe is finite, is the universe itself finite? A 2D mathematician might suggest that the universe looks like a rectangle. Does, then, the universe have a point beyond which travel is impossible? If so, it likely is different from all other points of the universe, by its construction, since it has a boundary. Imagine you are playing Asteroids. As you move your rocket ship towards the right — once it reaches the boundary, it will reappear on the left side. In mathematical language we say that “the left edge of the rectangle has been identified, point by point, with the right edge.”
Figure 1. A finite 2D world with no boundary. Source.
How to Build a Torus When You Are Two-Dimensional
In three dimensions, we can glue these boundaries. First, we bend the rectangle to produce a cylinder, while being careful to glue only the left and right edges together. Since we are operating in a ‘perfect’ thought-based mathematical world, the cylinder is malleable. Now we bend it to achieve this second gluing, which looks like a donut. We call this a torus.
Recall that our reader currently lives in a 2D universe. A 2D person would not be able to see this torus in multiple dimensions. However, one would understand the space perfectly well, since the space can be identified back to the initial rectangle with edges. In this case, the universe is clearly finite, and with no edges.
Another surface with similar properties would be a two-dimensional sphere. A ladybug traveling on the sphere will notice that, locally, their world looks like a flat plane, and that the surface has no edges. On a small scale, if the mathematician ladybug splits the plane into triangles, the sum of the angles in each triangle will be 180 degrees. This is, in fact, known as a defining feature of euclidean geometry. However, on a larger scale, things start breaking.
A very large triangle drawn on the surface of the sphere has an angle sum far exceeding 180 degrees. Consider the triangle formed by the North Pole and two points on the equator of the sphere. The angle at each point on the equator is 90 degrees, so the total sum exceeds 180 degrees. Therefore, we conclude that — on a global scale — we can not apply euclidean geometry to a sphere. We call this hyperbolic geometry.
Figure 3. Angles add up to more than 180 degrees. Source.
Shape and Connection
To conclude this exercise of imagination, we must say that there is a wonderful connection between the shape (topology) of a surface and the type of geometry it inherits. This relationship is given, mathematically, by the Gauss-Bonnet equation:
\(\kappa A = 2 \pi \xi\)
Where \(\kappa \) is the Gaussian curvature of a surface, A is the element of area of the surface, \(\pi \) is our good old friend 3.1415…, and \(\xi\) is the Euler characteristic of the surface. In other words, the “geometry” is on the left side of the equation and the “topology” is on the right side of the equation; and the equal sign shows how deeply connected they are.
Studying various shapes and how they connect is not only interesting, but also important. Particularly, if a two-dimensional entity can deduce what sort of global geometry reigns their immediate world, they may be able to deduce the possible shapes for the things in their larger universe.
To Really Conclude…
Although geometry processing has been around for a few decades, its philosophy is not yet well-studied. So far, this area shares themes with the broader philosophy of geometry, in certain topics:
a priori vs empiricist knowledge,
the evaluation of space and its elements,
geometric human intuition,
construction of geometries, and
the study of shape and connection.
Given the computational aspect of geometry processing, there are additional important philosophical questions that arise, however they go beyond the scope of this post. It is difficult to determine whether such questions are best understood from the perspective of philosophy of geometry or philosophy of computer science, or both.
As geometry processing continues to impact our technological world, the philosophy of geometry processing will grow too. In the meantime, we should continue critically questioning our programming, our modeling, and — most importantly — the ethical implications of our work. We must also ensure that we have a functional understanding of the parts that constitute geometry processing. At the same time, it must be said that the discipline is more than just the sum of its parts.
Architects often want to create visually striking structures that involve curved materials. The conventional way to do this is to pre-bend materials to form the shape of the desired curve. However, this generates large manufacturing and transportation costs. A proposed alternative to industrial bending is to elastically bend flat beams, where the desired curve is encoded into the stiffness profile of the beam. Wider sections of the beam will have a greater stiffness and bend less than narrower sections of the beam. Hence, it is possible to design an algorithm that enables a designer or architect to plan curved structures that can be assembled by bending flat beams. This is a topic currently being explored by Bernd Bickel and his student Christian Hafner (our project mentor).
A model created using the concept of active bending (source)
For the flat beam to remain bent as the desired curve, we must ensure that the beam assumes this form at its equilibrium point. In Lagrangian mechanics, this occurs when energy is minimised, since this implies that there is no other configuration of the system which would result in lower overall energy and thus be optimal instead. Two different questions arise from this formulation. First, given the stiffness profile \(K\), what deformed shapes \(\gamma\) can be generated? Second, given a curve \(\gamma\), what stiffness profiles \(K\) can be generated? In answering the first problem we will find the curve that will result from bending a beam with a given width profile. However, we are more interested in finding the stiffness profile of a beam which will result in a curved shape \(\gamma\) of our choice. Hence, we wish to solve the second problem.
Below, we will give a full formulation of our main problem, and discuss how we transformed this into MATLAB code and created a user interface. We will conclude by discussing a more generalised case involving joint curves.
Problem formulation
To recap, the problem we want to solve is, given a curve \(\gamma \), what stiffness profiles \(K\) can be generated? For each point \(s\) on our curve \(\gamma: (0,\ell) \to \mathbb{R}^2\), we have values for the position \(\gamma(s)\), turning angle \(\alpha(s)\), and signed curvature \(\alpha'(s)\). In addition, we define the energy of the beam system to be \[W[\alpha] := \int_0^\ell K(s)\alpha'(s)^2 ds,\] where \(K(s)\) is the stiffness at point \(s\). At the equilibrium state of the system, \(W[\alpha]\) will take a minimum value. Intuitively, this implies that, in general, regions of larger curvature will have a lower stiffness. However, it is not true that two different points of equal curvature will have equal stiffness, since there are other factors at play.
Now, before finding the \(K\) that minimises \(W[\alpha]\), we must set additional constraints \[\alpha(0) = \alpha_0, \quad \alpha(\ell) = \alpha_\ell\] and \[\gamma(0) = \gamma_0, \quad \gamma(\ell) = \gamma_\ell.\] In words, we are fixing the turning angles (tangents) and positions of the two boundary points. These constraints are necessary, since they dictate the position in space where the curve begins and ends, as well as the initial and final directions the curve moves in.
We have now formulated a problem that can be solved using variational calculus. Without going into detail, we find that stationary points of this function are given by the equation \[K \kappa = a + \langle b, \gamma \rangle,\] where \(\kappa\) is the signed curvature (previously \(\alpha’\)), and \(a \in \mathbb{R}\) and \(b \in \mathbb{R}^2\) are constants to be found. However, a stiffness profile cannot be generated for all curves \(\gamma\). More specifically, it was shown by Bernd Bickel and Christian Hafner that a solution exists if and only if there exists a line \(L\) that intersects \(\gamma\) exactly at its inflections. With this information in hand, we can begin to create a linear program that computes the stiffness function \(K\).
The top four curves can be created using active bending, but the bottom four cannot (source)
Creating a linear programme
In order to find the stiffness \(K\), we need to solve for the constants \(a\) and \(b\) in the equation \( K \kappa = a + \langle b, \gamma \rangle \), which can be discretised to \[ K(s_i) = \frac{a + \langle b, \gamma(s_i) \rangle }{\kappa(s_i)}.\] It is possible that there is more than one solution to \(K\), so we want some way to determine the “best” stiffness profile. If we think back to our original problem, we want to ensure that our beam is maximally uniform, since this is good for structural integrity. Hence, we solve for \(a\) and \(b\) in the above inequality in such a way that the ratio between maximum and minimum stiffness is minimised. To do this we set the minimum stiffness to an arbitrary value, for example 1, and then constrain \(K\) between 1 and \(M = \text{min}_i K_i\), thus obtaining the inequality \[1 \leq \frac{a + \langle b, \gamma(s_i) \rangle }{\kappa(s_i)} \leq M.\] This can be solved using MATLAB’s linprog function (read the documentation here). In this case, the variables we want to solve for are \(a \in \mathbb{R}\), \(b \in \mathbb{R}^2\), and \(M\) (a scalar we want to minimise). So, using the linprog documentation \(x\) is the vector \((a,b_1,b_2,M)\) and \(f\) is the vector \((0,0,0,1)\), since \(M\) is the only variable we want to minimise. Since linprog only deals with inequalities of the form \(A \cdot x \leq b\), we can split the above inequality into two and write it in terms of the elements of \(x\), like so: \[-\Bigg(\frac{1}{\kappa(s_i)}a + \frac{\gamma_x(s{i})}{\kappa(s_i)}b_1 + \frac{\gamma_y(s_{i})}{\kappa(s_i)}b_2 + 0\Bigg) \leq -1,\] \[\frac{1}{\kappa(s_i)}a + \frac{\gamma_x(s_{i})}{\kappa(s_i)}b_1 + \frac{\gamma_y(s_{i})}{\kappa(s_i)}b_2 -M \leq 0.\] These two equations are of a form that can be easily written into linprog to obtain the values of \(a\), \(b\) and \(M\) and hence \(K\).
User Interface
Once we solved the linear program outlined above, we created a user interface in MATLAB that would allow users to draw and edit a spline curve and see the corresponding elastic strip created in real time. Custom splines can be imported as a .txt file in the following format
or alternatively, the file that is already in the folder can be used. Users can then run the user interface and edit the spline in real time. To add points, simply shift-click, and a new point will be added at the midpoint between the selected point and the next point. The user can right-click to delete a point, and left-click and drag to move points around. If there are zero or one control points remaining, then the user can add a new point where their mouse cursor is by shift-clicking. The number of control points must be greater than or equal to the degree plus one for the spline to be formed. There are certain cases in which the linear program cannot be solved. In these cases the elastic strip is not plotted, and the user must move the control points around until it is possible to create a strip.
Demonstration of the user interface
Here is the link to the GitHub repository, for those who want to try the user interface out.
Joints Between Two or More Strips
The current version of our code is able to generate elastic strips for any (feasible) spline curve generated by the user. However, an as of yet unsolved problem is the feasibility condition for a pair of elastic strips with joints. Solving this problem would allow us compute a pair of stiffness functions \(K_1\) and \(K_2\) that yield elastic strips that can be connected via slots at the fixed joint locations. With a bit of maths we were able to derive the equilibrium equations that would produce such stiffness functions. Suppose the two spline curves are given by \(\gamma_{1}: [0,\ell_{1}] \to \mathbb{R}^2\) and \(\gamma_{2}: [0,\ell_{2}] \to \mathbb{R}^2\) and that they intersect at exactly one point such that \(\gamma_{1}(t_a) = \gamma_{2}(t_b)\). Then we must solve the following two equations: \[ K_{1}(t) \kappa_{1}(t) = \left \{ \begin{array}{ll} a_1 + \langle b_1, \gamma_{1}(t) \rangle + \langle c, \gamma_{1}(t) \rangle & 0 \leq t \leq t_a \\ a_1 + \langle b_1, \gamma_{1}(t) \rangle + \langle c, \gamma_{1}(t_a) \rangle & t_a < t \leq \ell_{1} \end{array}\right. \] \[ K_{2}(t) \kappa_{2}(t) = \left \{ \begin{array}{ll} a_2 + \langle b_2, \gamma_{2}(t) \rangle – \langle c, \gamma_{2}(t) \rangle & 0 \leq t \leq t_b \\ _2 + \langle b_2, \gamma_{2}(t) \rangle – \langle c, \gamma_{2}(t_b) \rangle & t_b < t \leq \ell_{2} \end{array}\right. \]
Similar to the one-spline case, these can be translated into a linear program problem which can be solved. Due to the time constraint of this project, we were not able to implement this into our code and build an associated user interface for creating multiple splines. Furthermore, above we have only discussed the situation with two curves and one joint. Adding more joints would increase the number of unknowns and add more sections to the above piece-wise-defined functions. Lastly, we are still lacking a geometric interpretation for the above equations. All of these issues regarding the extension to two or more spline curves would serve as great inspiration for further research!
Last week I was working on the project “Robust computation of the Hausdorff distance between triangle meshes” under Dr. Leonardo Sacht’s supervision with TA Erik Amezquita and SGI fellows Bryce Van Ross and Deniz Ozbay.
If we have triangle meshes \( A, B \), then \( h(A,B) = \max_{p \in A}d(p,B) \) is called the Hausdorff distance from \( A\) to \( B \), where \(d\) is Euclidean distance. In general case, this function is not symmetric, so the final metric is defined as \(H(a,b) = \max(h(A,B), h(B,A))\).
The Hausdorff distance is very significant in Geometry Processing on the grounds that it may be used for determining the difference between two meshes. The method that we are studying is called the “branch-and-bound” method. The main idea is to calculate the common lower bound for distance from the whole mesh \( A\)to mesh \(B\) and individual upper bounds for distances from every triangle mesh \(T_{A}\) of mesh \( A \) to mesh \(B\). If the upper bounds of some triangles is smaller than the lower bound, then we throw them away and consider the remaining subdivided ones.
My task was to code the function that returns the upper bounds for distances from every triangle mesh \(T_{A}\) of mesh \( A \) to mesh \(B\). We are going to use the distances between vertices of triangle mesh \(T_{A}\) and triangle inequality to find the upper bound: \[u( T_{A}, B) = \max_{j = 1}^{3}(\max(|v_{j} – v_{j + 1}|, |v_{j} – v_{j + 2}|) + \min_{b \in B}(|v_{j} – b|))\] where \(v_{1}, v_{2}, v_{3}\) are the vertices of triangle mesh \(T_{A}\) .
This is how the implemented algorithm looks like in Matlab:
I really enjoyed working on the project this week with my team, TA and supervisor. I am looking forward to continuing work and improving what we have already done.
For the past two weeks, we have been doing research on how to optimize making 2D cuts onto a surface such that we get an optimal parameterization of the surface onto a disk. This question more generally is asking how we can take a triangulated surface and perform operations on it such that we see a representation of this surface on a disk that we can apply texture mapping or other various changes to.
Firstly let us explain why we want to make these cuts. We make these cuts so that we can use existing methods of disk parameterization from a surface with a boundary onto a disk. We can also think about making a 2D surface from our 3D surface using cuts in a more intuitive way. We can think about it like how a net works, that we know it is true that if we took enough cuts it would be true that we get a locally distortion free 2D representation of the surface; however doing this would both be computationally intensive and not very useful for performing tasks that we want, such as texture mapping, and it does not give an intuitive understanding of what the general surface is or what each section of the 2D map corresponds with each section of the surface. Thus the question becomes how do we take a minimal amount of cuts such that we have a 2D parameterization of the surface that minimizes distortion such that we can still use it for various texture mapping applications.
Disk Parameterization:
First we worked on mapping 3D meshes with existing boundaries to a disk. We implemented existing methods of disk conformal parameterization using three different weight functions. We start by mapping the mesh boundary to the circle, preserving the distance between chord edges. Placing the interior points is then a matter of solving a system of linear equations. In our first weighting method, each edge is weighted the same so in the resulting disk, interior points are the centroid of their neighbors. This gives one linear equation for each interior point, and since the boundary points are fixed, we can solve for the exact locations of each point on the disk.
Weighing each edge uniformly results in unnecessary distortion. Our second and third methods use the angles in the faces sharing the edge to add weights to the equations in the first method. The second method uses harmonic weights, with formula \(w_{i,j}=\frac{1}{2}\left( cot(i,j) + cot(i,j)\right)\) where \(\alpha_{i,j}\) and \(\beta_{i,j}\) are the angles opposite the edge connecting the \(i\) and \(j^{th}\) vertices. These weights reduce distortion, but obtuse angles can result in negative weights, so the resulting map is not guaranteed to be bijective. The third method, mean value weights, resolves this by instead using \(w_{i,j}=\frac{tan(\frac{\gamma_{i,j}}{2}) + tan(\frac{\delta_{i,j}}{2})}{2|| v_i – v_j ||}\) where \(\gamma_{i,j}\) and \(\delta_{i,j}\) are the angles adjacent to the edge connecting the \(i\) and \(j^{th}\) vertices. In our future disk parameterizations, we use these last two methods.
Methods 1, 2, and 3 from left to right
Virtual Boundary:
A downside to disk parametrization is that because the boundary has a fixed mapping to the circle, the resulting parameterization has a lot of distortion near the perimeter. There are a variety of free-boundary parameterization methods, but we focused on the virtual boundary method. By adding one or more layers of faces to the existing boundary, we create a new, virtual boundary that is mapped to the circle while the real boundary is allowed a less constrained shape. We can then ignore the added faces and use the resulting parameterization. The more layers in the virtual boundary, the more “free” the resulting boundary. For example, the mountain range mesh below has a boundary that is more rectangular in shape than circular. As we add more layers to the virtual boundary parameterization, we can see the outline become more square as well.
1 Layer
5 Layers
10 Layers
However, it is not always necessary to add multiple layers, but instead simply increase the size of the virtual boundary, this is due to how the virtual boundary is created following the method used in “Parametrization of Triangular Meshes with Virtual Boundaries” by Yunjin Lee, Hyoung Seok Kim, and Seungyong Lee. We create the virtual boundary by listing the vertices in the existing boundary as \(v_1,v_2,\dots,v_n\) and label the vertices in the virtual boundary \(v’_1,v’_2,\dots,v’_{2n}\) where \(v_i\) corresponds to \(v’_{2i}\), we then make sure the distance between \(v_i , v_{i+1} \) proportional to the distance between \(v_{2i} , v_{2i+2}\) and assigning \(v_{2i+1}\) to the midpoint between the two, to do this we simply say that \(v_{2i}=av_i\). Thus we can increase the size of the virtual boundary simply by increasing the a value, this allows a quicker method of getting a virtual boundary that gives the same results as running multiple layers of virtual boundaries.
a = 1.05
a = 1.5
Many meshes that we want to parameterize in 2D do not have boundaries. In order to do this, we make a cut along existing edges in the mesh to create our own boundary. A cut can be visualized as exactly how it sounds—taking a pair of scissors and slicing through the mesh. It is created from a path of edges by duplicating the non-endpoint vertices in the path, and connecting them to create a loop of boundary edges. This connected loop will naturally form by reassigning the faces on one side of the path to contain the new vertices rather than the old duplicated vertices.
The trouble comes with determining which side of the cut each face is on. We created our own cutting algorithm, and solved this problem by using the orientation of the faces with edges on the path. Given a directed path, the faces with edges on one side of the path will be oriented in the direction of the path, while the faces on the other side will be oriented in the opposite direction. Thus we can easily separate the faces sharing edges along the cut. A more difficult task is separating the faces that only share one vertex with the cut path. To solve this quickly, we store the directions of edges coming from each point on the path that are confirmed to be on one side of the path. We then can determine which side of the path a face is based on how close its edges’ direction is to these edges. There are conceivable shares which will cause this to fail, but in practice this can quickly implement the vast majority of cuts.
We also tried a different algorithm to determine which side of the cut each face is on, making use of geometric properties of the surface, in specific the normal of the surface. We compared the normal of the face to the normal of the edge on the cut. We took the plane defined by the normal on the cut, and compared if the normal of the face projected onto the plane was on the left or right of the plane, when considering the y-axis of the plane to be the edge. This result showed initial success, and continuous to work along small cuts, or cuts along areas of low curvature, however when dealing with areas of large curvature the cut has a large issue, that is when the surface has high curvature it is possible for the edge to be behind the surface with respect to where the edge is, which gives the normal on the other side than we desire, and thus gives it to be on the right when it is on the left, and other such issues, for this reason we have discarded this approach, however we are open to the possibility that using some operation on the normals it might be possible to use this approach.
Euclidean Max Cut:
Now let us discuss our main method of how to choose a cut to make, we considered existing methods, namely the method of using a minimal spanning tree as used by Alla Sheffer in “Spanning Tree Seams for Reducing Parametrization Distortion of Triangulated Surfaces,” seemed to have the opposite desire than we want, that is, to minimize the size of the cut. It is clear simply from the boundaries we were provided as well as intuition that the larger the cut the smaller the distortion would be, as if we had a large enough cut there would be no distortion, however we also know that artifacts exist along the cuts, thus we created a new method. This method was a way of maximizing the distance of the cut while minimizing the artifacts. This was maximizing the Euclidean distance between the end points of the path of the cut while making the cut that takes the minimum geodesic distance.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
As we can see however this creates unintuitive cuts, and cuts that do not have our desired purpose, for this we considered a new method of cuts.
Medial Mesh Max Cut:
To deal with this issue, we implemented a new method using the medial mesh, that is we consider a simplified mesh of the object making use of the midpoints of the faces within the interior of the surface. We did not do implementation of finding this medial mesh ourselves, rather using the code created in “Q-MAT: Computing Medial Axis Transform by Quadratic Error Minimization” by Pan Li, Bin Wang, Feng Sun, Xiaohu Guo, Caiming Zhang and Wenping Wang. We then after obtaining the Medial Mesh, create the minimal spheres on each vertex, that is the smallest sphere with the center at a vertex in the medial mesh, such that it intersects with a vertex on the surface. We then created a graph with vertices representing the spheres, and creating an edge if the spheres overlap, we then calculated the maximum path existing on the graph, with the edge lengths equal to the distance between the center of the spheres. We then identified the vertices on the surface such that they were the closest to the points on each sphere such that the Euclidean distance between the two end point spheres is maximized, we use these as the two points to create a cut on and we create a cut with the minimum geodesic distance between the two points. This was hopefully to avoid areas of high curvature, and thus create a better cut still of a reasonably large size.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
Of note is that areas that are not on the cut area directly are “tucked in” so to into the areas that we can barely see. The one below is in the left area with a high density of vertices. The rest follow from left to right the areas of high vertices.
The surface, and the cut we are making upon it.
A color function we add to the surface in order to see where different areas are mapped onto a 2D Parameterization.
The 2D Parameterization with Virtual Boundary.
We again get areas of high curvature being sort of tucked into the surface.
The bottom left area that shows the face.
The right middle area showing the back legs and tail.
How Multiple Cuts work:
Thus far we have only discussed and shown single cuts along simple paths. It is also possible, and often preferred, to make multiple cuts on intersecting paths. Making these cuts sequentially, each new cut is equivalent to making a cut on a shape that already has a boundary. If the cut path has one endpoint on the boundary and the rest of the points on the interior, then the resulting boundary is still a simple closed path, as is necessary for our 2D parameterization. The method is the same as before, except the endpoint on the existing boundary is duplicated in addition to the other points. Using multiple cuts increases the freedom of our cut paths and allows us to cut along branching pieces of a mesh.
We can extend our previous Euclidean maximum cut method to multiple cuts in an iterative process. After making our first cut as usual, we repeat the same method but restricting one of the endpoints to be on the previous cuts. Our second cut selects the point that is farthest from our new boundary and chooses the shortest path between them. This is repeated until there are the desired number of cuts. Below are some examples of multiple cuts using this method. One possible automated method of stopping would be to calculate the standard deviation of the distance from the boundary and continue to create cuts until it is below some threshold.
The surface, and the cuts we are making upon it.
The 2D Parameterization and a color function that we take from 2D parameterization to the surface.
The surface, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
The surface, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
[example of cutting through diagonals and where we don’t want, downside of this method]
We can also use the sphere approximation to make multiple cuts. The spheres’ centers indicate important spots in the mesh, like joints and protrusions. We decided to use a minimum spanning tree to determine a sequence of cuts that connects each of these areas. First we project each sphere center to a vertex on the mesh. Then we make a weighted graph out of these points where an edge between two points is weighted as their geodesic distance. This graph gives a minimum spanning tree that can be broken into a sequence of cut paths that connects each point from our sphere approximation. The paths sometimes require simple adjustments to remove overlap, and then they are ready for our cutting algorithm.
The surface, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
The surface, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
The surface, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
The surface from the front, and the cut we are making upon it.
The surface from the back, and the cuts we are making upon it.
The 2D Parameterization a color function that we take from 2D parameterization to the surface.
Conclusion:
This research was very interesting, but we did not get the sort of success that we would like: While we did get surfaces with minimal distortion, specifically those with multiple cuts, we did not get the desired usage that we would like. Thus it is necessary to state what future direction should be taken if this research would be continued. There are three main areas for future work.
The first is making cuts along lines that preserve symmetry, that is if the surface has symmetry along some axis, or some feature, the 2D parameterization also has these symmetries. This would be useful both for minimizing the distortion and would provide intuitive results for texture mapping. These lines of symmetry, however, would be incredibly difficult to calculate, so we would like for them to be usually user specified, which brings us to the second area: how the inclusion of landmark areas would be useful, that is specifying areas where the cut cannot travel through. This allows us to intuitively avoid areas of high curvature that are important to the surface and how it could be applied to different surfaces as useful for an artistic application.
We could also try to do this in an automated way by using a skeletal approximation of the surface, using either a simplified medial mesh or a method presented in the SGP 2021 Graduate School Course Shape Approximations & Applications: using the sphere mesh proxy simplification of a surface. This method considers the edges between the spheres as the “bones” and makes a cut on each bone section such that it goes over the entirety of the bone, but avoids areas of high curvature by making minimal geodesic cuts along each bone; the method connects the bones by paths the minimize the total geodesic distance. A secondary way to use this is to make cuts along the “bones” but additionally taking a cut perpendicular to the bone at each point where the sphere begins. To think of this on a 3D model of a cat, for example, we would take a cut along the spine, and then take a perpendicular cut along the cat. To see if either of these methods work and if so which is better, further research would be needed.
AKA: Takeaways from Weeks 1-2 of High-order Directional Field Design Research
Author: Bryce Van Ross
It’s incredible how much one can learn in a month, and I’m looking forward to learning more (especially theoretically). In that same spirit, I highlight some takeaways from research made in my first SGI project. This project was guided by research mentor Dr. Amir Vaxman and TA Klara Mundilova, where I worked with fellow SGI student Jonathan Mousley.
Question(s): Usually we think of vectors and computations like divergence and curl as interrelated, and they are. But can we determine something more nuanced about these properties with respect to some vector field if we encounter a complex (i.e. having multiple vertices/edges) triangular mesh? Yes, we can. But it depends on your choice of approach of partitioning your mesh. For the sake of my research, we focus on the face-based representation.
Face-based representations of vector fields can then be broken down into vertex-based and edge-based approaches, per face per triangle. This means we are working with vector fields on faces that are gradients of (piecewise linear) functions that are either defined on the vertices or on the midpoints of edges. Depending on the choice of approach, then your computations are different. But which way is better and what are the consequences?
Answer(s): This is a natural question. Vertex-based (for certain reasons) seems to yield better approximations, which lead to better attempts at mimicry of continuity. In this sense, vertex-based computations are considered conforming, whereas edge-based computations are deemed nonconforming (w.r.t. continuity). Suppose we wanted to express a given vector field u in terms of familiar computations. Naturally, we would prefer to use vertex-based computations. However, we must remember that degrees of freedom (D.O.F.) must be maintained. Surprisingly, using purely vertex-based computations (or, purely edge-based computations) are in violation of D.O.F. More surprisingly, we find that our only solution is to use a mix of both the conforming and nonconforming terms. So, even though the gradients are distinct, both are equally valuable in terms of a reduction of u. So, there is a need to incorporate both \(G_v\) (the vertex-based gradient) and \(G_e\) (the edge-based gradient). This mixture could be complicated, but isn’t…it only requires the sum of 3 terms. The first term includes \(G_v\) and computes the divergence but is curl-free. For the second term, including \(G_e\), it computes the curl yet is div-free. The last term, referred to as \(h\), is both divergence-free and curl-free. Note: that \(G_v\) and \(G_e\) can be interchanged w.r.t. the first two terms if such equation is multiplied by the rotation matrix \(J\). Ultimately, u (or any vector space) has non-trivial representation (a.k.a. there’s more than meets the eye).
Currently, we’re finishing the first half of the 2-week Robust computation of the Hausdorff distance between triangle meshes project. This research is lead by mentor Dr. Leonardo Sacht, TA Erik Amézquita, and in my immediate team, I work with fellow students Deniz Ozbay and Talant Talipov. The below is a brief summary of what’s happened so far:
A mathematical visualization of the Hausdorff distance of two meshes
Applications: Primarily computer vision, computer graphics, digital fabrication, 3D-printing, and modeling. For example, in computer vision, it is often desirable to identify a best-candidate target relative to some initial template. In reference to the set of points within the template, the Hausdorff distance can be computed for each potential target. The target with the minimum Hausdorff distance would qualify as being the best fit, ideally being a close approximation to the template object.
Motivation: Objects are geometrically complex. There are different ways to compare objects to each other via a range of geometry processing techniques and geometric properties. Distance is often a common metric of comparison. But what type of distance should we use, which distances are favorable, and why? These are important questions.
Pitfalls of using other types of distance for triangular meshes
For our research, we focus on computing the Hausdorff distance \(h\). “Hausdorff” may seem familiar to you if you know topology. There, a (topological) space is considered Hausdorff if any two elements can be separated into disjoint (open) sets. The key idea here is the separation property with respect to points.
In geometry processing, this idea is extended (in some sense) to the separation of triangle meshes. The Hausdorff distance \(h\) is fundamentally a maximum distance among desirable distances between 2 meshes. These desirable distances are minimum distances of all possible vectors resulting from points taken from the first mesh to the second mesh. But why is \(h\) significant? If \(h\) converges to zero (the smallest possible distance), then this implies that our meshes, and therefore the objects themselves, are very similar. This, like most things in math, implies within some epsilon, representing marginal change such as a slight deformation, rotation, translation, compression, or stretch. If \(h\) is large, then this implies that the two objects are dissimilar. Intuitively, this is due to a lack of ideal correspondence from triangle to triangle. In short, \(h\) serves as a means of computing the similarity between digital objects in terms of maximally separating the meshes’ points according to their minimum distances.
Tasks: To compute \(h\), we find the maximizer, the point (in the first mesh) corresponding to the computation of \(h\). This point is found via an algorithmic process called the branch and bound technique. Sparing the details, the result of applying this technique will provide a (very small) region where the minimizer is claimed to exist, after a series of triangle subdivisions and deletions. There are different ways to implement this technique. Our collective goal this week was to achieve accurate \(h\) given any two meshes. Once this is ensured, Week 2 would focus on making our code more robust (efficient/fast). This work can be simplified into three primary tasks: encoding the lower bound and upper bound of possible values of \(h\), and the subdivision method. My work focused on the lower bound, for which the other two functions are dependent.
Accomplishments: By Day 2, we started with writing pseudocode for the lower bound, using only the vertices of the first mesh and computing their distances with respect to vertices of the second mesh. This was computed per face, and the minimum was found. Although this method was correct, it didn’t account for either the edge or interior cases of a given triangle. Looping through an edge would be immediately doable, but the interior case would be more challenging. Thankfully, Dr. Sacht had us search through gptoolbox for such functionality that would account for all three cases. Once finding this function, we were able to reduce my code to two lines! The irony is although this was valid, the referenced function was basically an empty shell. The function was actually calling a C++ function that would have to compiled and linked against other libraries and due to our time constraints we decided it was best to address this issue in a future time. Ultimately, we ended up having to write from scratch once more. In searching for similar point-triangle distance algorithms, we initially found an approach using normal vectors, which would offer projective power of the former mesh vertices relative to the latter mesh. The computations were erroneous and misrepresentative. Since then, our team has been using a combinatorial plane-vector approach to find the lower bounds per vertex.
Hopefully we’ll finish soon… we’re excited for the next steps!
Isogeometric analysis (IGA) is an analysis technique that combines computer aided design (CAD) with traditional finite element analysis (FEA). It uses the same spline basis functions to construct the geometry and the solution space, which is beneficial as traditional FEA requires geometric approximation that can lead to inaccuracies [1].
For our week 3 project, my group and I used IGA simulation software that has been developed for modelling material transport in neurons. The software solves Navier-Stokes equations to obtain the velocity field and models the transport process by reaction-diffusion-transport equations. For our purposes, we used the solver to simulate material transport in complex neurons and heat transfer processes in various geometries including a simple block model and a rod model [2].
Our project was split into two main parts: geometric modeling and analysis. For geometric modelling, we used two open source software packages (HexGen and Hex2Spline [3]) for the construction of geometries, and we then used the aforementioned IGA software for simulation purposes. I was in charge of using the IGA software to run the simulations and visualizing the results with Paraview (a visualization application).
In this post, I’ll demonstrate our results using the simple model the team created and also share some of the results I got using the IGA software on some additional models.
Here is the example model we used:
After the geometric modeling stage, I received 3 input files that contained the control mesh, an initial velocity field and the simulation parameters (particle concentrations, diffusion coefficient, velocities of material, etc.). These files were used as input for the IGA solver [2], which consists of four stages: spline construction, mesh partition for parallel computing, solving Navier-Stokes equations, and final transport simulation.
Due to the high computational needs of the last two stages, we set up the simulation environment at the Pittsburgh Supercomputer Center (PSC). We ran the simulation by connecting to the remote supercomputer via an ssh client. Finally, after getting the results, I was able to visualize them using Paraview.
In the above example, the color bar shows the transport velocity magnitude. I also added streamlines that show the direction of flow.
Here is an animation of our results that depicts the movement of the particles:
I also ran the solver on a rod model (that Angran Li, our mentor, provided):
The model shows the flow of heat through the rod as it is input from the right end of the model. Again, the colorbar represents the magnitude of the velocity of particles at that point.
Finally, I also ran the IGA solver on a neuron model (found on NeuroMorpho.org and edited with the help of Angran).
The geometry of the neuron (visualised using Paraview)
Images show results we got by running the IGA solver on the neuron geometryThe transport of material through the neuron over time
From using the remote computer to editing the neurons for our analysis, I learned a ton of new techniques during the week. I enjoyed learning about the IGA solver and its applications ranging from neuroscience to engineering. I’d like to end this post by thanking Angran his support and for bearing with me through my 1000 questions—it was a pleasure learning with you! 🙂
References:
[1] T. J. Hughes, J. A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.Computer Methods in Applied Mechanics and Engineering, 194(39-41), 4135-4195, 2005.
[2] A. Li, X. Chai, G. Yang, Y. J. Zhang. An Isogeometric Analysis Computational Platform for Material Transport Simulations in Complex Neurite Networks.Molecular & Cellular Biomechanics, 16(2):123-140, 2019. GitHub link: https://github.com/CMU-CBML/NeuronTransportIGA
[3] Y. Yu, X. Wei, A. Li, J. G. Liu, J. He, Y. J. Zhang. HexGen and Hex2Spline: Polycube-Based Hexahedral Mesh Generation and Unstructured Spline Construction for Isogeometric Analysis Framework in LS-DYNA. Springer INdAM Serie: Proceedings of INdAM Workshop “Geometric Challenges in Isogeometric Analysis”. Rome, Italy. Jan 27-31, 2020. GitHub link: https://github.com/CMU-CBML/HexGen_Hex2Spline
Written by Deniz Ozbay, Tal Rastopchin, and Alexander Rougellis
Surface Deformation and Remeshing
Representing a curved surface with a mesh of polygons can prove to be difficult because, given that meshes are discrete surfaces, there are times when small changes to the triangulation of the surface can cause big changes to the overall geometry and measured quantities (as can be seen with the example of the Schwarz Lantern). Sometimes we want to deform surfaces, and the deformation of the surface can result in a “bad” triangulation. Representing these surfaces with such deformations is done by optimizing these deformations \(\delta S\) on the given surface by solving \(\min_{\delta S} f(S + \delta S)\) where \(S\) is that given surface. When the real valued function \(f: S \rightarrow \mathbb{R}\) gives us surface area, our optimization will minimize the surface area, which is also known as mean curvature flow. Minimizing surfaces can be done with gradient descent (given a small triangle mesh) and Newton’s Method (given a large mesh). After minimizing, if the surface is determined to be “bad”, then remeshing is needed using remeshing operations such as edge flip, edge split, and edge collapse (to name a few basic operations that could be used).
Mean Curvature Flow
Given that one of the goals of this project is to optimize functions on meshes, our first task was to put together a simple implementation of mean curvature flow in MATLAB. Professor Etienne Vouga introduced mean curvature flow as a deformation of a surface that minimizes surface area. He explained that if we have a function for the surface area of a mesh, as well as the gradient, we can use an optimization method like gradient descent in order to compute the deformation induced by the mean curvature flow. After our introduction to geometry processing during the tutorial week we knew that we could use the doublarea function to write a function that returned the surface area of a mesh. However, computing the gradient of this function is tricky—what would even be the domain of the surface area function? If the doublearea function relies on computing triangle areas using the cross product, and we are summing the result over all triangles in a mesh, is the domain some sort of “collection of triangles” that we are summing over?
To answer this question, Professor Etienne Vouga pointed us to the paper “Can Mean-Curvature Flow be Modified to be Non-singular?” and explained that we could express the gradient of the surface area function as the Laplacian applied to the \(x\), \(y\), and \(z\) vertex coordinate columns. In particular, the paper explained that “informally, mean-curvature flow can be thought of as a flow that pushes a point on a surface towards the average position of its neighbors.” The paper specifically explains that when \(M\) is a two dimensional manifold, \(\Phi_t : M \rightarrow \mathbb{R}^3\) is a smooth family of immersions of the manifold \(M\), and \(g_t( \cdot, \cdot)\) is a metric induced by the immersion at time \(t\), we have that \(\Phi_t\) is a solution to the mean-curvature flow if
where \(\Delta_t\) is the Laplace-Beltrami operator defined with respect to the metric \(g_t\). If we interpret the Laplacian as a local averaging operator, this equation exactly captures the idea that mean-curvature flow is a deformation that pushes a point on a surface towards the average position of its neighbors.
For our implementation, Professor Vouga explained that instead of worrying about the metric induced by the flow we could just compute the Laplacian matrix for the first step and use it throughout the entire simulation. One thing we learned was that sometimes for computation and derivation it can be easier to look at the same function, like surface area, from a bunch of different perspectives. For example, once we knew that the gradient of the surface area was the Laplacian of the \(V\) matrix, we could compute the Hessian by reshaping the \(V\) matrix into a column vector and reshaping the Laplacian matrix into a larger block matrix with 3 copies of the Laplacian matrix along the diagonal. Computing the Hessian in this way allowed us to also try to implement a Newton’s method approach to minimize surface area.
After a lot of discussion, programming, testing, and playing around with initial conditions, we finally got our MATLAB implementation of mean curvature flow to work.
A figure displaying the open cylinder mesh used as a seed surface for the mean curvature flow.
A figure displaying the catenoid mesh resulting from the mean curvature flow.
The figure on the left is a cylinder that was our initial surface and the figure on the right is the resulting minimal surface produced by the mean curvature flow. It was really cool to see the end result—the resulting minimal surfaces can be very pretty. We’re excited to learn more about optimization on surfaces as well as how remeshing could potentially improve this simulation process.
Current and Future Work
Currently, we are working on implementing Newton’s method and the gradient descent method in the same program as remeshing. Our algorithm runs Newton’s method as long as it decreases the surface area. If a Newton step does not decrease the surface area, the algorithm then switches to the gradient descent method. Then it remeshes the triangles using edge flipping to obtain a Delaunay triangulation. To build this code, we made use of numerous different structures, like an edge face matrix, edge vertex matrix and a vertex face matrix. Then, we were able to implement the edge flipping condition. At this step we first find the neighboring two faces of the edge we check to flip. We check if the angles at the third vertices of the corresponding two faces is bigger than \(\pi\) radians, and if so, the edge is flipped to get the remeshing closer to a Delaunay triangulation. This allows us to optimize the surface and do the remeshing at the same time, which we think is very cool!
While implementing our algorithm, we did lots of debugging and got some funny shapes. Of course, one of the best parts about trying to make the code work was working together as a team. We had a lot of fun trying to decipher what went wrong when we got surfaces compressing to lines then to funny looking disks. Our next step is to work on coming up with algorithms for edge split and edge collapse, so we can fully implement remeshing on the shape. We are really excited to see what the final product will look like, to try lots of cool surfaces and see what the optimized remeshing for each will be!
