Find the inner and outer contour of a set of points

[eivindtn]

Hi again @marcomusy. I have a set of points that I want to find the inner and outer contour of these points. This means the inner and outer contour of the red points below. The points order aren't organized. This is so I can compute the deviation between the outer and inner contour to the blue lines. I have tried to reconstruct a surface of the points with recoSurface(pts), and tried to find the boundaries of the mesh, but I didnt suceed since the mesh was so coarse for my points. Also creating the surface with delauney2d. I have tried to reconstruct the same issue with some code below. In advance, thank you!

from vedo import *

cyl1  = Cylinder(pos=(0,0,0), r=2, height=4, axis=(1,0,0), alpha=.5, cap=0, res=100).triangulate()
cyl2 = Cylinder(pos=(0,0,2), r=1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
cyl3 = Cylinder(pos=(0,0,2), r=1.1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
intersect_1 = cyl1.intersectWith(cyl2).join(reset=True).c('b')
intersect_2 = cyl1.intersectWith(cyl3).join(reset=True).c('b')
#show(cyl1,cyl2,cyl3, intersect_1, intersect_2).close()

#Trying to cut out the the surface between the two intersect lines
surf = cyl1.clone()
surf.cutWithMesh(cyl3, invert= True) #These two lines doesn't work for me, to cut out the section between cyl2 and cyl3? Have I done it wrong?
surf.cutWithMesh(cyl2, invert= True) # I tried using the cutWithCylinder also, instead of mesh, but the cut did end up with the same as cutWithMesh. An empty mesh.
show(surf,cyl2, cyl3).close()

#when figuring what I have done wrong with cutWithMesh, extract the points of the surf and find the countour. Maybe randomize the order of the points.
pts = v.Points(surf.points())

#find a way to find the inner and outer contour?

[marcomusy]

Hi It doesn't work because there are no vertices on the surface of the cylinder1 (only faces!). So I make i cylinder by extruding a closed line made from a circle (i might have extruded the circle but it would be capped). Check out:

from vedo import *

# cyl1 = Cylinder(pos=(0,0,0), r=2, height=4, axis=(1,0,0), alpha=.5, cap=0, res=100).triangulate() # not good
cyl1 = Line(Circle(r=2), closed=True).extrude(4, res=50).triangulate().shift(0,0,-2).rotateY(90)

cyl2 = Cylinder(pos=(0,0,2), r=1,   height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
cyl3 = Cylinder(pos=(0,0,2), r=1.1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
intersect_1 = cyl1.intersectWith(cyl2).join(reset=True).c('b').lw(3)
intersect_2 = cyl1.intersectWith(cyl3).join(reset=True).c('v').lw(3)
show(cyl1,cyl2,cyl3, intersect_1, intersect_2, axes=1).close()

#Trying to cut out the the surface between the two intersect lines
surf = cyl1.clone()
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.1, invert=False) # much faster than cutWithMesh
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.0, invert=True)
surf.subdivide(3, method=1).lw(1) # make mesh denser in points
bnds = surf.boundaries().c("red") # extract boundaries

show(surf,cyl2, cyl3, bnds, axes=1).close()

[eivindtn]

Okay! Thank you !My problem is that the extracted surface, the variable surf above. I only have the points, and need the boundaries. I'm not sure if this is right to say but I have only have the surf.points() in my issue and not having the mesh. So I need to create mesh from the points or extract the boundaries in another way? I tried to recreate the problem below using your code:

from vedo import *

# cyl1 = Cylinder(pos=(0,0,0), r=2, height=4, axis=(1,0,0), alpha=.5, cap=0, res=100).triangulate() # not good
cyl1 = Line(Circle(r=2), closed=True).extrude(4, res=50).triangulate().shift(0,0,-2).rotateY(90)

cyl2 = Cylinder(pos=(0,0,2), r=1,   height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
cyl3 = Cylinder(pos=(0,0,2), r=1.1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
intersect_1 = cyl1.intersectWith(cyl2).join(reset=True).c('b').lw(3)
intersect_2 = cyl1.intersectWith(cyl3).join(reset=True).c('v').lw(3)
#show(cyl1,cyl2,cyl3, intersect_1, intersect_2, axes=1).close()

#Trying to cut out the the surface between the two intersect lines
surf = cyl1.clone()
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.1, invert=False) # much faster than cutWithMesh
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.0, invert=True)
surf.subdivide(3, method=1).lw(1) # make mesh denser in points
bnds = surf.boundaries().c("red") # extract boundaries

#extract only the one region of surf points
pts = surf.points()
pts_2 = []
for p in pts:
    cpt = intersect_1.closestPoint(p)
    if mag2(p-cpt)<5:
        pts_2.append(p)

# I have only the points, pts_2
np.random.shuffle(pts_2) #shuffle for disorganized array
#somehow extract the boundias from pts_2?
pts_2 = Points(pts_2)
surf_2 = recoSurface(pts_2)
bnds_2  = surf_2.boundaries().c('g')
#this is the issue I have, the boundaries doesnt match with the contour on the points. Is there another way? 
show(pts_2,bnds_2, axes=1).close()

[marcomusy]

If the cloud is a generic point cloud then this becomes a hard problem with no general solution I guess. You can play with creating a "tube" around the line to erode the mesh and then use these new boundaries.

import numpy as np
from vedo import *

cyl1 = Line(Circle(r=2), closed=True).extrude(4, res=50).triangulate().shift(0,0,-2).rotateY(90)

cyl2 = Cylinder(pos=(0,0,2), r=1,   height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
cyl3 = Cylinder(pos=(0,0,2), r=1.1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
intersect_1 = cyl1.intersectWith(cyl2).join(reset=True).c('b').lw(3)
intersect_2 = cyl1.intersectWith(cyl3).join(reset=True).c('v').lw(3)

#Trying to cut out the the surface between the two intersect lines
surf = cyl1.clone()
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.1, invert=False)
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.0, invert=True)
surf.cutWithPlane(normal='z')
surf.subdivide(method=1).lw(1) # make mesh denser in points

# show(surf,cyl2, cyl3, axes=1).close()

#extract only the one region of surf points
pts = surf.points(copy=True)
np.random.shuffle(pts) #shuffle for disorganized array

#somehow extract the boundias from pts_2?
pts_2  = Points(pts)
surf_2 = recoSurface(pts, radius=0.04, dims=(200,200,50)).c('y')
# surf_2 = ConvexHull(pts).alpha(1) # this doesn't seem to work well in this case
bnds_2 = surf_2.boundaries().c('g')

modl = bnds_2.implicitModeller(distance=0.03, res=(150,150,50))
surf_2.cutWithMesh(modl)

show(pts_2, surf_2, bnds_2, axes=1).close()

Or you may try to identify the closest points to the center etc... but I think there is no "magic" solutions.

[eivindtn]

Nice ! I see ! This seems to work for my issue, after playing with some parameters. I can get the outer and inner boundaries of the points :) The code below isn't so fast, is it because of the implicitModeller?

1. Is it possible to extract the inner and the outer boundaries points without iterate with the closestpoint? I saw you did cut with a plane with the example above.
2. I want to plot the distance between the boundaries lines and the intersecion lines. So compute the distance between intersect_1 and the inner_bnds, and intersect_2 and the outer_bnds.
3. Is it possible to find point correspondece to compute the distance between the points. Cast a ray 360 degrees from the center point or something? Get average distance or an average transformation ?

Below is a started example. If you have a better suggestion or want me to look more into some of your examples! Feel free to say!

import numpy as np
from vedo import *

cyl1 = Line(Circle(r=2), closed=True).extrude(4, res=50).triangulate().shift(0,0,-2).rotateY(90)

cyl2 = Cylinder(pos=(0,0,2), r=1,   height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
cyl3 = Cylinder(pos=(0,0,2), r=1.1, height=3, axis=(0,0.3,1), alpha=.5, cap=0, res=100).triangulate()
intersect_1 = cyl1.intersectWith(cyl2).join(reset=True).c('b').lw(3)
intersect_2 = cyl1.intersectWith(cyl3).join(reset=True).c('v').lw(3)

#Trying to cut out the the surface between the two intersect lines
surf = cyl1.clone()
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.1, invert=False)
surf.cutWithCylinder(center=(0,0,2), axis=(0,0.3,1), r=1.0, invert=True)
surf.cutWithPlane(normal='z')
surf.subdivide(method=1).lw(1) # make mesh denser in points

# show(surf,cyl2, cyl3, axes=1).close()

#extract only the one region of surf points
pts = surf.points(copy=True)
np.random.shuffle(pts) #shuffle for disorganized array

#somehow extract the boundias from pts_2?
pts_2  = Points(pts)
surf_2 = recoSurface(pts, radius=0.04, dims=(200,200,50)).c('y')
# surf_2 = ConvexHull(pts).alpha(1) # this doesn't seem to work well in this case
bnds_2 = surf_2.boundaries().c('g')

modl = bnds_2.implicitModeller(distance=0.03, res=(150,150,50))
surf_2.cutWithMesh(modl)
bnds_3 = surf_2.boundaries().c('y')
bnds_3.shift(0.03, 0.02, 0.01) # to create a distance between

#extact the inner and outer boundaries, it is a better way for this?
inner_bnds = []
for p in bnds_3.points():
    cpt = intersect_1.closestPoint(p)
    if mag2(p-cpt)<0.003:
        inner_bnds.append(p)
outer_bnds = []
for p in bnds_3.points():
    cpt = intersect_2.closestPoint(p)
    if mag2(p-cpt)<0.003:
        outer_bnds.append(p)
outer_bnds = Points(outer_bnds).c('y')
inner_bnds=Points(inner_bnds).c('r')

#Smooth the points 
outer_bnds_s = outer_bnds.clone().smoothMLS1D(f=0.4).c('y')
inner_bnds_s = inner_bnds.clone().smoothMLS1D(f=0.4).c('g')

# here I want to compute a distance between the inner_bnds and intersect_1, and outer_bnds and intersect_2.
show(inner_bnds_s, outer_bnds_s, intersect_1, intersect_2, axes=1).close()

[marcomusy]

The code below isn't so fast, is it because of the implicitModeller?

Yes, and recoSurface() is also slow.

1. maybe you can create another cylinder to cut off inner/outer ring?
2. you may use distanceToMesh() (see examples/basic/distance2mesh.py)
3. i guess you have no real point to point correspondence so you need to figure out a rule that makes sense for your specific problem.

the rest of the code looks fine to me!

[eivindtn]

Thank you ! Sorry for the late reply!

1. Yes, that could be a possibility!
2. I did look into that example, and i'm adding it in for visualizing.
3. True, I don't have that. I ended up using ICP to find the deviation within transformation between the surface surf and the surf_2 points. I did the ICP with rigid = True. Is it possible to find the difference in scaling after ICP? I did compare surf.area() with the surf_2.area() . But this doesen't seem right because the recosurface() gives me a coarse mesh.

[marcomusy]

3) you can access the transformation object this way: surf_2.transform.GetMatrix() # vtkMatrix4x4

[eivindtn]

3. True, I don't have that. I ended up using ICP to find the deviation within transformation between the surface surf and the surf_2 points. I did the ICP with rigid = True. Is it possible to find the difference in scaling after ICP?

Yes, I have the transformation from surf_2.transform.GetMatrix(), but it is possible to get a measure of scaling between the two point clouds?

Also, does vedo have the possibility to convert a vtkMatrix4x4 to x,y, z, roll, pitz, yaw?

[marcomusy]

in principle you should extract the scaling from the matrix.

Also, does vedo have the possibility to convert a vtkMatrix4x4 to x,y, z, roll, pitz, yaw?

No!

[eivindtn]

in principle you should extract the scaling from the matrix.

Okey, I thought this matrix gave me the rotation and translation in terms to move the one pointcloud to another? Since the rigid = True is set.

[marcomusy]

..well it gives everything including the scaling, if you set rigid=True then this scaling should be one.

[eivindtn]

I understand. As always, gratitude for the fast response!