История об оптимизация функции alpha_composite в Pillow 2.0

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

С тех пор я уже смог написать для своих нужд более оптимальную реализацию, чем приведена в конце той статьи. И оказалось, что эта реализация быстрее alpha_composite() из новой версии Pillow, написанной на Си. Конечно, мне это польстило, но я все-таки решил попытаться улучшить реализацию из Pillow.

Для начала нужен тестовый стенд.

from PIL import Image from timeit import repeat from image import paste_composite  im1 = Image.open('in1.png') im1.load()  im2 = Image.open('in2.png') im2.load()  print repeat(lambda: paste_composite(im1.copy(), im2), number=100) print repeat(lambda: Image.alpha_composite(im1, im2), number=100)  out1 = im1.copy() paste_composite(out1, im2) out1.save('out1.png')  out2 = Image.alpha_composite(im1, im2) out2.save('out2.png') 

Скрипт берет лежащие в текущей папке файлы in1.png и in2.png и смешивает их сто раз сначала функцией paste_composite(), потом alpha_composite(). Я взял файлы из предыдущей статьи (раз, два). Время работы paste_composite() у меня составило в среднем 235мс. Время работы оригинальной alpha_composite() составило 400мс.

Почему alpha_composite() такая медленная, долго гадать не пришлось. Глядя на реализацию, видно, что функция производит все вычисления на числах с плавающей точкой. При работе с графикой числа с плавающей точкой позволяют достаточно просто описать алгоритм, но производительность у них всегда меньше, чем у целых чисел. Первое, что я сделал, перевел все на работу с целыми числами:

typedef struct {     UINT8 r;     UINT8 g;     UINT8 b;     UINT8 a; } rgba8;  Imaging ImagingAlphaComposite(Imaging imDst, Imaging imSrc) {     Imaging imOut = ImagingNew(imDst->mode, imDst->xsize, imDst->ysize);      for (int y = 0; y < imDst->ysize; y++) {         rgba8* dst = (rgba8*) imDst->image[y];         rgba8* src = (rgba8*) imSrc->image[y];         rgba8* out = (rgba8*) imOut->image[y];          for (int x = 0; x < imDst->xsize; x ++) {             if (src->a == 0) {                 *out = *dst;             } else {                 UINT16 blend = dst->a * (255 - src->a);                 UINT8 outa = src->a + blend / 255;                 UINT8 coef1 = src->a * 255 / outa;                 UINT8 coef2 = blend / outa;                  out->r = (src->r * coef1 + dst->r * coef2) / 255;                 out->g = (src->g * coef1 + dst->g * coef2) / 255;                 out->b = (src->b * coef1 + dst->b * coef2) / 255;                 out->a = outa;             }             dst++; src++; out++;         }     }     return imOut; } 

И сразу же получил прирост в 5,5 раз. Тест отрабатывал за 75 мс.

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

#!/usr/bin/env python import sys from PIL import Image, ImageMath  def chanel_diff(c1, c2):     c = ImageMath.eval('127 + c1 - c2', c1=c1, c2=c2).convert('L')     return c.point(lambda c: c if 126 <= c <= 128 else 127 + (c - 127) * 10)  im1 = Image.open(sys.argv[1]) im2 = Image.open(sys.argv[2])  diff = map(chanel_diff, im1.split(), im2.split())  Image.merge('RGB', diff[:-1]).save('diff.png', optimie=True) diff[-1].convert('RGB').save('diff.alpha.png', optimie=True) 

На вход скрипт принимает 2 имени файла с картинками. Между каждой парой каналов этих картинок находятся отличия: сначала просто ищется разница и складывается с серым фоном (значение 127), чтобы можно было увидеть отклонения в обе стороны. Затем все отклонения, сильнее чем в одно значение, усиливаются в 10 раз, чтобы они были заметны визуально. Результат сравнения записывается в 2 файла: rgb каналы в diff.png, альфа-канал в виде оттенков серого в diff.alpha.png.

Оказалось, что разница между результатом alpha_composite() и фотошопом достаточно существенная:

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

UINT16 blend = dst->a * (255 - src->a); UINT8 outa = src->a + (blend + 127) / 255; UINT8 coef1 = src->a * 255 / outa; UINT8 coef2 = blend / outa;  out->r = (src->r * coef1 + dst->r * coef2 + 127) / 255; out->g = (src->g * coef1 + dst->g * coef2 + 127) / 255; out->b = (src->b * coef1 + dst->b * coef2 + 127) / 255; out->a = outa; 

Время работы увеличилось до 94 мс. Но и отличий от референса стало намного меньше. Впрочем, они не пропали совсем.

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

Очевидно, что полученные неточности — результат недостаточной точности коэффициентов при умножении. Можно хранить дополнительные значащие биты для каждой переменной, что поднимает точность при делении. Переменная blend вычисляется без деления, с ней ничего делать не нужно. Переменная outa получается делением blend на 255. При расчете сдвигаем значения blend и src->a на 4 бита вверх. В результате в outa будет 12 значащих бит. В обоих коэффициентах происходит деление на outa, которая теперь больше на 4 бита. Плюс сами коэффициенты хотелось бы получить на 4 бита больше. В результате делимое сдвигаем уже на 8 бит.

UINT16 blend = dst->a * (255 - src->a);  // 16 bit max UINT16 outa = (src->a << 4) + ((blend << 4) + 127) / 255;  // 12 UINT16 coef1 = ((src->a * 255) << 8) / outa;  // 12 UINT16 coef2 = (blend << 8) / outa;  // 12  out->r = ((src->r * coef1 + dst->r * coef2 + 0x7ff) / 255) >> 4; out->g = ((src->g * coef1 + dst->g * coef2 + 0x7ff) / 255) >> 4; out->b = ((src->b * coef1 + dst->b * coef2 + 0x7ff) / 255) >> 4; out->a = (outa + 0x7) >> 4; 

Само собой, на эти же 4 бита нужно сдвинуть результат всех вычислений, прежде чем помещать его в пиксели изображения. Такой вариант работает еще чуть медленнее: 108 мс, но зато уже практически не имеет пикселей, отличающихся более чем на одно значение от эталона. Результат по качеству достигнут, можно опять подумать о производительности.

Очевидно, что такой код не самый оптимальный с точки зрения процессора. Деление — одна из самых сложных задач, а здесь на каждый пиксель выполняется целых 6 делений. Оказалось, что избавиться от большинства из них достаточно просто. В той же PIL в файле Paste.c описан способ быстрого деления с округлением:
a + 127 / 255 ≈ ((a + 128) + (a + 128) >> 8) >> 8

В результате одно деление заменяется на 2 сдвига и одно сложение, что быстрее на большинстве процессоров. Сгруппировав один из сдвигов с уже имеющимся сдвигом на 4, я получил такой код:

UINT16 blend = dst->a * (255 - src->a); UINT16 outa = (src->a << 4) + (((blend << 4) + (blend >> 4) + 0x80) >> 8); UINT16 coef1 = (((src->a << 8) - src->a) << 8) / outa;  // 12 UINT16 coef2 = (blend << 8) / outa;  // 12  UINT32 tmpr = src->r * coef1 + dst->r * coef2 + 0x800; out->r = ((tmpr >> 8) + tmpr) >> 12; UINT32 tmpg = src->g * coef1 + dst->g * coef2 + 0x800; out->g = ((tmpg >> 8) + tmpg) >> 12; UINT32 tmpb = src->b * coef1 + dst->b * coef2 + 0x800; out->b = ((tmpb >> 8) + tmpb) >> 12; out->a = (outa + 0x7) >> 4; 

Он выполняется уже за 65 мс и не вносит изменений в результат по сравнению с предыдущим вариантом.

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