Super-super fast Point Inclusion test for Cylinder/Breps?

Hi there, this is probably similar to RadovanG’s solution above, but in VB.net.
It seems to run about 10x faster than the in-built point inclusion component, but my component also generates the point geometry. Just doing a list of true/false would be even faster.
Cheers,
Steven
170714 Point Inclusion.gh (14.6 KB)

Hi Steven,
Thanks for giving this a shot. I hacked my own mesh setup and inserted into your solution for comparison.

I noticed that the Pt inclusions aren’t exactly inside the cylinder, instead they are centered at the bottom center point of the cylinder (so that the lower half of the points are below the cylinder disc) otherwise we could have compared the numbers. Here’s a screenshot from your and my setup (the exact same settings, starting from the same mesh)

Fig.1. Your component (displaced “target” volume though) :

Fig 2. My component. Notice that I only publish the hit number counts in this component, which of course reduces time a bit (you publish one list with all the points inside). Due to the displacement we do not get the same number of inclusions (and I can’t show the hits, but I can make a new version to show that as well):


Edit: Here I added the pointlist (6848 inclusions) but that didn’t do much differrence (9ms instead of 8ms) :

// Rolf

Hi Rolf,

yep, I noticed the cylinder didn’t align. This appears to be an issue with the TryGetCylinder method in RhinoCommon. I’ll post the bug separately and see what they come up with. Would be easy enough to fix by shifting the cylinder up by half its axis length.

Regarding performance, I’ve had a bit more of a look at why my component takes to long to process all the points. I couldn’t understand why transforming all the points and checking Z height would be faster than just checking against the planes aligned with the ends of the cylinder.
It looks like just passing the 50,000 points into the component is taking 99% of the time. You’ve avoided this with your component by passing in a mesh and extracting the points inside the component.

Using the stopwatch() method my component takes 12ms to process all of the points.

Cheers,
Steven

Hi Steven, 12ms that is a very good result.

One significant difference between compiled components and scripted ones is that script components does some kind of caching of the ports, while compiled components doesn’t. Of course it takes some time to load also a compiled component, but not so much.

My own last trick to reduce the execution time of my compiled component another ~50% (from ~9 ms) was to add an extra in-port named Invalidate (Default = True). The default value means normal function, while when set to False the component skips reading the pointlist inport, reusing the list which was already loaded in the last run.

Avoiding to invalidate an indata ports - in this case the point list - is useful when changing only the other parameters like position and inclination but not the number of points (which over time usually remains the same for a mesh). With Invalidate set to False my last score is ~7ms (while changing the other parameters with sliders). One component variant which only need to present the number of points, runs at 5ms.

If the internal point list is = 0 the property reads the inport regardless of the Invalidate setting. Like so :

    Private m_Pts As New List(Of Point3d)
    Protected ReadOnly Property InPts() As List(Of Point3d)
        Get
            '// Update the Point cloud only if it hasn't been loaded yet.
            If InInvalidate OrElse (m_Pts.Count = 0) Then
                m_Pts.Clear()
                If m_Pts.Count = 0 Then
                    If Not m_DA.GetDataList(PortNo(InPrt.Pts), m_Pts) Then
                        Print("# Error - Invalid Point cloud input.")
                    End If
                End If
            End If
            Return m_Pts
        End Get
    End Property

Fig 1. OK, this one is extreme (for evaluating “score” only) but anyway (+2 ms for an output list with these points):

I don’t know if it’s possible to tweak a script component like this, but if you can, try this to see what you end up with.

// Rolf

Hey Rolf,

One significant difference between compiled components and scripted ones is that script components does some kind of caching of the ports, while compiled components doesn’t. Of course it takes some time to load also a compiled component, but not so much.

There was a discussion about it in the Grasshopper forum.
Scriptcomponents copy all the input in order to prevent weird errors in case some inexperienced scripter changes referenced values. Usually copying the input doesn’t take much time, but when doing things time critical, this becomes a problem.
However, you can hack it, and input data by reference. So instead of passing a list of points, you need to pass the component itself. Therefore you need to connect the component to an arbitrary input with tree access. Then you need to search inside the component for the connected component at the right spot and simply retrieve its data.
I made an example in c#, but you can easily translate into vb, since it is very simple.

HackToPassByRef.gh (7.7 KB)

1 Like

I did not finish to comment. Here I hopefully made it more clear.HackToPassByRef_commented.gh (10.2 KB)

1 Like

Thank you for this hint. I notice a slight difference in performance between the two… :

Wow.

BTW, due to all the ill-meant advice on this forum (…) I have started to use also C#. Hope to become fluent in C# as well. :head_bandage:

Edit:

Yes, the comments were helpful. Its easy for a newbie in the .Net world to get lost in the strange syntax. Your comments made it crystal clear.

This trick makes even the scripts components useful where performance matters and that’s worth a lot of thanks. < both thumbs up >

// Rolf

Messed up Port Index Numbers


BTW, after writing a number of components and messing up the port number indexes, I tried an Enum approach which auto-numbers the ports in a similar manner as the registering of the ports does. A property "`PortNo()`" returns the correct index number based on the auto-indexed Enum **names**.

Code example which applies for the component pictured far below. Syntax to retrive a port index based on its port name “Pts” :

 PortNo(InPrt.Pts)

When reading the inport “Pts

 Dim m_Pts As List(OF Point3d)
 DA.ReadDataList(PortNo(InPrt.Pts), m_Pts)

When writing the total number of points to the out-port “PtsInsideCnt

 DA.ReadDataList(PortNo(OutPrt.PtsInsideCnt), Cnt)

The code behind:

C#.NET

    /// <summary>
    /// Defines port names, which will be autoindexed based on the order in which 
    /// they are defined in the Enum.  Use PortNo(InPrt.EnumName) to get the name 
    /// index of an In-port automatically (prefix with OutPrt for out-ports).
    /// </summary>
    public enum InPrt
    {
        C
        Dir
        Rad
        Pts
        OnOff
        Invalidate
    }

    public enum OutPrt
    {
        PtsInsideCnt
        CylBrep
    }

    private int PortNo(InPrt EnumName)
    {
        return Convert.ToInt32(EnumName);
    }
    private int PortNo(OutPrt EnumName)
    {
        return Convert.ToInt32(EnumName);
    }

VB.NET

Same as above.

    Enum InPrt
        C
        Dir
        Rad
        Pts
        OnOff
        Invalidate
    End Enum

    Enum OutPrt
        PtsInsideCnt
        CylBrep
    End Enum

    Private Overloads Function PortNo(EnumName As InPrt) As Integer
        Return Convert.ToInt32(EnumName)
    End Function

    Private Overloads Function PortNo(EnumName As OutPrt) As Integer
        Return Convert.ToInt32(EnumName)
    End Function

This approach reduces my tendency to mess up the port numbers.

// Rolf

Well in my opinion before doing “microtuning”, I often look at the algorithm in first instance. This thread is a perfect example when it comes to performance improvements and that is why I enjoyed observing it. I could even thing of a lot more improvements when it comes to minor optimisation, however the most critical part is the way things are calculated, and not necessary how well this implementation is made. This is also why the language doesn’t matter that much, especially between VB and C#, since both are compiled into the same basis. Generally seen, I think the similarities between the languages are greater in contrast to their differences. I’m not sure if you really need to switch over to C#, but I guess the best reason to do so is because of having the most common denominator inside the Rhino Universe, which opens a whole new world to already written scripts and better understanding of Rhinocommon. It is also “easier” to understand Java and C/C++ code which again is widely spread among other platforms…

I’m actually going for both languages, since I tend to “think in VB & Delphi code”, which I then can translate to C#, if for no other reason than that over time I will at least know how to read C# code, and benefit from existing code, just like you point out.

I’m old enough though, to not throw a stick into the spokes of a spinning wheel by trying to solve new problems using a new language when I can go right at it using what I already know, that’s just making things harder than they need to be. But translating a solution to see if I can get the same (logcial) result in C# is a good way to learn, I think.

Anyway, learning C# seems to be an experince that may have been illustrated in the strip below. May I suggest learning C# is somewhere between the Infantry and the Ranger? :slight_smile:

Doesn’t feel like flying though.

// Rolf

Rolf,

It seems you’re overcomplicating the problem. A cylinder is a simple thing: it has an axis and a radius. For point inclusion you just measure the distance of the test point to the central axis of the cylinder, it should be less than the radius.

Then exclude the semi-spherical caps at the end, by checking if the closest point on the extended axis to the test point is within the finite axis-line-segment or outside it.

Doing it like this takes 5 lines of code and computes the solution for 40-60.000 point in about 400ms on my macbook.

Jack

This is the exact solution @Sweris_Waddanders suggested 8 days ago.

not to forget the huge improvement by passing the point data by reference:

cylindertest2.gh (279.7 KB)

1 Like

[quote=“jtb, post:52, topic:45714”]
It seems you’re overcomplicating the problem… Doing it like this takes 5 lines of code and computes the solution for 40-60.000 point in about 400ms on my macbook.[/quote]
OK, I’m down on 5 ms. In total. My testcase was 46830 points (see pictures, and the mesh model being used is attached in one post above, I think).

You beat 5 ms, or even come close, and then post again. I give you an extra line in addition to your initial five. :slight_smile:

My solution isn’t very complicated at all. I think if was … off the top of my head … 7 lines. Based on suggestions here in this thread. Which truly is an amazing thread.

Also, my component was made in VS, not the script component. But @TomTom has shown that it’s very possible to come close to my 5ms also with a script component. See his posts.

// Rolf

wow, 5 ms. Awesome.
Okay, its hard to tell by not having the same computer specs, because the machine used by me is actually very slow. But, just out of curiosity which algorithm did you choose?

Let me come back on the algorithm. I don’t think it was the (potentially) fastest one, but when down at 5ms it seems like I wouldn’t put more effort into exactly this one…

I’m having an old i7 3930K ona Gigabyte LGA 2011 mother board. And I installed 32GB (from 16) yesterday, so that’s after the screeenshots in this thread. System info says:

// Rolf

…I should have read the whole thread before replying…

I read your initial post and reacted a bit too quick, because the solution seemed so obvious.

I also missed that, as David Rutten pointed out, Sweris_Waddanders had already suggested the same solution.

5ms is very impressive indeed!

I read also read “fast forward” at times… :slight_smile:

BTW, notice that I have presented two different results of relevance, one is 7ms which involves the output of the points found inside the cylinder. But since my needs often are only knowing the number of points inside, I also have a component variant which presents only the number, not the point list itself. That component runs at 5ms (pictured below). This result requires the invalidation of the inport point list to be inhibited, to be “preserved” between runs, that is.

if also the point list is presented the time goes up a little, to 7ms. Add another 3 or 4 ms if also the outside points are added to an outport. Producing a bool-filter list for the Dispatch component to work with is about ~10ms, but total time is longer since the execution time of the Dispatch component has to be added in order to get the final point lists.

Fig.1. Here only the total number of inclusion points is presented in the outport, which takes 5ms. But if presenting the point list it takes 7ms in total :

(All the test runs I have presented on this thread is based on the exact same centroid point bsed on the mesh itself, a centroid point from which the same size/inclination cylinder is generated and for which point inclusion test are performed. In my case I have use for repeated recalculations due to change in the input parameters Centroid, Direction and Radius while the point cloud, the number of mesh vertices, remains constant)

// Rolf

I put together some code which:

  • is based purely on geometry and algebra

  • and does not create instances of Point3D or Vector3D types

  • uses just pure double types in computing

  • so it is supposed to be fast
    All explanation about alghorithm is embended in the code as comments.
    On my PC ( Intel® Xeon® Processor E3-1241 v3 3.50 GHz, 16GB RAM, Win10 Pro) speed of execution for 47.000 points is not mesurable in milliseconds (0 miliseconds).
    With 500.000 of points time execution becomes 5-6 miliseconds.
    Here is the code:

      public static bool[] IsPointInsideCylinder(Cylinder cylinder, Point3d[] points, double tolerance)
      {
          bool[] pointStatus = new bool[points.Length];
          //pointStatus[i] = true => points[i] is inside of cylinder
          Vector3d N = cylinder.Axis;
          N.Unitize(); // N should be already unit size vector, but for any case...
          Point3d C1 = cylinder.Center;
          //we find base of cylinder = circle 
          Circle baseCircle = cylinder.CircleAt(cylinder.Height1);
          double radiusWithTolerance = baseCircle.Radius + tolerance;
          double SquareRadiusWithTolerance = radiusWithTolerance * radiusWithTolerance;
    
          // Any point T(T.X; T.Y; T.Z) on cylinder.Axis has parametrical formula:
          //   T = C1+ t * N  where t is parameter (t is real number)
          // In that case coordiantes of T point can be writen as: 
          // (1)
          //       T.X = C1.X + t*N.X
          //       T.Y = C1.Y + t*N.Y    
          //       T.Z = C1.Z + t*N.Z
    
          // So, for point C1 parameter is t_C1 = 0
          double t_C1 = 0;
          // and for point C2 parameter is t_C2 = cylinder.TotalHeight,  Point3d C2 = C1 + t_C2 * N
          double t_C2 = cylinder.TotalHeight;
    
          // Let P(P.X; P.Y; P.Z) be any point in space
          //  and T(T.X; T.Y; T.Z) be projection of point P to axis, so T= C1 + t * N for one and just one "t".
          // Because T is projection of P on axis we know that vector P-T is orthogonal to axis,
          //  so it must be (P-T)*N = 0  (DOT PRODUCT of two orthogonal vectors is 0)
          // (2)
          //      (P-T) * N =  0   =>  we will use algebraic definition of dot product to calcualte it: 
          //        ( P.X -T.X ) * N.X  +  (P.Y - T.Y) * N.Y  +  (P.Z - P.Z) * N.Z = 0
          //   and if we substitute T.X, T.Y and T.Z in previous formula with presentation from (1) and do all calculus we get
          //      t = (N.X * (P.X - C1.X) + N.Y * (P.Y - C1.Y) + N.Z * (P.Z - C1.Z)) / (N.X*N.X + N.Y*N.Y + N.Z*N.Z)
          //   becasue N is unit vector value of (N.X*N.X + N.Y*N.Y + N.Z*N.Z) = N.SquareLength = 1
          //   so our parameter t for point T becomes:
          //  (3)
          //      t = (N.X * (P.X - C1.X) + N.Y * (P.Y - C1.Y) + N.Z * (P.Z - C1.Z))
          //
          //  And now we can use formula (1) to calucalte T.X, T.Y and T.Z  if we need it
          //
          //  So now we can calculate square distance of point P to point T by:
          //  (4)
          //      SquareDistancePT = ( P.X -T.X ) * ( P.X -T.X )  +  (P.Y - T.Y) * (P.Y - T.Y)  +  (P.Z - P.Z) * (P.Z - P.Z)
          //  
          //  Now we have everything we need to chcek if point P is inside cylinder.
          //  First condition: if distance between point P and axis is larger than Cylinder radius than point IS NOT inside cylinder,
          //                   this statement is equivalent for square distances so we can write it down as
          //      if (SquareDistancePT> SquareRadiusWithTolerance) { //point is not inside cylinder}
          //
          //   In case that SquareDistancePT<= SquareRadiusWithTolerance  point may be inside cylinder so
          //      we have to check if point T (projection of point P to axis) is inside cylinder.
          //   It is obvious that point T for which we have determined parameter t by (3) IS NOT inside cylinder if:
          //      if (t < t_C1 || t > t_C2) {//point is not inside cylinder} 
    
          //  Before we use parameters t_C1 and T_C2 in previuos formula we should add tolerance to them as follows:
          //     t_C1 = t_C1 - tolerance
          //     t_C2 = t_C2 + tolerance
          t_C1 = t_C1 - tolerance;
          t_C2 = t_C2 + tolerance;
          // So we have everything to calculate if an arbitrary point P is inside cylinder
          for (int i = 0; i < points.Length; i++)
          {
              double t = N.X * (points[i].X - C1.X)  +  N.Y * (points[i].Y - C1.Y)  +  N.Z * (points[i].Z - C1.Z);
              double TX = C1.X + t * N.X;
              double TY = C1.Y + t * N.Y;
              double TZ = C1.Z + t * N.Z;
              var SquareDistancePT = (points[i].X - TX) * (points[i].X - TX) + 
                                     (points[i].Y - TY) * (points[i].Y - TY) + 
                                     (points[i].Z - TZ) * (points[i].Z - TZ); 
    
              if (t < t_C1 || t > t_C2 || SquareDistancePT > SquareRadiusWithTolerance)
              {
                  pointStatus[i] = false; //point not inside cylinder
              }
              else
              {
                  pointStatus[i] = true;
              }
          }
          //
          return pointStatus;
      }
    

Radovan

4 Likes

Speechless. :smiley:

// Rolf