Реализация движущейся медианы с помощью scipy median_filter - PullRequest
0 голосов
/ 20 июня 2019

ПРИМЕЧАНИЕ: вопросы начинают становиться действительно большими, поэтому я переформулировал все в другом вопросе: Вычисление скользящей медианы с помощью scipy generic_filter и numpy median_filter дает различные выходные данные

Я хочу реализовать движущуюся медиану для одномерного вектора в Python. В основном, я хочу сделать последовательность вызовов numpy.median размера 5, но мне также нужна маска, так как для каждого вызова я должен избавиться от центрированного элемента и затем применить медиану:

    import numpy as np

    v = np.array([0, 1, 2, 3, 4]) # vector of size 5
    mask = np.array([True, True, False, True, True]) # mask that gets rid of the centered element
    print(np.median(v[mask])) #(1 + 3) / 2 = 2.0

Это то, что показано на изображении ниже. Я не могу присоединиться к этому сообщению (может быть, оно слишком большое?), Я прошу прощения за это.

https://www.noelshack.com/2019-25-4-1561022011-median1.png

Проблема в том, что вызовы numpy на таких маленьких массивах слишком медленные. Я знаю, что есть ось параметра, и я мог бы переупорядочить свой вектор как матрицу, чтобы каждая строка имела 5 элементов, и я мог применить медиану за один вызов. Но моя проблема заключается в создании этой новой матрицы, я бы предпочел использовать оригинальный вектор и сэкономить время. Я также предпочел бы вызывать существующие функции, поскольку они были бы гораздо более оптимизированы, чем то, что я мог бы сделать.

Я заметил, что в scipy есть эта функция: median_filter . Но я не могу понять его вывод. На изображении, указанном ниже, я пытаюсь объяснить, что я понимаю. Но, как вы можете видеть внизу изображения, есть разница между тем, что я ожидаю, и тем, что делает функция.

https://www.noelshack.com/2019-25-4-1561023099-median2.png

Мне кажется, я неправильно понимаю использование параметра footprint. Вот как я вызываю функцию, чтобы получить результат, показанный оранжевым вектором. Но с этими параметрами я ожидаю получить результат, показанный синим вектором.

    import numpy as np
    from scipy import ndimage, misc

    v = np.array([0., 1., 2., 3, 4., 5., 6., 7., 8., 9., 10.])
    print(ndimage.median_filter(v, footprint=np.array([1, 1, 0, 1, 1]), 
                                       output=np.float64, mode="mirror", origin = 0))

    #gives : [2. 2. 3. 4. 5. 6. 7. 8. 9. 9. 9.]
    #but I expect : [1.5 1.5 2.  3.  4.  5.  6.  7.  8.  8.5 8.5]

Ты хоть представляешь, что не так с моим звонком? Или как я мог получить движущуюся медиану, используя встроенные функции Python? Заранее спасибо.

ОБНОВЛЕНИЕ ________________________________________________________ ОБНОВЛЕНИЕ

Хорошо, коллега показал мне эту функцию: scipy generic_filter. И это действительно делает то, что я хочу, но не так быстро, как я хочу. В моей программе я должен сделать много медиан и хотел бы, чтобы этот шаг занимал как можно меньше времени. Ниже вы увидите, что при создании подматрицы для моих входных данных код выполняется быстрее, чем с generic_filter. Есть ли способ избежать создания подматрицы? Есть ли в Python способ создать матрицу представлений для моих входных данных, чтобы я мог сразу запустить numy медиану и избежать затрат на создание новой матрицы? Еще раз спасибо за вашу помощь.

# Note : How I deal with borders is not important for the question. 
# With my method, I truncate the filter. 
# generic_filter applies a mirror on the input data and uses the same filter.

import numpy as np
import scipy.ndimage as sc

v = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

def myMovingMedian5(IN) :
    # 5 masks, 4 for borders, 1 for all other values.
    mask0 = np.array([False, True, True])
    mask1 = np.array([True, False, True, True])
    mask2 = np.array([True, True, False, True, True])
    mask3 = np.array([True, True, False, True])
    mask4 = np.array([True, True, False])

    nR = IN.shape[0]

    # Generate a sub matrix to compute most of the medians (except on borders).
    # Can I avoid this step and use IN directly ?
    # Or is it possible to make a matrix of views of IN and avoid the creation of new data ?
    TMP = np.zeros((nR - 4, 4))
    indTMP = 0
    for i in range(2, nR - 2) :
        TMP[indTMP, 0:4] = (IN[i - 2:i + 3])[mask2]
        indTMP = indTMP + 1
    #TMP :
#     [[ 0.  1.  3.  4.]
#      [ 1.  2.  4.  5.]
#      [ 2.  3.  5.  6.]
#      [ 3.  4.  6.  7.]
#      [ 4.  5.  7.  8.]
#      [ 5.  6.  8.  9.]
#      [ 6.  7.  9. 10.]]

    # Allocate OUT matrix with 4 more elements for the borders
    OUT = np.zeros(nR)
    # Replace its center part by applying the median on each line of TMP
    OUT[2:nR - 2] = np.median(TMP, axis = 1)


    # Add remaining 4 medians on borders
    OUT[0] = np.median((IN[0:3])[mask0])
    OUT[1] = np.median((IN[0:4])[mask1])

    beforeLast = nR - 2
    OUT[beforeLast] = np.median((IN[beforeLast - 2:beforeLast + 2])[mask3])
    OUT[nR - 1] = np.median((IN[nR - 3:nR])[mask4])

    return OUT

print(myMovingMedian5(v))
%timeit myMovingMedian5(v)

print(sc.generic_filter(v, sc.median, footprint=np.array([1, 1, 0, 1, 1]), mode = "mirror", output=np.float64))
%timeit sc.generic_filter(v, sc.median, footprint=np.array([1, 1, 0, 1, 1]), mode = "mirror", output=np.float64)

отпечатков: myMovingMedian5:

[1.5 2. 2. 3. 4. 5. 6. 7. 8. 8. 8.5]

121 мкс ± 2,33 мкс на цикл (среднее ± стандартное отклонение из 7 циклов, по 10000 циклов каждый)

scipy generic_filter:

[1.5 1.5 2. 3. 4. 5. 6. 7. 8. 8.5 8.5]

310 мкс ± 1,5 мкс на цикл (среднее ± стандартное отклонение из 7 циклов, 1000 циклов каждый)

ЗАКЛЮЧИТЕЛЬНОЕ ОБНОВЛЕНИЕ _____________________________________________ ЗАКЛЮЧИТЕЛЬНОЕ ОБНОВЛЕНИЕ

Я проверил время scipy median_filter и, хотя он не дает правильного ответа, он определенно быстрее:

%timeit sc.median_filter(v, footprint=np.array([1, 1, 0, 1, 1]), output=np.float64, mode="mirror")

# [2. 2. 3. 4. 5. 6. 7. 8. 9. 9. 9.]
# 12.3 µs ± 62.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Но почему мои вызовы general_filter и median_filter не дают одинаковых выходных данных?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...