How to get the Intersection Surface and Vertices between Ellipsoid and Frustum and colour them?

Question

I want to get the intersection between an ellipsoid and a conical frustum. The frustum axis always lies along the line-of-sight between its small-radius basis and the ellipsoid centre. To do so I have first some utility functions

def get_r_dist(r1, div, dist):
    return (r1 + dist * div)

def point_in_conical_frustum(point, cf, direction, r1, div, h):
    v = point - cf
    proj_length = v.dot(direction)

    if proj_length < 0 or proj_length > h:
        return False

    r_dist = get_r_dist(r1, div, proj_length)
    perp_dist = (v - proj_length * direction).length
    return perp_dist <= r_dist

Then the construction of the ellipsoid and frustum

# ELLIPSOID
# =========
ce = mathutils.Vector((2.5, 2.5, 2.5))
l = mathutils.Vector((1.5, 1.0, 0.75))
nhor: int = 60
nver: int = 60
bpy.ops.mesh.primitive_uv_sphere_add(segments = nhor, ring_count = nver, radius = 1.0, location = ce)
ell = bpy.context.object
bpy.context.object.scale = (l.x, l.y, l.z)

# CONICAL FRUSTUM
# ================
cf = mathutils.Vector((0.0, 0.0, 0.0))
fru_direction = ce - cf
fru_direction.normalize()
r1: float = 0.15
h: float = (ce - cf).length
div: float = np.pi / 60.0
r2: float = r1 + h * div
bpy.ops.mesh.primitive_cone_add(vertices = 32, radius1 = r1, radius2 = r2, depth = h)
fru = bpy.context.object
up = mathutils.Vector((0.0, 0.0, 1.0))
rotation_quat = up.rotation_difference(fru_direction)
fru.rotation_mode = 'QUATERNION'
fru.rotation_quaternion = rotation_quat
fru.location = cf + fru_direction * (0.5 * fru.dimensions.z)
bpy.context.view_layer.objects.active = fru
bpy.context.scene.cursor.location = cf
bpy.ops.object.origin_set(type = 'ORIGIN_CURSOR', center = 'BOUNDS')
bpy.context.scene.cursor.location = cursor_origin

The intersection points are computed as

#INTERSECTION
#============
intersection_mat = bpy.data.materials.new(name = "IntersecMat")
intersection_mat.diffuse_color = (0, 0, 1.0, 1.0)
spheres_collection = bpy.data.collections.new("IntersecSphere")
bpy.context.scene.collection.children.link(spheres_collection)

i: int = 0
for vertex in ell.data.vertices:
    world_coords = ell.matrix_world @ vertex.co

    if point_in_conical_frustum(world_coords, cf, fru_direction, r1, div, h):
        bpy.ops.mesh.primitive_uv_sphere_add(radius = 0.015, location = world_coords)
        sphere_obj = bpy.context.object
        sphere_obj.name = f"IntersecSphere_{i}"
        spheres_collection.objects.link(sphere_obj)
        sphere_obj.data.materials.append(intersection_mat)
        i += 1

This works fine, but I am pretty sure Blender has some built-in capability to make the same thing as this script. Can you indicate me how? Moreover, I would also like to have the whole surface intersection coloured. How can I achieve this?

EDIT

I was able to get the mesh of the intersection by using this code

# Select both objects
bpy.ops.object.select_all(action = 'DESELECT')
fru.select_set(True)
ell.select_set(True)
bpy.context.view_layer.objects.active = fru

# Join the objects
bpy.ops.object.join()

# Go in 'edit' mode
bpy.ops.object.mode_set(mode = 'EDIT')
bpy.ops.mesh.select_mode(use_extend = False, use_expand = False, type = 'VERT')

# Select vertices based on coordinates
bpy.ops.mesh.select_all(action = 'DESELECT')
joined_obj = bpy.context.object
mesh_joined = joined_obj.data
bm_mesh_joined = bmesh.from_edit_mesh(mesh_joined)

tol = 1.0e-3
for vert_joined_mesh in bm_mesh_joined.verts:
    if np.abs(vert_joined_mesh.co.length - r1) <= tol:
        vert_joined_mesh.select = True

# Link vertices
bpy.ops.mesh.select_linked(delimit = set())

# Enter knife mode
bpy.ops.mesh.intersect()

# Re-select vertices based on coordinates
bpy.ops.mesh.select_all(action = 'DESELECT')
for vert_joined_mesh in bm_mesh_joined.verts:
    if np.abs(vert_joined_mesh.co.length - r1) <= tol:
        vert_joined_mesh.select = True

# Link vertices
bpy.ops.mesh.select_linked(delimit = set())

# Delete the selected vertices
bpy.ops.mesh.delete(type = 'VERT')
bmesh.from_edit_mesh(mesh_joined)
bm_mesh_joined.free()

How can I automatically (i.e. programmatically) colour the complete intersection area in, say, yellow, and get the coordinates of all the nodes inside that surface?

Answer 1

As suggested in the comments, a possible solution is to use Boolean modifiers, and the delete the right parts. To avoid any issue, first make a copy of each object

fru_copy = fru.copy()
fru_copy.data = fru_copy.data.copy()

ell_copy = ell.copy()
ell_copy.data = ell_copy.data.copy()

bpy.context.collection.objects.link(ell_copy)
bpy.context.collection.objects.link(fru_copy)

Then add the Boolean modifier and apply it

bpy.ops.object.select_all(action = 'DESELECT')
ell_copy.select_set(True)
fru_copy.select_set(True)

bpy.context.view_layer.objects.active = ell_copy
bpy.ops.object.modifier_add(type = 'BOOLEAN')
boolean_modifier = ell_copy.modifiers[-1]
boolean_modifier.operation = 'INTERSECT'
boolean_modifier.object = fru_copy

Enter in edit mode, select the vertices to delete, and delete them

bpy.context.view_layer.objects.active = ell_copy
bpy.ops.object.mode_set(mode = 'EDIT')
bpy.ops.mesh.select_all(action = 'DESELECT')
bpy.ops.mesh.select_mode(type = 'VERT')

bm = bmesh.from_edit_mesh(ell_copy.data)

tol = 1.0e-3 
normal = (obj2_copy.location - obj1_copy.location).normalized()
for vert in bm.verts:
    world_co = ell_copy.matrix_world @ vert.co
    if np.abs((world_co - ce).length - r2) <= tol:
            if np.abs(normal.dot(world_co - ce)) <= tol:
                vert.select = True

bpy.ops.mesh.delete(type = 'VERT')

Create a new material for the intersection

intersect_mat = bpy.data.materials.new(name = "IntersectMat")
intersect_mat.diffuse_color = (1.0, 1.0, 0, 1.0)
if ell_copy.data.materials:
    ell_copy.data.materials[0] = intersect_mat    
else:
    ell_copy.data.materials.append(intersect_mat)

Finally, set the intersection origin, update the mesh, and come back in object mode

bpy.ops.object.mode_set(mode = 'OBJECT')
bpy.context.view_layer.objects.active = obj2_copy
bpy.ops.object.origin_set(type = 'ORIGIN_GEOMETRY', center = 'MEDIAN')
bpy.ops.object.mode_set(mode = 'EDIT')
bmesh.update_edit_mesh(ell_copy.data)
bpy.context.view_layer.update()
bpy.data.objects.remove(fru_copy, do_unlink=True)
bpy.ops.object.mode_set(mode = 'OBJECT')
bpy.ops.mesh.select_all(action = 'DESELECT')

These steps should give the desired outcome.