Heat Method

Yes, I did have a look at that plugin, and for pure unsmoothed geodesic distance, I agree it probably does make sense to use an implicit method like the one in the CGAL library which that plugin uses.
I like to understand what’s going on and be able to change and play with everything though.
I also found its contouring component buggy unless the mesh resolution is quite fine (it seems to fail whenever there are more than 2 contour lines across any single face).

I don’t remember seeing your implementation of the Hoppe paper though, that’s nice.

There’s also some fun stuff that happens if instead of keeping the gradient fixed after the first heat distribution step, you keep alternating between the two, and recalculating the gradient as you go, so the direction is allowed to change. For this you don’t even need to specify any heat sources.
For example, try commenting out lines 100,127 and 129 in the example above and you get this:


changing how many iterations you do of each (scalar smoothing/matching scalar to gradient) before swapping also gives different results.
I was thinking it could also be interesting to try extracting the critical points of these scalar fields to form the Morse-Smale complex and use it to generate a quad mesh, like in



(the nice thing about spectral methods like this is that you can get an actual coarse quad patch layout, not just a valence-semi-regular mesh, like I wrote about here).
9 Likes

Finding critical points:
HeatMethod_Morse1.gh (156.2 KB)


13 Likes

Thanks Daniel,
This is really interesting stuff. With the updates, it almost gave me what I was looking for, using “true,false” for the two bools. The heat points are set to include every vertex on the lower naked edge of the mesh. The image below shows three examples: on the left, poisson is false, and in the center, true. I prefer the continuity of the left curves, but I need the spacing to be as consistent as possible (in reality I need to set a max threshold between any two points on successive contours to stay below).

With poisson true, the spacing is more consistent but the curves breakdown at the non-fixed boundary. If I set both internal bools to true, (and poisson true), I get the result on the right. The spacing is good near the boundaries (makes sense regarding the two bools) but the center breaks down.

I realize I am sort of applying the tool as a black box here, but its actually very close to working…any thoughts? The distance field alone may be enough for me, and I might be able to build intermediate curves in a second step with more consistent spacing, but perhaps that can just happen with some tweaks to your script?

1 Like

I think what might be needed here is to fix the values at the boundaries, but having one as a source and the other 2 as sinks. I’ll have a look.

So I tried a version of the above shape with the two “sink” arms of different length. With poisson false, I now realize that the values produced are normalized and there are the same number on each branch. It seems like I can just modify the separate contouring script to generate contours in a different way. Instead of a fixed number of contours, something that can adjust the spacing to keep the distance below a threshold. Hopefully it won’t require an iterative approach for each successive contour.

Does the fork geometry come from scanned objects?

currently no; its actually a subD model and then quadremesh.


Boundary_equalize.gh (85.2 KB)

Yes, I see what you mean about arms of different length. For this it won’t work to fix the temperature at the sources and sinks in advance, since the 2 sinks will need different temperatures.

So what I did instead was to allow specifying specific curves as equalized isolines, so the value at that boundary can change, but it stays the same along the length of that boundary.

You can get different results depending on which boundary you set as the heat source and which ones as isolines:

This could have some other uses too…

4 Likes

That certainly looks promising. It seems to have some sensitivity to the input mesh though. In the attached gh file, the first step didn’t maintain the rings around each branch / which I assume led to the isolines not equalizing (or maybe that is not related). The only difference I can see is the input mesh. Your mesh was triangular, very even edge lengths, and mine quad (which I realize internally gets converted to triangles in the first step). Or is their something else happening here?Boundary_equalize_not.gh (102.6 KB)

Hi Wes, here you go:
Boundary_equalize_release.gh (98.1 KB)

I don’t think the quality of the input mesh you are using here is any problem. The cotan weighting generally allows it to adapt to quite varied triangle shapes and sizes.

The issue instead was that for your shape, the contours of pure geodesic distance from the boundary at the base have quite different directions to the angled cuts you have on the branch ends.

If we stick strictly to geodesic distances, the contour directions are completely determined by the initial heat source. The aims of equal width strips, smoothness and alignment with all boundary curves conflict with each other - it is impossible to achieve all 3 exactly together.

So what I’ve done here is to set the branch ends as fixed value sinks for the first step, which gives us our gradient field perpendicular to the boundaries, then for the second step which equalizes distances between contours (as much as possible within these constraints), switch these boundary values to equalize along the length of each of the 3 boundary curves but not keeping them fixed. This lets the isocurves slide past the ends, so the shorter branch ends up with fewer rings than the long one to keep the distances even, but still keeping the rings aligned at the ends.

Another way could be to rearrange the script to allow ‘sinks’ of different temperatures with the possibility to change their values while iterating. This could allow you shift the location of that saddle point left or right by varying the relative values of the 2 branch ends.

Yet another way would be to change the first step completely, and instead of starting with scalar heat diffusion, work directly with the vector field, fixing directions and lengths at the boundary and smoothing vectors on the interior. Then use this vector field as the gradient for the next step. I have more thoughts on this, and something new I’m working on that I’m quite excited about, but I’ll save it for another post…

4 Likes

Thanks,
That all makes sense and it seems to work pretty cleanly now. Looking forward to hearing about the new project. In the meantime, I am working with this in a bigger project, and wondering how hard it would be to automate the convergence of each of the two steps. I remember in an earlier post you discussed using a linear solver vs being able to “see” it happen in GH. I don’t really need to see the process, just the result. Assuming the solver is a bit harder, for now I was just thinking to change the two outer for loops to while loops. It seems by comparing the value to the previous value, this gets smaller each step. I will also use a similar criteria on the second step. Any smarter way? I tested it and so far its quite fast.

Yes, iterating internally should work fine and be pretty fast.

Note that for the first heat diffusion step you don’t need to wait for the values to stop changing - all that we need from this step is the direction of the gradient, which stops changing pretty quickly even though the temperature values continue to diffuse. You probably don’t even need to worry too much about threshold criteria - just do 500 iterations and move on to the next step. For a mesh the size of the one in your example this part is very fast, especially when you don’t have to output and render any intermediate steps.

For the Poisson step you could iterate until the avg values in the
for(int i = 0;i < equalize.Length;i++)
loop are changing by less than some amount.

Though again, unless you’re really pushing to optimize the calculation time, I’d say just pick a number of iterations that works for the kind of mesh sizes you have and do that as a fixed number.
Here with 2000 iterations. It doesn’t change much beyond this if you keep iterating, and still takes less than a second on that mesh:
Boundary_equalize_release_static.gh (96.9 KB)

1 Like

Sounds good. Also want to ask about one more step. I need to a) sort the tree of contours at the end (I think I can do this brute force, but if there is some inherent info in the script that already implies this it might be faster), and b) instead of a “count” of contours, I need to generate them according to a maximum vector distance (IE no point on the next curve further than x from the prev curve). Probably easier to do this if they are sorted. I think I could do that also very brute force, by starting with a guess as to the iso step size, checking successive curve deviations, changing the step size, and generate all contours again. When doing this on Breps, its a bit more obvious to figure out where the “longest” iso curves are and divide by a step. In this case, I wondered if the original gradient you calculate could provide any info on the number of “contours” that will be formed?

Guess I should have paid more attention. It appears that they are sorted (though from cold to hot). The branches alternate, which is ideal. I think I get how it sorts, since we are counting up by “iso” in the textureCoordinates, we only create lines at the correct “contour” level. This might be the spot to test for distance compared to the previous (though will have to get rid of the parallel for I think).

Hello, many thanks for the very interesting discussion and code!

I was wondering if @DanielPiker you could give us some additional explanation on how the iterative Poisson step works on your code, and why you are using the edge gradient for it. I’m having a bit of trouble understanding its mechanism.
Could this also be done with a direct method, provided that other libraries were available?

In addition, I tried to code the heat diffusion step in python using the equalized boundaries as you describe above, but I get some local maxima appearing right below the boundaries that do not have a fixed value. I was wondering if there is some trick to avoid this. I attach a screenshot.

Many thanks!

Hi @yanna_m,

Thanks - it’s a fun topic. While this started out as an implementation of the heat method as in the paper, I think it has diverged a fair bit from it with these later examples.
If you just want to use or implement yourself exactly the method as originally described, there are better examples to look at than mine in the links on Keenan’s page.

What kind of application or further development of this do you have in mind?

Yes - a direct method would be the more usual way to solve the Poisson step, and is what I think is used in all the other implementations of this.
I just thought doing it iteratively is fun and flexible, and I think maybe will allow some interesting possibilities for further development that might not be so easy with direct methods. It’s something I’m still exploring.

The idea of my iterative method is this - once we have a scalar field with its gradient in the right directions, we need to adjust the scalar values so that the gradient magnitude becomes uniform, without changing its direction (although in some of the later examples above I was experimenting with allowing the direction to also keep changing during this part).
So it checks what the current change is between the scalar values of the vertices at the end of each edge, and the projection of the edge direction onto the gradient direction, then uses these to push the scalar values of the end vertices up or down by the amount required make the scalar difference per unit distance (in the gradient direction) some constant value (with the same value used for the whole mesh).
Of course at the first step, all the edges connected to a vertex want to change its scalar value by a different amount, but I just keep averaging and repeating, and it converges to a solution.

Regarding that artefact with the equalized boundaries - I don’t remember getting that in my version. I can’t say what the problem is without seeing the code.

1 Like

Many thanks for your explanations! I re-implemented this in python, and it works really nice (although it is of course slower).

Hello, DanielPiker, your work is brilliant, I think this branch isocurve can be used in 3d printing, but now I have a problem, after generate the isolines, which will be used for the printing path, how can I tell different branches for the print order, if you can give some advice, I will really appreciate, thank you very much.

Hello! Regarding creating such paths for 3D printing, you might be interested in this approach
https://dl.acm.org/doi/fullHtml/10.1145/3424630.3425408
Also implemented here

Thank you