如何在不规则物体内部生成随机点和格点？

Question

我有一个不规则的 3d 对象，想知道这个对象的表面。对象可以是凸型或非凸型。我可以应用任何方法（如行进立方体、表面轮廓或等值面）来获取该对象的表面。

所有这些方法都为我提供了基本上包含边和顶点的三角网格。

我的任务是在对象内部生成随机点和格点。

我应该如何检查我的点是在内部还是外部？

有什么建议吗？ 非常感谢。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

from skimage import measure, io
from skimage.draw import ellipsoid
import skimage as sk 
import random 

I=np.zeros((50,50,50),dtype=np.float)

for i in range(50):
  for j in range(50):
    for k in range(50):
        dist=np.linalg.norm([i,j,k]-O)
        if dist<8:
            I[i,j,k]=0.8#random.random()  
        dist=np.linalg.norm([i,j,k]-O2)
        if dist<16:
            I[i,j,k]=1#random.random()  

verts, faces, normals, values = measure.marching_cubes_lewiner(I,0.7)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts[faces])
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)
plt.show()

%now forget the above code and suppose i have only verts and
%faces information. Now how to generate random points inside this Data

Data=verts[faces]
???????

Answer 1

对于封闭形状内的随机点：

1. 选择样品的线性密度
2. 制作包围形状的边界框
3. 选择框上的入口点
4. 选择出口点，计算方向余弦（w_x、w_y、w_z）。沿射线查找形状内的所有线段
5. 从入口点开始射线
6. 进入第一段并将其设置为 p_start
7. 具有选定线性密度的指数分布的样本长度 s
8. 找到点 p_end = p_start + s (w_x, w_y, w_z)
9. 如果它在第一段，则存储它，并使 p_start = p_end。转到第 7 步。
10. 如果不是，则转到另一个段的开头，并将其设置为 p_start。转至第 7 步。如果没有剩余线段，则您已完成一条射线，请转至第 3 步并生成另一条射线。

生成一些预定义数量的光线，收集所有存储的点，然后就完成了

Answer 2

我正在分享我编写的代码。如果有人对类似问题感兴趣，它可能对其他人有用。这不是优化代码。随着网格间距值的减小计算时间增加。也取决于网格的三角形数量。欢迎任何有关优化或改进代码的建议。谢谢

    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    import numpy as np 
    #from mayavi import mlab 

    verts # numpy array of vertex (triangulated mesh)
    faces # numpy array of faces (triangulated mesh) 

    %This function is taken from here 
    %https://www.erikrotteveel.com/python/three-dimensional-ray-tracing-in-python/
    def ray_intersect_triangle(p0, p1, triangle):
        # Tests if a ray starting at point p0, in the direction
        # p1 - p0, will intersect with the triangle.
        #
        # arguments:
        # p0, p1: numpy.ndarray, both with shape (3,) for x, y, z.
        # triangle: numpy.ndarray, shaped (3,3), with each row
        #     representing a vertex and three columns for x, y, z.
        #
        # returns: 
        #    0.0 if ray does not intersect triangle, 
        #    1.0 if it will intersect the triangle,
        #    2.0 if starting point lies in the triangle.

        v0, v1, v2 = triangle
        u = v1 - v0
        v = v2 - v0
        normal = np.cross(u, v)

        b = np.inner(normal, p1 - p0)
        a = np.inner(normal, v0 - p0)

        # Here is the main difference with the code in the link.
        # Instead of returning if the ray is in the plane of the 
        # triangle, we set rI, the parameter at which the ray 
        # intersects the plane of the triangle, to zero so that 
        # we can later check if the starting point of the ray
        # lies on the triangle. This is important for checking 
        # if a point is inside a polygon or not.

        if (b == 0.0):
            # ray is parallel to the plane
            if a != 0.0:
                # ray is outside but parallel to the plane
                return 0
            else:
                # ray is parallel and lies in the plane
                rI = 0.0
        else:
            rI = a / b

        if rI < 0.0:
            return 0

        w = p0 + rI * (p1 - p0) - v0

        denom = np.inner(u, v) * np.inner(u, v) - \
            np.inner(u, u) * np.inner(v, v)

        si = (np.inner(u, v) * np.inner(w, v) - \
            np.inner(v, v) * np.inner(w, u)) / denom

        if (si < 0.0) | (si > 1.0):
            return 0

        ti = (np.inner(u, v) * np.inner(w, u) - \
            np.inner(u, u) * np.inner(w, v)) / denom

        if (ti < 0.0) | (si + ti > 1.0):
            return 0

        if (rI == 0.0):
            # point 0 lies ON the triangle. If checking for 
            # point inside polygon, return 2 so that the loop 
            # over triangles can stop, because it is on the 
            # polygon, thus inside.
            return 2

        return 1


    def bounding_box_of_mesh(triangle):  
        return  [np.min(triangle[:,0]), np.max(triangle[:,0]), np.min(triangle[:,1]), np.max(triangle[:,1]), np.min(triangle[:,2]), np.max(triangle[:,2])]

    def boundingboxoftriangle(triangle,x,y,z):  
        localbox= [np.min(triangle[:,0]), np.max(triangle[:,0]), np.min(triangle[:,1]), np.max(triangle[:,1]), np.min(triangle[:,2]), np.max(triangle[:,2])]
        #print 'local', localbox
        for i in range(1,len(x)):
            if (x[i-1] <= localbox[0] < x[i]):
                x_min=i-1
            if (x[i-1] < localbox[1] <= x[i]):
                x_max=i

        for i in range(1,len(y)):
            if (y[i-1] <= localbox[2] < y[i]):
                y_min=i-1
            if (y[i-1] < localbox[3] <= y[i]):
                y_max=i

        for i in range(1,len(z)):
            if (z[i-1] <= localbox[4] < z[i]):
                z_min=i-1
            if (z[i-1] < localbox[5] <= z[i]):
                z_max=i

        return [x_min, x_max, y_min, y_max, z_min, z_max]



    spacing=5 # grid spacing 
    boundary=bounding_box_of_mesh(verts)
    print boundary 
    x=np.arange(boundary[0]-2*spacing,boundary[1]+2*spacing,spacing)
    y=np.arange(boundary[2]-2*spacing,boundary[3]+2*spacing,spacing)
    z=np.arange(boundary[4]-2*spacing,boundary[5]+2*spacing,spacing)

    Grid=np.zeros((len(x),len(y),len(z)),dtype=np.int)
    print Grid.shape

    data=verts[faces]

    xarr=[]
    yarr=[]
    zarr=[]

    # actual number of grid is very high so checking every grid is 
    # inside or outside is inefficient. So, I am looking for only 
    # those grid which is near to mesh boundary. This will reduce 
    #the time and later on internal grid can be interpolate easily.  
    for i in range(len(data)):
        #print '\n', data[i]
        AABB=boundingboxoftriangle(data[i],x,y,z)  ## axis aligned bounding box 
        #print AABB
        for gx in range(AABB[0],AABB[1]+1):
            if gx not in xarr:
                xarr.append(gx)

        for gy in range(AABB[2],AABB[3]+1):
            if gy not in yarr:
                yarr.append(gy)

        for gz in range(AABB[4],AABB[5]+1):
            if gz not in zarr:
                zarr.append(gz)


    print len(xarr),len(yarr),len(zarr)



    center=np.array([np.mean(verts[:,0]), np.mean(verts[:,1]), np.mean(verts[:,2])]) 
    print center 

    fw=open('Grid_value_output_spacing__.dat','w')

    p1=center #np.array([0,0,0])
    for i in range(len(xarr)):
        for j in range(len(yarr)):
            for k in range(len(zarr)):
                p0=np.array([x[xarr[i]],y[yarr[j]],z[zarr[k]]])
                for go in range(len(data)):
                    value=ray_intersect_triangle(p0, p1, data[go])
                    if value>0:
                        Grid[i,j,k]=value
                        break
                fw.write(str(xarr[i])+'\t'+str(yarr[j])+'\t'+str(zarr[k])+'\t'+str(x[xarr[i]])+'\t'+str(y[yarr[j]])+'\t'+str(z[zarr[k]])+'\t'+str(Grid[i,j,k])+'\n') 
        print i

    fw.close()     
    #If the grid value is greater than 0 then it is inside the triangulated mesh. 
    #I am writing the value of only confusing grid near boundary. 
    #Deeper inside grid of mesh can be interpolate easily with above information. 
    #If grid spacing is very small then generating random points inside the 
    #mesh is equivalent to choosing the random grid.