http://ndevilla.free.fr/median/median/index.html и внизу статьи ссылки на разные реализации, типа http://ndevilla.free.fr/median/median/src/quickselect.c

+ http://ndevilla.free.fr/median/median.pdf

Но я так, мимопроходил, ничего не знаю по теме

Пиши лучше алгоритм адаптации. Он, насколько я понимаю, работает независимо и только определяет размер окна. А там видно будет, стоит ли оптимизировать медиану или нет. Может, даже если ты потратишь кучу времени и оптимизируешь её на пару секунд суммарно, это ничего не даст по сранению с минутами, которые требуются для алгортма адаптации.

quickselect.c

Это я раньше использовал. Значительно медленней работает, чем скользящий алгоритм (т.к. каждый раз надо заново формировать массив по окну).

Проблема именно в реализации адаптивного медианного фильтра. Уж очень не хочется велосипедить 100500-й велосипед.

Кстати, чего-то я забыл про openMP. Сейчас добавил, веселей получилось: 0.5, 1.3 и 2.7 секунд для тех же изображений.

Алгоритм адаптации простой: 1. делаем копию изображения (I) с медианным фильтром в виде «креста» 3х3 элемента (Imed) + еще два изображения — минимум (Imin) по «кресту» и максимум (Imax) 2.1. если Imed == Imin или Imed == Imax, с данным пикселем нужно работать отдельно, увеличивая размер окна (здесь уже скользящая медиана бессмысленна, можно обычным opt_medX); если размер достиг предела, а так ничего и не вышло, оставляем значение I 2.2. иначе если Imin < I < Imax, то в новое изображение заносим I, в противном случае заносим Imed

Вот я и думаю, что если велосипедить, то нужно попробовать 2 варианта: либо этот, либо же отдельно посчитать медианы по кресту 3х3, квадрату 3х3, квадрату 5х5 и квадрату 7х7 (вместе с их минимаксами), затем уже работать с готовыми изображениями (тоже задача параллелизуется нормально).

Мне кажется, тебе подойдёт «running median при помощи двух куч» в текущей точке. Вот здесь, вроде бы, именно оно. Ты ему даёшь поток чисел, которые захватывает всё более увеличивающееся окно, он быстро и эффективно даёт тебе текущее значение медианы без пересчёта заново с каждым увеличением окна. Минимумы-максимумы тоже считаешь в потоке, это несложно вроде.

ХЗ, сколько у тебя возможных вариантов окон, но считать все значения для всех возможных окон от самого большого до самого маленького — чудовищно, по-моему.

Почему именно адаптивный? Ты купился на красивую картинку из книжечки?

По поводу неадаптивного - все нормальные пацаны делают так: https://nomis80.org/ctmf.html (Median Filtering in Constant Time)

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

Максимум и минимум в окне - это стандартные морфологические операции, как их за константу реализовать см лептонику.

Вот здесь, вроде бы, именно оно.

Не посмотрел внимательно. Это, во-первых, то, что ты давал в оп-посте, и, во-вторых, не оно. Это не «текущая», а «скользящая» медиана только для К последних значений («текущая» она лишь первые К значений, но инициализируется долго, чтобы это на каждой новой точке делать). Её, наверное, можно переделать в чисто «текущую», но готового кода нет.

Почему именно адаптивный?

Потому что не размывает границы.

Идея в использовании гистограмм

И где я возьму терабайты оперативы для гистограммы в double?

см лептонику.

Чудовищно медленная. Я те же "шагающие квадраты" шустрей реализовал.

Черт!

Я совсем забыл, что даже адаптивная медиана «сжирает» локальные максимумы. Для нормальной фильтрации (чтобы можно было по этим данным потом фотометрию делать) нужно будет метод маленько усложнить. Скажем, сделать обычную медианную фильтрацию 3х3 и по ней посчитать градиент (или даже лапгаусса), да проверить, где есть явные выбросы. В точках выбросов изображение и восстанавливать. Причем, наверное, лучше даже B-сплайны натянуть (благо, их есть у меня). Заодно можно будет и фон вычислить да вычесть…

P.S. Пока что я не заморачивался с краями. Квазиадаптивная (я не расширяю окно, если медиана равна локальному экстремуму, позже исправлю это) медианная фильтрация на CPU для изображения 1к х 1к дает:
0.05с с окном 3х3
0.1с с окном 5х5
0.2с с окном 7х7
0.25с с окном 9х9
...
1с с окном 25х25
9с с окном 101х101

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

В общем, вот так сделал. Пока что края не обрабатываются → надо будет исправить это. Ну и еще недостаток — медианная адаптивная фильтрация неполная: т.к. на реальных изображениях необходимость увеличить размер окна возникает лишь иногда при фильтрации крестом 3x3, я просто сделал обработку конкретных «плохих» пикселей для данного случая окном 5х5 пикселей.

Про границы еще подумаю, как получше бы сделать.

а как же mc? говорим одно, а пользуем другое, так получается?

При чем здесь mc? Скриншоты из моего велосипеда сделаны. Цветные картинки — это превьюшки в окне «Open FITS file».

Почитал, подумал, и нарисовал такую схему для предотвращения изменений локальных максимумов:

(цикл по пикселям изображения: Ixy)
start: wsize = 0 // окно в виде креста 3х3 пикселя
0. (wsize < wmax) ? goto 1 : (x > 0 ? Ixy = Ix-1y && continue : continue) // значение пикселя слева или без изменения
1. calculate Imin, Imax
2. (Imin < Ixy < Imax) ? continue : goto 3
3. calculate Imed, sigmaI
4. (Imin < Imed < Imax) ? goto 5 : ++wsize && goto 0
5. (Imax - Imin < 3*sigmaI) ? continue : Ixy = Imed && continue

Т.о., вычисление медианы понадобится лишь для обработки локальных экстремумов. И то, если диапазон пикселей в локальной окрестности входит в три сигмы, значение остается нетронутым. Если же выбивается, то значение считается выбросом и заменяется на медиану. Увеличение размера окна необходимо для обработки областей с перекопами. В моем случае эту обработку лучше не делать, поэтому пункт 0 выбрасывается, а пункт 4 принимает вид
4. (Imin < Imed < Imax) ? goto 5 : (x > 0 ? Ixy = Ix-1y && continue : continue)