О степенной функции, OpenMP и прочих граблях

Есть вот такое вот известное выражение, переложеное из википедии:

В том смысле, что у sRGB внизу - линейный участок, а дальше степенная функция.

Если его напрямую запрограммировать (вот прямо с powf()), то в зависимости от компилятора получается от 7 Mpix/sec (Visual C++ 2010, один поток) до 100Mpix/sec (Intel Parallel Studio XE, включен OpenMP). Пиксель - это 3 компонента по 4 байта (float), то бишь гигабайт с небольшим в секунду - это максимум.

Как-то медленно.

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

Табличный лукап

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

Линейная интерполяция

Можно брать из таблицы два значения и линейно интеполировать, прямо как велит нам Klaus Post. Для таблицы из 512 записей получается быстро, но неточно (200-210 Mpix/sec при максимальной ошибке около 0.025), а для таблицы в 64k записей - медленнее и тоже не очень точно (150Mpix/sec, максимальная ошибка 3-4x10^-4). Растить таблицу больше 64к записей (512kb памяти с учетом удвоения ради выравнивания, как велит нам Клаус) не хочется совсем.

Кубические сплайны

Кубические сплайны оказались отличным компромиссом: для таблицы в 4k записей (64 килобайта памяти) получаем 80 мегапикселей в секунду для реалиазации на C++, 120 Mpix/sec для SSE-варианта и ошибку не более чем 3.5x10^-6. Увеличение размера таблицы не нужно, уменьшение тоже на скорость уже почти не влияет, а на точность - влияет очень сильно.

OpenMP

Все цифры выше (кроме одной, в самом начале, где это оговорено) - это в один поток. OpenMP любит нас, мы ее, режем входной поток на куски по несколько десятков мегабайт, #pragma omp parallel for, запускаем, получаем 360-400Mpix/sec для SSE-реализации сплайнов на компиляторе Intel и четырех ядрах.

Правда вылезают интересные грабли, ради которых, собственно, пост и написан (для собственной памяти):

1. Режим округления у SSE сбивается при инициализации потока. То есть _MM_SET_ROUNDING_MODE надо звать уже внутри OpenMP-блока. Убил бы нахрен!
2. У Visual Studio 2010 с OpenMP все не так радужно. OpenMP есть, считает правильно, только вот OpenMP-версия не в 3-4 раза быстрее однопоточной, а в полтора-два раза медленнее.

Заметим, что 4-4.8 гигабайта в секунду - это быстрее, чем залить картинку в видеокарту, посчитать (допустим, мгновенно) и вылить оттуда результат. Я к тому, что на внешнем сопроцессоре мелкими кусками обрабатывать смысла нет, туда надо нести сразу много кода.

Да, а "в 4 раза быстрее простой реализации" - очень похоже на константу, не первый раз ее уже вижу.