3D streamlines from custom point/vector data (Python or SPM Vector Field)

[Alejandra_Rodriguez]

Hi everyone,

Vectors.gh (187.6 KB)

I’m working with a custom list of 3D points and vectors, where each point has an associated vector. The magnitude of the vectors is important, as it represents intensity — in areas where vectors are strong, there shouldn’t be streamlines starting or ending.

My goal: Generate accurate 3D streamlines (or isolines) that closely follow the vector directions and magnitudes. The result should visualize the flow field in space.

I’ve tried two workflows:

1. Python script inside Grasshopper:

• Based on nearest neighbor search for each step.
• Follows the field but doesn’t feel accurate in direction or spacing (see image below).
• Magnitude handling is limited.

2. SPM Vector Field Integration component:

• Inputs: points as sample positions, vectors as directions, and seed points.
• I’m not sure how to get good results in 3D – I only get short lines or incorrect paths.
• Can’t figure out how to make the streamlines follow magnitude properly or generate clean 3D isolines.

What I’m looking for:

• Help on how to structure the SPM input properly for this case.
• Advice on choosing/placing seed points to get better streamlines.
• Feedback or suggestions on improving the Python method, or if another plugin/tool would be more suitable.

Any help or suggestions would be super appreciated!

(Images of Python result and SPM result attached)

[laurent_delrieu]

Hello
I am not sure on your assumption. In a vector field magnitude could give acceleration, velocity, position to a particle. But you don’t really follow magnitude, except if you have something like a Drag… Also as you choose vector to be the axis of a plane, Magnitude is always 1 !

What I can advise you is to make a Vector Field using you points and vectors.
Here I use a component from Nautilus (next release) that make an interpolation using Invers Distance Weighting

Inverse distance weighting

Then using this field you can use Legacy component from Grasshopper to draw the streamlines. I uses 2 directions of vector to have the 2 sides of streamline.

For SPM, I used to have it but I didn’t manage to install it on Rh8 !

Here the derived class I use for the field

    public sealed class InverseDistanceWeightingField : GH_FieldElement
    {
        List<Point3d> lst_pts;
        List<Vector3d> lst_vectors;
        double power;
        
        public InverseDistanceWeightingField(List<Point3d> arg_lst_pts, List<Vector3d> arg_lst_vectors, double arg_power)
        {
            lst_pts = arg_lst_pts;
            lst_vectors = arg_lst_vectors;
            power = arg_power;
            if (arg_lst_pts.Count != arg_lst_vectors.Count) throw new Exception("FieldFromVectors number of points must be the same as number of vectors");                       
        }

        public override Guid ElementGuid
        {
            get { return Guid.NewGuid(); }
        }

        public override BoundingBox BoundingBox
        {
            get { return new BoundingBox(new Point3d(-10,-10,-10), new Point3d(10,10,10)); }
        }

        public override IGH_FieldElement Duplicate()
        {
            return new InverseDistanceWeightingField(lst_pts, lst_vectors, power);
        }


        public override Vector3d Force(Point3d pt)
        {
            return LDLIB.InverseDistanceWeighting.InverseDistanceWeighting_Vectors(lst_pts, lst_vectors, power, pt);
        }
    }

In Nautilus there is a component that can clean the curves and suppress curves too close from each other.

Perhaps it could be better to have the Zero vector in the other points of the grid ! This could improve the interpolation ?

And yes with some limitation of the points for interpolations, it is better

[Alejandra_Rodriguez]

It looks interesting, can you please share the script? @laurent_delrieu

[laurent_delrieu]

Here is the script
it must work with last Nautilus release 1.9.9.7, just available on Package Manager
Vectors_LD.gh (120.3 KB)

[diff-arch]

I think you’re already on the right track with your Python script.
I’d say stick with it, instead of relying on another plugin that may or may not do what you want to in the end.

As Laurent already pointed out, you’re inputting your vectors as z-axes of planes, which by nature all have magnitude 1.0. This conflicts with your idea of using the vector magnitudes in the simulation.
On top of that, you’re normalizing the input vectors in generate_streamline() and scaling them with the step variable, before doing your semi-implicit Euler integration:

vec.Unitize()
next_pt = current_pt + vec * step

Currently, this is superfluous since you’re already dealing with normal vectors, but it’s also detrimental to the idea of using their magnitudes. You’re essentially scaling all your vectors to the same magnitude of step.
Maybe try something like this instead, once you have more appropriate data:

scale = 1.0 / step
next_pt = current_pt + vec * scale

This makes your various vector magnitudes inversely proportional to scale.

Your simulation is currently also abysmally slow. This is mainly due to the unoptimized way you do your closest neighbor/point searches:

min_dist = float("inf")
closest_vec = None
for i, p in enumerate(pts):
    dist = pt.DistanceTo(p)
    if dist < min_dist:
        min_dist = dist
        closest_vec = vecs[i]

This brute-force approach if fine for a couple of points, but slows your whole process down a lot if your dealing with a couple of hundred or even thousands points.
Generally, to speed thing up, you can use a quadtree, k-D tree, or R-tree.
Simply put, theses tree structures partition space hierarchically to quickly eliminate large regions that can’t contain the nearest point, reducing the number of comparisons needed.
In Rhino, I would use one of the rg.RTree methods for closest neighbor searches, maybe rg.RTree.Point3dKNeighbors:

closest_index = list(rg.RTree.Point3dKNeighbors(pts, [pt], 1))[0][0]
closest_pt = pts[closest_index]
closest_vec = vecs[closest_index]

Currently the Python 3 component is also considerably slower than the IronPython one.
So, if you aren’t relying on any external Python dependencies (e.g. numpy, etc.), I’d switch to the latter to shave off even a couple more seconds.

[AndersDeleuran]

@diff-arch

Generally, to speed thing up, you can use a quadtree, k-D tree, or R-tree.

The Point3dList closest point methods (1,2,3) can also be quite performant. They probably can’t match trees for large cases, but are a nice simple alternative for a general fast-ish closest point method. Edit: Here’s a quick example, had to crank up the Points count to 10000 to even trigger the profiler:

250421_ClosestPoint_GHPython_00.gh (11 KB)

[AndersDeleuran]

Out if curiosity, I just had a closer look at this: If one subtracts the cost of inputting many items (~20ms for 95000 points here) Point3dList.ClosestPointInList is the substantially faster option. Presumably building the RTree takes some time, but may be worth it in other cases (e.g. finding N closest points, multiple searches):

Also note that not using type hints is another substantial performance improvement in this case. Where the run times increase drastically by using Point3d type hints instead:

250422_ClosestPoint_GHPython_00.gh (12.3 KB)

Anywho, some food for thought

Edit: One last finding, the performance benefit of not using type hints appears to be substantially less with the new script editor CPython component:

250421_ClosestPoint_GHPython_01.gh (19.7 KB)

[diff-arch]

Interesting indeed!

The R-tree insertion usually has a time complexity of O(n). The Point3dList must also use either an R-tree or some other spacial partitioning algorithm, otherwise it couldn’t perform quick nearest neighbor searches like this (in my opinion). Maybe they’ve multithreaded parts of it or something.
I’ve looked it up in the documentation, but haven’t found further details on what’s happening under the hood.

The penality for setting a type hint, I didn’t know of. I guess it must be casting all the input data that might be Grasshopper point instances or something to RhinoCommon point instances?

[AndersDeleuran]

@diff-arch

Maybe they’ve multithreaded parts of it or something… I’ve looked it up in the documentation, but haven’t found further details on what’s happening under the hood.

Yeah it would be interesting hearing from McNeel which algorithm/implementation it uses. Its always been fast, but I don’t recall it being this performant the last time I benchmarked solutions, which admittedly is like ten years ago.

@diff-arch

I guess it must be casting all the input data that might be Grasshopper point instances or something to RhinoCommon point instances?

Indeed, the explicit casting can be costly. There’s an old (also ten years, this Tuesday!) topic that goes into details over on the old Grasshopper forum:

GHPython Some Bottlenecks (Bugs?)

[diff-arch]

In another thread a couple of weeks ago, I wrote that I would like for Grasshopper to become more integrated into Rhino, instead of keeping it in this weird plug-in limbo status (in my opinion).
Eventhough I rather meant back then that the separate window should optionally become a viewport instead (like in Houdini or Blender), but this is also a good example.
If Grasshopper and Rhino would use the same data types, within Grasshopper data could probably be passed by reference, instead of always by value, and of course also without type casting. I would imagine that a user could even set a type hint for how they want the data to be treated.

[Alejandra_Rodriguez]

Thank you for sharing your workflows — it’s definitely very helpful for this specific case, especially since the vectors are quite well organized. However, I have some doubts about how to approach the method when the vectors are not so orderly and point in different directions.

I tried splitting the vectors into vertical and horizontal groups to generate two separate sets of isolines, I also used dot product to flip the vectors but the outcome is not good, I guess there is something to do with were the lines start and end, maybe some algorithms could be included to do that, like to be continious until certain point and then when its not possible anymore to end the line. The idea is that I need 3D isolines that make sense spatially — these vectors come from a FEA analysis and represent principal stress directions. What I’m looking for are isolines of principal stress that I can use to orient an orthotropic material. Continuity is also essential.

You can think of the geometry as being composed of strips — so I can’t have the abrupt changes in direction like those shown in the attached image. That’s why I came up with the idea of filtering the vectors into groups before generating the isolines.

Basically, I’m looking for a result similar to the isolines from Karamba, but extended into 3D. Do you think this could be done with your method, and adapted to different shapes like columns and beams with varying forms?

I’ve attached an image and the script with internalized data for reference.
forum.gh (172.2 KB)

Does someone have an idea of how to solve this? @AndersDeleuran @diff-arch @laurent_delrieu
Thanks again for your insights!