Mask certain coordinates (geometry with hole)

[DeepaMahm]

Hi @marcomusy ,

I am trying to visualize a geometry with a hole in the center.

I could specify the coordinates of the points in a square domain.

pts = [[x, y] for x, y in zip(points_2d[0, :], points_2d[1, :])]
pcloud = Points(pts, c='k', alpha=1, r=5)
vplt.show(pcloud)

But I want to mask certain points in the domain in vedo. i.e. display a hole in the center (please see the figure below, plotted generated using a different library).

e.g. np.sqrt(points_2d[0, :] 2 + points_2d[1, :] 2) < R

I would like to plot the below in vedo to visualize the results of a simulation by mapping values to coordinates.

Could you please let me know how to mask points in the circle (to visualize a hole in the square domain)?

[marcomusy]

Hi @DeepaMahm, try

from vedo import *
points_2d = np.random.rand(10000, 2)-0.5
pcloud = Points(points_2d, c='k', r=5)
pcloud.cut_with_cylinder([0, 0, -1], [0, 0, 1], r=0.2, invert=True)
show(pcloud, axes=1)

[DeepaMahm]

Hi,

Thanks a lot, this works great.

Could you please let me know if we can get the indices of the points that are removed?

There are values associated with each of these points and I would like to color these points based on the values and visualize. But I need to remove the values corresponding to the points removed.

[marcomusy]

Hi, what about this

from vedo import *
r = 0.2
points_2d = np.random.rand(10000, 2)-0.5
pcloud1 = Points(points_2d, c='k', r=5).rotate_x(30)
pcloud2 = pcloud1.clone().cut_with_cylinder(r=r, invert=True)
cyl = Cylinder(r=r, height=1).alpha(0.2)
d = pcloud1.distance_to(cyl, signed=True)
show(pcloud2, cyl, f"inside points: {len(d[d<0])}", axes=1)

[DeepaMahm]

Hello,

Thank you for the suggestion.

I'm sorry, I think I didn't explain the doubt properly in my previous post.

pcloud2 = pcloud1.clone().cut_with_cylinder(r=r, invert=False) removes the points within the circle of radius 0.2 when invert=False.

Now I want to find the incides of points pcloud2 in pcloud1 . i.e., I'm looking for a function like mask = pcloud1.clone().find_indicies(pcloud2).

If we could get the mask_indices, we can remove those items in values. (The length of value array and pcloud1 are 1000 and the length of pcloud2 is <1000 after using cut_with_cylinder )

Currently, I get the error

[vedo.pointcloud] ERROR: in cmap(), nr. of input points scalars 1000 != 888 ...skip coloring. when I try

pts1 = pcloud2.clone().point_size(15).cmap('bwr', list(values), vmin=value_min, vmax=value_max)
    pts1.add_scalarbar3d(
        title='Values',
        title_font='VictorMono',
        label_font='VictorMono',
    )
    show(pts1, pcloud2, axes=True)

Please see the complete code below.

from vedo import *

r = 0.2

points_2d = np.random.rand(1000, 2)-0.5
values = np.random.rand(1000, 1)
value_min = min(values)
value_max = max(values)

pcloud1 = Points(points_2d, c='w', r=5)
pcloud2 = pcloud1.clone().cut_with_cylinder(r=r, invert=False)

# ------------------------------------------------------------------------------------------------------------------
# mask = pcloud1.clone().find_indicies(pcloud2)
# mask = np.ones(pts.npoints, dtype=bool)
# mask[ids] = False
# values = values[mask]

# ------------------------------------------------------------------------------------------------------------------

if False:
    cyl = Cylinder(r=r, height=1).alpha(0.2)
    d = pcloud1.distance_to(cyl, signed=True)
    show(pcloud2, cyl, f"inside points: {len(d[d<0])}", axes=1)

pts1 = pcloud1.clone().point_size(15).cmap('bwr', list(values), vmin=value_min, vmax=value_max)
pts1.add_scalarbar3d(
    title='Values',
    title_font='VictorMono',
    label_font='VictorMono',
)
show(pts1, pcloud2, axes=True)

[DeepaMahm]

If I may post another doubt,

Can we get a surface plot like the following?

I am currently coloring the points based on the values, but not sure how to interpolate and get a surface plot.

Thanks a lot

[marcomusy]

The array d<0 is the numpy array mask of booleans in my example! Consider this:

from vedo import *

settings.default_font = "VictorMono"

r = 0.2

points_2d = np.random.rand(10000, 2) - 0.5
values = np.random.rand(10000)
vmin, vmax = values.min(), values.max()

pcloud1 = Points(points_2d, c='k', r=5).rotate_x(30)
pcloud2 = pcloud1.clone().cut_with_cylinder(r=r, invert=True)

cyl = Cylinder(r=r, height=1, res=360).alpha(0.2)
dists = pcloud1.distance_to(cyl, signed=True)
mask = dists < 0
print("The boolean mask is", mask)

pcloud1.pointdata['values'] = values
pcloud1.pointdata['MASK'] = mask.astype(np.uint8)

pts1 = pcloud1.clone().point_size(5)
pts1.cut_with_scalar(0.5, 'MASK')

pts1.cmap('bwr', 'values', vmin=vmin, vmax=vmax).add_scalarbar3d(title='values')
# pts1.cmap('RdYlBu', 'MASK').add_scalarbar3d(title='MASK')
pts1.scalarbar.rotate_x(90)

show(pts1, cyl, axes=True)

then to interpolate:

# ...
# pts1.scalarbar.rotate_x(90)

grid = Grid(res=[100,100]).rotate_x(30)
grid.interpolate_data_from(pcloud1, n=3)
grid.cut_with_cylinder(r=r, invert=True)
grid.cmap('bwr', 'values', vmin=vmin, vmax=vmax).wireframe(False).lw(0)
grid.add_scalarbar3d(title='interpolated values').scalarbar.rotate_x(90)

show(cyl, grid, axes=True)

[DeepaMahm]

Awesome! thanks a ton.

Sorry, I missed this dists = pcloud1.distance_to(cyl, signed=True).

[DeepaMahm]

Hi @marcomusy , When I try to get the boundary points for this geometry with hole in the center,

I could get the boundary points on the square domain. However, I am not sure how to identify the boundary points that lie on the circle ( hole) in the center.

Could you please help?

[marcomusy]

I guess #1061 is answering this ?

[DeepaMahm]

Hello @marcomusy ,

Thank you . I tried,

from vedo import *

if __name__ == '__main__':
    npts = 10000
    length = 10
    radius = 1

    points_2d = (np.random.rand(npts, 2) - 0.5) * length  # generate random numbers between [-5, 5]]
    pcloud1 = Points(points_2d, c='k', r=5)
    pcloud2 = pcloud1.clone().cut_with_cylinder(r=radius, invert=True)

    cyl = Cylinder(pos=(0, 0, 0), r=radius, height=1, res=360).alpha(0.2)
    dists = pcloud1.distance_to(cyl, signed=True)
    mask = dists < 0
    print("The boolean mask is", mask)

    pcloud1.pointdata['MASK'] = mask.astype(np.uint8)
    pts_domain = pcloud1.clone().point_size(5)
    pts_domain.cut_with_scalar(0.5, 'MASK')

    # pts = pcloud2
    msh = pcloud2.generate_delaunay2d(mode='xy').c('w').lc('o').lw(1)
    msh.add_ids()
    msh.keep_cell_types(["triangle"])  # to remove lines

    # method 1
    msh.compute_normals().clean().linewidth(0.1)
    pids = msh.boundaries(return_point_ids=True)
    boundary_pts = Points(msh.points()[pids], r=10, c='red5')

    # method 2
    bmsh = msh.boundaries()
    lines = [b.lw(3) for b in bmsh.split()]

    show([pcloud2, lines, boundary_pts])

I am unsure how to get the boundary points on the inner circle. Could you please have a look?

Got a little busy with thesis submission :) I hope you are doing great!

[marcomusy]

Consider the following:

from vedo import *

npts = 10000
length = 10
radius = 3

# generate random numbers between [-5, 5]]
points_2d = (np.random.rand(npts, 2) - 0.5) * length 
pcloud1 = Points(points_2d, c='k', r=5)

cyl = Cylinder(r=radius, height=3, cap=False, res=360).alpha(0.2)
pcloud1.distance_to(cyl, signed=False)
pcloud1.cmap("viridis", vmin=0, vmax=radius).add_scalarbar()

pcloud2 = pcloud1.clone().threshold("Distance", below=0.05).ps(10)

show([[pcloud1, cyl], [pcloud2, cyl]], N=2, axes=1).close()