Styjun

I have an irregular 3d object and want to know the surface of this object. The object can be both convex or non convex type. I can get the surface of this object applying any method like marching cube, surface contour, or isosurface.

All this methods give me triangulated mesh which is basically contains edges and vertex.

My task is to generate random and lattice points inside the object.

How should i check whether my point is inside or outside?

Any suggestion?
Thanks a lot.

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]
???????

For random points inside the closed shape:

1. Select linear density of samples


2. Make bounding box enclosing the shape


3. Select entry point on the box


4. Select exit point, compute direction cosines (w_x, w_y, w_z). Find all segments inside the shape along the ray


5. Start the ray from entry point


6. Get to first segment and and set it to p_start


7. Sample length s from exponential distribution with selected linear density


8. Find point p_end = p_start + s (w_x, w_y, w_z)


9. If it is in the first segment, store it, and make p_start = p_end. Go to step 7.


10. If it is not, go to the start of another segment, and set it to p_start. Go to step 7. If there is no segment left, you're done with one ray, go to step 3 and generate another ray.

Generate some predefined number of rays, collect all stored points, and you're done

I am sharing the code which I have written. It might be useful for others if anybody is interested for similar kind of problem. This is not the optimize code. As grid spacing value decrease computation time increase. Also depends upon the number of triangle of mesh. Any suggestion for optimizing or improve the code is welcome. Thanks

 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.