В пилотной части я рассказал о задаче как можно подробнее. Рассказ получился долгим и беспредметным — в нем не было ни одной строчки кода. Но без понимания задачи очень сложно заниматься оптимизацией. Конечно, некоторые техники можно применять, имея на руках только код. Например, кешировать вычисления, сокращать ветвления. Но мне кажется, что некоторые вещи без понимания задачи просто никогда не сделать. Это и отличает человека от оптимизирующего компилятора. Поэтому ручная оптимизация все еще играет огромную роль: у компилятора есть только код, а у человека есть понимание задачи. Компилятор не может принять решение, что значение "4" достаточно случайно, а человек может.
Напомню, что речь пойдет об оптимизации операции ресайза изображения методом сверток в реально существующей библиотеке Pillow. Я буду рассказывать о тех изменениях, что я делал несколько лет назад. Но это не будет повторение слово-в-слово: оптимизации будут описаны в порядке, удобном для повествования. Для этих статей я сделал в репозитории отдельную ветку от версии 2.6.2 — именно с этого момента и будет идти повествование.
Тестирование
Если вы хотите не просто читать, но и экспериментировать самостоятельно, вам пригодится пакет с тестами pillow-perf.
# Установка нужных для компиляции и тестов пакетов
$ sudo apt-get install -y build-essential ccache python-dev libjpeg-dev
$ git clone -b opt/scalar https://github.com/uploadcare/pillow-simd.git
$ git clone --depth 10 https://github.com/python-pillow/pillow-perf.git
$ cd ./pillow-simd/
# Коммит, с которого все начинается
$ git checkout bf1df9a
# Собираем и устанавливаем Pillow
$ CC="ccache cc" python ./setup.py develop
# Наконец-то запускаем тест
$ ../pillow-perf/testsuite/run.py scale -n 3Т. к. Pillow состоит из множества модулей и не умеет компилироваться инкрементально, для существенного ускорения повторных сборок используется утилита ccache. С помощью pillow-perf можно тестировать много операций, но нас интересует scale. -n 3 задает кол-во запусков операции. Пока код медленный, можно взять число поменьше, чтобы не уснуть. На старте производительность такая:
Scale 2560?1600 RGB image
to 320x200 bil 0.08927 s 45.88 Mpx/s
to 320x200 bic 0.13073 s 31.33 Mpx/s
to 320x200 lzs 0.16436 s 24.92 Mpx/s
to 2048x1280 bil 0.40833 s 10.03 Mpx/s
to 2048x1280 bic 0.45507 s 9.00 Mpx/s
to 2048x1280 lzs 0.52855 s 7.75 Mpx/s
to 5478x3424 bil 1.49024 s 2.75 Mpx/s
to 5478x3424 bic 1.84503 s 2.22 Mpx/s
to 5478x3424 lzs 2.04901 s 2.00 Mpx/sРезультат для коммита bf1df9a.
Эти результаты немного отличаются от полученных в официальном бенчмарке для версии 2.6. Тому есть несколько причин:
- В официальном бенчмарке используется 64-битная Ubuntu 16.04 с GCC 5.3. Я же буду использоваться 32-битной Ubuntu 14.04 с GCC 4.8, на которой делал все эти оптимизации впервые. В конце статьи станет понятно, почему.
- В статьях я начинаю рассказ с коммита, в котором исправлен баг, который не относится к оптимизации, но влияет на производительность.
Структура кода
Большая часть кода, который нас интересует, находится в файле Antialias.c, в функции ImagingStretch. Код этой функции можно разделить на три части:
// Пролог
if (imIn->xsize == imOut->xsize) {
// Вертикальный ресайз
} else {
// Горизонтальный ресайз
}Как я говорил ранее, ресайз изображения свертками можно сделать в два прохода: на первом меняется только ширина изображения, на втором — высота, или наоборот. Функция ImagingStretch за один вызов как раз может сделать либо то, либо другое. Тут можно увидеть, что она действительно вызывается дважды во время каждого ресайза. В функции выполняется общий пролог, а дальше, в зависимости от параметров, та или иная операция. Довольно необычный подход к избавлению от повторяющегося кода, коим в данном случае является пролог.
Внутри оба прохода выглядят примерно одинаково, с поправкой на меняющееся направление обработки. Для краткости приведу только один:
for (yy = 0; yy < imOut->ysize; yy++) {
// Подсчёт коэффициентов
if (imIn->image8) {
// Цикл для одного канала 8 бит
} else {
switch(imIn->type) {
case IMAGING_TYPE_UINT8:
// Цикл для нескольких каналов по 8 бит
case IMAGING_TYPE_INT32:
// Цикл для одного канала 32 бита
case IMAGING_TYPE_FLOAT32:
// Цикл для одного канала float
}
}Тут происходит ветвление на несколько форматов представления пикселей, поддерживаемых Pillow: одноканальные 8 бит (оттенки серого), многоканальные 8 бит (RGB, RGBA, LA, CMYK, некоторые другие), одноканальные 32 бита и float. Нас будет интересовать тело цикла для нескольких каналов по 8 бит, так как это наиболее часто встречающийся формат изображений.
Оптимизация 1: эффективно используем кэш
Хоть я и сказал выше, что два прохода похожи, между ними есть очевидное различие. Посмотрите на вертикальный проход:
for (yy = 0; yy < imOut->ysize; yy++) {
// Подсчёт коэффициентов
for (xx = 0; xx < imOut->xsize*4; xx++) {
// Свертка столбца пикселей
// Сохранение пикселя в imOut->image8[yy][xx]
}
}На горизонтальный проход:
for (xx = 0; xx < imOut->xsize; xx++) {
// Подсчёт коэффициентов
for (yy = 0; yy < imOut->ysize; yy++) {
// Свертка строки пикселей
// Сохранение пикселя в imOut->image8[yy][xx]
}
}У вертикального прохода во внутреннем цикле итерируются столбцы конечного изображения, а у горизонтального — строки. Горизонтальный проход представляет серьезную проблему для кеша процессора. На каждом шаге внутреннего цикла происходит обращение на одну строку ниже, а значит запрашивается значение из памяти, лежащее далеко от значения, нужного на предыдущем шаге. Это плохо при небольшом размере свертки. Дело в том, что на современных процессорах линейка кеша, которую процессор может запросить из оперативной памяти, всегда равна 64 байтам. Это значит, что если в свертке участвует меньше 16 пикселей, то часть данных гоняется впустую из оперативной памяти в кеш. А теперь представьте, что циклы поменялись местами и следующий пиксель сворачивается не строкой ниже, а следующий в этой же строке. Тогда большая часть нужных пикселей уже была бы в кеше.
Второй негативный фактор такой организации кода проявляется при длинной строке для свертки (т. е. при сильном уменьшении). Дело в том, что для соседних сверток исходные пиксели довольно значительно пересекаются, и было бы хорошо, если бы эти данные оставались в кеше. Но когда мы движемся сверху вниз, данные для старых сверток постепенно начинают вытесняться из кэша данными для новых. В результате, когда проходит полный внутренний цикл и начинается следующий шаг внешнего, в кеше уже не оказывается верхних строк, они все вытеснены нижними, их нужно снова запрашивать из памяти. А когда дело доходит до нижних, в кеше уже все вытеснено верхними. Получается круговорот, в результате которого в кеше никогда нет нужных данных.
Почему же обход сделан так? В псевдокоде выше видно, что второй строчкой в обоих случаях идет подсчет коэффициентов для свертки. Для вертикального прохода коэффициенты зависят только от строки конечного изображения (от значения yy), а для горизонтального от текущего столбца (от значения xx). Т. е. в горизонтальном проходе нельзя просто поменять местами два цикла — расчет коэффициентов должен быть внутри цикла по xx. Если же начать считать коэффициенты внутри внутреннего цикла, это убьёт всю производительность. Особенно когда для расчета коэффициентов применяется фильтр Ланцоша, в котором есть тригонометрические функции.
Вычислять коэффициенты на каждом шаге нельзя, но тем не менее они одинаковые для всех пикселей одного столбца. Значит можно просчитать все коэффициенты заранее для всех столбцов, а во внутреннем цикле брать уже посчитанное. Давайте так и сделаем.
В коде есть выделение памяти для коэффициентов:
k = malloc(kmax * sizeof(float));Теперь нам понадобится массив таких массивов. Но можно упростить — выделить плоский кусок памяти и эмулировать двумерность адресацией.
kk = malloc(imOut->xsize * kmax * sizeof(float));Еще понадобится где-то хранить xmin и xmax, которые тоже зависят от xx. Под них тоже сделаем массив, чтобы не пересчитывать.
xbounds = malloc(imOut->xsize * 2 * sizeof(float));Также внутри цикла используется некое значение ww, которое нужно для нормализации значения свертки. ww = 1 / ?k[x]. Его можно вообще не хранить, и нормализовать сами коэффициенты, а не результат свёртки. Т. е. после того, как мы посчитали коэффициенты, нужно пройтись по ним еще раз и поделить их на их же сумму, в результате чего сумма всех коэффициентов станет равна 1:
k = &kk[xx * kmax];
for (x = (int) xmin; x < (int) xmax; x++) {
float w = filterp->filter((x - center + 0.5) * ss);
k[x - (int) xmin] = w;
ww = ww + w;
}
for (x = (int) xmin; x < (int) xmax; x++) {
k[x - (int) xmin] /= ww;
}Вот теперь можно наконец-то развернуть обход пикселей на 90°:
// Подсчёт всех коэффициентов
for (yy = 0; yy < imOut->ysize; yy++) {
for (xx = 0; xx < imOut->xsize; xx++) {
k = &kk[xx * kmax];
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
// Свертка строки пикселей
// Сохранение пикселя в imOut->image8[yy][xx]
}
}Scale 2560?1600 RGB image
to 320x200 bil 0.04759 s 86.08 Mpx/s 87.6 %
to 320x200 bic 0.08970 s 45.66 Mpx/s 45.7 %
to 320x200 lzs 0.11604 s 35.30 Mpx/s 41.6 %
to 2048x1280 bil 0.24501 s 16.72 Mpx/s 66.7 %
to 2048x1280 bic 0.30398 s 13.47 Mpx/s 49.7 %
to 2048x1280 lzs 0.37300 s 10.98 Mpx/s 41.7 %
to 5478x3424 bil 1.06362 s 3.85 Mpx/s 40.1 %
to 5478x3424 bic 1.32330 s 3.10 Mpx/s 39.4 %
to 5478x3424 lzs 1.56232 s 2.62 Mpx/s 31.2 %Результат для коммита d35755c.
В четвертом столбце показано ускорение относительно предыдущего варианта, а под таблицей ссылка на коммит, где хорошо видно сделанные изменения.
Оптимизация 2: Ограничение выходных значений
В коде мы в нескольких местах видим такую конструкцию:
if (ss < 0.5)
imOut->image[yy][xx*4+b] = (UINT8) 0;
else if (ss >= 255.0)
imOut->image[yy][xx*4+b] = (UINT8) 255;
else
imOut->image[yy][xx*4+b] = (UINT8) ss;Это ограничение значения пикселя в рамках [0, 255], если результат вычислений выходит за пределы 8 бит. Действительно, сумма всех положительных коэффициентов свертки может быть больше единицы, а сумма всех отрицательных может быть меньше нуля. А значит при определенных исходном изображении можем получить переполнение. Это переполнение является результатом компенсации резких перепадов яркости и не является ошибкой.
Взглянем на код. В нем одна входная переменная ss и ровно одна выходная imOut->image[yy], значение которой присваивается в нескольких местах. Плохо тут то, что сравниваются числа с плавающей точкой. Было бы быстрее преобразовать все в целые и потом уже сравнивать, раз уж все равно в конце нам нужен целый результат. Итого получается вот такая функция:
static inline UINT8
clip8(float in) {
int out = (int) in;
if (out >= 255)
return 255;
if (out <= 0)
return 0;
return (UINT8) out;
}Использование:
imOut->image[yy][xx*4+b] = clip8(ss);Это тоже дает прирост производительности, хоть и небольшой.
Scale 2560?1600 RGB image
to 320x200 bil 0.04644 s 88.20 Mpx/s 2.5 %
to 320x200 bic 0.08157 s 50.21 Mpx/s 10.0 %
to 320x200 lzs 0.11131 s 36.80 Mpx/s 4.2 %
to 2048x1280 bil 0.22348 s 18.33 Mpx/s 9.6 %
to 2048x1280 bic 0.28599 s 14.32 Mpx/s 6.3 %
to 2048x1280 lzs 0.35462 s 11.55 Mpx/s 5.2 %
to 5478x3424 bil 0.94587 s 4.33 Mpx/s 12.4 %
to 5478x3424 bic 1.18599 s 3.45 Mpx/s 11.6 %
to 5478x3424 lzs 1.45088 s 2.82 Mpx/s 7.7 %Результат для коммита 54d3b9d.
Как можно заметить, эта оптимизация дала больший эффект для фильтров с меньшим окном и с большим выходным разрешением (единственное исключение это 320x200 Bilinear, но я не берусь сказать, почему). И действительно, чем меньше окно фильтра и больше конечное разрешение, тем больший вклад в производительность делает обрезание значений, которое мы и оптимизировали.
Оптимизация 3: Разворот циклов c постоянным количеством итераций
Если еще раз присмотреться к горизонтальному шагу, там можно насчитать целых 4 вложенных цикла:
for (yy = 0; yy < imOut->ysize; yy++) {
// ...
for (xx = 0; xx < imOut->xsize; xx++) {
// ...
for (b = 0; b < imIn->bands; b++) {
// ...
for (x = (int) xmin; x < (int) xmax; x++) {
ss = ss + (UINT8) imIn->image[yy][x*4+b] * k[x - (int) xmin];
}
}
}
}Итерируется каждый ряд и каждый столбец выходного изображения (т. е. каждый пиксель), а внутри итерируется каждый пиксель исходного изображения, который нужно свернуть. А вот что такое b? b — это итерирование по каналам изображения. Очевидно, что количество каналов на протяжении работы функции не меняется и никогда не превышает 4 (ввиду способа хранения изображения в Pillow). Следовательно, может быть всего 4 возможных случая. А с учетом того, что одноканальные 8-битные изображения хранятся иным способом, то случаев вовсе три. Соответственно можно сделать три отдельных внутренних цикла: для двух, трех и четырех каналов. И сделать ветвление на нужное количество каналов. Я покажу только код трехканального случая, чтобы не занимать слишком много места.
for (xx = 0; xx < imOut->xsize; xx++) {
if (imIn->bands == 4) {
// Тело для 4 каналов
} else if (imIn->bands == 3) {
ss0 = 0.0;
ss1 = 0.0;
ss2 = 0.0;
for (x = (int) xmin; x < (int) xmax; x++) {
ss0 = ss0 + (UINT8) imIn->image[yy][x*4+0] * k[x - (int) xmin];
ss1 = ss1 + (UINT8) imIn->image[yy][x*4+1] * k[x - (int) xmin];
ss2 = ss2 + (UINT8) imIn->image[yy][x*4+2] * k[x - (int) xmin];
}
ss0 = ss0 * ww + 0.5;
ss1 = ss1 * ww + 0.5;
ss2 = ss2 * ww + 0.5;
imOut->image[yy][xx*4+0] = clip8(ss0);
imOut->image[yy][xx*4+1] = clip8(ss1);
imOut->image[yy][xx*4+2] = clip8(ss2);
} else {
// Тело для двух и одного канала
}
}Можно на этом не останавливаться и переместить ветвление еще на один уровень выше, до цикла по xx:
if (imIn->bands == 4) {
for (xx = 0; xx < imOut->xsize; xx++) {
// Тело для 4 каналов
}
} else if (imIn->bands == 3) {
for (xx = 0; xx < imOut->xsize; xx++) {
// Тело для 3 каналов
}
} else {
for (xx = 0; xx < imOut->xsize; xx++) {
// Тело для двух и одного канала
}
}Scale 2560?1600 RGB image
to 320x200 bil 0.03885 s 105.43 Mpx/s 19.5 %
to 320x200 bic 0.05923 s 69.15 Mpx/s 37.7 %
to 320x200 lzs 0.09176 s 44.64 Mpx/s 21.3 %
to 2048x1280 bil 0.19679 s 20.81 Mpx/s 13.6 %
to 2048x1280 bic 0.24257 s 16.89 Mpx/s 17.9 %
to 2048x1280 lzs 0.30501 s 13.43 Mpx/s 16.3 %
to 5478x3424 bil 0.88552 s 4.63 Mpx/s 6.8 %
to 5478x3424 bic 1.08753 s 3.77 Mpx/s 9.1 %
to 5478x3424 lzs 1.32788 s 3.08 Mpx/s 9.3 %Результат для коммита 95a9e30.
Что-то похожее можно сделать и для вертикального прохода. Там сейчас такой код:
for (xx = 0; xx < imOut->xsize*4; xx++) {
/* FIXME: skip over unused pixels */
ss = 0.0;
for (y = (int) ymin; y < (int) ymax; y++)
ss = ss + (UINT8) imIn->image[y][xx] * k[y-(int) ymin];
ss = ss * ww + 0.5;
imOut->image[yy][xx] = clip8(ss);
}Тут нет отдельного итерирования по каналам, вместо этого xx итерируется по ширине, умноженной на 4. Т. е. xx проходится по каждому каналу вне зависимости от их количества в изображении. FIXME в комментарии как раз говорит о том, что это нужно исправить. Исправляется точно так же — ветвлением кода для разного количества каналов в исходном изображении. Приводить код здесь не буду, ссылка на коммит ниже.
Scale 2560?1600 RGB image
to 320x200 bil 0.03336 s 122.80 Mpx/s 16.5 %
to 320x200 bic 0.05439 s 75.31 Mpx/s 8.9 %
to 320x200 lzs 0.08317 s 49.25 Mpx/s 10.3 %
to 2048x1280 bil 0.16310 s 25.11 Mpx/s 20.7 %
to 2048x1280 bic 0.19669 s 20.82 Mpx/s 23.3 %
to 2048x1280 lzs 0.24614 s 16.64 Mpx/s 23.9 %
to 5478x3424 bil 0.65588 s 6.25 Mpx/s 35.0 %
to 5478x3424 bic 0.80276 s 5.10 Mpx/s 35.5 %
to 5478x3424 lzs 0.96007 s 4.27 Mpx/s 38.3 %Результат для коммита f227c35.
Как можно заметить, горизонтальный проход дал больше производительности для уменьшения картинки, тогда как вертикальный — для увеличения.
Оптимизация 4: Целочисленные счетчики
for (y = (int) ymin; y < (int) ymax; y++) {
ss0 = ss0 + (UINT8) imIn->image[y][xx*4+0] * k[y-(int) ymin];
ss1 = ss1 + (UINT8) imIn->image[y][xx*4+1] * k[y-(int) ymin];
ss2 = ss2 + (UINT8) imIn->image[y][xx*4+2] * k[y-(int) ymin];
}Если посмотреть в самый внутренний цикл, можно увидеть, что переменные ymin и ymax объявлены как float, но приводятся к целому на каждом шаге. Более того, за пределами цикла для присвоения им значений используются функции floor и ceil. Т. е. по факту в переменных всегда хранятся целые числа, но объявлены они почему-то как float. То же самое касается xmin и xmax. Меняем и замеряем.
Scale 2560?1600 RGB image
to 320x200 bil 0.03009 s 136.10 Mpx/s 10.9 %
to 320x200 bic 0.05187 s 78.97 Mpx/s 4.9 %
to 320x200 lzs 0.08113 s 50.49 Mpx/s 2.5 %
to 2048x1280 bil 0.14017 s 29.22 Mpx/s 16.4 %
to 2048x1280 bic 0.17750 s 23.08 Mpx/s 10.8 %
to 2048x1280 lzs 0.22597 s 18.13 Mpx/s 8.9 %
to 5478x3424 bil 0.58726 s 6.97 Mpx/s 11.7 %
to 5478x3424 bic 0.74648 s 5.49 Mpx/s 7.5 %
to 5478x3424 lzs 0.90867 s 4.51 Mpx/s 5.7 %Результат для коммита 57e8925.
Финал первого акта и босс
Признаюсь, я был очень рад полученным результатам. Удалось разогнать код в среднем в 2,5 раза. Причем для получения этого ускорения пользователю библиотеки не нужно было ставить дополнительное оборудование, ресайз как и раньше идет на одном ядре того же процессора, что и до этого. Нужно было только обновить версию Pillow до версии 2.7.
Но до релиза 2.7 оставалось еще какое-то время, а мне не терпелось проверить новый код на том сервере, на котором он должен был работать. Я перенес код, скомпилировал и сначала подумал, что что-то напутал:
Scale 2560?1600 RGB image
320x200 bil 0.08056 s 50.84 Mpx/s
320x200 bic 0.16054 s 25.51 Mpx/s
320x200 lzs 0.24116 s 16.98 Mpx/s
2048x1280 bil 0.18300 s 22.38 Mpx/s
2048x1280 bic 0.31103 s 13.17 Mpx/s
2048x1280 lzs 0.43999 s 9.31 Mpx/s
5478x3424 bil 0.75046 s 5.46 Mpx/s
5478x3424 bic 1.22468 s 3.34 Mpx/s
5478x3424 lzs 1.70451 s 2.40 Mpx/sРезультат для коммита 57e8925. Получен на другой машине и не участвует в сравнении.
Лолчто? Результаты почти такие же, как до оптимизации. Я 10 раз все перепроверил, поставил принтов, чтобы убедиться, что работает нужный код. Это не был какой-то сайд-эффект от Pillow или окружения, разница воспроизводилась даже на минимальном примере в 30 строчек. Я задал вопрос на Stack Overflow, и в конце концов удалось найти явную закономерность: код выполнялся медленно, если был скомпилирован с GCC под платформу 64 бита. И именно в этом было различие локальной Убунты и на сервере: локально была 32-битная.
Ну слава Муру, я не сошел с ума, это реальный баг в компиляторе. Причем баг был исправлен в GCC 4.9, но GCC 4.8 входил в актуальную на тот момент Ubuntu 14.04 LTS, т. е. скорее всего был установлен у большинства пользователей библиотеки. Игнорировать это было нельзя: что толку от оптимизации, если она не работает у большинства, в том числе и на продакшене, для которого она делалась. Я обновил вопрос на SO и кинул клич в твитере. На него пришел Вячеслав Егоров, один из разработчиков движка V8 и гений оптимизации, который помог докопаться до сути и нашел решение.
Чтобы понять суть проблемы, нужно углубиться в историю процессоров и их текущую архитектуру. Когда-то давно процессоры x86 не умели работать с числами с плавающей точкой, для них был придуман сопроцессор с набором команд x87. Он исполнял инструкции из того же потока, что и центральный процессор, но устанавливался на материнскую плату как отдельное устройство. Довольно быстро сопроцессор стали встраивать в центральный, и физически это стало одним устройством. Долго ли коротко ли, в третьем Пентиуме появился набор инструкций SSE (Streaming SIMD Extensions). Кстати, про SIMD инструкции будет вторая часть цикла статей. Несмотря на название, SSE содержал не только SIMD-команды для работы с числами с плавающей точкой, но и их эквиваленты для скалярных вычислений. Т. е. SSE содержал набор инструкций, дублирующий набор x87, но закодированный иначе и с немного отличающимся поведением.
Тем не менее, компиляторы не спешили генерировать SSE-код для вычислений с плавающей точкой, а продолжали использовать устаревший набор x87. Ведь наличие SSE в процессоре никто не гарантировал, в отличие от x87, который был встроен с незапамятных времен. Все изменилось с приходом 64-битного режима работы процессора. В 64-битом режиме стало обязательным наличие набора инструкций SSE2. Т. е. если вы пишете 64-битную программу для x86, вам доступны как минимум инструкции SSE2. Чем и пользуются компиляторы, генерируя для вычислений с плавающей точкой в 64-битном режиме SSE-инструкции. Напомню, что это никак не связано с векторизацией, речь об обычных скалярных вычислениях.
Именно это и происходит в нашем случае: используются разные наборы инструкций для 32-битного и 64-битного режима. Но пока это не объясняет, почему более современный SSE-код оказывается в разы медленнее устаревшего набора x87. Для объяснения этого феномена нужно разобраться, как именно процессор выполняет инструкции.
Когда-то процессоры действительно «выполняли» инструкции. Они брали инструкцию, декодировали, выполняли её полностью, клали результат куда сказано. Процессоры были довольно тупы. Современные процессоры намного умнее и сложнее, они состоят из десятков разных подсистем. Даже на одном ядре, без всякой параллельной обработки, процессор в один момент времени, в один такт, выполняет несколько инструкций. Просто выполнение происходит на разных стадиях: какая-то инструкция еще только декодируется, какая-то запрашивает что-то из кеша, какая-то передана в арифметический блок. Каждая подсистема процессора занимается своей частью. Это и называется конвейер.
На картинке разные подсистемы обозначены разным цветом. Видно, что хотя команды требует от 4 до 5 тактов на выполнение, благодаря конвейеру каждый такт выбирается одна новая команда и одна завершает свое выполнение.
Конвейер работает тем эффективнее, чем равномернее он заполнен и чем меньше подсистем простаивает. В процессоре даже есть подсистемы, которые планируют оптимальное заполнение конвейера: меняют местами инструкции, дробят одну инструкцию на несколько, объединяют несколько в одну.
Одна из вещей, которая способна заставить конвейер работать не оптимально — зависимость по данным. Это когда для выполнения какой-то инструкции необходим результат выполнения другой инструкции и многие подсистемы простаивают, ожидая когда будет выполнена предыдущая инструкция.
На рисунке видно, что инструкция 2 ждет завершения выполнения инструкции 1. В это время большинство подсистем простаивает. Выполнение всей цепочки команд замедлилось на три такта.
А теперь, вооружившись всеми этими знаниями, давайте посмотрим на псевдокод инструкции, конвертирующей целое число в число с плавающей точкой:
Instruction: cvtsi2ss xmm, r32
dst[31:0] := Convert_Int32_To_FP32(b[31:0])
dst[127:32] := a[127:32]В младшие 32 бита записывается результат конвертации. Дальше хоть и написано, что в dst переносятся биты из какого-то регистра a, по сигнатуре самой инструкции видно, что она работает только с одним xmm регистром, а значит dst и a — это один и тот же регистр, и старшие 96 бит не переносятся, а остаются нетронутыми. Вот тут и происходит зависимость по данным. Инструкция сделана таким образом, что гарантирует сохранность старших битов, а значит для её выполнения нужно дождаться результата от всех остальных операций, работающих с этим регистром. Вот только в реальности мы не пользуемся старшими битами в этих регистрах, все вычисления происходят только в младшем 32-битном float. Нам все равно, какие значения будут в старших битах, на результат работы это не повлияет. Такая зависимость по данным называется ложной.
К счастью, такую зависимость вполне можно разорвать. Компилятор должен непосредственно до инструкции конвертации cvtsi2ss вставить очистку регистра, то есть инструкцию xorps. У меня язык не поворачивается назвать этот фикс интуитивным и даже логичным. Скорее всего, это и не фикс вовсе, а хак, сделанный на уровне декодера команд, который заменяет xorps + cvtsi2ss на какую-то внутреннюю инструкцию со следующим псевдокодом:
dst[31:0] := Convert_Int32_To_FP32(b[31:0])
dst[127:32] := 0К сожалению, фикс для GCC 4.8 достаточно уродливый, он состоит из ассемблерной вставки и кода препроцессора, который проверяет уместность фикса. Приводить его не буду, ссылки на коммит, как обычно, под результатами тестов. Фикс полностью излечивает 64-битный код.
Scale 2560?1600 RGB image
320x200 bil 0.02447 s 167.42 Mpx/s
320x200 bic 0.04624 s 88.58 Mpx/s
320x200 lzs 0.07142 s 57.35 Mpx/s
2048x1280 bil 0.08656 s 47.32 Mpx/s
2048x1280 bic 0.12079 s 33.91 Mpx/s
2048x1280 lzs 0.16484 s 24.85 Mpx/s
5478x3424 bil 0.38566 s 10.62 Mpx/s
5478x3424 bic 0.52408 s 7.82 Mpx/s
5478x3424 lzs 0.65726 s 6.23 Mpx/s Результат для коммита 81fc88e. Получен на другой машине и не участвует в сравнении.
Ситуация, описанная тут, совсем не редка. Код, который переводит целые числа в числа с плавающей точкой, а потом производит достаточно небольшие вычисления, встречается почти в каждой программе. Тот же ImageMagick тоже подвержен этой проблеме: 64-битные версии, скомпилированные GCC 4.9 примерно на 40% быстрее скомпилированных более ранними версиями компилятора. Как по мне, это довольно серьезный прокол в системе команд SSE.
Комментарии (62)
0serg
22.02.2017 08:53Для уважаемого homm: DistortNeo предлагает примерно такой вариант ресайза:
import cv2 extra_smoothing_factor = 0.3 # more = less aliasing, but more blur im = cv2.imread(r"d:\pixar.jpg") downsampling_factor = im.shape[0] / 512 im = cv2.GaussianBlur(im, (0,0), sigmaX=downsampling_factor*extra_smoothing_factor) im = cv2.resize(im, (512, 214), interpolation=cv2.INTER_LINEAR) cv2.imwrite(r"d:\pixar.cv2.jpg", im)
Это вполне жизнеспособный подход и он должен быть довольно быстрым, хотя и менее качественным. Заменой гауссова размытия на более качественный low-pass filter kernel можно добиться лучшего результата, но из коробки в cv2 я сходу альтернатив не нашёл
И encyclopedist предложил отличную оптимизацию для общего ресайза свертками, подобный подход даст хорошую прибавку к скорости
homm
22.02.2017 12:23DistortNeo предлагает примерно такой вариант ресайза
Конечно же нет. DistortNeo предлагал «всё делать за один проход без использования промежуточного изображения». У вас тут добавляется еще два прохода для блюра. Но давайте все равно посмотрим на ваше предложение.
Это вполне жизнеспособный подход и он должен быть довольно быстрым
Он по определению не может быть более быстрым. Начну с того, что непонятно как вы задали extra_smoothing_factor = 0.3. При этом значении результат получается ну очень aliased. Так-то и в свертках можно сильно уменьшить радиус фильтра и получить то же самое по качеству и выигрышу в производительности. Так что давайте возьмем значение, когда результат более-менее похоже на бикубик.
extra_smoothing_factor = 0.8
Смотрите сами, GaussianBlur это тоже свертка. Радиус её равен 0.8 ? W / w ? ?2.5, где ?2.5 — коэффициент, на который opencv домножает сигму, чтобы определить границы GaussianBlur. Мне не известно точное значение, известно что для качественного блюра оно должно быть от 2 до 3.
Т.е. в вашем варианте происходит W?H сверток с радиусом W / w ? 2. Если делать как я, то понадобится w?h сверток с радиусом W / w ? 2. Ну очевидно же, что при уменьшении W?H >> w?h. Что и подтверждается практикой:
In [38]: im = Image.open('../pixar.jpg') ...: im.load() ...: %time im.resize((512, 214), Image.BICUBIC) ...: Wall time: 34.1 ms In [39]: extra_smoothing_factor = .8 ...: im = cv2.imread(r"../pixar.jpg") ...: %time im = cv2.GaussianBlur(im, (0,0), sigmaX=downsampling_factor*extra_smoothing_factor) ...: %time im = cv2.resize(im, (512, 214), interpolation=cv2.INTER_LINEAR) ...: cv2.imwrite(r"../pixar.cv2.03.png", im) ...: Wall time: 582 ms Wall time: 1.39 msDistortNeo
22.02.2017 12:51Оптимальный уровень размытия Гауссом:
extra_smoothing_factor = 0.3 * sqrt(s * s - 1), где s — коэффициент уменьшения. Уменьшаем в 2 раза — значит, надо размывать с 0.7 — 0.75. Уменьшаем в 8 раз — будет сигма 2.38 и размер фильтра 2 * 8 + 1 = 17 пикселей. Не забываем также, что при больших коэффициентах можно использовать быструю аппроксимацию фильтра Гаусса. А если он хреново реализован в OpenCV — что ж, страдать теперь? Берём C++ и реализовываем.
При большом коэффиценте
sникуда не деться — будет 2 прохода, и самый затратный — это антиалиасинг. Но отличие вот в чём: в вашем случае коэффициенты фильтра будут разные для каждого пикселя. В моём же случае это полноценная свёртка — т.е. код получается легко векторизуемым. Впрочем да, ваш код тоже можно эффективно векторизовать, если предварительно посчитанные коэффициенты размещать в памяти немного по-другому.homm
22.02.2017 13:18Оптимальный уровень размытия
Оптимальный для чего? Я же говорю, что результат получается очень aliased.
DistortNeo
22.02.2017 13:26Оптимальный для чего? Я же говорю, что результат получается очень aliased.
Оптимальный для уменьшения фотографических изображений (не абстрактной графики).
Вы, видимо, даже не попробовали написанное. Ну ладно, я помогу.
Вот исходник и результат уменьшения в 2, 4 и 8 раз, описанный мной.
homm
22.02.2017 13:49Понял, в чем проблема. В opencv координаты перевернуты и
im.shape[0]— это высота изображения. Картинка сильно горизонтальная, поэтомуdownsampling_factorиз начального комментария 0serg был сильно занижен. 0.3 действительно нормальное значение.
Тем не менее sigmaX получается такая же, производительностью я правильно намерил. Так что все равно получается в 15 раз медленнее. В OpenCV блюр не аппроксимация, но реализован очень хорошо и при сигме 3.9 работает со скоростью аппроксимации, а то и быстрее.
DistortNeo
22.02.2017 14:18Так что все равно получается в 15 раз медленнее.
… при использовании OpenCVТеперь смотрим ваши результаты:
im.resize((512, 214), Image.BICUBIC) -> Wall time: 34.1 ms
cv2.resize(im, (512, 214), interpolation=cv2.INTER_LINEAR) -> Wall time: 1.39 ms
Если делать низкочастотную фильтрацию быстрее, чем за 30 мс, то получим явное преимущество в скорости работы. И это преимущество раскрывается при больших коэффициентах уменьшения (8 и выше). Ну а то, что фильтр Гаусса в OpenCV такой медленный — это камень в огород OpenCV,
0serg
22.02.2017 13:37+1Конечно же нет. DistortNeo предлагал «всё делать за один проход без использования промежуточного изображения». У вас тут добавляется еще два прохода для блюра. Но давайте все равно посмотрим на ваше предложение.
Он свои мысли плохо излагает, но идея у него была именно эта.
Он по определению не может быть более быстрым. Начну с того, что непонятно как вы задали extra_smoothing_factor = 0.3. При этом значении результат получается ну очень aliased.
Чисто эмпирически, для Гаусса нет хорошей четкой границы частотного отсечения. По хорошему там вообще нормальный ФНЧ нужен, просто «из коробки» ничего лучше Гаусса OpenCV вроде не предлагает
Т.е. в вашем варианте происходит W?H сверток с радиусом W / w ? 2. Если делать как я, то понадобится w?h сверток с радиусом W / w ? 2. Ну очевидно же, что при уменьшении W?H >> w?h. Что и подтверждается практикой:
В этом плане да, согласен. Вообще говоря лишние пиксели которые не нужны потом в расчетах там по хорошему бы просто не нужно считать, но с точки зрения имплементации это серьезная проблема в которую подход DistortNeo так или иначе упрется.
zuborg
22.02.2017 14:14Я правильно понимаю, что усредняются значения цветов пикселей?
Если в исходном изображении будет шахматная сетка из пикселей со значениями (0,0,0) и (255,255,255), то каким будет результат уменьшения в два раза? Серый однородный фон? Серый цвет будет (127,127,127)? Или (128,128,128)? Или же (188,188,188)?0serg
22.02.2017 14:25Да, уже обсуждалось
zuborg
22.02.2017 15:10Ну, в таком случае просто некорректно называть эту процудуру «ресайзом».
Итоговое изображение не будет соответствовать исходному. Особенно весело будет, есть ресайзнутое таким образом изображение будет ресайзнуто ещё раз или несколько.
Как минимум, надо предупреждать пользователей об этом, не все ж в курсе нелинейности sRGB пространства и подобных фокусов.ToSHiC
22.02.2017 21:19У вас глаза так же работают, если будете смотреть издалека на шахматную сетку — будете видеть равномерный серый цвет. Можете проверить сами на любой газете :) Так что в этом смысле ресайз очень даже честно работает.
DistortNeo
22.02.2017 21:49Игнор гаммы обычно не приводит к негативным последствиям, зато позволяет не тратить ресурсы на дополнительные преобразования гаммы. Глядя на изображение, вы визуально не поймёте, было ли оно уменьшено с учётом гаммы, или нет. Только на специальных паттернах можно заподозрить неладное.
zuborg
22.02.2017 23:39Чем больше разница в значениях пикселов, которые сольются в один итоговый, тем больше будет отклонение значения итогового пиксела от истинного. На более-менее плавных градиентах да, разница будет не заметной. Но для изображений с большим кол-вом контрастных деталей, а также на больших степенях сжатия уже будут видны артефакты. Причем будет не только уменьшаться средняя яркость, но и искажаться цвета.
Прямое преобразование (из sRGB в линейное) — это прямой lookup по таблице в 256 элементов — не сказать что сильно большая потеря в скорости. Обратное преобразование уже сложнее, но все равно это максимум бинарный поиск по этой же таблице, без всяких умножений и логарифмов.DistortNeo
22.02.2017 23:53+1Прямое преобразование (из sRGB в линейное) — это прямой lookup по таблице в 256 элементов — не сказать что сильно большая потеря в скорости. Обратное преобразование уже сложнее, но все равно это максимум бинарный поиск по этой же таблице, без всяких умножений и логарифмов.
Есть маленький нюанс: векторизовать такой код уже не получится. А значит, можно просесть в скорости на порядок только на этом.
При использовании векторных операций, скорее всего, быстрее выйдет честно (или приближённо — небольшая ошибка роли не играет) возводить число в степень.
zuborg
22.02.2017 23:27будете видеть равномерный серый цвет
Я не против серого цвета, я против того, что его значение вычисляется некорректно.
popov654
01.03.2017 10:45Неужели было столь проблемно предупредить пользователей о необходимости обновить версию компилятора? Да и в конце концов, кому очень нужен быстрый ресайз прямо здесь и сейчас — может и 32-битную ОС для этого поставить.
Но всё равно очень круто :)homm
01.03.2017 19:33Часто компилятор в Линуксе — часть операционной системы и обновляется только с ней. Релизы Убунты LTS, например выходят раз в 2 года и актуальны года 4.
DistortNeo
Алгоритм можно сделать ещё более быстрым, если вместо двух проходов всё делать за один проход без использования промежуточного изображения. Особенно сильно разница будет заметна при параллельно-векторной реализации, когда производительность будет упираться в память.
Если на пальцах: пусть нам надо построить строку j с помощью билинейной интерполяции. Тогда для её построения нам достаточно только двух строк из исходного изображений (для бикубической — четыре). Сначала делается интерполяция по вертикали — это простая сумма двух строк из исходного изображения с определёнными коэффициентами. А сделать интерполяцию строки по горизонтали проблем нет.
Ну и сильно проще становится, когда коэффициент увеличения фиксированный, например 2 или 4.
homm
Это верно только при увеличении. Можно сделать частный случай.
DistortNeo
Это почему же?
А вообще, для уменьшения не стоит использовать использовать ни билинейную, ни бикубическую интерполяцию в чистом виде — алиасинг попрёт. Либо комбинировать с предварительной фильтрацией (размытием), либо строить адаптивное ядро.
homm
https://habrahabr.ru/post/321744/
encyclopedist
Не важно сколько именно, главное что нужно конечное количество строк. Смотрите например как я в своё время сделал для RawTherapee
ToSHiC
Вроде же в итоге надо прочитать/записать ровно столько же пикселов в памяти, за счёт чего происходит ускорение? Кэш процессора более эффективно используется?
encyclopedist
Да, идея была именно в этом. Но на самом деле надо мерить конечно же, это давно было.
DistortNeo
Интерполяция с помощью линейных фильтров выглядит следующим образом:
Т.е. строится непрерывное изображение, затем просто берутся значения в новых узлах.
При этом не имеет никакого значения — увеличивается ли картинка или уменьшается.
homm
Вы сказали, что «для её построения нам достаточно только двух строк», что не так и ваша формула с суммой от -? до +? это подтверждает.
DistortNeo
Видимо, я переборщил с матаном (из диссера выдрал). Вот более понятная формула:
Суть такова: есть дискретное изображение u, заданное в целочисленных узлах (i, j), и ядро интерполяции K(x, y) = K(x)K(y) — примеры ядер в Вашей предыдущей публикации есть. Далее делается свёртка дискретного изображения (сумма дельта-функций) с непрерывным ядром — получается непрерывное изображение u(x, y). Далее просто спокойно считаем значения там, где они нам нужны.
Так как в случае билинейной интерполяции K(x) равно 0 вне интервала (-1, 1), то в сумме будет использоваться не более двух точек с ненулевыми коэффициентами по каждой из координат. В случае бикубической — не более 4 точек.
homm
А куда из вашей формулы делись hx и hy?
DistortNeo
Я их приравнял единице. Это шаг сетки. ui, j = u(ihx, ihy).
homm
Почему вы приравняли их 1?
K(x) и K(y) от целых чисел всегда равно нулю (см графики). По вашей последней формуле x,y и i,j — целые, поэтому K(x) и K(y) не равны нулю только в 0 и ваша последняя формула полностью эквивалентна u(x, y) = ui,j
DistortNeo
(i, j) — целые, (x, y) — вещественные
А приравнял единице, чтобы считать проще было. Просто изначально задача была более общая: математическая модель построения пиксельного изображения на матрице камеры.
homm
Классно вы. Так вот, если вы их не просто выкинете, а сделаете сетку u(x, y) в масштабе конечного изображения, вы и получите метод сверток, в котором размер фильтра K зависит от масштаба и соответственно для уменьшения нужно будет больше, чем 2 строки и 2 ряда. Именно об этом методе и идет речь в статьях.
Ваше исходное утверждение «можно сделать быстрее, если делать за один проход» может оказаться верным, надо проверять. Но ваша фраза «нам достаточно только двух строк из исходного изображений» не верна.
Вот что бывает, когда «приравнял единице, чтобы считать проще было».
DistortNeo
В том-то и дело, что не зависит.
Было изображение ui, j, хотим увеличить/уменьшить в s раз и получить vi, j. Считаем: vi, j = u(i / s, j / s). Так как координаты (i / s, j / s) вещественные, интерполируем по формуле выше.
Вот чтиво — изучайте B. Generalized interpolation.
homm
Я не очень понимаю, какую цель вы преследуете. Вы меня пытаетесь убедить, что у меня считается не так, как я думаю? Или что так как у меня, считать неправильно? Выше вы написали:
Если вы считаете то, что вы описываете (приравнял единице) «чистым видом», то я её и не использую как раз. У меня как раз таки «адаптивное ядро», как вы советуете.
DistortNeo
Нет, я просто пытаюсь убедить, что для построения одной строки результирующего изображения достаточно доступа к двум строкам исходного. Доказываю через формулу. И это будет работать и для увеличения, и для уменьшения, потому что интерполяция от этого не меняется.
Разве это не логично? Если мы хотим найти значение промежуточной точки с помощью билинейной интерполяции, нам нужно только 4 окрестные точки (квадрат 2х2), и ничего больше. Если с помощью бикубической, то 4х4.
homm
Вы можете строить результирующее изображение, используя хоть одну строку исходного изображения, хоть ноль. Пожалуйста, кто вам запретит. Но результат будет вот такой. Для ваших целей это может быть приемлемо, для моих нет.
Я не знаю, зачем вы это хотите. Дело ваше, опять же. В первой статье я описал, чего хочу я и почему.
3dcryx
Я не могу понять, что имеется ввиду что нужно 4 или 16 точек для интерполяции. Т.е если у нас картинка 1024x1024 и мы хотим уменьшить ее до 64x64, разве не нужно разбить исходное изображение на квадратики 16x16 и вычилить коэффициенты интерполяционного многочлена по всем точкам (МНК там какимнить), возможно даже учтя некоторые точки соседних квадратов?
homm
Когда исходное изображение разбивается на квадараты и берется среднее, это называется supersampling.
DistortNeo
Первый раз слышу такой термин. В научной среде он не используется совсем, это называется downscaling using box filter.
DistortNeo
Если вы под адаптивным ядром при уменьшении понимаете растяжение ядра сразу на много пикселей, т.е. при уменьшении в 10 раз шапочка линейного фильтра накрывает 21 пиксель, то да, такой способ вполне себе имеет право на жизнь.
Но такое уменьшение будет довольно далеко от реальных моделей, где в качестве функций рассеяния точки берётся функция Гаусса или прямоугольник.
Поэтому я бы советовал при уменьшении делать сначала фильтрацию Гаусса, затем применять неадаптивный фильтр.
encyclopedist
Функцию Гаусса или прямогольник берут только в учебниках, чтобы считать было проще. Как показала история и исследования, для изображений подходит функция типа Ланцоша.
0serg
Это не «исследования» показали, это чистая математика.
Все ядра для интерполяции — это вариации на тему фильтра низких частот. Идеальный ФНЧ — это sinc(x), но он требует бесконечно широкого окна. В реальных фильтрах же для удобства вычислений используют разные его приближения, Ланцош из которых является одним из наиболее точных
encyclopedist
Я собственно это и имел ввиду.
DistortNeo
А ещё он требует выполнения условий теоремы Котельникова. Но на практике изображений с ограниченной частотой, за редким исключением, просто не существует. Поэтому если увеличивать изображение с помощью sinc, то полезет рингинг (эффект Гиббса).
На деле, какой фильтр не выбери, на выходе будет комбинация между алиасингом, размытием и рингингом. Если бороться с одним эффектом, вылезеть другой, оптимального фильтра нет. Ланцош — это просто хороший компромисс между этими эффектами.
DistortNeo
Возьмите любую современную научную публикацию по суперразрешению (желательно многокадровому). В модели понижения разрешения используют в подавляющем большинстве случаев Гаусса, совсем редко — прямоугольный фильтр.
Есть научно обоснованная математическая модель, а есть простые фильтры, которые применяют, когда математическую модель по какой-то причине применить нельзя. Противоречий здесь нет.
0serg
Вы мешаете в кучу две совершенно разные вещи. Пожалуйста, не надо этого делать. Не порождайте у меня желания разгромить Ваш диссер (а то уже честно говоря захотелось)
DistortNeo
Как раз наоборот: я их не мешаю, а разделяю. Мешают их те, кто пытается две вещи делать одновременно: интерполяцию одновременно с подавлением алиасинга.
0serg
В этом объединении нет ничего странного или плохого, это просто удобно и вполне эффективно. Вы же зачем-то мешаете PSF и SR с обычным downsampling-ом. Подробнее написал отдельно ниже.
homm
Спасибо большое. За сим предлагаю закончить.
0serg
DistortNeo, у Вас концептуально верная формула, но Вы совершенно неверно ее применяете
При уменьшении изображения в N раз частота Найквиста уменьшается тоже в N раз.
Дабы избежать алиасинга Вам нужно в качестве ядра взять фильтр низких частот который отрежет частоты выше частоты Найквиста. Очевидно что для разных N и разных частот Найквиста этот фильтр будет разным. И простейший способ получить такой фильтр — это взять ФНЧ для частоты Найквиста в исходном изображении, а затем растянуть его в пространственном домене в N раз (это сожмет его в частотном в N раз же, что нам для соответствия новой частоте Найквиста и требуется). Поэтому при сжатии в 2 раза ядро «правильной билинейной» интерполяции будет уже задано на [-2,2] а не [-1,1]. При работе же с билинейкой [-1,1] в таком сценарии Вы получите замечательный алиасинг т.к. частоты между «старой» и «новой» частотами Найквиста не будут подавлены. Вот при апсемплинге такой проблемы нет и там ядро всегда будет определяться частотой Найквиста в исходном сигнале, там ядро действительно получается константным.
Функции рассеяния, кстати, тут вообще не при чем. Вы ее применили один раз, при оцифровке исходного изображения и все, дальше можете забыть: она навсегда стала частью оцифрованного сигнала. Исключением является только тот случай когда Вы последующей обработкой пытаетесь ее компенсировать.
homm
Спасибо вам большое! Как вы видите, меня больше интересует прикладная сторона и моих знаний недостаточно чтобы объяснить почему так правильно, а так неправильно. Но теперь у меня будет место, куда я смогу ссылаться.
0serg
Давно хочу написать статью, но увы… А там много всего интересного есть. К примеру когда мы апскейлим изображение в методиках super-resolution, то оказывается что при наличии достаточного количества немного отличающихся картинок можно легко обратить алиасинг. Поэтому там «самый правильный» способ апскейлинга — это тот который алиасинг не трогает, пространство между пикселями в растянутом изображении просто заполняется нолями, эдакая сеточка из пикселей получается :).
DistortNeo
У вас просто однобокий взгляд на проблему. Судя по всему, вы занимаетесь чисто прикладными задачами, скорее всего — компьютерной графикой, визуализацией, либо не очень сложной обработкой.
Эффект Даннинга — Крюгера имеет место быть, причём в отношении нас обоих.
Оба варианта правильные, просто применяются в совершенно различных задачах.
DistortNeo
Да, я это прекрасно знаю. Но я предпочитаю разделять эти две операции. При уменьшении я сначала делаю фильтрацию, затем decimation по моей формуле. Так получается быстрее, чем при объединении этих операций, и вот почему: при интерполяции коэффициенты у каждого пикселя получаются свои (особенно если уменьшение идёт в нецелое число раз), поэтому пересчитывать эти коэффициенты желательно по-минимуму.
0serg
Такой подход имеет право на жизнь, но он менее точен и больше портит изображение т.к. ФНЧ эффективно применяется дважды. И это не единственный способ, Вам описывают совершенно корректную и более точную методику. Которая при правильной реализации к тому же может быть даже более быстрой чем Ваша ибо у Вас там два прохода по памяти получится, а в таком подходе нужен лишь один.
DistortNeo
Зависит от коэффициента уменьшения:
Ну а количество проходов по памяти, естественно, должно быть по-минимуму. В идеале, всё должно делаться за один проход, просто реализация будет более трудоёмкой.
0serg
Для целых коэффициентов это так, но там и этап фильтрации получится абсолютно идентичным ресайзу свертками. Для остальных — возможно, но это надо проверять.
DistortNeo
Ну и отдельно про функцию рассеяния точки: иногда она используется.
Задачу уменьшения изображения можно сформулировать так: по изображению U, полученному с камеры с определённым разрешением и определённой функцией рассеяния точки (PSF), построить изображение V, как если бы оно было получено с такой же камеры, но в S раз меньшим разрешением и пропорционально растянутой PSF.
Соответственно, здесь низкочастотный фильтр, применяемый для уменьшения изображения, напрямую связан с PSF камеры. А на PSF камеры влияет дефокус/апертура и собственно сенсор. Если мы хотим получить изображение, максимальное приближенное к реальному, PSF придётся использовать. А если просто уменьшить картинку — да что угодно.
0serg
Задача уменьшения изображения — это просто downsampling исходного сигнала, полностью самостоятельная задача не имеющая никакого отношения к PSF(т.к. downsamplится сигнал уже свернутый с PSF).
SR-методы опираются на универсальную модель формирования цифрового изображения «свернули непрерывное изображение с PSF, сделали выборку по дискретной сетке (дискретизацию)». И все они пытаются так или иначе восстановить непрерывное изображение, хотя и формируют лишь его представление с другой PSF и/или частотой дискретизации. Оба этапа (свертка и дискретизация) вносят определенные искажения. Свертку можно убрать deconvolution и более сложными image-aware аналогами, для этого (если PSF >> размера пикселя) в принципе достаточно даже одного изображения. Дискретизация вносит алиасинг (мешающий деконволюции), поэтому выше разрешения исходного изображения (не прибегая к уловкам с «угадыванием» изображения) так не подняться, но если у нас есть набор картинок, то алиасинг можно подавить и получить новую, более мелкую дискретизацию (картинку с субпиксельным разрешением)
Так вот: гаусс и вообще PSF и дискретизация — разные вещи дающие разные эффекты с которыми борются по разному. А такой задачи ресайза как Вы описываете и в помине никто не ставит.
DistortNeo
Если коэффициент уменьшения небольшой (1-1.5) или картинка просто аффинно деформируется, то пиксельный PSF будет играть весьма большую роль, т.к. его частотные характеристики будут сильно искажать характеристики интерполяционного фильтра.
Такая задача существует, просто не используется в чисто прикладных задачах. Видимо, для расширения кругозора местных обитателей и разрешения холивара, которого на самом деле нет, я в ближайшее время напишу статью о том, какие вообще существуют подходы для работы с изображениями и где они используются.
Наилучшее качество при уменьшении достигается, кстати, не Ланцошом, а специальными методами обработки, которые оперируют с изображением не как с абстрактным сигналом, а как с изображением, на котором есть определённые информация, специфичная именно для изображений, о качестве которой стоит задуматься, например:
Content aware image downsampling.
0serg
Пиксельный PSF — это просто часть общего PSF, его нет нужды рассматривать каким-то специальным образом.
По остальному — Вы в какую-то специфичную задачу закапываетесь. Там есть разные подходы конечно, но я могу точно сказать что на основании специфики этой Вашей задачи Вы много фигни разной выше в комментариях написали о другой задаче (ресемплинге изображений).
Content-aware методы это отдельный разговор, да. Я в комментариях это вроде везде пытался оговорить отдельно.