I couldn't figure out why yet, but sometimes it works quite well and sometimes it's weird:

Generally, I like the way you implemented it, but have to do a more detailed review when it's working in the case above.

You can use 0.005f as default radius. We do the same in drawing code (see e.g. draw/intern/draw_curves.cc and cycles/blender/curves.cpp).

I couldn't figure out why yet, but sometimes it works quite well and sometimes it's weird

My mistake, i had tbb disabled. It seems to solve the constraints correctly but for some reason the changes are not applied to the curve positions.

• Since this can make combing significantly slower, there should be a way to turn this off. Maybe another icon next to the symmetry icons?
• I found that reducing the parallel batch size improves performance for me noticeably. I reduced it from 256 to 64 and that made combing feel much smoother in my test file with ~700 curves.
• I did manage to get curves into the surface sometimes, can you describe under what circumstances that can happen?
• Do you have thoughts on how this could be made more efficient in the future? Since most time is spend on the bvh tree lookup, there are probably two main options: Use a faster bvh tree implementation (hopefully we can use embree at some point), and, more importantly, do fewer bvh tree lookups.
• What is the expected behavior for curves that are on the inside of the mesh, and are combed to the outside?
• I find a it a bit surprising that contact points detected first. My possibly naive assumption would be that one would have to find the closest point on the surface in every iteration.
• Why is it necessary to store multiple contact points on the same ray for every point? Shouldn't it be enough to only get the closest one?

Since this can make combing significantly slower, there should be a way to turn this off. Maybe another icon next to the symmetry icons?

No objection.

I found that reducing the parallel batch size improves performance for me noticeably. I reduced it from 256 to 64 and that made combing feel much smoother in my test file with ~700 curves.

Interesting. I picked the number somewhat arbitrarily, but i suppose there is not very much overhead at all to choosing a finer grain size. Should probably also do a bunch of tests with different extremes, like lots of curves, lots of control points, high detail mesh, and see which has the biggest impact.

I did manage to get curves into the surface sometimes, can you describe under what circumstances that can happen?

Yeah i've seen that too, need to dig into this. My current working hypothesis:
The collision is detected by a raycast from the original position to the unconstraint result. This can cut through meshes, which is currently allowed (segment collision could be possible but is out of scope atm).

Do you have thoughts on how this could be made more efficient in the future? Since most time is spend on the bvh tree lookup, there are probably two main options: Use a faster bvh tree implementation (hopefully we can use embree at some point), and, more importantly, do fewer bvh tree lookups.

Yes, the collision detection is likely much more expensive than the actual constraint solving. I don't know enough about the details of the BVH trees to make concrete suggestions at this point. One thing i noticed is that the curve tools all seem to use BVHs with branching factor of 2, perhaps a larger factor could be beneficial (more children per node but shallower tree overall). It might also be possible to build an optimized raycast function to return the N closest points, as a middle ground between the closest point or all possible intersections. (see last point on why detecting N contacts is useful)

What is the expected behavior for curves that are on the inside of the mesh, and are combed to the outside?

If a curve starts inside the mesh it should not collide until it moves outside.

I find a it a bit surprising that contact points detected first. My possibly naive assumption would be that one would have to find the closest point on the surface in every iteration.

Yes, this is not ideal atm. New contacts should be detected a few times for each step as an outer loop (that is part of the issue above to detect backward faces). The inner loop is then the actual solving of constraints, using the contacts detected in the outer loop.

An implicit solver would replace the inner loop by constructing a matrix and solving all the constraints simultaneously, while Gauss-Seidel forgoes the matrix construction and pays with slower convergence/less accuracy. It solves one constraint at a time, but that can violate other constraints in turn. So the solver loops over the constraints multiple times until the overall error is small enough.

Why is it necessary to store multiple contact points on the same ray for every point? Shouldn't it be enough to only get the closest one?

This becomes important in cases like concave meshes where a single contact is not enough to move the point to the outside of the mesh. When one of the constraints is solved the other one is still violated and requires another iteration. If only one contact is stored at a time this can mean a lot of unnecessary lookups. If a reasonable set of contacts (whatever that means) can be found in one lookup then the inner loop can find a solution without having to repeat BVH lookups over and over.

I've tested setting both substeps and iterations of the solver to 1. This results in a single iteration like the previous method, and the performance is comparable.
The main problem is that the new solver works a bit differently: For any segment except the first it will move both points instead of just the 2nd point. This is a more physically accurate approach, but it requires more iterations to conserve segment length. So i would like to provide both methods and use the faster "follow-the-leader" approach by default for editing. For simulating hair physics the more balanced approach is preferable. There is a more in-depth discussion of the problem and possible mitigations for numerical stiffness in this paper https://www.researchgate.net/publication/235259170_Fast_Simulation_of_Inextensible_Hair_and_Fur

Combining the single-iteration FTL method with collisions could be difficult. We could adopt the cloth-sculpt approach and accept that collision can violate length constraints somewhat. That just requires a persistent data layer for segment lengths because we can't simply measure it at the start of the stroke.

Storing the segment lengths might be a good idea in any case to avoid drift: I would expect curves to get a little bit longer with each stroke due to numerical error when accuracy is traded for performance. I can already see this with the more detailed curves in the Einar test file: 5 solver iterations is not enough to fully converge and the hairs (~15 points each) get a bit more stretched out with each stroke. Changing the order of constraint evaluation might improve the speed of convergence too, the GS method can be susceptible to order of constraints.

I'm doing some profiling to investigate threading. One thing i noticed is that thread utilization is quite bad in case of the plain outer for loop over the TLS here. The individual thread-local blocks can have wildly different sizes, because which thread a curve gets assigned to depends on its position relative to the brush. Each thread picks a set of 256 curves at a time and filters it by distance. Some threads end up filtering out almost all curves, and some accept hundreds of curves.

This is ok for the initial combing, but doing a straight loop over these and solving each block individually causes very uneven loads. The solver in turn parallelizes over curves, but if the initial set is only a few curves to begin with most threads remain idle. Putting an outer parallel_for loop allows tbb to schedule the next block of curves earlier, but i think it would help even more to recombine the filtered curves.

EDIT
Here's a graph to compare the effect of the simple for loop (blue), the naive parallel_for loop with varying task sizes (red) and a parallel_for after combining all the thread-local changed_curves into a single array. Interestingly the naive parallel loop is actually worse than the sequential loop, but after combining arrays the thread load is much more consistent.

/* Simple sequential loop */
for (auto curves : changed_curves) {
  self_->constraint_solver_.step_curves(
      *curves_orig_, surface, transforms_, start_positions_, VArray<int>::ForSpan(curves));
};

/* Parallel loop with varying task size */
threading::parallel_for_each(changed_curves, [&](const Vector<int> &changed_curves) {
  self_->constraint_solver_.step_curves(
      *curves_orig_, surface, transforms_, start_positions_, VArray<int>::ForSpan(changed_curves));
});

/* Combine TLS curves into a single array (solver parallelizes internally) */
int totcurves = 0;
for (auto curves : changed_curves) {
  totcurves += curves.size();
}
Array<int> all_changed_curves(totcurves);
totcurves = 0;
for (auto curves : changed_curves) {
  all_changed_curves.as_mutable_span().slice(totcurves, curves.size()).copy_from(curves);
  totcurves += curves.size();
};
self_->constraint_solver_.step_curves(
    *curves_orig_, surface, transforms_, start_positions_, VArray<int>::ForSpan(all_changed_curves));

It seems like the solver could receive a single IndexMask of curves to affect. If the brush uses separate stores of changed curves per thread, it should be straightforward to convert from that to a single IndexMask. Plus it means the rest of the brush probably has similar utilization problems anyway.

BLI_index_mask_ops.hh is a somewhat relevant header.

Testing collision detection for different mesh sizes and BVH branching factors: As expected the detection time grows more or less logarithmically. The majority of time is spent in the broadphase bounds checks, with a tree depth of log_b(n), where b is the BVH branching factor. The brushes currently use a branching factor of 2. The sweet spot seems to be 4 (which makes sense for surface meshes), but the overall effect of branching factors isn't very large.

For improving collision performance i think we can only expect minor improvements by optimizing the BVH. To make collision a viable feature in editing and realtime animation there are a few options that come to mind:

1. Use simpler meshes. The "Einar" beard has about 22k curves with 440k control points and the surface mesh has 140k triangles. This kind of high detail mesh isn't needed for reasonable collision accuracy, it's common practice to have a lower-detail mesh for collisions. Big question is how that could be integrated into the curve editing workflow, without collision and visual mesh getting desynced and without too much mental overhead for artists.
2. Use simpler curves. Probably not an option for initial hair styling (we want all the detail there), but could be useful for later animation and physics. Calculating collisions on guide curves or softbody cages would be a lot cheaper than colliding every single hair, especially when self-collision is a goal.
3. Avoid the BVH: A possible option is to use a SDF volume grid, which allows very fast lookups and should provide enough accuracy (mentioned in "Fast Simulation of Inextensible Hair and Fur"). The cell size should be based on the max. collision distance per substep, much like the BVH lookup distance is now. This should be cached as much as possible.

The collision detection performance increases substantially with more substeps. The reason is that the number of near-phase triangle checks decreases with the size of the search radius, which is in turn determined by position delta per substep ("velocity"). There is a certain optimum when the search radius gets smaller than the average mesh triangle, beyond which additional substeps start increasing the overall cost.

Here i'm comparing a simple torus (20000 faces) with the "Einar" test file (combing the beard). For Einar the optimum seems to be around 60 substeps, at which point the performance becomes almost usable for editing.

The optimum depends on the "average" size of triangles on the mesh, which could make it difficult to calculate reliably. Since the cost rises only slowly with additional substeps it may be preferable to choose a default that is too large vs. one that is too small. It might be possible to adjust substeps automatically, but this could backfire as well, e.g. if triangle density is uneven. I'd like to compare this to the SDF method, which is independent from mesh topology.

FYI, I'm currently merging master in this patch and check if we can make the initial functionality master ready for Blender 3.5. I might strip it down a bit for that, while keeping the general API in place to allow for more complex/better implementations in the future.

The main thing that changed compared to when you were working on this patch originally is that we now want to support workflows where the user is just grooming guides and children hair is generated. Hence, performance is a bit less of a problem.

In D15451#469159, @Jacques Lucke (JacquesLucke) wrote:

FYI, I'm currently merging master in this patch and check if we can make the initial functionality master ready for Blender 3.5. I might strip it down a bit for that, while keeping the general API in place to allow for more complex/better implementations in the future.

Great! There is room for improvements but it should work as a first step.

The main thing that changed compared to when you were working on this patch originally is that we now want to support workflows where the user is just grooming guides and children hair is generated. Hence, performance is a bit less of a problem.

Yes, i've been thinking about formalizing the concept of "physics proxies". We can't really expect realtime collision handling with 100k+ points and dense (10k+) meshes, even if the BVH tree gets an order of magnitude faster. Guide hairs are one way to "proxify" hair physics and decouple the physics performance from visual detail. Cage meshes or "embedding" in specialized softbodies would be another.

Hi! I just wanted to throw in a comment regarding surface collision when grooming. It would be very useful to have the option of defining a custom collision collection.

The character could for instance have a coat, scarf or similar that needs to be included in the collision. This would also enable the user to voxelize/remesh the collision geometry for optimization purposes.

Another common approach is to use a "hair cap" for optimizing. A hair cap is a single sided mesh that is only used for attaching hair and is usually skinned to the armature. The head geometry can potentially be very dense (if it requires complex shape keys) so it's nice to use a hair cap, since it does not require the same polygon density.