Графический формат JPEG уменьшает размер изображений без особо заметной для глаза потери качества — упрощая тем самым их хранение и передачу. Студенты из БГУИР — Артём Подгайский, Сергей Буйвид, Юрий Наскевич и Дмитрий Степанчук — в рамках Зимней школы RISC-V YADRO изучили работу декодера JPEG для архитектуры RISC-V, нашли пути для его оптимизации и далее расскажут о своем проекте.

Для начала рассмотрим этапы преобразования изображений. Основной — это дискретное косинус-преобразование Фурье (ДКП), в результате которого растровые данные преобразуются в сумму базисных сигналов.

На иллюстрации левый верхний угол — это низкочастотные, а правый нижний — высокочастотные сигналы. Каждый блок 8x8 кодируемого изображения представлен на рисунке в виде суммы базисных блоков с некоторым коэффициентом. Поскольку человеческий глаз более чувствителен к низкочастотным сигналам, то высокочастотные можно отбросить в результате квантования. Дальнейший порядок операций таков:

1. Из заголовка файла извлекается таблица Хаффмана, которая используется для декодирования битового представления закодированного изображения.

2. Выполняется извлечение коэффициентов ДКП и их деквантование с использованием таблицы квантования, извлеченной из заголовка файла.

3. Блоки изображения 8х8 восстанавливаются при помощи обратного ДКП.

4. Полученные компоненты YCbCr преобразуются в RGB.

5. Полученное RGB изображение сохраняется в память в BBM-формате (bitmap).

Вот так выглядит структура исходной реализации декодера:

├── 3rdparty
│ ├── CMakeLists.txt
│ └── googletest
│
└─ . .
├── CMakeLists.txt
├── include
│ ├── sjpg_bit_stream.h
│ ├── sjpg_decoder.h
│ ├── sjpg_huffman_table.h
│ ├── sjpg_log.h
│ ├── sjpg_segments.h
│ └── TMP1_sjpg_decoder.h
├── README.md
├── resources
│ └── lena.jpg
├── tests
│ ├── CMakeLists.txt
│ ├── test_bit_stream.cpp
│ ├── test_huffman_table.cpp
│ └── test_jpeg_decoder.cpp

Мы оценили его работу на RISC-V с помощью профайлера и получили крайне неудовлетворительный результат:

Почти 90% времени программы занимает ДКП, а в нем — вычисление косинусов:

std::vector<float> idctNaive (const std::vector<int16_t> &data) {
  std::vector<float> result(8 * 8, 0.0f);
    for (auto y = 0; y < 8; ++y) {
      for (auto x = 0; x < 8; ++x) {
        auto sum = 0.0f;
        for (auto u = 0; u < 8; ++u) {
          for (auto v = 0; v < 8; ++v) {
            float cu = (u == 0) ? 1.0f / std::sqrt(2.0f) : 1.0f;
            float cv = (v == 0) ? 1.0f / std::sqrt(2.0f) : 1.0f;
            float t0 = cu * std::cos((2 * x + 1) * u * M_PI / 16.0);
            float t1 = cv * std::cos((2 * y + 1) * v * M_PI / 16.0);
            // colum major
            auto data_value = data[u * 8 + v];
            sum += (data_value * t0 * t1);
          }
        }
        sum *= 0.25;
        result[x * 8 + y] = sum;
      }
    }
    return result;
}

Для ускорения удобней всего провести мемоизацию, то есть создать таблицу косинусов, из которой потом можно брать готовые значения. Постараемся максимально оптимизировать ДКП. Сначала мы создали возможность выбора функции DKP, добавив в структуру sjpg_idct.h:

template <typename T>
  using idct_call_t =
    std::function<std::vector<T>
    (const std::vector<int16_t> &)>;

В sjpg_decoder.h прописали:

void setIdctCall (idct_call_t <idct_num_t>_idctCall){
  idctCall = _idctCall
}

Так выглядит таблица косинусов в sjpg_idct.h:

static constexpr float idctCoefficients [8][8] = {


  { 0.707106781f, 0.980785280f, 0.923879533f, 0.831469612f, 0.707106781f, 0.555570233f, 0.382683432f, 0.195090322f },
  { 0.707106781f, 0.831469612f, 0.382683432f, 0.195090322f, 0.707106781f, 0.980785280f, 0.923879533f, 0.555570233f },
  { 0.707106781f, 0.555570233f, 0.382683432f, 0.980785280f, 0.707106781f, 0.195090322f, 0.923879533f, 0.831469612f },
  { 0.707106781f, 0.195090322f, 0.923879533f, 0.555570233f, 0.707106781f, 0.831469612f, 0.382683432f, 0.980785280f },
  { 0.707106781f, 0.195090322f, 0.923879533f, 0.555570233f, 0.707106781f, 0.831469612f, 0.382683432f, 0.980785280f },
  { 0.707106781f, 0.555570233f, 0.382683432f, 0.980785280f, 0.707106781f, 0.195090322f, 0.923879533f, 0.831469612f },
  { 0.707106781f, 0.831469612f, 0.382683432f, 0.195090322f, 0.707106781f, 0.980785280f, 0.923879533f, 0.555570233f },
  { 0.707106781f, 0.980785280f, 0.923879533f, 0.831469612f, 0.707106781f, 0.555570233f, 0.382683432f, 0.195090322f },


};

Пора начать исследование производительности. Вот конфигурация тестового стенда:

После внедрения мемоизации оценили прогресс:

Производительность декодера выросла многократно, но этого еще недостаточно, чтобы ДКП работало, например, в промышленных библиотеках. Здесь нужна более интегральная оптимизация. Рассмотрим одномерное дискретное косинусное преобразование как матричное преобразование:

Если матрицу косинусного преобразования преобразовать в произведение нескольких разреженных матриц, то мы сможем значительно уменьшить количество операций, которые будут использоваться для ДКП. Эта идея была предложена Кристофером Лёффлером в статье 1989 года Practical fast 1-D DCT algorithms with 11 multiplications.

Наша матрица ДКП факторизуется во множество разреженных матриц. Их можно представить в коде не алгоритмическим матричным умножением, а разбить на конкретный набор операций для каждой разреженной матрицы. Таким образом, из матрицы ДКП можно получить матрицу обратного ДКП простым взятием инверсной матрицы.

Почти очевидно, что составляющие нашей матрицы ДКП являются ортогональными, а значит, обратные матрицы совпадают с транспонированными. Такая матричная схема не очень удобна, поэтому Лёффлер в статье использовал схему с алгоритмическими конструкциями. Она делает реализацию более наглядной.

Преодолев опечатки и недосказанности статьи, мы это реализовали. Следующий шаг оптимизации — сократить количество умножений. Вот операция вращения из статьи Лёффлера:

Она занимает четыре умножения и два сложения, но ее можно свести к трем умножениям и трем сложениям. В общем случае умножение занимает больше времени, чем сложение, поэтому так мы ускорим вычисление. Используя формулы суммы синусов и косинусов, мы можем преобразовать выражение

в выражение

Оценим производительность снова:

Оптимизация Лёффлера позволила еще раз многократно увеличить производительность ДКП. На Lichee Pi -O2 ДКП теперь занимает чуть больше одной миллисекунды. Это уже позволяет использовать решение в промышленной библиотеке, поэтому одна из реализаций ДКП в известной open source-библиотеке libjpeg — это и есть схема Лёффлера.

Изначально она у нас была реализована с помощью чисел с плавающей точкой, но это имеет свои недостатки:

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

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

• В конце алгоритма необходимо умножить все значения на 1/8, что добавляет дополнительные 64 умножения — на целых числах их можно заменить битовыми сдвигами.

• Главная проблема — мы зависим от наличия FPU, который может быть далеко не везде.

Поэтому мы сделали реализацию на целых числах:

#define PRECISION 8


constexpr int32_t FIX(double x){
  return ((x) * (1 << PRECISION) + 0.5);
}


#define FIX_SQRT2 FIX(1.4142135623730951)
#define FIX_SIN_C1 -FIX(0.19509032201612825)
#define FIX_COS_C1 FIX(0.9807852804032304)
#define FIX_SIN_C3 -FIX(0.5555702330196022)
#define FIX_COS_C3 FIX(0.8314696123025452)
#define FIX_SQRT2_SIN_C6 -FIX(1.3065629648763766)
#define FIX_SQRT2_COS_C6 FIX(0.5411961001461971)


#define DESCALE(x, n) (((x) + (1 << ((n) 1))) >> (n))

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

ptr [0] = DESCALE((s2_0 + s2_7), (PRECISION + 3));
ptr [7] = DESCALE((s2_0 - s2_7), (PRECISION + 3));
ptr [1] = DESCALE((s2_1 + s2_6), (PRECISION + 3));
ptr [6] = DESCALE((s2_1 - s2_6), (PRECISION + 3));
ptr [2] = DESCALE((s2_2 + s2_5), (PRECISION + 3));
ptr [5] = DESCALE((s2_2 - s2_5), (PRECISION + 3));
ptr [3] = DESCALE((s2_3 + s2_4), (PRECISION + 3));
ptr [4] = DESCALE((s2_3 - s2_4), (PRECISION + 3));

Оценим производительность после целочисленной реализации схемы Лёффлера:

Почему производительность повысилась только во втором случае? Сначала мы думали, причина в том, что наше тестовое ядро на FPGA было 32-битным, а реализация использовала 64-битные числа. Но потом пришли к тому, что если инструкции RISC-V занимают 32 бита, то для записи в них константы места остается меньше. Загрузка 32-битной константы в регистр занимает две инструкции, а 64–битной — еще больше.

С 32-битной реализацией схемы Лёффлера нам удалось повысить производительность в обоих baremetal-тулчейнах более чем в два раза. На Lichee Pi победа осталась за реализацией с плавающей точкой:

Вот итоговые результаты в более удобном виде:

Самая простая имплементация с мемоизацией работает во много раз дольше, чем схема Лёффлера. Ее различные вариации — с 64-битной, 32-битной арифметикой или с плавающей точкой — могут быть чуть быстрее или медленнее, в зависимости от конкретной конфигурации.

Николай Петровский

куратор проекта, к. т. н., доцент кафедры ЭВС БГУИР

Проект простейшего JPEG-декодера, выполненный студентами БГУИР, имеет практическую и теоретическую значимость для области компрессии изображений. Рассматриваемый декодер реализует часть стандарта JPEG без режима прогрессивного сжатия. Референсная реализация выполняет все этапы реального декодера — от энтропийного сжатия, переупорядочивания коэффициентов и обратного дискретно-косинусного преобразования второго типа для цветного изображения.

Сложность проекта состояла в том, что наивная реализация ДКП-преобразования с применением чисел с плавающей точкой выполнялась на софт-процессоре RISC-V несколько часов для тестового изображения разрешением 512х512. Однако после реализации на основе фиксированной запятой и аппроксимации схемой Лёффлера производительность увеличилась в десятки раз.

Участники Зимней Школы RISC-V в БГУИР работали с софт-процессорным ядром RISC-V, реализованным на отладочной FPGA-плате. Это позволило им не только изучить архитектуру RISC-V, но и получить практический опыт работы с мультипроцессорной асимметричной системой.

Максим Вашкевич

куратор проекта, д. т. н., доцент, профессор кафедры ЭВС БГУИР

Еще одной целью проекта было показать студентам, как на первый взгляд «сухие» математические формулы становятся незаменимыми помощниками в решении сложных инженерных задач. Например, в оригинальной статье Лёффлера приведен быстрый алгоритм для прямого восьмиточечного ДКП, тогда как в декодере JPEG требуется реализация обратного ДКП. Поначалу это вызвало затруднения, однако мы подсказали ребятам, что матрица ДКП является ортогональной. А это значит, что обратная к ней матрица совпадает с транспонированной матрицей прямого ДКП.

Затем мы разобрали еще один важный момент: оказалось, что быстрый алгоритм ДКП можно представить в виде произведения разреженных матриц. Следовательно, обратный алгоритм может быть выражен как произведение тех же матриц, но в обратном порядке и в транспонированной форме.

Другие проекты Зимней школы RISC-V сезона 2024–2025:

• Измерение латентности и пропускной способности, оптимизация приложений на C#

• Программный генератор случайных числовых последовательностей на RISC-V с использованием PUF в DRAM

• Оптимизация функций компьютерного зрения (библиотека OpenCV) для RISC-V