3D-поиск пиков с использованием CWT

Question

Я пытаюсь обнаружить пики в массиве 3d (гистограмма), используя метод, подобный scipi.signal.find_peaks_cwt . И обсуждается в этой статье .

Источник : Пан Ду, Уоррен А. Киббе, Саймон М. Лин. * Биоинформатика, том 22, выпуск 17, 1 сентября 2006 года, страницы 2059–2065, https://doi.org/10.1093/bioinformatics/btl355

Я использую разность гауссовского приближения для вычисления ядра вейвлета. Гауссианы имеют круглую форму и используют ковариационную матрицу шкалы, чтобы определить их ширину (cov = np.identity(3) / s).

В части, связанной с поиском хребта, я испытываю трудности. Я создаю стек трехмерных массивов, содержащих свертку между входом (гистограммой) и каждым из ядер вейвлета. В моей реализации метод скоринга передает стек выходных данных от самых узких до самых широких вейвлетов, пытаясь «обнаружить» локальные максимумы при разных разрешениях.

Код

def get_gaussian_kernel(size, s=1.0):
    if size == 1:
        return np.array([[[1]]])
    assert(size > 1)
    cov = np.identity(3) / s
    rv = multivariate_normal((0.5, 0.5, 0.5), cov)
    grid = np.mgrid[:size, :size, :size].T
    kernel = rv.pdf(grid / (size - 1))
    return kernel


def marr_dog_func(size, s=64):
    '''
    Marr D. , Hildreth E. and Brenner Sydney
    Theory of edge detection207Proc. R. Soc. Lond. B.
    https://doi.org/10.1098/rspb.1980.0020
    '''

    # Marr wavelet approximated by Difference of Gaussians
    k1 = get_gaussian_kernel(size, s)
    k2 = get_gaussian_kernel(size, 1.6 * s)
    return k2 - k1


def cwt_peak_detect_3d(h, widths, tolerance=1.0, func=None):
    if not func:
        func = marr_dog_func
    stack = []
    for w in widths:
        kernel = func(w)  # Wavelet kernel
        h2 = convolve(h, kernel, mode='same')  # Convolve with wavelet
        stack.append(h2)  # Add to stack

    def score_point(idx):
        score = 0
        peak = max(0, stack[0][idx])
        for s in stack[1:]:
            v = max(0, s[idx])
            if v < tolerance * peak:
                break
            score += 1
            peak = v
        return score

    scored = np.zeros_like(h)
    for idx in np.ndindex(h.shape):
        scored[idx] = score_point(idx)
    return scored


def peak_find_3d(pixels, n_bins=49, bias=0):
    bin_range = np.arange(-1, n_bins + 2) / n_bins
    bins = (bin_range,) * 3
    H, edges = np.histogramdd(pixels, bins)

    widths = [
        (n_bins // 16) * 2 + 1,
        (n_bins // 8) * 2 + 1,
        (n_bins // 4) * 2 + 1,
        (n_bins // 2) * 2 + 1,
    ]

    h2 = cwt_peak_detect_3d(H, widths)

    coords, values = zip(*list(np.ndenumerate(h2)))
    normalized_coords = np.array(coords) / (h2.shape[0] - 1)
    stuff = sorted(zip(normalized_coords, values), key=lambda x: x[1],
                   reverse=True)
    top_coords, top_values = zip(*stuff)
    return np.array(top_coords), top_values

Демо: центр фрагмент функции 3D-вейвлета для различной ширины

import matplotlib.pyplot as plt

kernel_size = 127
widths = [
    kernel_size // 2 ,
    kernel_size // 4 ,
    kernel_size // 8 ,
    kernel_size // 16 ,
]

kernels = [marr_dog_func(w) for w in widths]
fig, ax = plt.subplots(1, len(kernels))
ax_flat = ax.flatten()
for a, k in zip(ax_flat, kernels):
    print(k.shape)
    center_x = k.shape[0] // 2 + 1
    k2d = k[center_x]
    a.imshow(k2d, interpolation='nearest', cmap=plt.cm.gray,
             vmin=-32, vmax=32)
plt.show()

Выход

Демонстрация: 3D-график, показывающий вывод

from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
import matplotlib.pyplot as plt

np.random.seed(546)

cdict = {
    'red': [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]],
    'green': [[0.0, 0.5, 0.5], [1.0, 0.5, 0.5]],
    'blue': [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]],
    'alpha': [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]
}
cmap = LinearSegmentedColormap('foo',segmentdata=cdict)

s1 = np.random.random((10, 3))
#s2 = np.random.random((1000, 3))
stuff = np.vstack((s1, s1, s1, s1))
n_bins = 19
top_cooords, top_values = peak_find_3d(stuff, n_bins)
top_cooords = np.array(list(top_cooords))
top_values = top_values
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(top_cooords[:, 0], top_cooords[:, 1], top_cooords[:, 2],
           c=top_values, cmap=cmap, vmin=0, vmax=5)
plt.show()

Ввод гистограмма

Выходная гистограмма

В моем случае вывод очень шумный для тривиального ввода и совсем не то, что я ищу. Неожиданно различия в пиках между строками преобразования имеют порядок ошибки точности с плавающей запятой. Мои widths слишком редки? Является ли мой алгоритм поиска гребней проблемой или есть проблема в том, как я вычисляю CWT?