То есть у тебя доп. массив? Это зря... Постарайся от него уйти, например юзать массив char** указателей на существующие строки.

У меня обычное изображение типа float* (одномерный массив).

Дополнительный массив используется для вычисления медианы. Без него не обойтись, т.к. во-первых, алгоритм quick_select переставляет данные в массиве; а во-вторых, если вместо промежуточного массива float* использовать массив float**, т.е. указателей на пиксели изображения, то, сами понимаете, функция будет выполняться еще медленнее, т.к. разыменование указателя - довольно дорогая операция.

Дополнительный массив используется для вычисления медианы.

Это понятно.

Без него не обойтись,

4.2. У тебя основной затуп в копировании в этот доп. массив. И копирование и жрет кучу времени.

т.к. во-первых, алгоритм quick_select переставляет данные в массиве;

Поменяй quick_select, чтобы он принимал float** (сори, что написал char), который ссылается на куски входного массива. Так у тебя не будет копирования, а будет работа только с входным массивом. Оверхед на разыменование - ничтожен по сравнению с копированием.

4.2. У тебя основной затуп в копировании в этот доп. массив.

По сравнению с длительностью quick_select, это копирование - копейки.

Поменяй quick_select, чтобы он принимал float** (сори, что написал char), который ссылается на куски входного массива

Нельзя: quick_select изменяет входные данные, поэтому входной массив в этом алгоритме копируется в новую область памяти. Можно сделать, чтобы quick_select менял местами не данные, а указатели на них (по времени, вроде бы, никакой разницы не будет). Дополнять массив указателей, по идее, будет так же затратно, как дополнять массив указателей. Так что, выгоды от этого явно не получится.

> 127секунд для 3Кх3К с фильтром 32х32

исходя из гигагерцовой частоты (навскидку так)

[code]

1e+9*127/(3000*3000*32*32)

13.780381944444445 [/code]

13 тактов на обрабатываемое значение. И чего Вы хотите от этой жизни? Ну в разы можно ускорить... и то не факт. Возьмите окно поменьше, или другой фильтр;-)

Хотя (мысли вслух) гипотетически при движении по одной строке можно попытаться как то использовать результат сортировки пред. отсчета... т.е. нужно хранить данные окна в таком хитром виде, что бы можно было быстро выкинуть ненужный столбец и добавить новый. По идее добавление небольшого числа значений к уже отсортированному списку куда быстрей, чем сортировка списка целиком.

В общем да, вот у нас есть двусвязный список значений в окне, есть массив по столбцам, в котором прописаны адреса элементов списка Вначале строки сортируем список и находим его середину. При сортировке меняются тока связи списка, сами элементы лежат где лежали, поэтому массив со столбцами адресов остается актулальным. Дальше когда идем вдоль строки, на каждом след отсчете выносим ненужные элементы (адреса их мы знаем), только следим что бы не потерять середину, потом докидываем новый столбец и в процессе докидывания смотрим куда уедет середина. Число операций сокращается в ширину окна раз.

http://en.wikipedia.org/wiki/Selection_algorithm, см секцию Linear general selection algorithm - Median of Medians algorithm. Оно паралелится хорошо, даж в книжечке The Art of Concurrency by Clay Breshears его рассматривали.

Но это всё медленно, это всё O(n log(n) ). O(1) получится только с гистограммами, так как ты и нашёл в примере. Поэтому есть смысл подумать о том, нужен ли тебе этот float.

Идея алгоритма в котором «чёрт ногу сломает» весьма проста: для нахождения медианы используется гистограмма (идеологически напоминает алгоритм lsb radix sort). А для пущего ускорения - поддерживается 2 гистограммы: для значений пикселов (256 столбиков) и для значений старших 4-бит этих пикселов (16 столбиков). Меньшая гистограмма используется для ускоренного поиска медианы (те шаг 16 а не 1):

class BlockMedian : public std::unary_function < BlockParam, unsigned char > {
  ptrdiff_t pitch;
public:
  BlockMedian( ptrdiff_t s ) throw() : pitch(s) {};
  unsigned char operator()( BlockParam const & param ) const throw() {
    unsigned char const* curr = param.First();
    unsigned char const* last = param.Last();
    ptrdiff_t w = param.Width(), i, c, 
    fHistogram  [256] = {0}, // full, fine histogram
    cHistogram  [16]  = {0}; // coarse histogram
    ptrdiff_t pos = ptrdiff_t ( last - curr ) * w / (pitch << 1) + 1;
    
    while ( curr != last) {
      for (i = w; i; ) {
        c = curr[--i]; // a kind of magic :)
        ++(fHistogram [c     ]);
        ++(cHistogram [c >> 4]);
      }
      curr += pitch;
    }
    
    for ( w = c = i = 0; (i < 0xff) && (c < pos); )
      if (cHistogram[w]) { c += fHistogram[i]; ++i; }
      else { i += 16; ++w; }
    return i;
  }
};

Оптимизация под процессоры состоит в ускорении в подсчёте гистограмм, используя предыдущие значения. Те при скольжении окна фильтра из гистограммы вычитается гистограмма последнего столбика окна и добавляется гистограмма следующего столбика. Ну и аналогично при скольжении окна вниз на строку.

Вся мега идея в том, чтоб на каждый столбик изображения иметь по гистограмме значений интенсивности в рамках текущего положения фильтра, и этим всем добром при помощи SIMD инструкций оперировать.

Да, это красиво, но будет работать тока для счетного кол-ва значений. Можно конечно округлять и строить гистограмму с шириной столбца скоко-то-там-сотых. Но не факт что это будет быстрей того, чего я предлагаю;-)

13 тактов на обрабатываемое значение

Не совсем так: 1e+9*127/(3000*3000) \approx 14000 тактов на один пиксель, что очень медленно.

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

Вот-вот. Так обычно медианную фильтрацию и ускоряют. Вот только данные-то при упрощенном поиске медианы меняются местами - «выдрать» старые из полуотсортированного массива и «воткнуть» туда новые не так то и просто (накладные расходы сожрут всю пользу).

Основные методы упрощенной фильтрации пользуются гистограммой (что сводит задачу к сложности О(кол-во пикселей в изображении)): намного быстрее извлечь ненужные уже данные из гистограммы, вставить туда нужные и найти середину, чем сортировать. Однако, у меня float (т.к. разрядность АЦП зачастую неизвестна) => с гистограммой не выйдет.

> Не совсем так: 1e+9*127/(3000*3000)

Ну Вы блин даете... я говорю про число тактов на одну чиселку проходящую через процессор, а Вы о своих хотениях. Алгоритм который у Вас сейчас реализован кардинально ускорить можно тока на CUDA (ИМНО должно параллелится идеально), 13 тактов на флот это очень приличный результат.

Гистограмма даст кардинальный выигрыш тока если число точек в окне много больше чем число точек в гистограмме, иначе шило на мыло выходит.

Вариант с «выдиранием» старых данных я Вам выше расписал.

Да, вспомнил ещё вот что: http://stereopsis.com/radix.html

Те можно таки зная конкретное представление float использовать гистограмные трюки и тут!

На CUDA этот алгоритм ускорится не более чем в N/n раз (N - кол-во ядер GPU, n - кол-во ядер CPU), хотя реальные накладные расходы сделают чуть ли не N/(2n).

(но на CUDA я это еще перепишу - я пока общие алгоритмы составляю)

Я-то думал, может есть еще какие-то варианты...

// ваш не подходит, т.к. вы опять забыли про накладные расходы на разыменование указателей + придется содержать еще и массив указателей на эти указатели, чтобы отслеживать данные после сортировки

Хотя на CUDA все будет хучже чем кажется поскоку в сортировке много условий...

Про мой вариант - а Вы попробуйте... все эти доп массивы фигня, порядка 10^3 элементов... число операций уменьшается на порядок, вот что важно.

вы опять забыли про накладные расходы на разыменование указателей

пг'остите что влезаю, но я не верю что Ваш код настолько оптимален чтобы считать «накладные расходы на разыменование указатей»

//чую, чую как пахнет преждевеременная оптимизация

Можете сами посмотреть. Конечно, до оптимизации там далеко, но все-таки, если пользоваться указателями на указатели, то мне придется:

1. при первоначальной инициализации заполнить массив float *tmp[y*w1+x] = &image[y*w+x], чтобы указатели в нем последовательно указывали на пиксели в первоначальном квадрате фильтра
2. заполнить массив указателей на эти указатели float *tmp1[n] = tmp[n]. Этот массив мы будем сортировать
3. провести первый вызов quick_select над массивом tmp1. Извлечь значение медианы.
4. заменить первые h1 указателей в tmp, чтобы они указывали на следующий столбец. «автоматически» изменится содержимое массива tmp1.
5. перейти к п.3.

(w,h - размер изображения; w1,h1 - размер фильтра)

Просто я сомневаюсь, что все эти манипуляции будут быстрее.

// могу, конечно, попробовать проверить - но, думаю, это будет лишняя трата времени

У тебя основная проблема в том, что 3kx3k картинка, состоящая из float'ов, занимает 36 мб - это ни в какой кэш не влезет. Использование нескольких потоков только ухудшает ситуацию с кэшем. Поэтому, в первую очередь, нужно работать над улучшением data locality в алгоритме.

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

Очень полезно твой мешок пикселов размещать на huge pages. Если ядро более-мене свежее (.32+), то у mmap есть флаг MAP_HUGETLB, иначе libhugetlbfs в руки.

Спасибо. Может, и применю ваш совет.

Просто я - совсем не программист, и что такое «huge pages» понятия не имею.

> y*sizex+x0 facepalm.png

*tmp++ = *ptr++;

memcpy?

да и if на каждый пиксел... а вы такты считаете.

memcpy

Не имеет смысла: у меня 32-разрядная система, так что memcpy и мое будет работать абсолютно одинаково.

да и if на каждый пиксел... а вы такты считаете.

А и правда: функцию на 2 можно разделить, спасибо. Правда, сильно производительности это не поднимат.

Просто я - совсем не программист, и что такое «huge pages» понятия не имею.

Штатный размер странички на x86 - 4 кб. Для каждой страницы процессор делает трансляцию её адреса из физического в виртуальный. В кэшах L2 и L3 данные хранятся с физическими адресами, а в L1 - с виртуальными. Каждый раз, когда данные перемещаются между L1 и L2/3, процессер должен вычислить соответствие между физическим и виртуальным адресом. Результат трансляции хранится в кэше TLB, потому что вычисление достаточно накладное. Размер TLB ограничен, поэтому использование большого количества 4к страниц выносит TLB, и процессор вынужден перетранслировать адреса опять.

У процессора есть возможность использовать страницы с бОльшей гранулярностью - 2 или 4 мб. При этом для одной такой страницы всё равно используется один слот в TLB.

А как сделать malloc с флагом HUGE_TABLES (если не лезть в дебри сис.вызовов)?

А как сделать malloc с флагом HUGE_TABLES (если не лезть в дебри сис.вызовов)?

Используй mmap вместо malloc.

> 3kx3k картинка, состоящая из float'ов, занимает 36 мб

Ну и что? Алгоритм то однопроходный, локальность и так высокая. Ище раз - 13 ТАКТОВ НА ЧИСЛО, ЭТО ПРАКТ ТЕОРЕТИЧЕСКИЙ ПРЕДЕЛ ПРОИЗВОДИТЕЛЬНОСТИ. Когда бывает затык с памятью, цифры гораздо хуже...

Использование нескольких потоков только ухудшает ситуацию с кэшем.

Не в этом случае. Гляньте в исходном посте данные профилирования - переход на 4ре потока ускоряет счет больше чем в 3 раза, а когда проблема с кэшем то обычно замедление по числу потоков. И вообще 36Мб... это фигня.

Поэтому, в первую очередь, нужно работать над улучшением data locality в алгоритме.

а чтобы стало понятно про что это гугли про data oriented design

Ну и что? Алгоритм то однопроходный, локальность и так высокая.

Да нифига она не высокая, данные для вычисления медианы для одной точки лежат даже не в одной странице.

Ище раз - 13 ТАКТОВ НА ЧИСЛО, ЭТО ПРАКТ ТЕОРЕТИЧЕСКИЙ ПРЕДЕЛ ПРОИЗВОДИТЕЛЬНОСТИ.

ОБОСНУЙ ЦИФРАМИ.

> Да нифига она не высокая, данные для вычисления медианы для одной точки лежат даже не в одной странице.

И что? Для непосредственной работы нужно число страниц по числу строк в окне фильтрации, это в кэш точно лезет и этого хватает на 1000 точек, так?

ОБОСНУЙ ЦИФРАМИ.

Обосновать ЧТО? 13 тактов обоснованно выше, что 13 тактов это теоретич. предел обосновать? Там на одно число при сортировке log(N) проверок, это раз. Переход на 4ре потока дает выигрыш почти в 4ре раза, а если бы подсистема памяти не справлялась, то давал бы замедление, это два.

> Для непосредственной работы нужно число страниц по числу строк в окне фильтрации,

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

Короче по ПСП памяти там огроменный запас, и ковырять надо именно алгоритм с т.з. уменьшения числа операций. Банально флопсов не хватает... вообще редкая ситуация;-)

Не путайте, используется не алгоритм quick_sort, а алгоритм quick_select. Он не выполняет полную сортировку и предназначен только для поиска медианы.

Я думаю, что можно будет расширить алгоритм: т.е. для малых размеров окна (3x3, 5x5, 7x7 и 9x9) использовать самый быстрый способ на макросах; для размеров от 10x10 до 128x128 - текущий алгоритм; а для бóльших размеров попробовать реализовать алгоритм с плавающей гистограммой (т.е. выделять в гистограмме ограниченное кол-во точек и заполнять/освобождать их по мере пробега окна; гистограмму придется реализовать в виде списка указателей - можно попробовать, вдруг реально скорость получится выше). Хотя, конечно, фильтрация медианой с размером окна больше 16х16 - явление довольно редкое.

Если вдруг понадобятся непрямоугольные окна, всё становится в разы веселей ;)

В этом случае работает этот же самый метод, только выбор координат чуть более сложный. Естественно, немного увеличивается время вычислений. Но вряд ли глазом будет заметно.

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

Ок, quick_select. Он таки все значения окна проходит, так? Сколько операций на значение?

ИМНО оптимальным будет что то вроде бинарного дерева, элементы которого лежат в векторе. После первичного заполнения просто сносите первые n элементов вектора и втыкаете на их место новые n (не забывая сохранять связность и балансировку дерева), где n - размер окна.

И определитесь таки с типовым размером окна - вроде как все эти radix алгоритмы дают профит тока когда у вас число строк в окне много больше чем разрядность float?

Ок, quick_select. Он таки все значения окна проходит, так? Сколько операций на значение?

А черт его знает. Но это быстрее, чем quick_sort.

ИМНО оптимальным будет что то вроде бинарного дерева, элементы которого лежат в векторе. После первичного заполнения просто сносите первые n элементов вектора и втыкаете на их место новые n (не забывая сохранять связность и балансировку дерева), где n - размер окна.

Сейчас попробую с двумя массивами указателей (один - пополняемый, второй - сортируемый).

И определитесь таки с типовым размером окна - вроде как все эти radix алгоритмы дают профит тока когда у вас число строк в окне много больше чем разрядность float?

До radix я еще не дошел, но тоже буду изучать.

> А черт его знает. Но это быстрее, чем quick_sort.

Я думаю что типа log(N/2), разница небольшая.

Не верю своим глазам!

Но использование массива указателей и массива указателей на указатели (второй - сортируется) дало прирост производительности для больших N:

 *              45секунд для 3Кх3К с фильтром 41x41
 *              27секунд для 3Кх3К с фильтром 32х32
 *              7.3секунд для 3Кх3К с фильтром 16x16
 *              1.5секунд 3Кх3К с фильтром 6x6
 *              1.1секунда 3Кх3К с фильтром 5х5 (med_opt)

В 4ре раза, неожиданно, да?;-) Существует огромное кол-во мифов об оптимизации в головах программистов...

Попробуйте таки помудрить с бинарным деревом. И чего то я не понял, почему у вас сортируется ** ...

В 4ре раза, неожиданно, да?;-)

Не в 4раза, а где-то на 30% (например, было 40с, стало 27).

С бинарным деревом неохота заморачиваться. Хочу еще попробовать метод с двухуровневой гистограммой. Он, правда, иногда будет двухпроходным (если полоса детализирующей гистограммы не будет совпадать со значением грубой), но, вполне возможно, что на больших N будет быстрее.

И чего то я не понял, почему у вас сортируется ** ...

**arr - указатели на пиксели внутри окна фильтра; ***sel - указатели на **arr. Сортировка меняет местами указатели sel, не изменяя массива arr, что позволяе при продвижении на строку вниз заменять только часть arr. А т.к. он уже частично подготовлен, получаем, как вы и говорили, выигрыш в производительности.

Сейчас попробую сократить один *, может, быстрее будет...

// кстати, я free(sel) забыл - добавил уже.

Да, с тройным указателем я намудрил... Подсократил на 1 указатель, получил прирост производительности еще в 10% :)

Закоммитил

> Не в 4раза, а где-то на 30% (например, было 40с, стало 27).

А, так оно у Вас в ч потока... так бы и писали. Лучше тестите сначала один поток, потом 4ре, а то неясно какой профит от потоков, какой от изменения алгоритма. Ну и 30-40% - это не то ускорение, за которое стоит так истово бороться.

Ну, все равно приятно, что чуть быстрее. Правда, на больших размерах фильтра (256х256) получилось почти в 2 раза медленнее. Сижу, придумываю, как бы так с гистограммой...

мм, чтобы не плодить тем, а не приходилось ли кому нагуглить какого-нить нелобового риалтаймового алгоритма для векторного медианного фильтра?

матчасть: http://demonstrations.wolfram.com/VectorMedianFilter/ на пальцах: векторная медиана в метрическом пространстве есть точка, сумма растояний от которой до остальных минимальна.

А возьмите Вы линейный фильтр... идеологически похоже, рабоать будет на порядки быстрее?

Это же изображение! Сигналы нельзя фильтровать средним арифметическим.

Я в этом плохо разбираюсь. А если линейный фильтр с каким то пральным шаблоном?

Любые линейные фильтры можно применять только после медианной фильтрации. Если нет желания уменьшать разрешение, обычно делают несколько экспозиций, а потом находят их общую медиану.

Ну а дифференциальные фильтры работают довольно быстро: например, модуль градиента (два Собеля + вычисление sqrt(sob1^2+sob2^2) ) занимает чуть больше двух секунд (это на CPU).

Ладно, промучился я с гистограммами и понял, что ну их нафиг. И так скорость на небольших размерах окна достаточно адекватна, а большие размеры используются крайне редко.