Suppose you have a mesh with ugly triangles and you want to make it nicer. Have these beautiful equilateral triangles. Or, you want to deform a mesh by dragging on some control points. Or, maybe you want to compress your meshdata on a server, and only send parts of it, but display all of it on the client.

How would you go about achieving this?

One solution that is a basis for lots of usefull methods for working with 3d/2d/whatever meshes are least-squares meshes. Let's have a look at them!

First, let's generate us a 2d mesh:

grid = generate_grid()
grid = disturb_mesh(grid)

plot_mesh(grid)

print("Vertices: {}, Edges: {}".format(len(grid["vertices"]), len(grid["edges"])))

Vertices: 36, Edges: 85

This mesh can be represented with two matrices:

• $V$ which represents the positions of the vertices.
• $E$ which represents the edges.

$V = \begin{pmatrix} \vec{v_1} \ \vec{v_2} \ \vdots \end{pmatrix} = \begin{pmatrix} x_1 & y_1 \ x_2 & y_2 \ \vdots \end{pmatrix}$

$E = \begin{bmatrix} 0 & 0 & 1 & \dots \ 0 & 0 & 0 & \dots \ 1 & 0 & 0 & \dots \ \vdots \end{bmatrix}$

Hereby an entry $E_{ij}$ in Matrix $E$ at row $i$ and column $j$ means "Vertex $i$ is connected to vertex $j$ ".

Fairing a Mesh

Now, to actually do some mesh fairing, we need to define some anchor vertices. These are the ones we want to stay put, as much as possible. We will go outwards from them, trying to figure out the positions of all the other vertices.

These anchors can later be used as control points to deform a mesh. In this case, let's use the border vertices:

border_verts, degrees = get_border_verts(grid)
plot_mesh(grid, highlight_verts=border_verts)

Imagine the mesh getting deflated – like a net, lying on the ground. The information about connections is still there, but no information about where in space specific vertices should be. So we basically know $E$ , but have no idea about $V$ . Now w/ the anchor points, we "rehydrate" the net. We take in positional information, that then gets propagated to all the vertices that have had their positions removed.

Let's see it in action, and then think about how this was calculated.

grid0 = grid
grid1 = laplace_mesh(grid0, border_verts)
grid2 = laplace_mesh(grid1, border_verts)
grid3 = laplace_mesh(grid2, border_verts)

grids = [
    grid0,
    grid1,
    grid2,
    grid3
]

for i, agrid in enumerate(grids):
    if i > 0:
        plt.figure()
    plot_mesh(agrid, highlight_verts=border_verts)
    plt.title("Fairing-Operation applied {} times".format(i))

We can see just how quickly the algorithm distributes the vertices in a fair way. The triangles are nice to look at! (Also note that our anchors got moved, too. We will get to that a bit later).

Now to prove to you, that this indeed only needs the locations of the anchor points, let's set all the other vertices to $v_i = \begin{pmatrix}0\0 \end{pmatrix}$

grid_deflated = mesh_deflate_except(grid, border_verts)

plot_mesh(grid_deflated)

And now, let's apply the operation:

grid_inflated = laplace_mesh(grid_deflated, border_verts)

plot_mesh(grid_inflated)

Isn't that crazy? The entire internal structure of the mesh has been reconstructed, based on a few known positions. And smoothly too.

So, how did we do it?

The mechanism

First, let's think of what it means for a mesh to have a "nice", "fair" structure. If you're a 3D-modeler, you're probably thinking about uniform triangles, with similar interior angles, something like this:

This also happens to be a "nice" structure from a computational standpoint, saving you from numerical instabilities. So we'll go ahead with this.

Now, since we're working with vertex positions and edges, we need to define this fairness in terms of vertex positions (since we don't wanna change the connections). Let's derive this from a picture of the perfectly fair mesh:

Just look at it for a while, until you can formulate the fairness you see in terms of positions. How is the central vertex dependent on it's neighbours?

Did you find it?

It's really simple, and doesn't need much maths to describe it: The central vertex lies in the middle of it's neighbours. It's position is the mean of the position of its neighbours (we call these directly-connected neighbours the "1-ring neighbourhood").

So let this be our fairing condition:

$\vec{v_i} = \frac{1}{N_i} \sum^{neighbours(v_i)}{\vec{v_j}}$

with $N_i$ being the number of 1-ring neighbours of $v_i$ , and $v_j$ being vertices drawn from the 1-ring neighbourhood around $v_i$ .

To compute this efficiently, we need to rewrite this in matrix-vector form. So let's first expand this naively, and then simplify it (with $j$ being the neighbours of $v_1$ , $k$ being the neighbours of $v_2$ and so on):

$\begin{bmatrix} \vec{v_1} \ \vec{v_2} \ \vdots \end{bmatrix} = \begin{bmatrix} \frac{1}{N_1} \sum{\vec{v_j}} \ \frac{1}{N_2} \sum{\vec{v_k}} \ \vdots \end{bmatrix}$

We simplify the left side and factor out the $N$ -terms:

$V = \begin{bmatrix} \frac{1}{N_1} & 0 & \dots \ 0 & \frac{1}{N_2} & \dots \ \vdots & \vdots & \ddots \ \end{bmatrix} \begin{bmatrix} \sum{\vec{v_j}}\ \sum{\vec{v_k}}\ \vdots \end{bmatrix}$

How can we replace that ugly sum-vector on the right? Well, Matrix-Vector notation lends itself very nicely to sums. The vector specifies which elements we want to sum, and the matrix (or row-vector) defines in what combination we want to sum them:

$\begin{bmatrix} 1 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} v_{1}\ v_{2}\ v_{3}\ v_{4} \end{bmatrix} =v_{1} +v_{4}$

So we only need to know the combination matrix to get the appropriate sums... Hmmm... We already have that!

Remember at the top of the page? $E$ ? Our matrix that defines which vertex is connected to which? Exactly what we need!

So we can replace the sums with $EV$ :

$\begin{bmatrix} \sum{\vec{v_j}}\ \sum{\vec{v_k}}\ \vdots \end{bmatrix}= \underbrace{\begin{bmatrix} 0 & 1 & 1 & \dotsc \ 1 & 0 & 1 & \dotsc \ 1 & 1 & 0 & \dotsc \ \vdots & \vdots & \vdots & \ddots \end{bmatrix}}$E \underbrace{\begin{bmatrix} v\ v_{2}\ v_{3}\ \vdots \end{bmatrix}}_V

Thus we receive:

$V=\underbrace{\begin{bmatrix} \frac{1}{N_{1}} & & & \ & \frac{1}{N_{2}} & & \ & & \frac{1}{N_{3}} & \ & & & \ddots \end{bmatrix}}${N^{-1}}\underbrace{\begin{bmatrix} 0 & 1 & 1 & \dotsc \ 1 & 0 & 1 & \dotsc \ 1 & 1 & 0 & \dotsc \ \vdots & \vdots & \vdots & \ddots \end{bmatrix}} v_{1}\ v_{2}\ v_{3}\ \vdots \end{bmatrix}}_{V} }\underbrace{\begin{bmatrix

Now to solve this, we can reformulate this into a usual system of linear equations: $Ax = b$ with $A = (I - N^{-1}E)$ , $x=V$ and $b=0$ (with $I$ being the identity matrix):

$(I - N^{-1}E)V=0$

But of course this doesn't get us far, the obvious answer here is $0$ . After all, we're just solving for the position of vertices that lie in the center of their neighbours. This is certainly true if all vertices are at $(0,0)$ .

So how do we get this more interesting?

We need to add our previous knowledge about existing vertex positions. Some constraints that don't allow the simple solution of zero.

This is where our anchor points come in. These points will add prior knowledge, or wishes for positions respectively. I will name our anchors

$U\ =\ \begin{bmatrix} u_{1}\ u_{2}\ u_{3}\ \vdots \end{bmatrix}$

Let's say we want to keep the positions of vertex 1 and 2. This leads us to the following constraint:

$\begin{bmatrix} 1 & 0 & 0 & \dotsc \ 0 & 1 & 0 & \dotsc \end{bmatrix}\begin{bmatrix} v_{1}\ v_{2}\ v_{3}\ \vdots \end{bmatrix} =\begin{bmatrix} u_{1}\ u_{2}\ \end{bmatrix}$

Let's call this new matrix of constraints $F$ and the right side as $b_k$ . We can just tack it onto $A$ (and $b$ respectively) and suddenly our system solves for real values!

$A =\begin{bmatrix} A\F \end{bmatrix}, b=\begin{bmatrix} 0\ b_{k} \end{bmatrix}$

This system now has more constraints than unkowns and as such can be solved with the least square approach (Normal equation $A^TAx=A^Tb$ ). This however doesn't conserve the positions exactly, as it tries to satisfy all constraints well, loosing precision in the meantime.

Now, let's get a little more control over this. Let's add some preference for either fairness or precision: A weight that gets multiplied with the positional constraints $F$ and $b_k$ .

We'll measure the $L_2$ distance between the positions of the original mesh and the smoothed one and compute the mean over that. This way we can put a number on how well the positions are kept.

viridis = plt.get_cmap("viridis")

def f(w):
    laplaced_mesh = laplace_mesh(grid_deflated, border_verts, anchor_weight=w)

    errors = []
    for i, vertLaplace in enumerate(laplaced_mesh["vertices"]):
        if i not in border_verts:
            errors.append(0)
            continue

        vertLaplace = vertLaplace
        vertOriginal = grid_deflated["vertices"][i]

        errors.append(
            np.linalg.norm(vertLaplace - vertOriginal)
        )

    plot_mesh(
        laplaced_mesh,
        cmap=viridis,
        c=errors)
    plt.colorbar()
    plt.title("w={}, Mean Error: {:.2}".format(w, np.mean(errors)))

plot_mesh(grid)

for i in range(1,4):
    plt.figure()
    f(i)

As you can see, the more we increase the weights of the anchor points, the more precise their position is kept. At the same time, the mesh isn't faired as much, since it stays closer to the original.

This, in a nutshell is, what least squares meshes are. An approximation of some anchor vertices, combined with a propagating constraint that moves all unconstrained vertices. And this propagation condition (also the fairing condition) we used is usually called the laplace-operator and is useful for lots of things!

Now to go back to the beginning, you've seen how it can be used to fair meshes. But what about the other two? Compressing meshes and deforming meshes?

The compression works by not having to send/store the positions of all non-anchor points. There are clever algorithms to chose the best anchor points to be able to regain the original mesh as closely as possible.

The deformation works by not using existing positions as anchorpoints, but by letting the user move anchor points. As the change in position gets propagated from the anchor to the rest of the mesh, it gets deformed.

You can download the original .ipynb and the code of this post here :).

cya,

max