Обнаружение пиков в шумном 2d массиве

Ответ 1

Я добавляю этот ответ, потому что это решение, которое я использовал. Это комбинация комментария Би Рико здесь (30 мая в 18:54) и ответ, заданный в этом вопросе: Найти пик 2d гистограммы.

Как выясняется, используя алгоритм обнаружения пика из этого вопроса Распознавание пиков в 2D-массиве только усложняет ситуацию. После применения фильтра Гаусса к изображению все, что нужно сделать, это запросить максимальный бит (как указал Би Рико), а затем получить максимум в координатах.

Поэтому вместо того, чтобы использовать функцию обнаружения-пиков, как я сделал выше, я просто добавляю следующий код после получения гистограммы Gaussian 2D:

# Get 2D histogram.
H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
# Get Gaussian filtered 2D histogram.
H1 = gaussian_filter(H, 2, mode='nearest')
# Get center of maximum in bin coordinates.
x_cent_bin, y_cent_bin = np.unravel_index(H1.argmax(), H1.shape)
# Get center in x,y coordinates.
x_cent_coor , y_cent_coord = np.average(xedges[x_cent_bin:x_cent_bin + 2]), np.average(yedges[y_cent_g:y_cent_g + 2])

Ответ 2

Отвечая на последнюю часть вашего вопроса, всегда у вас есть точки на изображении, вы можете найти их координаты, выполнив в некотором порядке локальные максимумы изображения. Если ваши данные не являются точечным источником, вы можете применить маску к каждому пику, чтобы избежать максимального значения максимума при выполнении будущего поиска. Я предлагаю следующий код:

import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
import copy

def get_std(image):
    return np.std(image)

def get_max(image,sigma,alpha=20,size=10):
    i_out = []
    j_out = []
    image_temp = copy.deepcopy(image)
    while True:
        k = np.argmax(image_temp)
        j,i = np.unravel_index(k, image_temp.shape)
        if(image_temp[j,i] >= alpha*sigma):
            i_out.append(i)
            j_out.append(j)
            x = np.arange(i-size, i+size)
            y = np.arange(j-size, j+size)
            xv,yv = np.meshgrid(x,y)
            image_temp[yv.clip(0,image_temp.shape[0]-1),
                                   xv.clip(0,image_temp.shape[1]-1) ] = 0
            print xv
        else:
            break
    return i_out,j_out

#reading the image   
image = mpimg.imread('ggd4.jpg')
#computing the standard deviation of the image
sigma = get_std(image)
#getting the peaks
i,j = get_max(image[:,:,0],sigma, alpha=10, size=10)

#let see the results
plt.imshow(image, origin='lower')
plt.plot(i,j,'ro', markersize=10, alpha=0.5)
plt.show()

Изображение ggd4 для теста можно загрузить с помощью:

http://www.ipac.caltech.edu/2mass/gallery/spr99/ggd4.jpg

Первая часть - получить некоторую информацию о шуме на изображении. Я сделал это, вычислив стандартное отклонение полного изображения (на самом деле лучше выбрать небольшой прямоугольник без сигнала). Это говорит нам, сколько шума присутствует в изображении. Идея получить пики состоит в том, чтобы запросить последовательные максимумы, которые превышают определенный порог (скажем, 3, 4, 5, 10 или 20 раз). Это то, что делает функция get_max. Он выполняет поиск максимумов до тех пор, пока один из них не станет ниже порога, налагаемого шумом. Чтобы избежать нахождения одного и того же максимума много раз, необходимо удалить пики из изображения. В общем случае форма маски для этого сильно зависит от проблемы, которую нужно решить. для случая звезд должно быть хорошо удалить звезду с помощью функции Гаусса или что-то подобное. Я выбрал для простоты квадратную функцию, а размер функции (в пикселях) - это переменная "размер". Я думаю, что из этого примера каждый может улучшить код, добавив более общие вещи.

EDIT:

Исходное изображение выглядит так:

Пока изображение после идентификации светящихся точек выглядит следующим образом:

Ответ 3

Слишком много n00b для, чтобы прокомментировать Алехандро в другом месте. Я бы немного уточнил его код, чтобы использовать предварительно выделенный массив numpy для вывода:

def get_max(image,sigma,alpha=3,size=10):
    from copy import deepcopy
    import numpy as np
    # preallocate a lot of peak storage
    k_arr = np.zeros((10000,2))
    image_temp = deepcopy(image)
    peak_ct=0
    while True:
        k = np.argmax(image_temp)
        j,i = np.unravel_index(k, image_temp.shape)
        if(image_temp[j,i] >= alpha*sigma):
            k_arr[peak_ct]=[j,i]
            # this is the part that masks already-found peaks.
            x = np.arange(i-size, i+size)
            y = np.arange(j-size, j+size)
            xv,yv = np.meshgrid(x,y)
            # the clip here handles edge cases where the peak is near the 
            #    image edge
            image_temp[yv.clip(0,image_temp.shape[0]-1),
                               xv.clip(0,image_temp.shape[1]-1) ] = 0
            peak_ct+=1
        else:
            break
    # trim the output for only what we've actually found
    return k_arr[:peak_ct]

При профилировании этого и кода Алехандро, используя его пример изображения, этот код примерно на 33% быстрее (0,03 секунды для кода Алехандро, 0,02 секунды для моего.) Я ожидаю на изображениях с большим количеством пиков, это будет еще быстрее - добавление вывод в список будет медленнее и медленнее для большего количества пиков.

Ответ 4

Я думаю, что первый шаг, который здесь требуется, - выразить значения в H через стандартное отклонение поля:

import numpy as np
H = H / np.std(H)

Теперь вы можете поместить порог на значения этого H. Если шум считается гауссовым, выбирая порог 3, вы можете быть абсолютно уверены (99.7%), что этот пиксель может быть связан с реальным пиком и нет шум. См. здесь.

Теперь возможен дальнейший отбор. Мне не совсем ясно, что именно вы хотите найти. Вы хотите точное расположение пиковых значений? Или вы хотите, чтобы одно место для кластера пиков находилось в середине этого кластера?
В любом случае, начиная с этой точки со всеми значениями пикселей, выраженными в стандартных отклонениях поля, вы должны иметь возможность получить то, что хотите. Если вы хотите найти кластеры, вы можете выполнить поиск ближайшего соседа в сетчатых точках > 3-сигма и установить порог на расстоянии. То есть только соединяйте их, когда они достаточно близко друг к другу. Если подключено несколько точек сетки, вы можете определить это как группу/кластер и вычислить некоторый (сигма-взвешенный?) Центр кластера. Надеюсь, мой первый вклад в Stackoverflow вам пригодится!

Ответ 5

Как я это сделаю:

1) нормализуют H между 0 и 1.

2) выберите пороговое значение, как предполагает tcaswell. Это может быть от 0,9 до 0,99, например

3) используйте маскированные массивы, чтобы сохранить только координаты x, y с H выше порога:

import numpy.ma as ma
x_masked=ma.masked_array(x, mask= H < thresold)
y_masked=ma.masked_array(y, mask= H < thresold)

4) теперь вы можете усреднять по маске в замаскированных координатах, с весом чего-то вроде (H-threshold) ^ 2 или любой другой силой, большей или равной одной, в зависимости от вашего вкуса/тестов.

Комментарий: 1) Это не является надежным по отношению к типу пиков, которые у вас есть, так как вам, возможно, придется адаптировать thresold. Это второстепенная проблема; 2) Это НЕ работает с двумя пиками, как есть, и даст неправильные результаты, если второй пик выше порогового значения.

Тем не менее, он всегда даст вам ответ без сбоев (с плюсами и минусами вещи).