Tengo un conjunto de valores de datos para una función 3D escalar, organizados como entradas x,y,z en una matriz de forma (n,3) y los valores de función f(x,y,z) en una matriz de forma (n,) .
EDITAR: Por ejemplo, considere la siguiente función simple
data = np.array([np.arange(n)]*3).T F = np.linalg.norm(data,axis=1)**2Me gustaría convolucionar esta función con un núcleo esférico para realizar un suavizado 3D. La forma más fácil que encontré para realizar esto es mapear los valores de la función en una cuadrícula espacial 3D y luego aplicar una convolución 3D con el kernel que quiero.
Esto funciona bien, sin embargo, la parte que asigna la función 3D a la cuadrícula 3D es muy lenta, ya que no encontré la manera de hacerlo solo con NumPy. El siguiente código es mi implementación real, donde data es la matriz (n,3) que contiene las posiciones 3D en las que se evalúa la función, F es la matriz (n,) que contiene los valores correspondientes de la función y M es la (N,N,N) matriz que contiene la cuadrícula espacial 3D.
step = 0.1 # Create meshgrid xmin = data[:,0].min() xmax = data[:,0].max() ymin = data[:,1].min() ymax = data[:,1].max() zmin = data[:,2].min() zmax = data[:,2].max() x = np.linspace(xmin,xmax,int((xmax-xmin)/step)+1) y = np.linspace(ymin,ymax,int((ymax-ymin)/step)+1) z = np.linspace(zmin,zmax,int((zmax-zmin)/step)+1) # Build image M = np.zeros((len(x),len(y),len(z))) for l in range(len(data)): for i in range(len(x)-1): if x[i] < data[l,0] < x[i+1]: for j in range(len(y)-1): if y[j] < data[l,1] < y[j+1]: for k in range(len(z)-1): if z[k] < data[l,2] < z[k+1]: M[i,j,k] = F[l]¿Existe una forma más eficiente de llenar una cuadrícula espacial 3D con los valores de una función 3D?
Para cada elemento de datos, está escaneando píxeles de cuboide para verificar si está dentro. Hay una opción para omitir este escaneo. Podría calcular los índices correspondientes de estos píxeles usted mismo, por ejemplo:
data = np.array([[1, 2, 3], #14 (corner1) [4, 5, 6], #77 (corner2) [2.5, 3.5, 4.5], #38.75 (duplicated pixel) [2.9, 3.9, 4.9], #47.63 (duplicated pixel) [1.5, 2, 3]]) #15.25 (one step up from [1, 2, 3]) step = 0.5 data_idx = ((data - data.min(axis=0))//step).astype(int) M = np.zeros(np.max(data_idx, axis=0) + 1) x, y, z = data_idx.T M[x, y, z] = FTenga en cuenta que solo se asigna un valor de píxeles duplicados a M.
Todo lo que necesita es remodelar F[:, 3] (solo f(x, y, z)) en una cuadrícula. Es difícil ser más preciso sin datos de muestra:
Si los datos no están ordenados, debe ordenarlos:
F_sorted = F[np.lexsort((F[:,0], F[:,1], F[:,2]))] # sort by x, then y, then zElige solo f(x, y, z)
F_values = F_sorted[:, 3]Finalmente, remodele los datos en una cuadrícula:
M = F_sorted.reshape(N, N, N)Este método es más rápido que el original (aproximadamente 20 veces más rápido):
step = 0.1 mins = np.min(data, axis=0) maxs = np.max(data, axis=0) ranges = np.floor((maxs - mins) / step + 1).astype(int) indx = np.zeros(data.shape,dtype=int) for i in range(3): x = np.linspace(mins[i], maxs[i], ranges[i]) indx[:,i] = np.argmax(data[:,i,np.newaxis] <= (x[np.newaxis,:]), axis=1) -1 M = np.zeros(ranges) M[indx[:,0],indx[:,1],indx[:,2]] = F La primera parte configura las variables de cuadrícula requeridas. La función argmax proporciona una forma simple (y rápida) de encontrar el primer valor verdadero de la matriz transmitida. Esto produce un conjunto de índices para las direcciones x, y y z para cada uno de los valores de la función.
La matriz resultante M no es la misma que la producida por el código original, ya que el código original pierde datos. La lógica de y[j] < data[l,1] < y[j+1] donde y es un vector producido usando linspace significa que se perderán los valores mínimo y máximo para cada dirección ( data[l,1] podría ser igual a y[j] o y[j+1] !). Ejecútelo con un conjunto de datos de dos valores, cada uno con sus propias coordenadas y la matriz M será todo ceros.