How to find if a point is inside a mesh?

Question

I need a python script that would help me determine if a point is inside a cylindrical mesh.

I tried the .pointInside(point, selected_all=False) but it is deprecated in the newer versions of blender and google search yielded a couple of scripts (below) and they don't work for me.

Ray-casting method: (http://blenderartists.org/forum/showthread.php?228683-Point-in-mesh-scripts-not-working)

import bpy
import mathutils
def pointInsideMesh(point,ob):
 axes = [ mathutils.Vector((1,0,0)) ]
 outside = False
 for axis in axes:
    mat = ob.matrix_world
    mat.invert()
    orig = mat*point
    count = 0
    while True:
        location,normal,index = ob.ray_cast(orig,orig+axis*10000.0)
        if index == -1: break
        count += 1
        orig = location + axis*0.00001
    if count%2 == 0:
        outside = True
        break
return not outside
print(pointInsideMesh( mathutils.Vector((3,0,0)),bpy.context.active_object))

And one described here (http://blenderartists.org/forum/showthread.php?228683-Point-in-mesh-scripts-not-working)

## normals should be pointing out
import bpy
obj = bpy.context.object
cur = bpy.context.scene.cursor_location.copy()
cur = obj.matrix_world.inverted() * cur
cpom = obj.closest_point_on_mesh(cur)
vec = cur - cpom[0]
dot = cpom[1].dot(vec)
if dot < 0.0: print(dot, 'inside')
else: print(dot, 'outside')

Answer 1

let's look at the documentation for object.ray_cast(start, end).

The ray_cast function returns 3 values: (location, normal, index):

location, The hit location of this ray cast, float array of 3 items in [-inf, inf]

normal, The face normal at the ray cast hit location, float array of 3 items in [-inf, inf]

index, The face index, -1 when no intersection is found, int in [-inf, inf]

It will return the index of the first face encountered on the path between start and end Vectors.

• If the start vector is outside of the Object, and the face index is -1, you already know the point is not inside the object.
• But if it does return a face index, then you start counting how many consecutive faces it intersects by doing a ray_cast from the Vector of the most recent intersection (plus a small offset towards the destination to push it away from the most recent face) to the end point.
• When at some point the face index returns -1, you know there are no more faces between the checked point and the end point, then you add up the total number of intersections.
  • If that number is even, it went in and out, and is currently out.
  • If it's odd, it's still inside.

In code that might look something like this:

def is_inside(ray_origin, ray_destination, obj):

    # the matrix multiplations and inversions are only needed if you
    # have unapplied transforms, else they could be dropped. but it's handy
    # to have the algorithm take them into account, for generality.
    mat = obj.matrix_local.inverted()
    f = obj.ray_cast(mat * ray_origin, mat * ray_destination)
    loc, normal, face_idx = f

    if face_idx == -1:
        return False
    
    max_expected_intersections = 1000
    fudge_distance = 0.0001
    direction = (ray_destination - loc)
    dir_len = direction.length
    amount = fudge_distance / dir_len
    
    i = 1
    while (face_idx != -1):
        
        loc = loc.lerp(direction, amount)    
        f = obj.ray_cast(mat * loc, mat * ray_destination)
        loc, normal, face_idx = f
        print(face_idx)
        if face_idx == -1:
            break
        i += 1
        if i > max_expected_intersections:
            break

    return not ((i % 2) == 0)

Here a test blend using Sverchok Scripted Node with that algorithm.

caveat: The fudge distance is not very nicely calculated, if might help precision to repeat the algorithm from a few randomly picked points around the object, and take the most common return value.

edit: I just realized you can track the indices of intersected faces and adjust the fudge factor of the ray until the ray_cast no longer returns the index of a previously intersected face, letting it progress on..

Another approach

using obj.closest_point_on_mesh. Offered by Kosvor on sverchok issue tracker:

def is_inside(p, max_dist, obj):
    # max_dist = 1.84467e+19
    point, normal, face = obj.closest_point_on_mesh(p, max_dist)
    p2 = point-p
    v = p2.dot(normal)
    print(v)
    return not(v < 0.0)

this assumes all faces of the object are pointing outwards

Answer 2

A non node version

Here's a bmesh version I've been using for when the data isn't in a Mesh datatype.

from mathutils import Vector
from mathutils.bvhtree import BVHTree

def are_inside(points, bm):
    """
    input: 
        points
        - a list of vectors (can also be tuples/lists)
        bm
        - a manifold bmesh with verts and (edge/faces) for which the 
          normals are calculated already. (add bm.normal_update() otherwise)
    returns:
        a list
        - a mask lists with True if the point is inside the bmesh, False otherwise
    """

    rpoints = []
    addp = rpoints.append
    bvh = BVHTree.FromBMesh(bm, epsilon=0.0001)

    # return points on polygons
    for point in points:
        fco, normal, _, _ = bvh.find_nearest(point)
        p2 = fco - Vector(point)
        v = p2.dot(normal)
        addp(not v < 0.0)  # addp(v >= 0.0) ?

    return rpoints

Here I show a Vector Grid (points) and a Torus (a bmesh). The red dots are outside (False), and white dots are inside (True).

warning

• This doesn't produce desired results on low poly meshes.
• subdividing low poly meshes for the sake of inputting them into the algorithm, also doesn't work.

a better solution will come.

Answer 3

If false-inside positives are important to you, consider the following tweak of the dot-product method.

One issue with the dot-product method is that if the angle between the closest-face normal vector and point-closest-mesh-point vector is close to 90 degrees, rounding errors will result in some points that are outside of the mesh to be classified as inside.

Below method adds a tolerance (in degrees) parameter:

from math import pi, acos

def is_inside(target_pt_global, mesh_obj, tolerance=0.05):

    # Convert the point from global space to mesh local space
    target_pt_local = mesh_obj.matrix_world.inverted() * target_pt_global

    # Find the nearest point on the mesh and the nearest face normal
    _, pt_closest, face_normal, _ = mesh_obj.closest_point_on_mesh(target_pt_local)

    # Get the target-closest pt vector
    target_closest_pt_vec = (pt_closest - target_pt_local).normalized()

    # Compute the dot product = |a||b|*cos(angle)
    dot_prod = target_closest_pt_vec.dot(face_normal)

    # Get the angle between the normal and the target-closest-pt vector (from the dot prod)
    angle = acos(min(max(dot_prod, -1), 1)) * 180 / pi

    # Allow for some rounding error
    inside = angle < 90-tolerance

    return inside

The downside, however, is that due to the tolerance parameter, there might be some points that are inside the mesh that will be classified as outside (false negatives). Adjust tolerance based on your sensitivity to false-positives.

In my test with ~25K random points, tolerance of 0.02 eliminated all false-positives without any false-negatives.

Answer 4

def insideMesh(x,y,z, mesh):
    p = Vector((x,y,z))
    max_dist = 1.0e+20
    hit, point, normal, face = mesh.closest_point_on_mesh(p, max_dist)
    p2 = point-p
    v = p2.dot(normal)
    return not(v < 0.0)

E.G.

def add_cube(x,y,z,r):
    bpy.ops.mesh.primitive_cube_add(radius=r,location=(x, y, z))
    return bpy.context.active_object

#add cube with dimension = 10 on x,y,z
add_cube(0,0,0,5)
mesh =  bpy.data.objects['Cube']
print('inside: '+str(insideMesh(0,0,-5,mesh)))
#(0,0,-5) results true since mesh{x,y,z}  is bounded by  -r to r

don't forget

import bpy
from mathutils import Vector

Answer 5

i numpyfied the version of @zeffi

def numpy_apply_transforms(ob, co):
    m = np.array(ob.matrix_world)    
    mat = m[:3, :3].T 
    loc = m[:3, 3]
    return co @ mat + loc



def volume_boolean_data(vg_name, bool_name, reverse=True):

    _t = time.time()

    target_object = bpy.context.object
    bool_object   = bpy.data.objects[bool_name]

    #create mesh data
    bm = bmesh.new()
    bm.from_mesh(bool_object.data)
    bmesh.ops.transform(bm, matrix=bool_object.matrix_world, verts=bm.verts) #local to global coord
    bm.verts.ensure_lookup_table()
    bm.edges.ensure_lookup_table()
    bm.faces.ensure_lookup_table()
    bvh = BVHTree.FromBMesh(bm, epsilon=0.0001)

    _d = datetime.timedelta(seconds=time.time()-_t)
    print(f"after bmesh {_d}")

    #create array of points
    l = len(target_object.data.vertices)
    points = np.zeros((l * 3), dtype=np.float, )
    target_object.data.vertices.foreach_get('co',points)
    points.shape = (l, 3, )
    points = numpy_apply_transforms(target_object, points) #local to global coord

    _d = datetime.timedelta(seconds=time.time()-_t)
    print(f"after points apply transform {_d}")


    def check_if_inside(row):
        nonlocal bvh, reverse
        fco, normal, _, _ = bvh.find_nearest(row)
        p2 = fco - Vector(row)
        v  = p2.dot(normal)
        return v < 0.0

    def check_if_inside_reverse(row):
        nonlocal bvh, reverse
        fco, normal, _, _ = bvh.find_nearest(row)
        p2 = fco - Vector(row)
        v  = p2.dot(normal)
        return not v < 0.0

    if reverse: values = np.apply_along_axis(check_if_inside, 1,points)
    else: values = np.apply_along_axis(check_if_inside_reverse, 1,points)

    _d = datetime.timedelta(seconds=time.time()-_t)
    print(f"after along axis in {_d}")

    if vg_name not in o.vertex_groups: 
           vg = o.vertex_groups.new(name=vg_name)
    else:  vg = o.vertex_groups[vg_name]

    for i, v in enumerate(values):
        vg.add([i], v, 'REPLACE')

    _d = datetime.timedelta(seconds=time.time()-_t)
    print(f"finished in {_d}")

    return

Answer 6

def is_inside(point,ob):
    box_size = 0.15
    bpy.ops.mesh.primitive_cube_add(size=2, enter_editmode=False, align='WORLD', 
location=point, scale=(box_size,box_size,box_size))
    box=bpy.data.objects[bpy.context.active_object.name]
    try:
        bpy.ops.object.modifier_add(type='BOOLEAN')
        bpy.context.object.modifiers["Boolean"].operation = 'INTERSECT'
        bpy.context.object.modifiers["Boolean"].object =    ob
        bpy.ops.object.modifier_apply(modifier="Boolean")
        result = (box.dimensions.x >= box_size) and (box.dimensions.y >= box_size) and (box.dimensions.z >= box_size)
    return result
finally:
    bpy.data.objects.remove(box, do_unlink=True)

None of the other answers worked for me with an object that I imported from OBJ files, not sure if that was the cause. The object was plane that was then made into a mesh with solidify. The above code is my hack that works, but is slow.

I use the following code to drive the method: obj_object is the mesh.

The flipping of z with y is because I am targeting Table Top Simulator. delta = 0.5

nb_x = math.floor(obj_object.dimensions.x / delta)
nb_z = math.floor(obj_object.dimensions.z / delta)

half_x = obj_object.dimensions.x / 2
half_z = obj_object.dimensions.z / 2
x = -half_x
while x <= half_x :
    z = -half_z
    while z <= half_z :
        p = mathutils.Vector((x, z, 0))
        if is_inside(p, obj_object) :
            print("{x=%f, y=%f}," % (x,z))
        z = z + delta
    x = x + delta