I accelerated blurhash generation by 3̶6̶ 8̶7̶ 128 times
You can't teach an old dog new tricks, so I took up the old tricks. Blurhash is a compact way to represent a blurred image preview as an ASCII string. Developed by the Finnish company Wolt (similar to Delivery Club). I've been wanting to implement something like this for a long time to your APIso that any client can load content on their site more smoothly and gracefully. But no matter how much I looked at it, I was always haunted by the speed of the work; it was written too slowly and “head-on”. But now the time has come to finally figure out why it is working so slowly.
Main repository − https://github.com/woltapp/blurhashit contains a reference implementation in as many as four languages (C, Swift, Kotlin, TS). But I started my journey with libraries for Python.
0. Interface C <=> Python
First, let's figure out how the library generally works. The sources showthat interaction with the library is implemented through the cffi interface, but implemented in a slightly atypical way: instead of a real interface to any system library, files are literally copied verbatim into the project decode.c And encode.c and there are no other dependencies. This approach has the right to life for small libraries with a simple interface, which in this case is blurhash.
C function create_hash_from_pixels accepts:
Number of hash components (in other words, its resolution)
x_componentsAndy_componentsSize of the image sent to it
widthAndheightalong with the image itselfrgb.Increment in bytes to move to the next line
bytes_per_row– in case the line takes up more space than is needed to accommodate it. Note to the owner: this option, if necessary, allows you to get a “crop” of the image for free. To do this, you just need to shift the starting pointer to the first pixel of the button image and specifybytes_per_rowcorresponding to not cropped.String for the result of the work
destination.
From Python's side, the call goes like this:
def encode(image, x_components, y_components):
if not isinstance(image, Image.Image):
image = Image.open(image)
if image.mode != 'RGB':
image = image.convert('RGB')
red_band = image.getdata(band=0)
green_band = image.getdata(band=1)
blue_band = image.getdata(band=2)
rgb_data = list(chain.from_iterable(zip(red_band, green_band, blue_band)))
width, height = image.size
image.close()
rgb = _ffi.new('uint8_t[]', rgb_data)
bytes_per_row = _ffi.cast('size_t', width * 3)
destination = _ffi.new("char[]", 2 + 4 + (9 * 9 - 1) * 2 + 1)
result = _lib.create_hash_from_pixels(
x_components, y_components,
width, height, rgb, bytes_per_row,
destination,
)Here you can already grab your hair: all the pixels of the image are iterated, then a list is created from them, which is already placed in a blue array rgb. Let's see how it works:
from blurhash import encode, Image im = Image.open('/Code/imgs/bologna2k.jpg').resize((360, 240), Image.BICUBIC) assert encode(im, 6, 4) == 'W7EUfLs: M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
%timeit encode(im, 6, 4)Я вставил assert, чтобы убедиться, что ничего не сломаю в процессе оптимизаций. Для тестов Уже тут нас ждет неприятный сюрприз:
ValueError: Operation on closed imageЧто это значит? А это тот самый image.close() на 11-й строке. Ну ок, пока что обойдем его, создавая копию изображения перед каждым вызовом функции:
...: assert encode(im.copy(), 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im.copy(), 6, 4)
126 ms ± 995 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)126 мс на крошечное изображение разрешением 360x240 пикселей — это конечно мощно. Хорошо, что на самом деле нам не нужно получать никакие отдельные каналы изображения, чтобы снова их соединить. Чтобы получить представление картинки в виде куска памяти в Pillow есть метод .tobytes():
- red_band = image.getdata(band=0)
- green_band = image.getdata(band=1)
- blue_band = image.getdata(band=2)
- rgb_data = list(chain.from_iterable(zip(red_band, green_band, blue_band)))
+ rgb_data = image.tobytes() ...: assert encode(im.copy(), 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im.copy(), 6, 4)
109 ms ± 487 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)Вот уже 15% прироста. Заодно избавляемся от лишних импортов из itertools и six (что?).
Но что же там с image.close()? На сама деле он не бесполезен, т.к. функция может принимать как готовое изображение, так и открывать его из файла. И кстати тут присутствует баг, что если открытое изображение было не в режиме RGB, то закроется не открытое изображение, а конвертированное (строка номер 5). Но вот закрывать готовое изображение явно не надо. Исправить это можно, задав контекст:
from contextlib import nullcontext
def encode(image, x_components, y_components):
if isinstance(image, Image.Image):
image_context = nullcontext()
else:
image = Image.open(image)
image_context = image
with image_context:
if image.mode != 'RGB':
image = image.convert('RGB')
rgb_data = image.tobytes()
width, height = image.size
rgb = _ffi.new('uint8_t[]', rgb_data) ...The rest of the code remains the same. But now you can remove the explicit call .copy() before the call encode()
In [2]: from blurhash import encode, Image ...: im = Image.open('/Code/imgs/bologna2k.jpg').resize((360, 240), Image.BICUBIC) ...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
108 ms ± 30 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)Большого ускорения это не дало, зато посмотрите насколько стали стабильнее результаты без лишнего копирования в памяти.
1. Кешируем всё что можно
Давайте наконец-то взглянем на encode.c. Тут всё работа (и 99.9% времени выполнения) происходит в функции multiplyBasisFunction. Она вызывается для каждого пикселя результирующего изображения (для каждого xComponents и yComponents.
static struct RGB multiplyBasisFunction(int xComponent, int yComponent, int width, int height, uint8_t *rgb, size_t bytesPerRow) {
struct RGB result = { 0, 0, 0 };
float normalisation = (xComponent == 0 && yComponent == 0) ? 1 : 2;
for(int y = 0; y < height; y++) {
for(int x = 0; x < width; x++) {
float basis = cosf(M_PI * xComponent * x / width) * cosf(M_PI * yComponent * y / height);
result.r += basis * sRGBToLinear(rgb[3 * x + 0 + y * bytesPerRow]); result.g += basis * sRGBToLinear(rgb[3 * x + 1 + y * bytesPerRow]); result.b += basis * sRGBToLinear(rgb[3 * x + 2 + y * bytesPerRow]); } } float scale = normalization / (width * height); result.r *= scale; result.g *= scale; result.b *= scale; return result; }What's going on here? Some coefficient is calculated basis and is multiplied by each pixel. Bah, yes it is cosine transform! That is, in essence, each blurhash is a miniature JPEG without headers and with a fixed quantization table. Interesting approach. Another interesting thing is that calculations are performed in linear color space (function sRGBToLinear). And this is the first candidate for serious optimization here. The fact is that the pixels in an image can only contain 256 different values. But the pixels themselves, even in our example, are 360 * 240 * 3 = 259,200, a thousand times more. That is, it would be much more profitable to compile a table of all possible responses of the function in advance sRGBToLinear and then contact her later. But you can go even further: the fact is that these values will not change, no matter how many times you call sRGBToLinear. This means that you can calculate this table on the first call encode and cache it forever. In memory, such a table will take only 256 * 4 = 1024 bytes.
float *sRGBToLinear_cache = NULL;
static void init_sRGBToLinear_cache() {
if (sRGBToLinear_cache != NULL) {
return;
}
sRGBToLinear_cache = (float *)malloc(sizeof(float) * 256);
for (int x = 0; x < 256; x++) {
sRGBToLinear_cache[x] = sRGBToLinear(x);
}
}
static struct RGB multiplyBasisFunction(int xComponent, int yComponent, int width, int height, uint8_t *rgb, size_t bytesPerRow) {
struct RGB result = { 0, 0, 0 };
init_sRGBToLinear_cache();
for(int y = 0; y < height; y++) {
for(int x = 0; x < width; x++) {
float basis = cosf(M_PI * xComponent * x / width) * cosf(M_PI * yComponent * y / height);
result.r += basis * sRGBToLinear_cache[rgb[3 * x + 0 + y * bytesPerRow]];
result.g += basis * sRGBToLinear_cache[rgb[3 * x + 1 + y * bytesPerRow]];
result.b += basis * sRGBToLinear_cache[rgb[3 * x + 2 + y * bytesPerRow]];
}
}
...
}We check:
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
14.2 ms ± 6.37 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)Воу-воу-воу! Мы только начали, а уже ускорили весь код в 8.87 раз! Впрочем, это только начало.
Взгляните, как вычисляется basis. Дело в том, что cosf — довольно тяжелая функция, она выполняется ни один такт (а скорее всего, даже не десять). Но значения косинусов зависят только от счетчиков, которые раз за разом проходят одни и те же значения. К счастью, компилятор достаточно умен, чтобы не считать cos(y) для каждого пикселя. Но вот с cos(x) он сам не в состоянии справится. Давайте ему поможем:
static struct RGB multiplyBasisFunction(int xComponent, int yComponent, int width, int height, uint8_t *rgb, size_t bytesPerRow) {
struct RGB result = { 0, 0, 0 };
float *cosx = (float *)malloc(sizeof(float) * width);
for(int x = 0; x < width; x++) {
cosx[x] = cosf(M_PI * xComponent * x / width); } for(int y = 0; y < height; y++) { float cozy = cosf(M_PI * yComponent * y / height); uint8_t *src = rgb + y * bytesPerRow; for(int x = 0; x < width; x++) { float basis = cozy * cosx[x]; result.r += basis * sRGBToLinear_cache[src[3 * x + 0]]; result.g += basis * sRGBToLinear_cache[src[3 * x + 1]]; result.b += basis * sRGBToLinear_cache[src[3 * x + 2]]; } } free(cosx); ...In [1]: from blurhash import encode, Image ...: im = Image.open('/Code/imgs/bologna2k.jpg').resize((360, 240), Image.BICUBIC) ...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
3.45 ms ± 3.25 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Обещанное ускорение в 36 раз практически на пустом месте. В принципе тут бы конечно стоило остановиться, такой скорости уже вполне хватит для внедрения этой библиотеки. Но дальше во мне проснулся спортивный интерес.
2. Делаем всё за один проход
Для чего нужная функция multiplyBasisFunction? Она считает ровно один коэффициент косинусного преобразования. Делает она это считая basis для каждого пикселя исходного изображения. Но вот в чем штука: для одинаковых xComponent и yComponent наборы посчитанных косинусов будет одинаковыми. А мы их все равно считаем каждый раз. Следовательно, можно считать их заранее и передавать в функцию уже готовые массив.
Сама функция при этом станет только проще:
static struct RGB multiplyBasisFunction(
int xComponent, int yComponent, int width, int height, uint8_t *rgb, size_t bytesPerRow,
float *cosx, float *cosy
) {
struct RGB result = { 0, 0, 0 };
for(int y = 0; y < height; y++) {
uint8_t *src = rgb + y * bytesPerRow;
for(int x = 0; x < width; x++) {
float basis = cosy[y] *cosx[x]; result.r += basis * sRGBToLinear_cache[src[3 * x + 0]]; result.g += basis * sRGBToLinear_cache[src[3 * x + 1]]; result.b += basis * sRGBToLinear_cache[src[3 * x + 2]]; } } ... }And here is the one calling her blurHashForPixels will change noticeably:
const char *blurHashForPixels(int xComponents, int yComponents, int width, int height, uint8_t *rgb, size_t bytesPerRow, char *destination) {
if(xComponents < 1 || xComponents > 9) return NULL;
if(yComponents < 1 || yComponents > 9) return NULL;
float factors[9 * 9][3];
init_sRGBToLinear_cache();
float *cosx = (float *)malloc(sizeof(float) * width * xComponents);
if (! cosx) return NULL;
float *cosy = (float *)malloc(sizeof(float) * height);
if (! cosy) {
free(cosx);
return NULL;
}
for(int x = 0; x < xComponents; x++) {
for(int i = 0; i < width; i++) {
cosx[x * width + i] = cosf(M_PI * x * i / width);
}
}
for(int y = 0; y < yComponents; y++) {
for(int i = 0; i < height; i++) {
cosy[i] = cosf(M_PI * y * i / height);
}
for(int x = 0; x < xComponents; x++) {
struct RGB factor = multiplyBasisFunction(x, y, width, height, rgb, bytesPerRow,
cosx + x * width, cosy);
factors[y * xComponents + x][0] = factor.r;
factors[y * xComponents + x][1] = factor.g;
factors[y * xComponents + x][2] = factor.b;
}
}
free(cosx);
free(cosy);Here we prepare all the coefficients in advance.
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
3.4 ms ± 3.65 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)К сожалению прироста от этого практически никакого, я и так уже очень серьезно сократил количество расчетов косинусов. Значит нужно придумать что-то ещё.
Немного подумав, можно заметить, что не только косинусы между разными вызовами multiplyBasisFunction не меняются, но и пиксели изображения. Сейчас мы проходим по изображению и достаем значения из sRGBToLinear_cache 6 * 4 = 24 раза для каждого канала. Какой в этом смысл? Правильно — никакого. Значит пришла пора для самого крупного изменения в коде и логике. Нужно сделать все вычисления за один проход.
Здесь возникает дилемма: если сделать в лоб и поместить обход xComponent и yComponent внутрь каждого пикселя, то получится 4 уровня вложенности циклов, причем внутренние циклы будут иметь очень небольшую дистанцию. Я не проверял, как это будет работать, но из опыта знаю, что код может сильно замедлиться из-за постоянных ошибок предсказателя ветвлений процессора.
Я решил попробовать другой подход: схлопнуть два внутренних цикла в один и не различать отдельно xComponent и yComponent, а сделать один плоский массив. Но тогда возникает вопрос: а как нам на каждом шаге узнать значения этих xComponent и yComponent, чтобы получить правильные индексы для косинусов? Целочисленно делить и брать модуль от деления счетчика — тоже не самая быстрая операция. Я решил, что намного проще будет просто продублировать нужные коэффициенты столько раз, сколько у нас xComponents и yComponents. То есть каждый массив cosX и cosY будет иметь размерность не width * xComponents, а width * factorsCount, где factorsCount = xComponents * yComponents. Это не увеличит количество расчетов, только незначительно увеличит расход памяти. В общем итоговый вид функции multiplyBasisFunction стал таким:
static void multiplyBasisFunction(
struct RGB *factors, int factorsCount, int width, int height, uint8_t *rgb, size_t bytesPerRow,
float *cosX, float *cosY
) {
for(int y = 0; y < height; y++) {
uint8_t *src = rgb + y * bytesPerRow;
float *cosYLocal = cosY + y * factorsCount;
for(int x = 0; x < width; x++) {
float pixel[3]; float *cosXLocal = cosX + x * factorsCount; pixel[0] = sRGBToLinear_cache[src[3 * x + 0]]; pixel[1] = sRGBToLinear_cache[src[3 * x + 1]]; pixel[2] = sRGBToLinear_cache[src[3 * x + 2]]; for (int i = 0; i < factorsCount; i++) { float basis = cosYLocal[i] *cosXLocal[i]; factors[i].r += basis * pixel[0]; factors[i].g += basis * pixel[1]; factors[i].b += basis * pixel[2]; } } } for (int i = 0; i < factorsCount; i++) { float normalization = (i == 0) ? 1 : 2; float scale = normalization / (width * height); factors[i].r *= scale; factors[i].g *= scale; factors[i].b *= scale; } } ...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
1.44 ms ± 950 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
New record 87.5 times!
3. SIMD (SSE + NEON)
Finally my favorite topic. Now it’s going to get hot, I thought.
I have always written SIMD only for x86, but here I had a goal to finally try to make full support for ARM (I ran all previous tests under the rosette on the M1 Pro).
Minimal changes were required for SSE. And by the way, this version is really SSE (the very first version), that is, it will run on the oldest processor, apparently starting from the Pentium III.
static void multiplyBasisFunction(
float factors[][4], int factorsCount, int width, int height, uint8_t *rgb, size_t bytesPerRow,
float *cosX, float *cosY
) {
for(int y = 0; y < height; y++) {
uint8_t *src = rgb + y * bytesPerRow;
float *cosYLocal = cosY + y * factorsCount;
int x = 0;
for(; x < width; x++) {
float *cosXLocal = cosX + x * factorsCount;
#ifdef __SSE__
__m128 pixel = _mm_set_ps(0, sRGBToLinear_cache[src[3 * (x+0) + 2]],
sRGBToLinear_cache[src[3 * (x+0) + 1]], sRGBToLinear_cache[src[3 * (x+0) + 0]]);
for (int i = 0; i < factorsCount; i++) {
__m128 basis = _mm_set1_ps(cosYLocal[i] * cosXLocal[i + 0 * factorsCount]);
__m128 factor = _mm_loadu_ps(factors[i]);
factor = _mm_add_ps(factor, _mm_mul_ps(basis, pixel));
_mm_storeu_ps(factors[i], factor);
}
#else
float pixel[4];
pixel[0] = sRGBToLinear_cache[src[3 * x + 0]];
pixel[1] = sRGBToLinear_cache[src[3 * x + 1]];
pixel[2] = sRGBToLinear_cache[src[3 * x + 2]];
for (int i = 0; i < factorsCount; i++) {
float basis = cosYLocal[i] * cosXLocal[i];
factors[i][0] += basis * pixel[0];
factors[i][1] += basis * pixel[1];
factors[i][2] += basis * pixel[2];
}
#endif
}
}
...
}I received my well-deserved increase:
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
973 μs ± 1.13 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)Это кстати обещанные 129 раз! Наивный же код, без интринсиков NEON и без rosetta, выполнялся на M1 pro c ещё большей скоростью (эти измерения не идут в общий зачет):
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
808 μs ± 454 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)Эквивалент кода под NEON получился вполне прямолинейным:
#elif defined(__aarch64__)
float32x4_t pixel = {sRGBToLinear_cache[src[3 * (x+0) + 0]], sRGBToLinear_cache[src[3 * (x+0) + 1]], sRGBToLinear_cache[src[3 * (x+0) + 2]]}; for (int i = 0; i < factorsCount; i++) { float32x4_t basis = vdupq_n_f32(cosYLocal[i] * cosXLocal[i + 0 * factorsCount]); float32x4_t factor = vld1q_f32(factors[i]); factor = vmlaq_f32(factor, basis, pixel); vst1q_f32(factors[i]factor); } #elseJoyful, I launched the benchmark:
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
962 μs ± 24.5 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
Опачки. Опачки. Код замедлился в 0.84 раза :-/ Даже несмотря на операцию сдвоенного умножения и сложения vmlaq_f32.
На самом деле к этому моменту я уже давно перешел в основной репозиторий blurhash и тестировал код на чистом Си, это в статье я продолжаю пользоваться blurhash-python. Поэтому мне не составила труда сдампить объектный файл encode.o ванильной версии и посмотреть, какие оптимизации применял компилятор, оказавшийся умнее меня. Оказалось, что компилятор не просто разворачивает внутренний цикл по factorsCount, он ещё и читал коэффициенты инструкциями ld4, которые загружали гораздо больше данных. (Кстати, тут я снова нарвался на баг в Докере с замедлением этих инструкций, если включена Rosetta).
Чтож, попытка оказалась провальной. Можно было бы конечно разворачивать циклы в своем коде тоже, подсматривая у компилятора. Но вместо того, чтобы пытаться соперничать с ним, я решил пойти другим путём, я решил посмотреть что будет, если попытаться ему помочь.
4. Развязываем компилятору руки
Первое, что можно было сделать — выровнять используемые структуры по 4 элемента, чтобы компилятор мог свободно распоряжаться ими как одним вектором. Для этого я отказался от struct RGB и перешел на просто массив из 4 элементов.
Дальше я взглянул на код SEE внимательнее и мне бросилось в глаза то, что было незаметно на скальном коде:
#ifdef __SSE__
__m128 pixel = _mm_set_ps(0, sRGBToLinear_cache[src[3 * (x+0) + 2]], sRGBToLinear_cache[src[3 * (x+0) + 1]], sRGBToLinear_cache[src[3 * (x+0) + 0]]); for (int i = 0; i < factorsCount; i++) { __m128 basis = _mm_set1_ps(cosYLocal[i] * cosXLocal[i + 0 * factorsCount]); __m128 factor = _mm_loadu_ps(factors[i]); factor = _mm_add_ps(factor, _mm_mul_ps(basis, pixel)); _mm_storeu_ps(factors[i]factor); } #elseFor each element factorsCounton each iteration of the inner loop I load and save factors[i]. I download and save, download and save. How else? We need to somehow summarize them. This is true, but what’s interesting is that for the next pixel we will again load and save the same values. It would be ideal to keep these values in registers, but this will not work, because firstly, the number factorsCount variable. Secondly, it can be quite large, up to 9 * 9, no registers will be enough. But what could be done is to download and save less often. Not for every pixel, but for every fourth, for example. This is not easy to do, but very simple:
static void multiplyBasisFunction(
float factors[][4], int factorsCount, int width, int height, uint8_t *rgb, size_t bytesPerRow,
float *cosX, float *cosY
) {
for(int y = 0; y < height; y++) {
uint8_t *src = rgb + y * bytesPerRow;
float *cosYLocal = cosY + y * factorsCount;
int x = 0;
for(; x < width - 3; x += 4) {
float *cosXLocal = cosX + x * factorsCount;
float pixel0[4] = {sRGBToLinear_cache[src[3 * (x+0) + 0]], sRGBToLinear_cache[src[3 * (x+0) + 1]], sRGBToLinear_cache[src[3 * (x+0) + 2]]};
float pixel1[4] = {sRGBToLinear_cache[src[3 * (x+1) + 0]], sRGBToLinear_cache[src[3 * (x+1) + 1]], sRGBToLinear_cache[src[3 * (x+1) + 2]]};
float pixel2[4] = {sRGBToLinear_cache[src[3 * (x+2) + 0]], sRGBToLinear_cache[src[3 * (x+2) + 1]], sRGBToLinear_cache[src[3 * (x+2) + 2]]};
float pixel3[4] = {sRGBToLinear_cache[src[3 * (x+3) + 0]], sRGBToLinear_cache[src[3 * (x+3) + 1]], sRGBToLinear_cache[src[3 * (x+3) + 2]]};
for (int i = 0; i < factorsCount; i++) {
float basis0 = cosYLocal[i] * cosXLocal[i + 0 * factorsCount];
float basis1 = cosYLocal[i] * cosXLocal[i + 1 * factorsCount];
float basis2 = cosYLocal[i] * cosXLocal[i + 2 * factorsCount];
float basis3 = cosYLocal[i] * cosXLocal[i + 3 * factorsCount];
factors[i][0] += basis0 * pixel0[0] + basis1 * pixel1[0] + basis2 * pixel2[0] + basis3 * pixel3[0];
factors[i][1] += basis0 * pixel0[1] + basis1 * pixel1[1] + basis2 * pixel2[1] + basis3 * pixel3[1];
factors[i][2] += basis0 * pixel0[2] + basis1 * pixel1[2] + basis2 * pixel2[2] + basis3 * pixel3[2];
}
}
for(; x < width; x++) {
float pixel[4];
float *cosXLocal = cosX + x * factorsCount;
pixel[0] = sRGBToLinear_cache[src[3 * x + 0]];
pixel[1] = sRGBToLinear_cache[src[3 * x + 1]];
pixel[2] = sRGBToLinear_cache[src[3 * x + 2]];
for (int i = 0; i < factorsCount; i++) {
float basis = cosYLocal[i] * cosXLocal[i];
factors[i][0] += basis * pixel[0];
factors[i][1] += basis * pixel[1];
factors[i][2] += basis * pixel[2];
}
}
}As a result, the compiler can now twist and optimize this code as it sees fit, and what’s most interesting is that it does an excellent job of it. The result for rosetta x86 is the same as from manual use of SSE internals:
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
975 μs ± 1.81 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)ARM results:
...: assert encode(im, 6, 4) == 'W7EUfLs:M{8yM{wJ$%MepHrE%LV[NaVtBUXSsBNbVuOrM|NaMeXm'
...: %timeit encode(im, 6, 4)
589 µs ± 2.27 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)Which is 37% faster than without the latest optimization.
Total
To be honest, I really wanted to get the coveted hundred. Are there ways to speed things up further? I think there is. First, you can try switching to integer calculations. At various points I tried to do this, and it gave some results, but the code became very complicated and when I opened the next optimization in the main branch, I was too lazy to change it to integers each time. In addition, I think that it may still be possible to keep the working part of the array factors in registers, making several passes over the image. I tried, but judging by the speed, the compiler did not appreciate my idea. You need to definitely switch to internals to control this.
Now I have two pull requests open, to the main library and to the Python library:
https://github.com/woltapp/blurhash/pull/256
https://github.com/woltapp/blurhash-python/pull/25
We haven't heard anything from the authors yet. If anyone has a desire, you can go through the rest of the implementations and roll out these or similar optimizations. In addition, in addition to the encoder, there is also a decoder, which I didn’t even look at. I am sure that it can also be greatly accelerated.
Optimization | x86 GCC 12.2.0 | ARM Clang 14.0.3 | ||
Without | 126 ms | 55.5ms | ||
Python interface | 108 | 1.17x | 47.7 | 1.16x |
Value cache | 3.45 | 36x | 2.65 | 21x |
One pass | 1.44 | 87x | 1.01 | 55x |
Compiler Help | 0.98 | 128x | 0.59 | 94x |
