Mesh Curvature confusion

[RIL (Rolf)]

OK, I admit it - curvature isn’t my ballpark.

What I want to do is to highligt “sharp” convex curvature, like when a surface with thickness folds around the edges from the outside to the inside surface. Simplified illustration below (not a mesh though) shows how only the ridge along the red line should be hightlighted:

The actual mesh I’d like to apply the curvature analysis to has holes, and I would like to highligt the “ridges” of the edges of those holes, preferably without gradients on the lesser curved areas/surfaces (the selected faces):

My problem is that I have not found any single curvature method that does this highlightning consistently enough (probably due to my lack of understanding & experience of curvature analysis), but I imagine that perhaps a combination of curvature methods would isolate exactly the desired extreme curvatures?

Isolating faces (or vertices) with extreme “positive cylindric” curvature seems to be what I’m after (while “double” curvatures like hyperboloid or ellipsoid are less relevant).

I’ve tried using the MeshCurvature plugin but I have failed to get a pure and consistent highlightinbg all the edges with steep change in curvature (sharp creases aside). Example result with very “shattered” or spotted highlightning:

I want the face indicies of these edges. Example mesh attached:
Forum - Edge Curvature.gh (862.6 KB)

Advice from curvature experts welcome.

// Rolf

[clement]

Hi @RIL,

wouldn’t it be possible to iterate over the topological mesh edges and analyze the connected 2 face normals ?

_
c.

[RIL (Rolf)]

It probably would, but I don’t really know what iteration algorithms to use for Rhino meshes. Halfedge meshes seems intuitive to traverse, but Rhino meshes is more like in the wild. I guess I would have to make some structs to keep track of “visited” edges and so on, but that seem like a whole project in itself before debugged and doing the right thing. Not lazy, and I really want to become a mesh master in due time, but for now I would prefer some advice that would get me out of trouble with something that already works.

// Rolf

[Petras_Vestartas]

I would also suggest to measure angle between connected faces. But then there would be also an issue what is max and min value?I think you need to measure angle between those faces

For instance looping through topology edge, skipping all naked ones and measuring normals:

 //Get connected face at each edge
  int[] faces =   mesh.TopologyEdges.GetConnectedFaces(edgeIndex);

if(faces.Length>1){

   //measure faces normals_
   Vector3d vec0 =  mesh.FaceNormals[faces[0]];
   Vector3d vec1 =  mesh.FaceNormals[faces[0]];

}

there is also method to get mesh faces edges but I think you know that:
mesh.TopologyEdges.GetEdgesForFace(0)

Other option would be unweld mesh by some angle, but on this noisy mesh it would not be clean

[RIL (Rolf)]

@Petras_Vestartas

what is max and min value?

I don’t think that would be very problematic. It would suffice to modify settings trying until the algorithm does its thing.

I would actually scrap faces on the edge-ridges, even some faces on the lesser curved surface could go with it (easy to restore with fill holes), and in so doing the mesh would more likely separate with Split Disjoint Meshes and along creases if using Unweld/Weld tricks etc.

@Petras_Vestartas, I’ll try the face normals thing and see how far I get. Thank you for the hint.

// Rolf

[DanielPiker]

There’s something in Kangaroo for this (I needed it to get angles when creating hinges from a mesh)
You can use the ‘HingePoints’ utility to get the points of the pair of triangles around each internal edge, and then take the angle between the normals.
If you also take this angle in a plane normal to that edge, you can get positive and negative values depending whether the crease is concave or convex:

MeshAngles.gh (33.2 KB)

[clement]

@RIL

It probably would, but I don’t really know what iteration algorithms to use for Rhino meshes.

I’ve thought of something simple as below, it can be threaded easily and you can get either the faces, vertices or edge lines…

import Rhino
import rhinoscriptsyntax as rs
import math

def DoSomething():
    
    mesh_id = rs.GetObject("Mesh", 32, True, False)
    if not mesh_id: return
    
    mesh = rs.coercemesh(mesh_id, True)
    mesh.FaceNormals.ComputeFaceNormals()
    
    lines = []
    
    for i in xrange(mesh.TopologyEdges.Count):
        rc = mesh.TopologyEdges.GetConnectedFaces(i)
        if rc.Count == 2:
            n0 = mesh.FaceNormals[rc[0]]
            n1 = mesh.FaceNormals[rc[1]]
            
            if Rhino.Geometry.Vector3d.VectorAngle(n0, n1) >= math.radians(90.0):
                lines.append(mesh.TopologyEdges.EdgeLine(i))

DoSomething()

_
c.

[RIL (Rolf)]

Thank you @clement for your advice. I had already done something similar in C# but I found that there’s need for more sophisticated analysis of the adjacent faces before actually knowing how to deal with a certain extreme angle between two face normals.

But first, the code I already had (essentially the same thing):

  private void RunScript(Mesh M, double AngleMin, double AngleMax, ref object M2, ref object FI, ref object FF)
  {
    if (M == null)
      return;
    
    // Ensure existing face Normals
    var mesh = (Mesh) M.DuplicateMesh();
    mesh.Normals.ComputeNormals();

    var angle_min = AngleMin;
    var angle_max = AngleMax;
    var rad = 0.0;

    var face_indicies_ = new int[mesh.TopologyEdges.Count];
    var face_filter_ = new bool[mesh.TopologyEdges.Count];

    // connected faces per edge
    System.Threading.Tasks.Parallel.For(0, mesh.TopologyEdges.Count, e =>
      {
      int[] faces = null;
      {
        faces = mesh.TopologyEdges.GetConnectedFaces(e);

        // measure angle between faces normals, skipping 
        // naked edges (edges with less than two faces)

        if(faces.Length > 1) // non-naked
        {
          var fix = faces[0];
          var vec0 = mesh.FaceNormals[fix];
          var vec1 = mesh.FaceNormals[fix + 1];

          rad = Vector3d.VectorAngle(vec0, vec1);
          if (rad > angle_min && rad < angle_max)
          {
            face_indicies_[e] = fix;
            face_filter_[e] = true;
          }
        }
      }
      });

    // collect those face indices which was marked 
    // at the angle test above
    
    var final_face_indicies = new List<int>();
    for (var i = 0; i < face_filter_.Length; i++)
    {
      if (face_filter_[i])
        final_face_indicies.Add(face_indicies_[i]);
    }

    // ensure face indices are ordered with the greatest 
    // indexes first (reducing internal list end-to-start)
    
    var ordered = final_face_indicies.OrderByDescending(x => x);
    mesh.Faces.DeleteFaces(ordered, true);

    // Output
    M2 = mesh;
}

The problem with a simple comparison between only pairs of normals is that only from a pair of faces you don’t know on which side of an edge you’ld like to mark for highlight, deletion or for whatever. But if one knew that one of the faces is adjacent to at least two other faces, which do not form an extreme curvature, then the other face is the one you’ld like to mark as “extreme”. And so on.

So for this reason it seem like there’s need for a temp face-state class (or struct) in which it could be recorded whether a face has one or two adjacent “non-extreme” neighbours (and therefore mark it as to “armor” it from being regarded as “extreme” if also having an adjacent face with extreme angle), and only in a second pass one could then know to pick on the face on the “other side” of edges having an extreme angle between the face normals. Something like that. In other words, in only a single pass and without knowing what the “average curvature” is among the nearest faces, one can hardly know which one out of two adjacent faces with extreme angle is the deviant, so to speak.

// Rolf

[clement]

Hi Rolf, i understand the problem. Maybe this idea helps:

If you get the indices of the topological mesh vertices which are on the edge line, you could get the other vertices (or vertex) from the 2 faces connected to the edge (on the opposite site). Then from these vertices, get the connected faces and check their normal deviation for extreme angle as well. I guess to keep track of this is the hardest part as you wrote, since you are visiting faces multiple times.

I’ve done a similar script once to detect mesh errors by measuring the vertex normals to connected face normals. Great to find tiny foldovers and overlaps…
_
c.

[RIL (Rolf)]

@clement

rom these vertices, get the connected faces and check their normal deviation for extreme angle as well.

Good idea. I definitely will try something along those lines when I have the time.

// Rolf

[RIL (Rolf)]

@DanielPiker

You can use the ‘HingePoints’ utility to get the points of the pair of triangles around each internal edge, and then take the angle between the normals.

This was a neat strategy, but it took me a while to figure out what you actually said above. I had to try to illustrate the strategy to understand what you are actually saying.

Fig 1. Anyway, the middle edge of the three yellow edges in the picture is the current “hinge” edge, and the other two being used as vectors for calculating the crossproduct (normals) for the two faces on each side of the “hinge” edge (crossproduct = the two green vectors pointing downward). Finally, the plane (white) was defined by the middle “hinge” edge serving as a normal for the plane:

Cool strategy. I just wanted to share visually your strategy here. A useful and very fast way of calculating face angles & face normal directions whithout having to perform expensive lookups of faces based on indicies/indirections.

// Rolf

[DanielPiker]

Thanks Rolf, this does show it more clearly than what I’d written!

One thing I found when writing this component (and you might need to watch out for if scripting something similar) was that if you want to distinguish between a ‘mountain’ and a ‘valley’ fold on the surface, you need to know which of the 2 faces sharing an edge is on the left, and which is on the right (when looking along the direction of the edge, and with ‘up’ being outside the surface). Without this you can get the angle between the face normals, but not know whether they are folding towards or away from each other.

There is an overload that Steve added to GetConnectedFaces that helps with this:
https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_Collections_MeshTopologyEdgeList_GetConnectedFaces_1.htm
That boolean lets you know which side is which.

[RIL (Rolf)]

@DanielPiker

you can get the angle between the face normals, but not know whether they are folding towards or away from each other.

There is an overload that Steve added to GetConnectedFaces that helps with this:

Wow, that extra info will be very helpful!

// Rolf

[osuire]

Hi Daniel,

This tool is EXACTLY what I needed for the job I’m doing right now.
Many thanks.

[RIL (Rolf)]

== CODE CHALLENGE ==

@DanielPiker

You can use the ‘HingePoints’ utility to get the points of the pair of triangles around each internal edge, and then take the angle between the normals.

Sorry for bumping up this old thread, but I have use for part of your idea regarding the HingePoints you hinted about earlier (the part which extracts the edge-points which I separated into its own method “SolveHingeEdgeVertices”). However, I encountered a problem when trying to collect the points using Parallel.For (the order of the resulting points seems to get messed up causing problems later when calculating angles based on them).

If I use a regular for-loop in the code (below) it works fine to calculate angles later, but if I use Parallel.For as in the code below, then the angle-calculation fails. The Parallel.For loop runs fine (without any data races) but the order seems to matter for the later angle-calculations.

Q1: Since I need the speed of the Parallel.For loop, I wonder why it matters if the points are not collected in strict order, but since it obviously matters…
Q2: … can someone hint about a way to process the points in parallel resulting in the exact same point-order as with a regular for-loop?

Ok, that was a mouth-full but at least I tried…

  public Mesh mesh;
  public int[] HingeEdgeTVI_0;
  public int[] HingeEdgeTVI_1;
  public int[] HingeEdgeTVI_2;
  public int[] HingeEdgeTVI_3;

  public void SolveHingeEdgeVertices()
  {
    var edge_tvi_0 = new ConcurrentBag<int>();  // tvi: TopologyVertice Indices
    var edge_tvi_1 = new ConcurrentBag<int>();
    var edge_tvi_2 = new ConcurrentBag<int>();
    var edge_tvi_3 = new ConcurrentBag<int>();
    var TE = mesh.TopologyEdges;
    var TV = mesh.TopologyVertices;

    //for (int tei = 0; tei < TE.Count; tei++) // angle calculation succeeds later
    Parallel.For(0, TE.Count, tei =>           // angle calculation fails later 
      {
      int[] connectedFaces = TE.GetConnectedFaces(tei);
      if (connectedFaces.Length == 2)
      {
        IndexPair tev_pair = TE.GetTopologyVertices(tei);
        edge_tvi_0.Add(TV.MeshVertexIndices(tev_pair.I)[0]);
        edge_tvi_1.Add(TV.MeshVertexIndices(tev_pair.J)[0]);
        int[] connected_tv = TV.ConnectedTopologyVertices(tev_pair.I, true);
        for (int j = 0; j < connected_tv.Length; j++)
        {
          if (connected_tv[j] == tev_pair.J)
          {
            var conn_count = connected_tv.Length;
            int tvi_0 = connected_tv[(j - 1 + conn_count) % conn_count];
            int tvi_1 = connected_tv[(j + 1 + conn_count) % conn_count];
            edge_tvi_2.Add(TV.MeshVertexIndices(tvi_0)[0]);
            edge_tvi_3.Add(TV.MeshVertexIndices(tvi_1)[0]);
          }
        }
      }
      }); // parallel.For

    HingeEdgeTVI_0 = (int[]) edge_tvi_0.ToArray();
    HingeEdgeTVI_1 = (int[]) edge_tvi_1.ToArray();
    HingeEdgeTVI_2 = (int[]) edge_tvi_2.ToArray();
    HingeEdgeTVI_3 = (int[]) edge_tvi_3.ToArray();
    //DEBUG
    RhinoApp.Write("Made it through! {0},{1},", HingeEdgeTVI_0.Length, 
        HingeEdgeTVI_1.Length);
    RhinoApp.WriteLine("{0},{1}", HingeEdgeTVI_2.Length, 
        HingeEdgeTVI_3.Length);
  }

// Rolf

[RIL (Rolf)]

@RIL

int[] connected_tv = TV.ConnectedTopologyVertices(tev_pair.I, true); // true==sort

Hm, perhaps the sort method (“true” option above) is the reason why the point order gets messed up when using parallel? Different threads trying to sort the same vertices/edges…

@stevebaer; any chance that the sorting can be done thread-safe? (perhaps not possible, but perhaps blocking before sorting and returning copies of the vertices to avoid conflicting threads…?).

Or, making copies of the mesh and direct processing partitions to the different meshes… thinking out loud.

Any other ideas about how to speed this up while preserving the point-order?

// Rolf

[PeterFotiadis (Pfotiad0)]

@RIL

Any other ideas

Is any (Mammon oriented) thingy related with this?

Anyway: Imagine that there’s a Tree with things (named modules) and you want to do some ops using other things as “tools” (sampled in public collections) and then put the results into a public Tree (named pieces). This is a simple/entry-level thread safe way to skin the cat:

Now … replace the modules with a public MTV collection (plus paths with MTV.Count) and give it a go.

NOTE: Not always // is the fastest bunny and NOT always the fastest bunny is the // bunny.

[RIL (Rolf)]

@PeterFotiadis

Is any (Mammon oriented) thingy related with this?

I wish it was.

I found that I can sort the edges before entering the parallel for loop, and remove the sort option inside the loop, then it produces consistent results. Feels g

  TV.SortEdges();  // <-- Yay!
  Parallel.For(0, TE.Count, tei =>
  {
      //...
      int[] connected_tv = TV.ConnectedTopologyVertices(tev_pair.I /*,true*/); // NO sorting here!
      //...
  }

One heavy test case went down from 7.0 sec to 5.1 sec when going Parallel.For. Problem solved.

// Rolf
// Rolf

[PeterFotiadis (Pfotiad0)]

Well … a happy end, then … but given the opportunity: make a habit of yours: NEVER attempt a // solution without a thread safe policy - some day (or some night) you’ll try something way more complex/challenging … and then Murphy’s Law blah, blah.

[RIL (Rolf)]

@PeterFotiadis

NEVER attempt a // solution without a thread safe policy

I tend to go for a code approach which by design doesn’t provoke any conflicts. Locks are not free, not threading either, so parallel is not always faster, and with locks it’s even less likely to be faster.

The only conflict candidate in this case was the sort option. Otherwise fool proof, if using

1. only local variables inside the loop
2. Fixed size arrays (in this case ConcurrentBags, slow, so I changed to Arrays)
3. (and no dynamic stuff like the sort option).
4. Test for consistency (optional… )

// Rolf