Rellenando huecos en una matriz numpy

Solo quiero interpolar, en los términos más simples posibles, un conjunto de datos 3D. Interpolación lineal, vecino más cercano, todo eso sería suficiente (esto es para comenzar un algoritmo, por lo que no se requiere una estimación precisa).

En las nuevas versiones de scipy, cosas como griddata serían útiles, pero actualmente solo tengo scipy 0.8. Así que tengo una matriz de “cubo” ( data[:,:,:] , (NixNjxNk)), y una matriz de indicadores ( flags[:,:,:,] , True o False ) del mismo tamaño. Quiero interpolar mis datos para los elementos de datos donde el elemento de bandera correspondiente es Falso, utilizando, por ejemplo, el punto de datos válido más cercano en los datos, o alguna combinación lineal de puntos “cercanos a”.

Puede haber grandes espacios en el conjunto de datos en al menos dos dimensiones. Aparte de codificar un algoritmo de vecino más cercano en toda regla utilizando kdtrees o similar, realmente no puedo encontrar un interpolador genérico de N-dimensional vecino más cercano.

Related of "Rellenando huecos en una matriz numpy"

Puede configurar un algoritmo de estilo de crecimiento de cristal desplazando una vista alternativamente a lo largo de cada eje, reemplazando solo los datos marcados con False pero con un True vecino. Esto da un resultado similar al “vecino más cercano” (pero no a una distancia euclidiana o de Manhattan; creo que podría ser el vecino más cercano si está contando píxeles, contando todos los píxeles de conexión con esquinas comunes) Esto debería ser bastante eficiente con NumPy ya que itera solo sobre iteraciones de convergencia y eje, no pequeños segmentos de los datos.

Crudo, rápido y estable. Creo que eso era lo que estabas buscando:

 import numpy as np # -- setup -- shape = (10,10,10) dim = len(shape) data = np.random.random(shape) flag = np.zeros(shape, dtype=bool) t_ct = int(data.size/5) flag.flat[np.random.randint(0, flag.size, t_ct)] = True # True flags the data # -- end setup -- slcs = [slice(None)]*dim while np.any(~flag): # as long as there are any False's in flag for i in range(dim): # do each axis # make slices to shift view one element along the axis slcs1 = slcs[:] slcs2 = slcs[:] slcs1[i] = slice(0, -1) slcs2[i] = slice(1, None) # replace from the right repmask = np.logical_and(~flag[slcs1], flag[slcs2]) data[slcs1][repmask] = data[slcs2][repmask] flag[slcs1][repmask] = True # replace from the left repmask = np.logical_and(~flag[slcs2], flag[slcs1]) data[slcs2][repmask] = data[slcs1][repmask] flag[slcs2][repmask] = True 

Para una buena medida, aquí hay una visualización (2D) de las zonas sembradas por los datos marcados originalmente como True .

introduzca la descripción de la imagen aquí

Usando scipy.ndimage, su problema se puede resolver con la interpolación del vecino más cercano en 2 líneas:

 from scipy import ndimage as nd indices = nd.distance_transform_edt(invalid_cell_mask, return_distances=False, return_indices=True) data = data[tuple(ind)] 

Ahora, en la forma de una función:

 import numpy as np from scipy import ndimage as nd def fill(data, invalid=None): """ Replace the value of invalid 'data' cells (indicated by 'invalid') by the value of the nearest valid data cell Input: data: numpy array of any dimension invalid: a binary array of same shape as 'data'. data value are replaced where invalid is True If None (default), use: invalid = np.isnan(data) Output: Return a filled array. """ if invalid is None: invalid = np.isnan(data) ind = nd.distance_transform_edt(invalid, return_distances=False, return_indices=True) return data[tuple(ind)] 

Ejemplo de uso:

 def test_fill(s,d): # s is size of one dimension, d is the number of dimension data = np.arange(s**d).reshape((s,)*d) seed = np.zeros(data.shape,dtype=bool) seed.flat[np.random.randint(0,seed.size,int(data.size/20**d))] = True return fill(data,-seed), seed import matplotlib.pyplot as plt data,seed = test_fill(500,2) data[nd.binary_dilation(seed,iterations=2)] = 0 # draw (dilated) seeds in black plt.imshow(np.mod(data,42)) # show cluster 

resultado: introduzca la descripción de la imagen aquí

Hace algún tiempo escribí este guión para mi doctorado: https://github.com/Technariumas/Inpainting

Un ejemplo: http://blog.technariumas.lt/post/117630308826/healing-holes-in-python-arrays

Lento, pero hace el trabajo. El kernel gaussiano es la mejor opción, solo compruebe los valores de tamaño / sigma.

Puede intentar abordar su problema como:

 # main ideas described in very high level pseudo code choose suitable base kernel shape and type (gaussian?) while true loop over your array (moving average manner) adapt your base kernel to current sparsity pattern set current value based on adapted kernel break if converged 

Esto puede implementarse de manera bastante directa (especialmente si el rendimiento no es una preocupación importante).

Obviamente, esto es solo heurística y necesita hacer algunos experimentos con sus datos reales para encontrar el esquema de adaptación adecuado. Cuando vea la adaptación del kernel como el peso del kernel, le recomendamos que lo haga en función de cómo se hayan propagado los valores. Por ejemplo, sus ponderaciones para los soportes originales son 1 y están relacionadas con la decadencia en la iteración en la que emergieron.

También la determinación de cuándo este proceso realmente ha convergido puede ser difícil. Dependiendo de la aplicación, puede ser razonable dejar algunas “regiones vacías” sin rellenar.

Actualización : Aquí hay una implementación muy simple a lo largo de las líneas *) descritas anteriormente:

 from numpy import any, asarray as asa, isnan, NaN, ones, seterr from numpy.lib.stride_tricks import as_strided as ast from scipy.stats import nanmean def _a2t(a): """Array to tuple.""" return tuple(a.tolist()) def _view(D, shape, strides): """View of flattened neighbourhood of D.""" V= ast(D, shape= shape, strides= strides) return V.reshape(V.shape[:len(D.shape)]+ (-1,)) def filler(A, n_shape, n_iter= 49): """Fill in NaNs from mean calculated from neighbour.""" # boundary conditions D= NaN* ones(_a2t(asa(A.shape)+ asa(n_shape)- 1), dtype= A.dtype) slc= tuple([slice(n/ 2, -(n/ 2)) for n in n_shape]) D[slc]= A # neighbourhood shape= _a2t(asa(D.shape)- asa(n_shape)+ 1)+ n_shape strides= D.strides* 2 # iterate until no NaNs, but not more than n iterations old= seterr(invalid= 'ignore') for k in xrange(n_iter): M= isnan(D[slc]) if not any(M): break D[slc][M]= nanmean(_view(D, shape, strides), -1)[M] seterr(**old) A[:]= D[slc] 

Y una simple demostración del filler(.) En acción, sería algo así como:

 In []: x= ones((3, 6, 99)) In []: x.sum(-1) Out[]: array([[ 99., 99., 99., 99., 99., 99.], [ 99., 99., 99., 99., 99., 99.], [ 99., 99., 99., 99., 99., 99.]]) In []: x= NaN* x In []: x[1, 2, 3]= 1 In []: x.sum(-1) Out[]: array([[ nan, nan, nan, nan, nan, nan], [ nan, nan, nan, nan, nan, nan], [ nan, nan, nan, nan, nan, nan]]) In []: filler(x, (3, 3, 5)) In []: x.sum(-1) Out[]: array([[ 99., 99., 99., 99., 99., 99.], [ 99., 99., 99., 99., 99., 99.], [ 99., 99., 99., 99., 99., 99.]]) 

*) Entonces aquí el nanmean(.) Solo se utiliza para demostrar la idea del proceso de adaptación. Sobre la base de esta demostración, debería ser bastante sencillo implementar un esquema de pesaje de adaptación y decaimiento más complejo. También tenga en cuenta que no se presta atención al rendimiento de ejecución real, pero aún así debería ser bueno (con formas de entrada razonables).

Tal vez lo que está buscando es un algoritmo de aprendizaje automático, como una neural network o una máquina de vectores de soporte.

Puede consultar esta página, que tiene algunos enlaces a paquetes SVM para python: http://web.media.mit.edu/~stefie10/technical/pythonml.html