9
$\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$
12
$\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

$\endgroup$
3
$\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.

$\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
1
$\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.

$\endgroup$
0
$\begingroup$

For blender 2.8 try this

def withinMesh(x,y,z,mesh):
    axes = [ mathutils.Vector((1,0,0)) , mathutils.Vector((0,1,0)), mathutils.Vector((0,0,1))  ]
    point = mathutils.Vector((x,y,z))
    outside = False
    mat = mesh.matrix_world.copy()
    mat.invert()
    for axis in axes:
        orig = mat @ point
        count = 0
        while True:
            result,location,normal,index = mesh.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

since ray-cast as per the doc now returns 4 array values

reminder

import bpy
import mathutils

$\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.