Seeing monitor

[ekvi]

Положение максимальной суммы 3х3 даст положение звезды

"Максимальная сумма на участке 16х16" пикселей дала координаты (858,513) - см 1ю илл. (rSpot=8 - радиус в пикселях вокруг максимально яркого пикселя).
Но в Максиме на этой координате - пустота. MaxImDL даёт координаты Полярной: (684,502), т.е. разница в 174 пикселя по абсциссе - 2я илл.
На 3й илл.- вырезки из 1й и 2й илл. с упомянутыми координатами центра Полярной на 2м снимке из серии снимков Полярной.
Поскольку по ординате смещение невелико, то цифры 513 и 502 - рядом. Но, даже, если учесть, что Максим читает файлы с конца (делаем флип + мирро), всё равно по абсциссе попадаем "пальцем в небо".

krussh верно указывает: нужно не суммировать, а брать МЕДИАННОЕ значение из массива яркостей пикселей, окружающих максимум. Вот этот алгоритм и опробуем.

Вот здесь https://yadi.sk/d/AMk5YgBUqp5j5 лежит болванка программы SeeingGauger на Паскале в Дельфи7, с исходниками. Для уменьшения веса архива в нём оставлено только 10 снимков из 300 (остальные можно скачать из архива moskow, пост #50).
Будет желание - поэкспериментируйте.

[krussh]

предлагаю взять стандартный тестовый кадр и для начала считать все для него. Например вот этот: http://optlab.astro.usu.ru/Python/seeing/test.fit
максимальный пиксель на кадре находится примерно в центре полярной.
координаты 573, 507. значение в пикселе 52.

координаты центра изображения
по Гауссу (krussh): 572.897 507.300
центроид (MaximDL): 572.913 507.319

[Gleb1964]

бегущий медианный фильтр это хороший метод подавления горячих пикселов.
А я предлагал, фактически, пробежать бегущим средним 3х3. Не знаю, правильно ли вы его реализовали.
Алгоритм такой - окошко 3х3 с весовыми коэффициентами (например, все 0.333 0.111) бежит по изображению с шагом 1пиксел. В каждой точке рассчитывается сумма перемножения коэффициентов окна на значения пикселов, эта сумма - новое значение пиксела в фильтрованном файле. Естественно, окошко на краях изображения отступает от края на 1пиксел - т.е. формат нового фильтрованного изображения будет на 2 пиксела меньше по горизонтали и по вертикали.
Можно бегать окнами и 5х5 и 7х7 и т.д., можно использовать разную систему весов в окошках, сумму весов хорошо нормировать к единице.
Так можно делать и бегущее среднее - все веса одинаковы, и фильтрацию гауссоидой - если сделать веса в форме гауссоиды, и высокочастотный фильтр, и фильтрирование типа Савитского-Галея (фитирование полиномиальной функции некого порядка), можно подчеркивающие определенные направления градиентов реализовать.
Фильтрирование гауссоидой даже немного лучше бегущего среднего, особенно, если размер фильтра сопряжен с ожидаемым размером изображения звезды.

для krussh: у меня нечем читать-конвертировать fit формат, подскажете, чем? или какой-нибудь формат попроще дадите в качестве примера

[krussh]

для работы с fits полезно использовать готовую библиотеку:
http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html
там же есть ссылка на некоторые программы для работы с файлами.

алгоритмы для поиска и фотометрии:
http://www.adass.org/adass/proceedings/adass98/mighellkj/

[ekvi]

для работы с fits полезно использовать готовую библиотеку:

Вадим, нижайший Вам поклон!
Теперь понятно, почему Вы выбрали Python'а.
С фитсами давно работаю, ~ с 2010го - с тех пор, как Вы снабдили меня Максимом.
Читать/писать - никогда не возникало проблем: утилиты на Паскале отлажены на 2х программах: AVIsum и FITsum.

[ekvi]

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

В SG применён алгоритм вычисления ЦТ. Но попробую и "бегущий с усредняющими весами".

[mo]

Наконец-то освободилась одна QHY5чб и искатель. На днях соберу из них комплекс, постоянно нацеленный на полярку.

Возникает вопрос в необходимой периодичности фотографирования полярки. Я без проблем смогу запустить фотографирование полярки всю ночь напролёт с максимальной скоростью, но объём данных будет огого. Будут ли эти данные иметь хоть какой-то смысл?

[krussh]

Наконец-то освободилась одна QHY5чб и искатель. На днях соберу из них комплекс, постоянно нацеленный на полярку.

Возникает вопрос в необходимой периодичности фотографирования полярки. Я без проблем смогу запустить фотографирование полярки всю ночь напролёт с максимальной скоростью, но объём данных будет огого. Будут ли эти данные иметь хоть какой-то смысл?

надо обрабатывать кадр как он появляется, сохранять положение и удалять его. частота съемки подбирается так, чтобы обработка успевала.
опять же можно все сделать на питоне или баше. проверяем папку, если есть кадр, то делаем обработку и удаляем. как только накопится 300-500 измерений, вычисляем seeing и пишем в базу.  при необходимости из нее вытаскиваем и строим графики.

[mo]

но время накопления нужно покороче, иначе вариация координат за время накопления трансформируется в размазывание изображения звезды

Формат (экспозиция, усиления и пр.) кадра уже урегулирован.
Оставался вопрос как часто нужен кадр. Ответ дан.

[ekvi]

попробуйте поменять знак на обратный

Многое для меня странно в этом "зазеркалье":
1. Обходишь кадр поперёк(т.е. высоту кадра заменяешь его шириной)- координаты Полярной приближаются к Макимовским.
2. ЦТ считаешь с весами - у конечного пятна распределение точек уже не ранжируется рядами, а приближается к случайному распределению.
3. Да, и от знака при суммировании моментов (яркость*расстояние до ЦТ) вид распределения также зависит (см. илл., на которой веса берутся, как 1/abs(r) или 1/sqr(r)).

[ekvi]

Возникает вопрос в необходимой периодичности фотографирования полярки

К сожалению, не умею обрабатывать видеопотоки, но мониторить небо нужно непрерывно, периодически обрабатывая содержимое потока - именно такая и нужна программа.
Частота же периода определяется, как и указал krussh, временем, необходимым программе на получение результата. В реальном времени получается, что программа постоянно должна переваривать содержимое потока и выдавать результат.

[ekvi]

опять же можно все сделать на питоне

- потому, что энтузиасты, живущие в альтернативном Линюксе, более пассионарны, и у них решены актуальные задачи.
Но для широкого пользователя программу придётся всё равно переводить под априорные стандарты Windows - не у каждого ЛА хватит упорства продираться сквозь архаику Пытона и своеобразие Линюкса.
Кстати, Питон позволяет сгенерировать EXE-файл. Это позволит остальным-прочим смертным, в том числе и Вашему покорному слуге, насладиться плодами Ваших трудов.

[krussh]

Многое для меня странно в этом "зазеркалье":

1) просто ищем максимальный пиксель на изображении. считаем, что звезда все же ярче шума.
2) берем область 2Nx2N пикселей вокруг найденного пикселя
3) проходим по всем пикселям (два вложенных цикла по ii и jj) и считаем
    а) сумму интенсивность*первый индекс: XSum = XSum + (ii * input[ii, jj]);
    б) сумму интенсивность*второй индекс: YSum = YSum + (jj * input[ii, jj]);
    в) сумму интенсивностей: MSum = MSum + input[ii, jj];
4) получаем координаты центроида как X = XSum / MSum и Y=YSum / MSum
можно ввести апертуру и считать суммы только в круге.
веса вводить не нужно.
необходимо помнить, что целые значения координаты соответствуют ЦЕНТРУ пикселя!

Но для широкого пользователя программу придётся всё равно переводить под априорные стандарты Windows - не у каждого ЛА хватит упорства продираться сквозь архаику Пытона и своеобразие Линюкса.

Питон работает в Windows не хуже, чем в Linux) За это я его и люблю. На какой версии Win Вы работаете? Я подберу пакеты и выложу на сервер.

Кстати, Питон позволяет сгенерировать EXE-файл. Это позволит остальным-прочим смертным, в том числе и Вашему покорному слуге, насладиться плодами Ваших трудов.

Наверное можно, слышал об этом. Но это будет не лучшее решение. для пользователей Win постараюсь написать программку как будет побольше свободного времени. Заготовка уже есть - это программа гидирования нашего 1.2-м телескопа.

[Gleb1964]

2) берем область 2Nx2N пикселей вокруг найденного пикселя
3) проходим по всем пикселям (два вложенных цикла по ii и jj) и считаем
    а) сумму интенсивность*первый индекс: XSum = XSum + (ii * input[ii, jj]);
    б) сумму интенсивность*второй индекс: YSum = YSum + (jj * input[ii, jj]);
    в) сумму интенсивностей: MSum = MSum + input[ii, jj];
4) получаем координаты центроида как X = XSum / MSum и Y=YSum / MSum

т.е., для простоты - квадратную область с изображением звезды можно суммированием по вертикали или по горизонтали преобразовать в два линейных варианта - среднее распределение сигнала вдоль оси X и Y, потом найти центры тяжести как в одномерном случае. Но для неквадратных областей тогда понадобятся веса.
Кстати, попробовал на предложенном примере - там звезда весьма резкая, \(\sigma \approx  1\), там можно очень просто максимум находить гауссоидой по трем точкам

[ekvi]

необходимо помнить, что целые значения координаты соответствуют ЦЕНТРУ пикселя!

Это постоянно имел в виду. И в РОС зрачковая функция так и считается, но в СинГэйдже этим как-то пренебрёг...
Всё остальное - элементарно и многажды перепробовано в различных вариациях.

"Питон позволяет сгенерировать EXE-файл" (об этом см. у Лутца в 1 и 2 рекламных главах).

Наверное можно, слышал об этом. Но это будет не лучшее решение

Да, потому, что экзе-файл всегда будет для уникальной ситуации.

На какой версии Win Вы работаете?

Windows XP, SP3, 32-битная, 1-ядерная.
На сайте Питона бывал неоднократно. К сожалению, многие утилиты, написанные питонцами и которыми Вы воспользовались при создании своей программы, работают только на 64-битных процессорах (и работа с фитс-файлами, и преобразования изображений, и БПФ и др).

Питон работает в Windows не хуже, чем в Linux

- т.к. сценарии (скрипты) - не что иное, как batch-файлы (например, AUTOEXEC.BAT - ), заимствованные линюксоидами из ДОС. Но, поскольку многие приверженцы Линюкса родились после 1990го, то о Windows они рассуждают как иностранцы о стране проживания, а о DOS и вообще не поминают.
Да, у сценария мощные рычаги: запуск сторонних программ, организация параллельных процессов и т.п. системные возможности, но кто их задействовал и ими воспользовался?! А ведь для рассматриваемой программы это не только актуально, но и реализуемо:
1. Запускаем Максима с требуемыми параметрами (обрезка снимков, поиск макс. пикселя и т.д.) и передаём полученные результаты в свою программу.
2. Операционная система - на автомате - загружает ядра процессора параллельными задачами, многократно сокращая время получения результата.
и т.д и т.п.

попробовал на предложенном примере - там звезда весьма резкая, \(\sigma \approx  1\), там можно очень просто максимум находить гауссоидой по трем точкам

krussh с самого начала подчёркивал, что требования его программы к снимкам специфические: горячие пиксели д.б. бледнее полезного сигнала. Это - недостаток настоящей версии его программы, т.к. обычно "горячие" пиксели всегда - самые яркие на снимке, особенно слабых объектов. В своей программе я пытаюсь учесть это.
Гауссоиду пробовал, но пока без эффекта.
В принципе, для данной программы конкретные координаты и не нужны, т.к. важны параметры распределения в суммарном пятне.
Но разборки эти нужны для дальнейшего использования разрабатываемых утилит.

[krussh]

Windows XP, SP3, 32-битная, 1-ядерная.
На сайте Питона бывал неоднократно. К сожалению, многие утилиты, написанные питонцами и которыми Вы воспользовались при создании своей программы, работают только на 64-битных процессорах (и работа с фитс-файлами, и преобразования изображений, и БПФ и др).

Python 2.7 для Windows x86, пакеты и инструкция по установке лежат по адресу http://optlab.astro.usu.ru/Python/Python_main/Python_27_x86/
в время установки требуется подключение к Интернет.

там же тестовая программка в которой используются все основные пакеты и тестовый файл.

[krussh]

.е., для простоты - квадратную область с изображением звезды можно суммированием по вертикали или по горизонтали преобразовать в два линейных варианта - среднее распределение сигнала вдоль оси X и Y, потом найти центры тяжести как в одномерном случае. Но для неквадратных областей тогда понадобятся веса.

вес 1 или 0 в зависимости от предполагаемого расстояния до ЦТ. можно вычислять расстояние от центра области и суммировать или же пропускать. в оригинальной версии так и считается:
на С#
//centering routine, return a shift of star from center of input matrix
        private PointF center(Int16[,] input, Int32 offset, float sky)
        {
            Int32 ii, jj;
            Double XSum = 0, YSum = 0, MSum = 0;
            float radius;

            for (ii = 0; ii < 2 * offset; ii++)
            {
                for (jj = 0; jj < 2 * offset; jj++)
                {
                    if ((input[ii, jj] - sky) > 0)
                    {
                        radius = (float)Math.Sqrt((ii - offset) * (ii - offset) + (jj - offset) * (jj - offset));
                        if (radius <= RAper)
                        {
                            XSum = XSum + (ii * (input[ii, jj] - sky));
                            YSum = YSum + (jj * (input[ii, jj] - sky));
                            MSum = MSum + (input[ii, jj] - sky);
                        }
                    }
                }
            }
            PointF result = new PointF((float)((XSum / MSum) - offset), (float)((YSum / MSum)) - offset);
            return result;
        }
на Python
##calculate center of mass
def centroid(R1, R2, R3, arr):
    total=0
    Ry=arr.shape[0]/2
    Rx=arr.shape[1]/2
   
    #mask
    X_index = np.arange(0, arr.shape[1], 1) ##index array
    Y_index = np.arange(0, arr.shape[0], 1) ##index array
    distance = np.sqrt(np.power(np.ones(arr.shape)*(X_index[None, :]-Rx), 2) + np.power(np.ones(arr.shape)*(Y_index[:, None]-Ry), 2)) ##distance array

    ##mean sky
    annulus_mask = np.copy(distance)
    annulus_mask[annulus_mask < R2]=0
    annulus_mask[annulus_mask > R3]=0
    annulus_mask[annulus_mask > 0 ]=1
    masked = arr*annulus_mask
    MSky=np.median(masked[np.nonzero(masked)])

    ##centroid
    aperture_mask = np.copy(distance)
    distance[distance <= R1]=1
    distance[distance > R1] =0
    masked = arr*distance
    total=np.sum(masked)
   
    X = np.sum(masked*X_index[None, :])/total
    Y = np.sum(masked*Y_index[:, None])/total
    return X-arr.shape[1]/2, Y-arr.shape[0]/2, MSky

[krussh]

1. Запускаем Максима с требуемыми параметрами (обрезка снимков, поиск макс. пикселя и т.д.) и передаём полученные результаты в свою программу.

лучше бы избегать платных программ. не то чтобы совесть мучает, но как-то не хорошо)

[krussh]

krussh с самого начала подчёркивал, что требования его программы к снимкам специфические: горячие пиксели д.б. бледнее полезного сигнала. Это - недостаток настоящей версии его программы, т.к. обычно "горячие" пиксели всегда - самые яркие на снимке, особенно слабых объектов. В своей программе я пытаюсь учесть это.

Для задачи seeing monitor это не существенно. Точно определить центроид все равно можно только для изображений с хорошим С/Ш. Ну и Полярная всегда яркая).
Если горячие пиксели начинают мешать, то лучше добавить простое вычитание темнового кадра, чем сложной математики (которая вообще то есть - алгоритм LACosmic, например).

[ekvi]

во время установки требуется подключение к Интернет.

- очень нежелательное требование, поскольку творить спокойней в келье, изолированной от внешних "флюидов".

там же тестовая программка в которой используются все основные пакеты и тестовый файл.

Никак не могу зацепить и скачать из Вашей папки в своём браузёре "Опера", который берет только архивные файлы.