12
$\begingroup$

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')
$\endgroup$
13
$\begingroup$

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.

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$
4
$\begingroup$

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

enter image description here

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ start with known point outside of mesh...ray_cast from outsize known point to test point. At each find, move epsilon....ray_cast again until no intersection. Odd number of intersections = inside, even number = outside. Can test each corner of bounding box to avoid low angle errors "slipping through." I plan to write this algorithm sometime soon :-) $\endgroup$ – patmo141 Nov 27 '17 at 2:01
3
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I am trying to figure out if a point is inside an object. It's throwing the following errors - for the first line inside the function ""Element-wise multiplication: not supported between 'Matrix' and 'Vector' types"" $\endgroup$ – Sourav Oct 13 at 16:31
  • $\begingroup$ @Sourav I can't test this right now, but try replacing '*' with '@' in the first line. It looks like a blender 2.8 issue (see blender.stackexchange.com/questions/129473/…) $\endgroup$ – Justas Oct 13 at 18:47
  • $\begingroup$ Yeah, the * has to be replaced with @. However, the script is producing tons of false positibe results for me. $\endgroup$ – Sourav Oct 13 at 18:49
1
$\begingroup$
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

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Been watching this one for the improvement.. but closest_point_on_mesh is quite expensive? If you're assuming the mesh is right-way-out and closed, won't shooting a ray down X do? $\endgroup$ – Robin Betts Jan 8 at 12:01
  • $\begingroup$ you 're right, couldn't get the coordinates to work correctly with the earlier ray cast approach so it had to be this $\endgroup$ – LiNKeR Jan 8 at 16:00
0
$\begingroup$

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

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.