填充 numpy 数组中的空白

Question

我只想用最简单的术语对 3D 数据集进行插值。线性插值，最近邻，所有这些就足够了（这是为了开始一些算法，所以不需要精确的估计）。

在新的 scipy 版本中，像 griddata 这样的东西会很有用，但目前我只有 scipy 0.8。所以我有一个“立方体”（data[:,:,:], (NixNjxNk)）数组和一个标志数组（flags[:,:,:,] code>、True 或 False）具有相同的大小。我想使用例如数据中最近的有效数据点或“接近”点的某种线性组合，将我的数据插入到标志的相应元素为 False 的数据元素中。

数据集中至少在两个维度上可能存在较大间隙。除了使用 kdtree 或类似的方法编写成熟的最近邻算法之外，我真的找不到通用的 N 维最近邻插值器。

Answer 1

使用 scipy.ndimage，您的问题可以通过两行最近邻插值来解决：

from scipy import ndimage as nd

indices = nd.distance_transform_edt(invalid_cell_mask, return_distances=False, return_indices=True)
data = data[tuple(ind)]

---

现在，以函数的形式：

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)]

使用示例：

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

结果：

Answer 2

您可以设置晶体生长式算法，沿每个轴交替移动视图，仅替换标记为 False 但具有 True 邻居的数据。这给出了类似“最近邻”的结果（但不是欧几里得距离或曼哈顿距离——我认为如果您正在计算像素，计算具有公共角的所有连接像素，它可能是最近邻）这对于 NumPy 来说应该相当有效因为它仅迭代轴和收敛迭代，而不是一小部分数据。

原始、快速、稳定。我认为这就是您所追求的：

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

为了更好地衡量，这里是由最初标记为 True 的数据播种的区域的可视化 (2D)。

Answer 3

不久前，我为博士学位编写了这个脚本： https://github.com/Technariumas/Inpainting

示例： http://blog.technariumas.lt/post/117630308826 /healing-holes-in-python-arrays

速度慢，但有效。高斯核是最好的选择，只需检查大小/西格玛值。

Answer 4

您可以尝试解决您的问题，例如：

# 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

这实际上可以以非常简单的方式实现（特别是如果性能不是最关心的问题）。

显然这只是启发法，您需要用实际数据进行一些实验才能找到合适的适应方案。当将内核适应视为内核重新权衡时，您可能希望根据值的传播方式来执行此操作。例如，原始支撑的权重是 1，它们的衰减与它们出现的迭代有关。

此外，确定该过程何时真正收敛可能是一个棘手的问题。根据应用的不同，最终保留一些“间隙区域”保持“未填充”可能是合理的。

更新：这是一个非常简单的实现，与上面描述的 *) 一致：

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]

以及 filler(.) 操作的简单演示，类似于：

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.]])

*)所以这里的nanmean(.)只是用来演示适应过程的想法。基于此演示，实施更复杂的适应和衰减权重方案应该非常简单。另请注意，没有关注实际的执行性能，但它仍然应该很好（具有合理的输入形状）。

Answer 5

也许您正在寻找的是机器学习算法，例如神经网络或支持向量机。

您可以查看此页面，其中包含一些 Python 的 SVM 包的链接：http: //web.media.mit.edu/~stefie10/technical/pythonml.html