GH: Mesh Masters: How trim away inner surface & "splinters"?

RIL · June 17, 2017, 10:53pm

I’m trying to clean up meshes of bones (from CT scans) which are full of “trash” inside which disturbs my further modelling.

For example I want to get rid of double surfaces (CT-scans results in inner surfaces separating the bone tissue from the marrow) and a bunch of irregular mesh splinters on the inside of the outmost surface. I want all faces inside cleaned up (deleted), including all polygons with non-manifold edges.

Q: What would be the most straight forward method (GH components to be used) to achieve such a clean-up?

Fig 1. This is an example of what can be found inside the outmost surface (A). B and C are the double (inner) surface and some stray facets, probably caused by bone marrow.

Example mesh to be cleaned up:
Forum Example Mesh Cleanup 000.3dm (3.2 MB)

Efter clean up all the holes to be closed.

Notice that this mesh is exploded, so it would be easy to simply delete all except the top surface, but in most cases the mesh is not separated like this one, so usually the mesh is one whole and the trick is to determine which facets are being part of the external surface, and delete the rest.

// Rolf

RicardoAmaral · June 18, 2017, 10:34am

Have you tried a boolean union?

RIL · June 18, 2017, 2:20pm

I was looking for a method that can be automated, which does not require manual selection of any different parts. Also, in some cases all these faces are one single mesh so there are no parts to boolean with. For this reason I need to delete individual mesh polygons being located inside the outer-most surface.

// Rolf

Helvetosaur · June 18, 2017, 9:11pm

Dunno… SplitDisjointMesh will split non-joined parts of a single mesh into separate pieces. If the “outer” is indeed one contiguous piece, you could select it, invert the selection and Delete. SplitDisjoint is also in GH.

–Mitch

RIL · June 18, 2017, 9:28pm

Yes, I have tried that, but one problem is that in some cases that mesh has “thickness” and thus there’s a joined inner surface as well. Using SplitDisjoint I get 8 list items, and the outer most surface is the 8th item in the list and looks like this ( that is, problem remains, since there’s that inner surface still inside) :

Challenging this

// Rolf

Helvetosaur · June 18, 2017, 9:30pm

Well, then there might not be a good collection of tools native in Rhino to do this - some mesh repair software like Magics have a “shrinkwrap” function, which is more or less what you want…

–Mitch

RIL · June 18, 2017, 9:47pm

Thanks for the hint. I will look into that.

But before I abandon Rhino on this sequence in the workflow, I would like to try some more tricks in order to be able to automate this inside Rhino/GH only.

In thing I’m trying right now is to use “pipes” starting from a centroid, stretching to each point in the point cloud, and cull the plints that are found inside those pipes (the end poins for the longest pipes ( = the outmost surface) will not be found inside any of those pipes…). OK, this is an evil idea but I will try it like so:

make a “centroid” from all points in the cloud (pictured)
draw lines from the centroid to all points (um, 36 thousand of them, but here culled down to only ~2000 points).
draw flat-capped pipes around the lines.
Boolean all the pipes ~~(arghh, no such thing to be found in GH…).~~
Filter away all points inside the pipes (PointsInBrep or perhaps MeshInclusion(?)

~~OK, so how can I automate §4 and §5 above? Can I extract the meshes of the pipes and boolean the mesh-pipes instead, or…?~~

Back to the drawing board.

Edit: It seems that there is a SolidUnion that can be used on the pipes, so I will try this workflow. Execution times expected to be “significant”, though…

// Rolf

AlW · June 18, 2017, 10:31pm

I think you’re on the right track with this approach. Possibly a cylindrical approach with slices would work too. Have you thought of seeing if there is a voxel library out there that would provide the functionality you need? If you choose a small enough voxel size it should be pretty straightforward to assign each point to a voxel. I’m assuming your issue is a common requirement in voxel image processing and probably well implemented in a library or two.

There’s also Paraview, but I’m not sure whether it has a function to do what you need, and no-one has written a 3dm reader for it. I have the impression from your postings here that you’re a pretty savvy programmer so perhaps a look at Paraview’s open source code would be of some help. Or the Kitware Visualization Toolkit (VTK) which Paraview uses.

RIL · June 18, 2017, 10:44pm

I’m not familiar with the voxel processing concept. But I guess I will have to look it up then…

It seems that the pipe->PointInBrep can work, but it takes aweful long time to process the union of the pipes. We’ll see where I end up with this.

But thanks for all your hints. I’m sure there is a solution for this, and all your hints are very helpful for a geometry-newbie like me.

Regarding my programming skills, well, it happens every now and then that some piece of code surprises me, but I would not over-estimate myself in the field of geometry… (or in any other field). And even if ikt would just so happen that I would get somthing to actually work, it needs to have some performance too.

In any case, I will look into every suggestion to find the best solution.

// Rolf

RIL · June 18, 2017, 11:52pm

I can confirm that it works, although the processing is slow.

In the following picture I reduced the number of pipes ( = reducing processing time for testing) and one can see from the picture that many points were “trapped” inside the pipes. I also unioned a Sphere (around the Centroid) in order to simplify the Boolean union of the Pipes converging at the Centroid (pipes being “cut off” by the sphere before being too much entangled with each other) :

Limitations
The Pipe radius needs to be a bit smaller than the distance between the points on the surface so that not also nearby surface points are included in the pipes.

A weakness with this method is that the inclination of the pipes near the “walls” at the the bottom part of the mesh) tends to include points located on the walls very closely above the end points of the pipes.

Due to the inclination of the pipes, half of some of the pipes (near the bottom edge) extends “through” the walls and thus traps legit points as well.

Reducing the pipe radius reduces the problem, and a small enough radius will elliminate the problem altogether, but a smaller radius will at the same time be more likely to miss points between the pipes. Points that should be removed.

Smarter people than me may see potential for enhancing this basic approach, reducing both the problem and the processing time(?)

Enhancement to try
Btw, while writing this I just came up with the idea to use two different radiuses for the pipes (using “Variable Pipe (VPipe)”. At the surfaces ends, the radius can be small enough to not trap nearby surface-points, small enough to not trap them even when the pieps are “leaning up agains the walls”. And the other end of the pipes could have the radius of the inner sphere. In this way each pipe will occupy much more volume and thus trap more “stray” points and…

Well, I’m sure of that more enhancements can be done. Performance is still a problem though.

// Rolf

AlW · June 19, 2017, 1:13pm

Voxel => volume pixel. ie, 3d pixel. Used in areas like MRI and CAT scan processing and visualization as well as a lot of other scientific 3d imaging.

Since you piqued my interest, I did a bit of internet browsing and came across this:

They’ve been bought by Autodesk, so will probably go the way of T-Splines, but for now it or the Grasshopper plugin version might help you out.

@bobmcneel, @brian: How did you let this get away from you?

RIL · June 19, 2017, 4:08pm

Thank you very much for this hint! I downloaded a copy. Will evaluate it in a few days (hoping for some rainy days… )

Before giving up on a native GH solution, is there a way to close/cap holes of a mesh in GH so that the MeshInclusion starts to work? I tested MeshInclusion on a GH generated/manipulated mesh Def (which reports no naked edges) but not one point is found inside although there are over 20 thousand of the,…).

Disappointing, to say the least.

What would a basic voxelation strategy look like in this case? (I’ve never been thinking along those lines, or cubes, so I would need a little push in the right direction to get going.

// Rolf

AlW · June 19, 2017, 8:03pm

Just because I know this stuff exists and have seen it used doesn’t mean I know anything about how to program it!

My brief review of Monoliths gallery videos and a very quick scan of their user manual indicated that it can:

read and write Rhino .3dm files
voxelize meshes (and maybe point clouds)
probably do what you want

A couple of other things to look intoner GH might be Cocoon, Aether, and Plankton library.

emilio · June 20, 2017, 9:32am

Hi Rolf

… Talking about performance, I tried your algorithm (nice idea BTW ) with some changes that … might make it faster … not sure though.

Here is the idea: splitting the point cloud into smaller ‘volumes’ and then apply your ‘pipes subtraction’ independently to them.
In order to make it faster, for any volume we only use a single pipe axis direction,
and for easier calculation, we orient the volume so that the pipes axes are vertical.
we will also use ‘square pipes’ … but I think cylindrical pipes could be used too … not sure about the difference in performance.

TAI.3dm (1.0 MB)

Here is, more or less, how it should work:

move the points to have the centroid in the origin
splt the point into 24 sets ( see the ‘walls’ that divide the volumes )
for each point set:
3A) orient it so that its ‘central axis’ lie on the Z axis
3B) purge the points by the ‘square pipes’
3C) orient the remaining point back to their original orientation

Unfortunately my Grasshopper is very far from fluent … so I’ll stick to plain Python scripts …

This is the test script I used to … build the idea

import Rhino
import rhinoscriptsyntax as rs
import scriptcontext

def xforms( xside, yside, zside, subindex, xfor, xbak ):
  '''
  xside, yside, zside : 1.0 OR -1.0
  subindex in ( 0, 1, 2 )
    calc transform, transform-back
    update xfor, xbak lists
  '''
  subside = [ xside, yside, zside ]
  vec0 = ( Rhino.Geometry.Vector3d.XAxis * xside +
      Rhino.Geometry.Vector3d.YAxis * yside +
      Rhino.Geometry.Vector3d.ZAxis * zside ) / 3.0
  vec0.Unitize()
  vec1 = [ Rhino.Geometry.Vector3d.XAxis * xside,
      Rhino.Geometry.Vector3d.YAxis * yside,
      Rhino.Geometry.Vector3d.ZAxis * zside ][ subindex ]
  vec2 = vec0 + vec1
  pla = Rhino.Geometry.Plane.WorldXY
  plb = Rhino.Geometry.Plane( Rhino.Geometry.Point3d.Origin, vec2 )
  xfa = Rhino.Geometry.Transform.PlaneToPlane( plb, pla )
  xfb = Rhino.Geometry.Transform.PlaneToPlane( pla, plb )
  if not xfa.IsValid:
    print 'XFA is not valid'
  if not xfb.IsValid:
    print 'XFB is not valid'
  xfor.append( xfa )
  xbak.append( xfb )

def purge( pset, xfor, xbak, side, deltaz ):
# ( unable to make this work )
# pset = list( xfor.TransformList( pset ) )
  for ix in range( len( pset ) ):
    pset[ ix ].Transform( xfor )
# rs.ObjectColor( rs.AddPoints( pset ), [ 255, 0, 0 ] )
  pset.sort( key = lambda x: x.Z, reverse = True )
# for ix in range( len( pset ) ):
#   print str( pset[ ix ] )
  hsid = side * 0.5
  for ix in range( len( pset ) ):
    pt = pset[ ix ]
    if not pt:
      continue
    xmi = pt.X - hsid
    xma = pt.X + hsid
    ymi = pt.Y - hsid
    yma = pt.Y + hsid
    zma = pt.Z - deltaz
    for iy in range( len( pset ) ):
      pp = pset[ iy ]
      if not pp:
        continue
      if ( ( pp.X > xmi ) and ( pp.X < xma ) and
          ( pp.Y > ymi ) and ( pp.Y < yma ) and
          ( pp.Z < zma ) ):
        pset[ iy ] = None
  res = []
  for pt in pset:
    if pt:
      res.append( pt )
  for ix in range( len( res ) ):
    res[ ix ].Transform( xbak )
  return res

def main():
# get points
  gids = rs.GetObjects( 'Points ?', filter = 1, preselect = True )
  pts = [ rs.PointCoordinates( gid ) for gid in gids ]
# get centroid
  cen = Rhino.Geometry.Point3d.Origin
  for pt in pts:
    cen += pt
  cen /= len( pts )
# move to origin
  for ix in range( len( pts ) ):
    pts[ ix ] -= cen
# split into 24 point sets
  # points
  pset = []
  for ix in range( 24 ):
    pset.append( [] )
  # transforms
  xfor = []
  # back transforms
  xbak = []
# build transforms
# 0 1 2
  xforms( 1.0, 1.0, 1.0, 0, xfor, xbak )
  xforms( 1.0, 1.0, 1.0, 1, xfor, xbak )
  xforms( 1.0, 1.0, 1.0, 2, xfor, xbak )
# 3 4 5
  xforms( -1.0, 1.0, 1.0, 0, xfor, xbak )
  xforms( -1.0, 1.0, 1.0, 1, xfor, xbak )
  xforms( -1.0, 1.0, 1.0, 2, xfor, xbak )
# 6 7 8
  xforms( -1.0, -1.0, 1.0, 0, xfor, xbak )
  xforms( -1.0, -1.0, 1.0, 1, xfor, xbak )
  xforms( -1.0, -1.0, 1.0, 2, xfor, xbak )
# 9 10 11
  xforms( 1.0, -1.0, 1.0, 0, xfor, xbak )
  xforms( 1.0, -1.0, 1.0, 1, xfor, xbak )
  xforms( 1.0, -1.0, 1.0, 2, xfor, xbak )
# 12 13 14
  xforms( 1.0, 1.0, -1.0, 0, xfor, xbak )
  xforms( 1.0, 1.0, -1.0, 1, xfor, xbak )
  xforms( 1.0, 1.0, -1.0, 2, xfor, xbak )
# 15 16 17
  xforms( -1.0, 1.0, -1.0, 0, xfor, xbak )
  xforms( -1.0, 1.0, -1.0, 1, xfor, xbak )
  xforms( -1.0, 1.0, -1.0, 2, xfor, xbak )
# 18 19 20
  xforms( -1.0, -1.0, -1.0, 0, xfor, xbak )
  xforms( -1.0, -1.0, -1.0, 1, xfor, xbak )
  xforms( -1.0, -1.0, -1.0, 2, xfor, xbak )
# 21 22 23
  xforms( 1.0, -1.0, -1.0, 0, xfor, xbak )
  xforms( 1.0, -1.0, -1.0, 1, xfor, xbak )
  xforms( 1.0, -1.0, -1.0, 2, xfor, xbak )
# build point sets
  for pt in pts:
    ax = abs( pt.X )
    ay = abs( pt.Y )
    az = abs( pt.Z )
# 0 1 2
    if ( pt.X > 0.0 ) and ( pt.Y > 0.0 ) and ( pt.Z > 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 0 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 0 + 1 ].append( pt )
      else:
        pset[ 0 + 2 ].append( pt )
# 3 4 5
    elif ( pt.X <= 0.0 ) and ( pt.Y > 0.0 ) and ( pt.Z > 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 3 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 3 + 1 ].append( pt )
      else:
        pset[ 3 + 2 ].append( pt )
# 6 7 8
    elif ( pt.X <= 0.0 ) and ( pt.Y <= 0.0 ) and ( pt.Z > 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 6 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 6 + 1 ].append( pt )
      else:
        pset[ 6 + 2 ].append( pt )
# 9 10 11
    elif ( pt.X > 0.0 ) and ( pt.Y <= 0.0 ) and ( pt.Z > 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 9 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 9 + 1 ].append( pt )
      else:
        pset[ 9 + 2 ].append( pt )
# 12 13 14
    elif ( pt.X > 0.0 ) and ( pt.Y > 0.0 ) and ( pt.Z <= 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 12 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 12 + 1 ].append( pt )
      else:
        pset[ 12 + 2 ].append( pt )
# 15 16 17
    elif ( pt.X <= 0.0 ) and ( pt.Y > 0.0 ) and ( pt.Z <= 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 15 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 15 + 1 ].append( pt )
      else:
        pset[ 15 + 2 ].append( pt )
# 18 19 20
    elif ( pt.X <= 0.0 ) and ( pt.Y <= 0.0 ) and ( pt.Z <= 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 18 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 18 + 1 ].append( pt )
      else:
        pset[ 18 + 2 ].append( pt )
# 21 22 23
    elif ( pt.X > 0.0 ) and ( pt.Y <= 0.0 ) and ( pt.Z <= 0.0 ):
      if ( ax > ay ) and ( ax > az ):
        pset[ 21 + 0 ].append( pt )
      elif ( ay > ax ) and ( ay > az ):
        pset[ 21 + 1 ].append( pt )
      else:
        pset[ 21 + 2 ].append( pt )
# debug
# for ix in range( 24 ):
#   print len( pset[ ix ] )
# rs.ObjectColor( rs.AddPoints( pset[ 0 ] ), [ 255, 0, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 1 ] ), [ 0, 255, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 2 ] ), [ 0, 0, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 3 ] ), [ 255, 255, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 4 ] ), [ 255, 0, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 5 ] ), [ 0, 255, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 6 ] ), [ 255, 0, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 7 ] ), [ 255, 127, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 8 ] ), [ 255, 127, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 9 ] ), [ 127, 0, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 10 ] ), [ 127, 127, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 11 ] ), [ 127, 0, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 12 ] ), [ 127, 255, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 13 ] ), [ 0, 255, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 14 ] ), [ 127, 255, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 15 ] ), [ 0, 127, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 16 ] ), [ 127, 127, 0 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 17 ] ), [ 0, 127, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 18 ] ), [ 127, 0, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 19 ] ), [ 0, 127, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 20 ] ), [ 127, 127, 255 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 21 ] ), [ 0, 0, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 22 ] ), [ 127, 0, 127 ] )
# rs.ObjectColor( rs.AddPoints( pset[ 23 ] ), [ 0, 127, 127 ] )

  ppts = []
  for ix in range( len( pset ) ):
    ppts.extend( purge( pset[ ix ], xfor[ ix ], xbak[ ix ], 15.0, 5.0 ) )
    print( '%d / 24' % ix )

# move points back
  for ix in range( len( ppts ) ):
    ppts[ ix ] += cen
# delete original points
# rs.DeleteObjects( gids )
# draw purged points
  rs.AddPoints( ppts )

main()

… for what it’s worth …

Regards

RIL · June 21, 2017, 1:32pm

Back from a trip to the capital, spent hours in traffic jams and… ah, never min d.

Thank you @AIW for the additional hints. I will definietly look into every thint. And @emilio … that was a nasty one! I’ll look into your approach as well! Impressing, as usual!

// Rolf