описание алгоритма.
Необходимо выделить точки в которых интенсивность превышает уровень шума в N раз по отношению к амплитуде шума.
В общем случае в изображении содержатся: ток смещения (bias), темновой ток, фон неба. Им соответствуют шумы: считывания, темновой шум, пуассоновский шум неба. Шумы суммируются. Кроме того, на некалиброванных кадрах есть горячие пиксели. Шум темнового тока и фона неба - пуассоновские, стандартное отклонение (сигма) = корню из среднего.
Решаем задачу для простого случая: изображения без вычитания дарков, на кадре нет туманностей = фон неба примерно постоянный по кадру. Этого вполне достаточно для большинства случаев. Относительно быстрый.
Применяем мягкий медианный фильтр для удаления горячих пикселей.
Data = filters.median_filter(Data, size=2, origin=0, mode= 'reflect')
Находим медианное среднее для кадра = приличная оценка среднего фона. Корень квадратный из этого значения дает сильно завышенную оценку уровня шума (так как мы не вычли ток смещения = bias)
##Poisson statistics
Bckgrnd_median = np.median(Data)
print('Median background=', Bckgrnd_median) #good background level estimate
RMS_median = np.sqrt(Bckgrnd_median)
print('Poisson noise=', RMS_median)
берем пиксели со значениями среднее+-3*оценка шума:
hist_l = Bckgrnd_median-3*RMS_median #low boundaries for histigram
hist_h = Bckgrnd_median+3*RMS_median #high boundaries for histigram
строим гистограмму, ищем пик гистограммы и вариацию (ширину распределения) исходя из предположения о нормальном (гауссовом) распределении.
Gauss = lambda x, A, B, C, D: A*np.exp(-(x-B)**2/(2*C*C))+D #model for Gauss distribution of background
p0 = np.array([max(counts), Bckgrnd_hist, 3, min(counts)])
#fit Gauss
try:
popt, pcov = curve_fit(Gauss, bins, counts, p0, ftol=0.1, maxfev=10000)
print('Background mode=', popt[1]) #fitted maximum of histogram
print('RMS =', popt[2]) #fitted RMS of histogram
Noise = popt[2]
## plt.plot(bins, counts)
## plt.plot(bins, Gauss(bins, *popt), color='green')
## plt.vlines(Bckgrnd_hist, 0, max(hist[0]), color='blue')
## plt.vlines(popt[1], 0, max(hist[0]), color='red')
## plt.show()
строчки в комментариях - визуализация результата
popt[1] - средний уровень фона (dias+темновой+небо)
popt[2] - амплитуда шума (сигма)
звездами считаем все, что имеет интенсивность выше среднего уровня фона+амплитуда шума*пороговое значение SN
Data_Peaks = Data - (popt[1] + SN*popt[2])
detected_peaks = Data_Peaks>0
дальше нужно/можно проверить размеры найденных деталей, интенсивность, положение (не слишком близко к краю)
работает не без проблем: на очень ярких звездах находит лучи, в тесных парах не разделяет звезды.
Сложности начинаются если на изображении много горячих пикселей и медианный фильтр с ними не справляется.
Особая [----] если фон неба сильно переменный (туманности). Тогда нужно оценивать фон неба локально в каждой точке. Можно сделать это используя медианный фильтр с шириной ядра больше, чем размеры звезд (это я использую иногда, медленно работает!). Можно попробовать минимальный (ищет минимум) фильтр. Можно оценивать в узлах сетки и интерполировать (так делает SExtractor и Photutils).