Find the points of a ndarray that fall inside another ndarray

I have created a NumPy structured array (called arr) with arcpy module doing:

arr = arcpy.da.FeatureClassToNumPyArray('MPtest','SHAPE@XYZ',explode_to_points=True)

The array looks like (only first rows are shown):

array([([309243.1420999998, 6143470.0293, 470.43829999999434],),
       ([309252.14080000017, 6143461.0857, 470.43829999999434],),
       ([309246.0201000003, 6143459.2754, 470.43829999999434],),
       ........................................................,
       ([309252.14080000017, 6143461.0857, 464.6000000000058],)], 
      dtype=[('SHAPE@XYZ', '<f8', (3,))])

It represents XYZ coordinates taken from the vertices of a 3d object ('MPtest') built in ArcGIS (a multipatch geometry).

I have another NumPy array (generated from reading a .las file using the laspy module), called point_cloud. This array looks like:

[((30922713, 614349673, 46679, 154, 17, 1, -10, 79, 5,  11570.850892),)
 ((30922712, 614349659, 46674, 112, 17, 1, -10, 78, 5,  11570.850897),)
 ((30922709, 614349645, 46663, 161, 17, 1, -10, 77, 5,  11570.850902),),
 ..................................................................,)],
[('point', [('X', '<i4'), ('Y', '<i4'), ('Z', '<i4'), ('intensity', '<u2'), ('flag_byte', 'u1'), ('raw_classification', 'u1'), ('scan_angle_rank', 'i1'), ('user_data', 'u1'), ('pt_src_id', '<u2'), ('gps_time', '<f8')])]

I'd like to be able to possibly get the indices of the points of point cloud which fall within arr. Is this even possible?

I've been trying to play around with functions like np.where, np.intersect1d, np.logical_and, and finally np.vstack, but so far, I wasn't able to do that. Also, I have a quite robust background in Python, but still NumPy is kind of new and very complex to my eyes (at least at first sight...).

Solution

Once you get the unstructured arrays (as I see that you already achieve in comments) you can use scipy.spatial.Delanuay as follows:

I just create a sample box and some points to clarify the example:

import numpy as np
from itertools import product
from scipy.spatial import Delaunay

arr = np.array(list(product([0,1], repeat=3)))
arr
array([[0, 0, 0],
       [0, 0, 1],
       [0, 1, 0],
       [0, 1, 1],
       [1, 0, 0],
       [1, 0, 1],
       [1, 1, 0],
       [1, 1, 1]])
point_cloud = np.array([[0.5,0.5,0.5], [2,2,2]])

You then create the Delanuay Triangulation:

dt = Delaunay(arr)

And find wich points of point_cloud are inside dt (The Delanuay triangulation of arr):

points_inside_arr = dt.find_simplex(point_cloud) >=0
points_inside_arr
array([ True, False], dtype=bool)

This results in a numpy boolean array indicating which of the points in point cloud are inside arr.