Discrete Laplacian

Kangaroo contains a Laplacian force component, which allows for cotangent weighting The literature I’ve been looking at uses this type of operator to calculate a Laplacian on a mesh, where each vertex is assigned the value of some function. In the Kangaroo example provided by @DanielPiker , its used to average spring forces at a vertex, Is there any way to adapt this to actually inputting function values instead of spring constants?

Keenan Crane and his colleagues have done some amazing, and beautiful things using this operator and it would be great to have something like this available in GH.


Yes, it should be possible to expose the cotan weighted Laplacian for general values.
Can you give an example of how you’d want it to work? or an example application?

The formula for the Laplacian cotangent weighting of a scalar function u is

image ,where i is the vertex being evaluated, and j are the vertices immediately surrounding it. The u’s are the scalar values at each vertex.

Right now, the Kangaroo component asks for the center point i, the neighbors, j , and I assume what’s being weighted are forces coming from elsewhere.There are a few plugins like Sandbox and Parakeet that provide the necessary mesh topology. Each point would also have to be assigned a scalar value as an input. The output would be the weighted value of the u’s at each node.

I know I’m probably making this sound simpler than it is and what’s more, Crane’s method ultimately turns all these values into large-scale, sparse, matrices for easier manipulation, something GH is ill equipped for. So, maybe more trouble than its worth.

I presume you are referring to the cotan weighting option for the Laplacian smoothing force in the old version of Kangaroo. (Kangaroo2 has instead the SoapFilm goal)
MeshMachine also used cotan weighted smoothing for the relocation of the vertices.

Here’s an example showing the effect of cotangent weighting on smoothing values per vertex (shown as colours). The mesh is deliberately uneven in triangle size.

cotanweight_example.gh (28.0 KB)
On the left is uniform weighted smoothing, and on the right, cotan weighted. You can see how it diffuses more evenly despite the variation in edge lengths.
The function cotanWeights in that script takes a plankton mesh and vertex index, and returns an array of weights corresponding to the ordered neighbouring vertices.

cotanweight_example2.gh (152.2 KB)
combining it with @laurent_delrieu’s mesh iso splitting script from:


much faster version:
cotanweight_example3.gh (154.5 KB)


I’m able to view @laurent_delrieu 's original definition but in his second example, and yours I get the following errors: I do have Plankton libraries in my Grasshopper/Libraries folder.

Can you make sure the assembly reference location you have set when you right click and Manage Assemblies is the same location as your libraries folder (not anything in downloads).

I guess this is “progress”. Now I’m getting this message:

That message can be ignored. Just press the reset button to start it.

It’s alive! Thanks Daniel. Now I just have to figure out what’s happening in its innards.

I just finished a semester teaching Keenan Crane’s Discrete Differential Geometry class, and it left me wanting to implement his whole machinery in Rhino/Grasshopper. Probably the worst way to do this will be in Python, but it’s the language I’m most comfortable with. Should still be possible, since there are some Sparse matrix solvers available through the math.net framework. (See the discussion here: Math.Net function in Python with Rhino for Mac?).

Probably won’t get to this until the Fall, but I’d be happy if someone else did it before me!

Well. my Python skills are still pretty rudimentary, but if you get around to making Crane’s work accessible to Grasshopper, you’ll have done humanity a great service – and me too! BTW, what programmatic resources were you using to teach this course, if not Python?

Hi David, that would be nice - more mesh processing tools are always good!
Happy to talk more if you need any help.

@akilli - Do you have an example of an application where there’s something you want to use the Laplacian for?

By the way, for all of you into mesh operators - some of the classes and talks for this year’s SGP conference are now available. These 2 are interesting:

here’s the full graduate summer school playlist:


Here an example of the geodesic heat method from Crane,Weischedel and Wardetzky.

Flowing heat from a point source using the cotan weighted Laplacian:

Taking the gradient of this:

Finally, normalising this vector field and calculating the corresponding scalar potential, and contouring it to get curves at the same geodesic distance from the point:

No sparse solver library required, just simple iteration :slight_smile:
(certainly using a matrix solver could make it much faster, but it’s quite nice to be able to watch it doing its thing)
I’ll post the files soon.


Wow, this is simply awesome!
I have been digging into Crane’s work recently during my research about automatic knitting pattern generation. All the geodesic and vector field operations would come in VERY handy in doing this, I’m especially intrigued about the stripe pattern generation. Have been trying to understand the concepts since then and been looking for gh-compatible and open-source implementations but haven’t found anything I could work with, yet.

It would be absolutely awesome if you could share the files - I would love to actually get a glimpse of understanding about how these mathematical papers are actually implemented. I’m not a mathematician and thus I’m lacking the proper knowledge how to read certain formulas, etc. :confused: Reading code (i.e. C# or Python) would make things a lot easier (I could imagine I’m not the only one for whom this is true).

PS: I think my main problems understanding these things are what a cotan weighted laplacian actually is, how to calculate it? Aswell as how to work with vector fields on meshes, creating them, normalising them, taking gradients, creating contours. All this seems some kind of black magic to me but I’m sure it actually isn’t. If anyone could point me in a direction where I can read up about these concepts in an understandable manner I would be really, really grateful! (really!)


Did you see this article about knitting. They use one of my script for geodesic distance calculation.

My script could be replaced by this plugin
But it will be very useful to have a rhinocommon implemented geodesic distance


No, didn’t see that actually! And I thought I had found almost anything regarding that topic :wink: But there’s always some things you don’t find, I guess!

My approach is quite similar. I don’t use geodesics for creating the contours (yet) but I also stumbled about your geodesic script when I was earlier in the process (very nice by the way!). My approach is currently based on the work of Mariana Popescu. As far as I can tell, a similar logic is also used in the paper you linked.

I think the Stripe Patterns approach by Crane could be a big improvement when dealing with knitting, as the stripe patterns basically imply the short rows needed for 3d-knitting.

I hope I have not gone too far off-topic with this. If so, my apologies!

For stripes pattern you can use directional reaction diffusion using direction from geodesic iso distance.
The only problem of my reaction diffusion is that it need a good triangular mesh. Because I didn’t use the cotan weigh …

If you look in the first example file I posted above - cotanweight_example.gh and inside the function:
double[] cotanWeights(PlanktonMesh P, int i)
you can see how these weights can be calculated. It looks at a single vertex of the mesh and the ring of connected neighbouring vertices, and returns an array of scalar values, one for each of these surrounding vertices.

One simple way of thinking about it is as a way of weighting the vertices so that they properly take into account the areas of the triangles.
Say we are looking at a vertex surrounded by a ring like this:
If we have a value on each of the surrounding vertices and we take a simple unweighted average of these values (by just adding up all their values and dividing by how many of them there are), then the right side would have a much bigger influence, because there are more vertices bunched together there, even though the total area of the triangles is fairly symmetrically distributed around the central vertex.
If we are trying to average some property over the surface area of the mesh, this simple uniform average causes unwanted distortion, depending on the mesh quality, as you can see with the way the colours on spread on the left in that first animation.
So instead we can calculate weightings to multiply the surrounding vertices by the area that ‘belongs’ to that connecting edge. Using the cotangent of the tips of the 2 adjacent triangles to the edge is just a geometric way of calculating this area:

(also, if you are confused why you don’t see any actual calls to trig functions in the code of my function - it’s because the cotangent of 2 vectors can be calculated efficiently by taking their dot product divided by their cross product)