Blur в Delphi. Часть I: Gaussian Blur

Размытие изображения (blur) — одна из базовых операций в обработке графики. Без неё не обходятся ни повышение резкости (Unsharp Mask), ни выделение краёв (Difference of Gaussians), ни шумоподавление, ни создание эффектов глубины резкости. Gaussian blur — эталон качества. Набор его свойств уникален и ни одно другое ядро его не повторяет. Именно поэтому он стал стандартом.

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

Содержание скрыть

Что такое Blur

Gaussian: Наивная реализация

Ядро свёртки

Строим колокол

Наивная 2D-свёртка: Считаем «в лоб»

Gaussian: Сепарабельная реализация

Сепарабельность: Что это такое и почему работает

Два быстрых прохода вместо одного долгого

Fixed-Point: Целочисленная арифметика

На чём основан фокус

Ядро должно быть целочисленным

Свёртка на целых числах

GBlur2: Легенда Delphi-форумов

Как определяется радиус (K.Size)

Почему это «вольная интерпретация»?

Выводы

Что ещё бросается в глаза:

Gaussian: Оптимизация, выжимаем максимум

Суть оптимизации

Транспонирование

Горизонтальный проход

Основная процедура

Системные средства: GDI+ и Direct2D

GDI+ Blur

Direct2D Gaussian Blur

Зачем включать в сравнение

Бенчмарк

Методика

Диаграммы

32-bit

64-bit

Разбор результатов

32-bit vs 64-bit: кому помогает разрядность?

Что из этого следует

А что насчёт рекурсивных (IIR) фильтров?

Основные алгоритмы

Почему они не вошли в эту статью?

Скачать

Что такое Blur

В общем виде размытие (blur) — это замена каждого пикселя изображения на взвешенное среднее его соседей. Чем больше радиус, тем больше соседей участвует и тем сильнее размытие.

Result = Σ(pixel[i] × weight[i]) / Σ(weight[i])

1	Result = Σ(pixel[i] × weight[i]) / Σ(weight[i])

Математически — это операция свёртки (convolution): на изображение накладывается матрица весов — ядро (kernel), и для каждого пикселя вычисляется сумма произведений значений соседних пикселей на соответствующие веса:

Pixel'[x, y] = Σ(i = -R..R) Σ(j = -R..R) Pixel[x+i, y+j] × Kernel[i, j]

1	Pixel'[x, y] = Σ(i = -R..R) Σ(j = -R..R) Pixel[x+i, y+j] × Kernel[i, j]

Где:

P'(x, y) — значение пикселя после фильтрации;
P(x + i, y + j) — значения исходных пикселей в окрестности R;
K(i, j) — коэффициент ядра (матрицы) свертки в позиции (i, j);
R — радиус ядра (так, что размер матрицы равен (2R+1) × (2R+1)).

Все алгоритмы размытия отличаются по сути только одним — формой ядра:

Тип	Форма ядра	Профиль
Box Blur	Все веса одинаковы	Прямоугольник
Gaussian Blur	Веса по кривой Гаусса	Колокол
Stack Blur	Веса линейно убывают от центра	Треугольник

Box — самый простой: все соседи равноправны, результат — простое среднее. Быстро, но на контрастных границах заметны прямоугольные артефакты.

Gaussian — эталон качества: ближние пиксели влияют сильнее, дальние — слабее, по гауссовой кривой. Размытие плавное и естественное, без артефактов.

Stack Blur — компромисс: треугольное ядро визуально близко к Gaussian blur, но вычисляется проще.

В этой статье мы начнём с Gaussian Blur — эталонного по качеству алгоритма — и шаг за шагом оптимизируем его от наивной 2D-свёртки до максимально быстрой реализации. Затем, во второй части, рассмотрим альтернативы: Box Blur и его трюк с тройным проходом, Stack Blur и даже хитрость с уменьшением-увеличением изображения. В финале — честный бенчмарк всех методов.

Почему именно Gaussian? Он сепарабелен (2D раскладывается в два 1D-прохода), изотропен (нет направленных артефактов), не создаёт ореолов (все веса положительные) и каскадируем (два размытия подряд с двумя малыми σ = одно с бОльшим σ). Оптический расфокус, рассеивание тепла, дым — всё это Gaussian. Поэтому глаз воспринимает такое размытие как «настоящее». Все остальные методы, упоминаемые выше, это быстрые приближения к нему.

Gaussian: Наивная реализация

Углубляться в математическую теорию не будем — наша задача практическая: изучить существующие алгоритмы размытия, сравнить их и выбрать оптимальные для дальнейшей работы. Но пару слов о ядре Гаусса сказать необходимо.

Гауссово размытие строится на знаменитой «колоколообразной» кривой нормального распределения. Для каждого соседнего пикселя на расстоянии (i, j) от центра вес вычисляется по формуле:

G(i, j) = (1 / 2πσ²) × exp(−(i² + j²) / 2σ²)

1	G(i, j) = (1 / 2πσ²) × exp(−(i² + j²) / 2σ²)

Где:

i, j — дискретные координаты в матрице (целые числа, где 0, 0 — центр).
σ (сигма) — параметр сглаживания.

Параметр сглаживания σ — среднеквадратичное отклонение нормального распределения — определяет «ширину колокола»: чем он больше, тем больше радиус размытия и тем сильнее эффект. На практике ядро ограничивают радиусом r = ⌈3σ⌉ — за этой границей веса настолько малы, что ими можно пренебречь. Например, при σ=2 получаем ядро 13×13 (радиус 6), при σ=5 — уже 31×31.

Наивная реализация делает ровно то, что написано в формуле: для каждого пикселя проходит по всему ядру (2r+1)×(2r+1) и суммирует взвешенные значения соседей. Это O(r²) операций на пиксель — при σ=5 это 961 умножение и сложение на каждый пиксель, при σ=50 почти 91 000. На изображении 600×600 при σ=50 это порядка 33 миллиардов операций (32 616 360 000, я зануда). Работает, качество эталонное, но скорость… мягко говоря, оставляет желать лучшего.

Ядро свёртки

Для наглядности посмотрим, как выглядит ядро Гаусса при σ=1 (радиус 3, матрица 7×7):

i \ j	-3	-2	-1	0	1	2	3
-3	0.0000	0.0002	0.0011	0.0018	0.0011	0.0002	0.0000
-2	0.0002	0.0029	0.0131	0.0216	0.0131	0.0029	0.0002
-1	0.0011	0.0131	0.0586	0.0966	0.0586	0.0131	0.0011
0	0.0018	0.0216	0.0966	0.1592	0.0966	0.0216	0.0018
1	0.0011	0.0131	0.0586	0.0966	0.0586	0.0131	0.0011
2	0.0002	0.0029	0.0131	0.0216	0.0131	0.0029	0.0002
3	0.0000	0.0002	0.0011	0.0018	0.0011	0.0002	0.0000

Особенности графика:

Центральный пик: Максимальное значение находится в точке (0, 0). Это «сердце» фильтра, которое определяет основной цвет пикселя при размытии.
Симметрия: График идеально симметричен относительно центра. Это означает, что размытие будет происходить равномерно во все стороны.
Экспоненциальный спад: Видно, как значения «прижимаются» к нулю уже на расстоянии 3-х пикселей от центра. Именно поэтому для σ=1 матрицы 7×7 (радиус 3) более чем достаточно — за её пределами веса становятся ничтожно малыми.

На графике представлена нормализованная функция Гаусса. В непрерывном случае объём под колоколом в точности равен 1. Наша дискретная матрица 7×7 покрывает ~99% этого объёма — оставшийся процент приходится на отсечённые хвосты за пределами радиуса 3σ.

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

Строим колокол

Функцию, которую мы обсуждали ранее, можно записать в развернутом (сепарабельном) виде как произведение двух 1D-фильтров:

Почему такая запись возможна, поговорим чуть позже. А пока основной вывод таков, что мы не будем хранить всё ядро целиком, а будем хранить ядро для 1D, чтобы затем перемножать соответствующие (по индексу) веса.

Для ядра заводим простую структуру:

type
  TKernel1D = record
    Radius: Integer;
    Weights: array of Double;
    // Weights[0] = центр, Weights[i] = смещение ±i
    // Симметричное ядро: храним только половину + центр
  end;

type

TKernel1D = record

Radius: Integer;

Weights: array of Double;

// Weights[0] = центр, Weights[i] = смещение ±i

// Симметричное ядро: храним только половину + центр

end;

Формируем ядро следующим образом:

function MakeGaussianKernel1D(Sigma: Double): TKernel1D;
var I: Integer; Sum, W: Double;
begin
  // Радиус = 3*Sigma, округляем вверх. Минимум 1
  Result.Radius := Max(1, Ceil(Sigma * 3));
  SetLength(Result.Weights, Result.Radius + 1);

  // Вычисляем веса
  Sum := 0;
  for I := 0 to Result.Radius do
  begin
    W := Exp(-(I * I) / (2 * Sigma * Sigma));
    Result.Weights[I] := W;
    if I = 0 then
      Sum := Sum + W     // центральный - один раз
    else
      Sum := Sum + 2*W;  // симметричные - дважды
  end;

  // Нормализуем, чтобы сумма = 1
  for I := 0 to Result.Radius do
    Result.Weights[I] := Result.Weights[I] / Sum;
end;

function MakeGaussianKernel1D(Sigma: Double): TKernel1D;

var I: Integer; Sum, W: Double;

begin

// Радиус = 3*Sigma, округляем вверх. Минимум 1

Result.Radius := Max(1, Ceil(Sigma * 3));

SetLength(Result.Weights, Result.Radius + 1);

// Вычисляем веса

Sum := 0;

for I := 0 to Result.Radius do

begin

W := Exp(-(I * I) / (2 * Sigma * Sigma));

Result.Weights[I] := W;

if I = 0 then

Sum := Sum + W // центральный - один раз

else

Sum := Sum + 2*W; // симметричные - дважды

end;

// Нормализуем, чтобы сумма = 1

for I := 0 to Result.Radius do

Result.Weights[I] := Result.Weights[I] / Sum;

end;

Пояснения:

Радиус Ceil(Sigma * 3) — это стандартное отсечение. По правилу трёх сигм вес гауссианы составляет менее 0.4% от центрального — этим можно пренебречь.
Массив весов Weights[0..Radius] хранит только правую половину колокола вместе с центром. Левая идентична: G(-x) = G(x).
Коэффициент 1/(2πσ²) из формулы опущен — он одинаков для всех элементов и сократится при нормализации.
Нормализация гарантирует, что сумма всех весов ядра равна 1. Без неё изображение станет темнее или светлее. Центральный вес входит в сумму один раз, остальные — дважды (слева и справа от центра).

Наивная 2D-свёртка: Считаем «в лоб»

Имея ядро, реализуем размытие буквально по определению: для каждого пикселя обходим окно (2r+1)×(2r+1) и суммируем взвешенные значения соседей. Двумерный вес получаем из нашего одномерного ядра: W(i, j) = K[Abs(i)] * K[Abs(j)]. Это корректно, потому что функция Гаусса сепарабельна (но это свойство подробно рассмотрим и воспользуемся им в следующем разделе, а пока считаем «в лоб»).

procedure BlurGaussian2D(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  K: TKernel1D;
  W, H, X, Y, KX, KY: Integer;
  SX, SY: Integer;
  Weight, TotalWeight: Double;
  SumR, SumG, SumB: Double;
  SrcCache, DstCache: TScanLineCache;
  SrcPx: TPixel32;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  // Строим 1D-ядро, потом используем как 2D = K[x] * K[y]
  K := MakeGaussianKernel1D(Sigma);

  Bitmap.PixelFormat := pf32bit;
  W := Bitmap.Width;
  H := Bitmap.Height;

  Result := TBitmap.Create;
  Result.PixelFormat := pf32bit;
  Result.SetSize(W, H);

  SrcCache := BuildScanLineCache(Bitmap);
  DstCache := BuildScanLineCache(Result);

  for Y := 0 to H - 1 do
  begin
    for X := 0 to W - 1 do
    begin
      SumR := 0; SumG := 0; SumB := 0;
      TotalWeight := 0;

      for KY := -K.Radius to K.Radius do
      begin
        SY := Clamp(Y + KY, 0, H - 1);

        for KX := -K.Radius to K.Radius do
        begin
          SX := Clamp(X + KX, 0, W - 1);

          // 2D-вес = 1D(KX) * 1D(KY) - сепарабельность,
          // но мы НАРОЧНО считаем «в лоб»
          Weight := K.Weights[Abs(KX)] * K.Weights[Abs(KY)];

          SrcPx := SrcCache[SY][SX];
          SumB := SumB + SrcPx.B * Weight;
          SumG := SumG + SrcPx.G * Weight;
          SumR := SumR + SrcPx.R * Weight;
          TotalWeight := TotalWeight + Weight;
        end;
      end;

      // TotalWeight = 1.0, но для пикселей у края может чуть отличаться
      with DstCache[Y][X] do
      begin
        B := Clamp(Round(SumB / TotalWeight), 0, 255);
        G := Clamp(Round(SumG / TotalWeight), 0, 255);
        R := Clamp(Round(SumR / TotalWeight), 0, 255);
        A := SrcCache[Y][X].A; // альфу не трогаем
      end;
    end;
  end;
end;

procedure BlurGaussian2D(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Sigma: Double;

K: TKernel1D;

W, H, X, Y, KX, KY: Integer;

SX, SY: Integer;

Weight, TotalWeight: Double;

SumR, SumG, SumB: Double;

SrcCache, DstCache: TScanLineCache;

SrcPx: TPixel32;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

// Строим 1D-ядро, потом используем как 2D = K[x] * K[y]

K := MakeGaussianKernel1D(Sigma);

Bitmap.PixelFormat := pf32bit;

W := Bitmap.Width;

H := Bitmap.Height;

Result := TBitmap.Create;

Result.PixelFormat := pf32bit;

Result.SetSize(W, H);

SrcCache := BuildScanLineCache(Bitmap);

DstCache := BuildScanLineCache(Result);

for Y := 0 to H - 1 do

begin

for X := 0 to W - 1 do

begin

SumR := 0; SumG := 0; SumB := 0;

TotalWeight := 0;

for KY := -K.Radius to K.Radius do

begin

SY := Clamp(Y + KY, 0, H - 1);

for KX := -K.Radius to K.Radius do

begin

SX := Clamp(X + KX, 0, W - 1);

// 2D-вес = 1D(KX) * 1D(KY) - сепарабельность,

// но мы НАРОЧНО считаем «в лоб»

Weight := K.Weights[Abs(KX)] * K.Weights[Abs(KY)];

SrcPx := SrcCache[SY][SX];

SumB := SumB + SrcPx.B * Weight;

SumG := SumG + SrcPx.G * Weight;

SumR := SumR + SrcPx.R * Weight;

TotalWeight := TotalWeight + Weight;

end;

// TotalWeight = 1.0, но для пикселей у края может чуть отличаться

with DstCache[Y][X] do

begin

B := Clamp(Round(SumB / TotalWeight), 0, 255);

G := Clamp(Round(SumG / TotalWeight), 0, 255);

R := Clamp(Round(SumR / TotalWeight), 0, 255);

A := SrcCache[Y][X].A; // альфу не трогаем

end;

Пару слов о вспомогательных типах. TScanLineCache — массив указателей на строки, заполненный один раз через ScanLine[Y]. Без него пришлось бы вызывать ScanLine на каждое обращение к пикселю, а внутри вложенного цикла это катастрофа. Мы подробно разбирали преимущества кэша ScanLine в предыдущей статье. TPixel32 — запись из четырёх байт (B, G, R, A), соответствующая формату pf32bit. Clamp — ограничение значения диапазоном, используется и для координат (обработка краёв), и для итоговых значений каналов (0…255).

type
  // Пиксель 32-бит: B, G, R, A - порядок как в TBitmap pf32bit
  PPixel32 = ^TPixel32;
  TPixel32 = packed record
    B, G, R, A: Byte;
  end;

  // Строка пикселей
  PPixel32Array = ^TPixel32Array;
  TPixel32Array = array[0..0] of TPixel32;

  // Кэш ScanLine - чтобы не вызывать ScanLine[Y] повторно
  TScanLineCache = array of PPixel32Array;

// Обрезка координаты до границ
function Clamp(Value, MinVal, MaxVal: Integer): Integer; inline;
begin
  if Value < MinVal then
    Result := MinVal
  else if Value > MaxVal then
    Result := MaxVal
  else
    Result := Value;
end;

type

// Пиксель 32-бит: B, G, R, A - порядок как в TBitmap pf32bit

PPixel32 = ^TPixel32;

TPixel32 = packed record

B, G, R, A: Byte;

end;

// Строка пикселей

PPixel32Array = ^TPixel32Array;

TPixel32Array = array[0..0] of TPixel32;

// Кэш ScanLine - чтобы не вызывать ScanLine[Y] повторно

TScanLineCache = array of PPixel32Array;

// Обрезка координаты до границ

function Clamp(Value, MinVal, MaxVal: Integer): Integer; inline;

begin

if Value < MinVal then

Result := MinVal

else if Value > MaxVal then

Result := MaxVal

else

Result := Value;

end;

В результате видим эталонное качество, идентичное любой академической реализации. Но скорость удручает: при σ=5 (ядро 31×31) на изображение 600×600 уходит порядка 3 секунд, а при σ=10 уже около 10 секунд. Сложность O(r²) на пиксель делает этот подход непригодным для практического использования. Ожидание результата при σ=50 тянется вечность… (см. рисунок выше). Картинки кликабельны.

К счастью, функция Гаусса обладает свойством, которое мы уже неявно использовали — сепарабельностью. Именно она превратит четыре вложенных цикла в два прохода по два, и позволит превратить O(r²) в O(r). Этим мы и воспользуемся в следующем разделе.

Gaussian: Сепарабельная реализация

Наивная 2D-свёртка работает, но при росте σ становится катастрофически медленной. Однако мы уже знаем, что двумерное ядро Гаусса раскладывается в произведение двух одномерных. Осталось воспользоваться этим на практике: вместо одного прохода с ядром (2r+1)×(2r+1) сделать два последовательных прохода с ядром (2r+1)×1 — сначала по горизонтали, затем по вертикали. Результат математически идентичен, а скорость совсем другая.

Сепарабельность: Что это такое и почему работает

Секрет таится в алгебре экспоненты. Раскроем показатель степени в формуле двумерного Гаусса:

G(i, j) = (1/2πσ²) × exp(−(i² + j²)/2σ²)

1	G(i, j) = (1/2πσ²) × exp(−(i² + j²)/2σ²)

Разделим суммы в показателе на степени:

G(i, j) = (1/2πσ²) × exp(- i²/2σ² - j²/2σ²)

1	G(i, j) = (1/2πσ²) × exp(- i²/2σ² - j²/2σ²)

Сумму i² + j² в показателе можно преобразовать, используя базовое свойство экспоненты: exp(a+b) = exp(a) × exp(b):

exp(−(i² + j²) / 2σ²) = exp(−i² / 2σ²) × exp(−j² / 2σ²)

1	exp(−(i² + j²) / 2σ²) = exp(−i² / 2σ²) × exp(−j² / 2σ²)

Таким образом получаем формулу:

Константу 1 / 2πσ² тоже раскладываем:

1 / 2πσ² = (1 / √(2π)σ) × (1 / √(2π)σ)

1	1 / 2πσ² = (1 / √(2π)σ) × (1 / √(2π)σ)

Подставляем обратно и получаем:

G(i, j) = (1/2πσ²) × exp(−(i² + j²)/2σ²) = 
(1/√(2π)σ) × (1/√(2π)σ) × exp(−i²/2σ²) × exp(−j²/2σ²) = 
(1/√(2π)σ) × exp(−i²/2σ²) × (1/√(2π)σ) × exp(−j²/2σ²) = 
G(i) × G(j)

G(i, j) = (1/2πσ²) × exp(−(i² + j²)/2σ²) =

(1/√(2π)σ) × (1/√(2π)σ) × exp(−i²/2σ²) × exp(−j²/2σ²) =

(1/√(2π)σ) × exp(−i²/2σ²) × (1/√(2π)σ) × exp(−j²/2σ²) =

G(i) × G(j)

Таким образом, получаем сепарабельность:

Где одномерное ядроK(t) определяется как:

K(t) = (1/√(2π)σ) × exp(−t²/2σ²)

1	K(t) = (1/√(2π)σ) × exp(−t²/2σ²)

На этом графике мы видим, как два независимых одномерных «колокола» (красный вдоль осиi и синий вдоль оси j) при пересечении буквально «вырезают» в пространстве двумерную поверхность.

Двумерное ядро математически точно раскладывается в произведение двух одномерных без потери качества. Именно поэтому в наивной реализации строка Weight := K.Weights[Abs(KX)] * K.Weights[Abs(KY)] давала точный результат.

А раз так — зачем обходить всё двумерное ядро за один проход? Можно сделать два последовательных одномерных прохода: сначала размыть каждую строку ядром K, затем каждый столбец тем же ядром. Вместо (2r+1)² операций на пиксель получаем 2×(2r+1) — при σ=10 это разница между 3721 и 122 операциями, в 30 раз меньше.

Посмотрим, как это выглядит в коде.

Два быстрых прохода вместо одного долгого

Идея проста: заводим промежуточный буфер Tmp. Первый проход размывает исходное изображение по горизонтали — каждый пиксель заменяется взвешенной суммой соседей вдоль строки. Второй проход берёт результат из Tmp и размывает его по вертикали. Два одномерных размытия подряд дают результат, идентичный одному двумерному.

procedure BlurGaussianSeparable(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  K: TKernel1D;
  W, H, X, Y, I: Integer;
  SX, SY: Integer;
  SumR, SumG, SumB, Weight: Double;
  SrcCache, TmpCache, DstCache: TScanLineCache;
  Tmp: TBitmap;
  SrcPx: TPixel32;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  K := MakeGaussianKernel1D(Sigma);

  Bitmap.PixelFormat := pf32bit;
  W := Bitmap.Width;
  H := Bitmap.Height;

  // Промежуточный буфер (после горизонтального прохода)
  Tmp := TBitmap.Create;
  try
    Tmp.PixelFormat := pf32bit;
    Tmp.SetSize(W, H);

    // Финальный результат (после вертикального прохода)
    Result := TBitmap.Create;
    Result.PixelFormat := pf32bit;
    Result.SetSize(W, H);

    SrcCache := BuildScanLineCache(Bitmap);
    TmpCache := BuildScanLineCache(Tmp);
    DstCache := BuildScanLineCache(Result);

    // Проход 1: горизонтальный
    for Y := 0 to H - 1 do
    begin
      for X := 0 to W - 1 do
      begin
        SumR := 0; SumG := 0; SumB := 0;

        for I := -K.Radius to K.Radius do
        begin
          SX := Clamp(X + I, 0, W - 1);
          Weight := K.Weights[Abs(I)];

          SrcPx := SrcCache[Y][SX];
          SumB := SumB + SrcPx.B * Weight;
          SumG := SumG + SrcPx.G * Weight;
          SumR := SumR + SrcPx.R * Weight;
        end;

        with TmpCache[Y][X] do
        begin
          B := Clamp(Round(SumB), 0, 255);
          G := Clamp(Round(SumG), 0, 255);
          R := Clamp(Round(SumR), 0, 255);
          A := SrcCache[Y][X].A;
        end;
      end;
    end;

    //  Проход 2: вертикальный
    for Y := 0 to H - 1 do
    begin
      for X := 0 to W - 1 do
      begin
        SumR := 0; SumG := 0; SumB := 0;

        for I := -K.Radius to K.Radius do
        begin
          SY := Clamp(Y + I, 0, H - 1);
          Weight := K.Weights[Abs(I)];

          SrcPx := TmpCache[SY][X];
          SumB := SumB + SrcPx.B * Weight;
          SumG := SumG + SrcPx.G * Weight;
          SumR := SumR + SrcPx.R * Weight;
        end;

        with DstCache[Y][X] do
        begin
          B := Clamp(Round(SumB), 0, 255);
          G := Clamp(Round(SumG), 0, 255);
          R := Clamp(Round(SumR), 0, 255);
          A := TmpCache[Y][X].A;
        end;
      end;
    end;

  finally
    Tmp.Free;
  end;
end;

procedure BlurGaussianSeparable(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Sigma: Double;

K: TKernel1D;

W, H, X, Y, I: Integer;

SX, SY: Integer;

SumR, SumG, SumB, Weight: Double;

SrcCache, TmpCache, DstCache: TScanLineCache;

Tmp: TBitmap;

SrcPx: TPixel32;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

K := MakeGaussianKernel1D(Sigma);

Bitmap.PixelFormat := pf32bit;

W := Bitmap.Width;

H := Bitmap.Height;

// Промежуточный буфер (после горизонтального прохода)

Tmp := TBitmap.Create;

try

Tmp.PixelFormat := pf32bit;

Tmp.SetSize(W, H);

// Финальный результат (после вертикального прохода)

Result := TBitmap.Create;

Result.PixelFormat := pf32bit;

Result.SetSize(W, H);

SrcCache := BuildScanLineCache(Bitmap);

TmpCache := BuildScanLineCache(Tmp);

DstCache := BuildScanLineCache(Result);

// Проход 1: горизонтальный

for Y := 0 to H - 1 do

begin

for X := 0 to W - 1 do

begin

SumR := 0; SumG := 0; SumB := 0;

for I := -K.Radius to K.Radius do

begin

SX := Clamp(X + I, 0, W - 1);

Weight := K.Weights[Abs(I)];

SrcPx := SrcCache[Y][SX];

SumB := SumB + SrcPx.B * Weight;

SumG := SumG + SrcPx.G * Weight;

SumR := SumR + SrcPx.R * Weight;

end;

with TmpCache[Y][X] do

begin

B := Clamp(Round(SumB), 0, 255);

G := Clamp(Round(SumG), 0, 255);

R := Clamp(Round(SumR), 0, 255);

A := SrcCache[Y][X].A;

end;

// Проход 2: вертикальный

for Y := 0 to H - 1 do

begin

for X := 0 to W - 1 do

begin

SumR := 0; SumG := 0; SumB := 0;

for I := -K.Radius to K.Radius do

begin

SY := Clamp(Y + I, 0, H - 1);

Weight := K.Weights[Abs(I)];

SrcPx := TmpCache[SY][X];

SumB := SumB + SrcPx.B * Weight;

SumG := SumG + SrcPx.G * Weight;

SumR := SumR + SrcPx.R * Weight;

end;

with DstCache[Y][X] do

begin

B := Clamp(Round(SumB), 0, 255);

G := Clamp(Round(SumG), 0, 255);

R := Clamp(Round(SumR), 0, 255);

A := TmpCache[Y][X].A;

end;

finally

Tmp.Free;

end;

Обратите внимание: оба прохода структурно идентичны — меняется только направление. В горизонтальном проходе Y фиксирован, сдвигается X + I; в вертикальном — наоборот, X фиксирован, сдвигается Y + I. Ядро K одно и то же, функция MakeGaussianKernel1D вызывается один раз.

Однако за оптимизацию взимается плата — дополнительный битмап Tmp размером с исходное изображение. Это несколько мегабайт для типичных размеров. Но это ничтожная плата за ускорение. Промежуточные значения округляются до byte после первого прохода, что вносит минимальную погрешность — на практике результат визуально неотличим от наивной 2D-свёртки, а попиксельная разница не превышает ±1.

Сравним с наивным результатом. Результаты подтверждают теорию: при σ=5 (219 против 2591 мс) ускорение примерно в 12 раз, при σ=50 (2015 против 245383 мс) — в 122 раза. Данные показывают порядок увеличения скорости на этих конкретных примерах. Более взвешенные данные увидим в бенчмарке. Сложность теперь линейна по радиусу, но мы по-прежнему работаем с Double, вызываем Clamp и Abs в горячем цикле, а промежуточные значения округляем до byte. Ещё есть куда стремиться. В следующем разделе избавимся от вещественной арифметики и посмотрим, какой выигрыш это даст.

Fixed-Point: Целочисленная арифметика

Алгоритм остаётся тем же — сепарабельная свёртка в два прохода. Меняется только способ вычислений: вместо Double используем Integer. Зачем? Целочисленные операции — умножение и сложение — на порядок дешевле вещественных, особенно на тех платформах, где нет аппаратного FPU. Но даже на современных десктопах разница ощутима: меньше тактов на операцию, лучше предсказание и конвейеризация.

На чём основан фокус

Идея fixed-point арифметики стара как мир: вместо дробного веса, например 0.135, храним целое число, сдвинутое влево на фиксированное количество бит. При сдвиге на 16 бит:

0.135 × 65536 = 8847

1	0.135 × 65536 = 8847

Теперь вместо Float * Float мы делаем Integer * Integer, а в конце сдвигаем результат обратно вправо на 16 бит (что эквивалентно делению на 65536, но выполняется за один такт). Точность при 16-битном сдвиге — порядка 0.0015%, для значений каналов 0…255 это даёт погрешность в тысячные доли, что совершенно незаметно.

Соответственно, всё, что нужно, это переписать формирование ядра, чтобы Weights стал массивом Integer, и заменить вещественное накопление сумм целочисленным. Посмотрим, как это выглядит.

Ядро должно быть целочисленным

Структура ядра меняется: вместо array of Double у нас будет array of Integer, а веса хранятся в формате fixed-point со сдвигом 16 бит:

const
  FIXED_SHIFT = 16;
  FIXED_ONE   = 1 shl FIXED_SHIFT; // 65536

type
  TKernel1DFixed = record
    Radius: Integer;
    Weights: array of Integer; // [0..Radius], фиксированная точка
  end;

const

FIXED_SHIFT = 16;

FIXED_ONE = 1 shl FIXED_SHIFT; // 65536

type

TKernel1DFixed = record

Radius: Integer;

Weights: array of Integer; // [0..Radius], фиксированная точка

end;

Формирование ядра чуть сложнее, чем в вещественном варианте, потому что появляется третий шаг, корректировка:

function MakeGaussianKernel1DFixed(Sigma: Double): TKernel1DFixed;
var
  I: Integer;
  Sum, W: Double;
  FloatWeights: array of Double;
  IntSum, Correction: Integer;
begin
  Result.Radius := Max(1, Ceil(Sigma * 3));
  SetLength(Result.Weights, Result.Radius + 1);
  SetLength(FloatWeights, Result.Radius + 1);

  // 1) Считаем float-веса и нормализуем
  Sum := 0;
  for I := 0 to Result.Radius do
  begin
    W := Exp(-(I * I) / (2 * Sigma * Sigma));
    FloatWeights[I] := W;
    if I = 0 then
      Sum := Sum + W
    else
      Sum := Sum + 2 * W;
  end;

  for I := 0 to Result.Radius do
    FloatWeights[I] := FloatWeights[I] / Sum;

  // 2) Переводим в fixed-point
  for I := 0 to Result.Radius do
    Result.Weights[I] := Round(FloatWeights[I] * FIXED_ONE);

  // 3) Корректируем, чтобы сумма = FIXED_ONE (65536)
  //    Иначе изображение чуть потемнеет/посветлеет
  IntSum := Result.Weights[0];
  for I := 1 to Result.Radius do
    IntSum := IntSum + 2 * Result.Weights[I];

  Correction := FIXED_ONE - IntSum;
  Result.Weights[0] := Result.Weights[0] + Correction;
end;

function MakeGaussianKernel1DFixed(Sigma: Double): TKernel1DFixed;

var

I: Integer;

Sum, W: Double;

FloatWeights: array of Double;

IntSum, Correction: Integer;

begin

Result.Radius := Max(1, Ceil(Sigma * 3));

SetLength(Result.Weights, Result.Radius + 1);

SetLength(FloatWeights, Result.Radius + 1);

// 1) Считаем float-веса и нормализуем

Sum := 0;

for I := 0 to Result.Radius do

begin

W := Exp(-(I * I) / (2 * Sigma * Sigma));

FloatWeights[I] := W;

if I = 0 then

Sum := Sum + W

else

Sum := Sum + 2 * W;

end;

for I := 0 to Result.Radius do

FloatWeights[I] := FloatWeights[I] / Sum;

// 2) Переводим в fixed-point

for I := 0 to Result.Radius do

Result.Weights[I] := Round(FloatWeights[I] * FIXED_ONE);

// 3) Корректируем, чтобы сумма = FIXED_ONE (65536)

// Иначе изображение чуть потемнеет/посветлеет

IntSum := Result.Weights[0];

for I := 1 to Result.Radius do

IntSum := IntSum + 2 * Result.Weights[I];

Correction := FIXED_ONE - IntSum;

Result.Weights[0] := Result.Weights[0] + Correction;

end;

Первые два шага очевидны: считаем гауссовы веса в Double, нормализуем, умножаем на 65536 и округляем. А вот третий шаг заслуживает пояснения.

При округлении каждого веса по отдельности накапливается ошибка: сумма целочисленных весов может оказаться 65534 или 65538 вместо ровных 65536. Казалось бы, мелочь — но при делении (сдвиге) на 65536 эта разница означает, что изображение систематически чуть темнеет или светлеет. Поэтому мы вычисляем разницу FIXED_ONE — IntSum и добавляем её к центральному весу. Он самый большой, так что коррекция в пару единиц на фоне значения примерно 14000 совершенно незаметна. После этого сумма всех весов точно равна 65536, и яркость изображения сохраняется.

Свёртка на целых числах

Структура процедуры один в один повторяет сепарабельную версию — два прохода, промежуточный буфер, тот же Clamp для краёв.

procedure BlurGaussianFixedPoint(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  K: TKernel1DFixed;
  W, H, X, Y, I: Integer;
  SX, SY: Integer;
  SumR, SumG, SumB: Integer;
  Weight: Integer;
  SrcCache, TmpCache, DstCache: TScanLineCache;
  Tmp: TBitmap;
  SrcPx: TPixel32;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  K := MakeGaussianKernel1DFixed(Sigma);

  Bitmap.PixelFormat := pf32bit;
  W := Bitmap.Width;
  H := Bitmap.Height;

  Tmp := TBitmap.Create;
  try
    Tmp.PixelFormat := pf32bit;
    Tmp.SetSize(W, H);

    Result := TBitmap.Create;
    Result.PixelFormat := pf32bit;
    Result.SetSize(W, H);

    SrcCache := BuildScanLineCache(Bitmap);
    TmpCache := BuildScanLineCache(Tmp);
    DstCache := BuildScanLineCache(Result);

    //  Проход 1: горизонтальный
    for Y := 0 to H - 1 do
    begin
      for X := 0 to W - 1 do
      begin
        SumR := 0; SumG := 0; SumB := 0;

        for I := -K.Radius to K.Radius do
        begin
          SX := Clamp(X + I, 0, W - 1);
          Weight := K.Weights[Abs(I)];

          SrcPx := SrcCache[Y][SX];
          SumB := SumB + SrcPx.B * Weight;
          SumG := SumG + SrcPx.G * Weight;
          SumR := SumR + SrcPx.R * Weight;
        end;

        with TmpCache[Y][X] do
        begin
          B := Clamp(SumB shr FIXED_SHIFT, 0, 255);
          G := Clamp(SumG shr FIXED_SHIFT, 0, 255);
          R := Clamp(SumR shr FIXED_SHIFT, 0, 255);
          A := SrcCache[Y][X].A;
        end;
      end;
    end;

    //  Проход 2: вертикальный
    for Y := 0 to H - 1 do
    begin
      for X := 0 to W - 1 do
      begin
        SumR := 0; SumG := 0; SumB := 0;

        for I := -K.Radius to K.Radius do
        begin
          SY := Clamp(Y + I, 0, H - 1);
          Weight := K.Weights[Abs(I)];

          SrcPx := TmpCache[SY][X];
          SumB := SumB + SrcPx.B * Weight;
          SumG := SumG + SrcPx.G * Weight;
          SumR := SumR + SrcPx.R * Weight;
        end;

        with DstCache[Y][X] do
        begin
          B := Clamp(SumB shr FIXED_SHIFT, 0, 255);
          G := Clamp(SumG shr FIXED_SHIFT, 0, 255);
          R := Clamp(SumR shr FIXED_SHIFT, 0, 255);
          A := TmpCache[Y][X].A;
        end;
      end;
    end;

  finally
    Tmp.Free;
  end;
end;

procedure BlurGaussianFixedPoint(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Sigma: Double;

K: TKernel1DFixed;

W, H, X, Y, I: Integer;

SX, SY: Integer;

SumR, SumG, SumB: Integer;

Weight: Integer;

SrcCache, TmpCache, DstCache: TScanLineCache;

Tmp: TBitmap;

SrcPx: TPixel32;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

K := MakeGaussianKernel1DFixed(Sigma);

Bitmap.PixelFormat := pf32bit;

W := Bitmap.Width;

H := Bitmap.Height;

Tmp := TBitmap.Create;

try

Tmp.PixelFormat := pf32bit;

Tmp.SetSize(W, H);

Result := TBitmap.Create;

Result.PixelFormat := pf32bit;

Result.SetSize(W, H);

SrcCache := BuildScanLineCache(Bitmap);

TmpCache := BuildScanLineCache(Tmp);

DstCache := BuildScanLineCache(Result);

// Проход 1: горизонтальный

for Y := 0 to H - 1 do

begin

for X := 0 to W - 1 do

begin

SumR := 0; SumG := 0; SumB := 0;

for I := -K.Radius to K.Radius do

begin

SX := Clamp(X + I, 0, W - 1);

Weight := K.Weights[Abs(I)];

SrcPx := SrcCache[Y][SX];

SumB := SumB + SrcPx.B * Weight;

SumG := SumG + SrcPx.G * Weight;

SumR := SumR + SrcPx.R * Weight;

end;

with TmpCache[Y][X] do

begin

B := Clamp(SumB shr FIXED_SHIFT, 0, 255);

G := Clamp(SumG shr FIXED_SHIFT, 0, 255);

R := Clamp(SumR shr FIXED_SHIFT, 0, 255);

A := SrcCache[Y][X].A;

end;

// Проход 2: вертикальный

for Y := 0 to H - 1 do

begin

for X := 0 to W - 1 do

begin

SumR := 0; SumG := 0; SumB := 0;

for I := -K.Radius to K.Radius do

begin

SY := Clamp(Y + I, 0, H - 1);

Weight := K.Weights[Abs(I)];

SrcPx := TmpCache[SY][X];

SumB := SumB + SrcPx.B * Weight;

SumG := SumG + SrcPx.G * Weight;

SumR := SumR + SrcPx.R * Weight;

end;

with DstCache[Y][X] do

begin

B := Clamp(SumB shr FIXED_SHIFT, 0, 255);

G := Clamp(SumG shr FIXED_SHIFT, 0, 255);

R := Clamp(SumR shr FIXED_SHIFT, 0, 255);

A := TmpCache[Y][X].A;

end;

finally

Tmp.Free;

end;

Если сравнить с предыдущей версией, изменений минимум:

SumR, SumG, SumB — теперь Integer, а не Double. Переполнения можно не бояться: максимальное значение канала 255, максимальная сумма весов 65536, итого 255 * 65536 = 16 711 680 — это далеко до границы Integer (2 147 483 647).
Round(SumB / TotalWeight) заменён на SumB shr FIXED_SHIFT — битовый сдвиг вправо на 16 вместо деления на 65536. Одна машинная инструкция вместо вещественного деления с округлением.
TotalWeight исчез — нормализация весов с коррекцией гарантирует, что сумма ровно 65536, делить больше не на что.

По сути, мы заменили самые дорогие операции в горячем цикле — умножение Double * Double и деление на умножение Integer × Integer и битовый сдвиг, не изменив ни структуру алгоритма, ни качество результата.

По сравнению с предыдущей реализацией, получили выигрыш примерно в 2 раза по скорости (при σ=5: 122 против 219 мс, при σ=50: 1192 против 2015 мс). Мы выжали из вещественной арифметики всё, перешли на целочисленную и получили ощутимый прирост. Но прежде чем двигаться дальше, стоит оглянуться. В Delphi-сообществе уже много лет ходит по рукам модуль GBlur2 — «проверенный и быстрый» Gaussian Blur. Его копируют из проекта в проект, на него ссылаются на форумах.

Поэтому давайте его тоже проанализируем, посмотрим, что у него под капотом и как он смотрится на фоне наших реализаций.

GBlur2: Легенда Delphi-форумов

Модуль GBlur2 кочует по Delphi-сообществам уже больше двадцати лет. Его копируют из проекта в проект, советуют на форумах как «быстрый Gaussian Blur». Заглянем внутрь.

Ключевая функция модуля: GBlur(Bitmap: TBitmap; Radius: Double). Обратите внимание: параметр называется Radius, хотя на самом деле это Sigma и настоящий радиус будет вычислен внутри. Да, этот параметр именно σ. Я разобрал алгоритм до винтика.

Инкапсулируем вызов процедуры GBlur2:

procedure BlurGBlur2(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var Sigma: Double;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  Result := TBitmap.Create;
  Result.Assign(Bitmap);
  Result.PixelFormat := pf24bit;
  // Вызов процедуры из GBlur2   
  GBlur(Result, Sigma);
end;

procedure BlurGBlur2(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var Sigma: Double;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

Result := TBitmap.Create;

Result.Assign(Bitmap);

Result.PixelFormat := pf24bit;

// Вызов процедуры из GBlur2

GBlur(Result, Sigma);

end;

Автор модуля весьма вольно интерпретировал создание ядра Gaussian Blur. Радиусом ядра мы считаем Sigma*3, в то время как здесь всё построено совершенно по другому. Blur2 реализует довольно нестандартный и «агрессивный» подход к определению рабочего радиуса. Вместо того, чтобы использовать фиксированное правило (вроде 3σ), автор модуля пытается динамически вычислить минимально необходимое окно, основываясь на точности представления цвета.

Как определяется радиус (K.Size)

Автор использует критерий «заметности» ошибки. Ключевая логика здесь:

delta := DataGranularity / (2 * MaxData); // Обычно 1 / (2 * 255) = 0.00196
...
while (temp < delta) and (KernelSize > 1) do
begin
  temp := temp + 2 * K.Weights[KernelSize];
  dec(KernelSize);
end;

delta := DataGranularity / (2 * MaxData); // Обычно 1 / (2 * 255) = 0.00196

...

while (temp < delta) and (KernelSize > 1) do

begin

temp := temp + 2 * K.Weights[KernelSize];

dec(KernelSize);

end;

Что это значит на практике:

Алгоритм начинает с максимально возможного размера ядра (MaxKernelSize).
Он суммирует веса «хвостов» распределения Гаусса (самые края ядра).
Как только суммарный вес отбрасываемых краев превышает порог delta (около 0.2% общей энергии фильтра), цикл останавливается.
Итог: K.Size — это минимальное расстояние от центра, при котором оставшаяся «масса» функции Гаусса еще оказывает влияние на итоговый цвет пикселя.

Почему это «вольная интерпретация»?

В классической математике мы сначала выбираем радиус, а потом строим под него ядро. Здесь же:

Сначала считается огромное ядро (на весь размер массива K.Weights).
Заполняется значениями исходя из сигмы σ (которая в коде названа radius).
А затем ядро обрезается (вычисляется K.Size), исходя из того, что некоторые значения стали слишком малы, чтобы хоть как-то влиять на байт цвета.

Выводы

«Настоящий» радиус в этом модуле — это динамическая величина.

Входной параметр radius — это строго σ(сигма).
Выходной K.Size — это эффективный радиус, рассчитанный так, чтобы не выполнять лишних умножений на почти нулевые коэффициенты.

Что ещё бросается в глаза:

pf24bit — модуль работает с 24-битными пикселями (3 байта, без альфа-канала). Наши реализации используют pf32bit, что удобнее для выравнивания и потенциально быстрее на современных процессорах.
Работает in-place — результат записывается прямо в переданный битмап. Поэтому мы делаем Result.Assign(Bitmap) перед вызовом.
Внутри — та же сепарабельная свёртка с вещественными весами. Никаких секретных оптимизаций: Double-накопители, Round, классический двухпроходный подход.

Итак, мы выяснили, что GBlur2 — это обычная сепарабельная свёртка с необычным (агрессивным) расчётом вещественного ядра. Модуль многим пригодился (и автору статьи тоже), критиковать его тут никто не собирается. Но всё же есть нюанс, о котором стоит напомнить:

const
  MaxKernelSize = 100;

1 2	const MaxKernelSize = 100;

Размер ядра жёстко ограничен константой. При большой σ (сигме) радиус упирается в этот потолок, и ядро обрезается — размытие перестаёт нарастать. Если вы передадите очень большую сигму, но MaxKernelSize в модуле ограничен (например, 50 или 100), то хвосты функции не успеют «затухнуть» до delta, и вы получите обрезание ядра по лимиту массива, что может привести к артефактам на краях.

Для бенчмарка мы расширим это ограничение до 300, чтобы замеры были честными, но в оригинальном модуле эта константа выставлена именно так.

Бенчмарк покажет, где именно GBlur2 оказывается в общем рейтинге. Забегая вперёд, скажу — на последнем месте, предыдущий fixed-point метод быстрее. А мы пока вернёмся к нашим оптимизациям, Gaussian ещё не сказал последнего слова. Окончательный оптимизированный Gaussian Blur — в следующем разделе.

Gaussian: Оптимизация, выжимаем максимум

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

Clamp для координат — на каждом шаге ядра мы проверяем, не вышел ли индекс за границу. Для пикселей в середине изображения (а это подавляющее большинство) проверка всегда проходит впустую — сосед гарантированно существует. Мы платим за граничные случаи, которые составляют доли процента.
Abs(I) для индекса ядра — ядро симметрично, и мы храним только правую половину. Значит, при отрицательных смещениях каждый раз вычисляем модуль. Мелочь, но в цикле, выполняющемся десятки миллионов раз, мелочи складываются.
Вертикальный проход бьёт по кэшу — при обходе столбца каждый следующий пиксель лежит на другой строке. Строки битмапа могут быть далеко друг от друга в памяти, и процессор постоянно промахивается мимо кэша. Горизонтальный проход читает память последовательно и работает заметно быстрее.
Два отдельных цикла по X и Y — структурно чистое решение, но между проходами мы полностью заканчиваем один и начинаем другой. Можно ли организовать данные так, чтобы вертикальный проход работал столь же дружелюбно к кэшу, как горизонтальный?

Все эти проблемы решаемы. Разобьём цикл на внутреннюю часть (где границы заведомо не нарушены) и краевые полосы (где нужен Clamp). Перестроим цикл так, чтобы индекс всегда был положительным, избавившись от Abs. А для вертикального прохода введём промежуточный буфер, организованный по столбцам. Посмотрим, что из этого получится.

Суть оптимизации

Оптимизация строится на трёх идеях: разделить цикл на края и середину, использовать симметрию ядра, а проблему кэша решить транспонированием. Начнём со вспомогательных процедур, затем — основная функция.

Комментарий к коду: Ядро формируется той же логикой, что и MakeGaussianKernel1DFixed. В коде демки оно продублировано как TFastKernel / BuildFastKernel — для удобства последующего переноса в итоговый модуль.

Транспонирование

Простая, но ключевая операция. Вместо того чтобы идти по столбцам (с постоянными кэш-промахами), мы переворачиваем изображение и снова идём по строкам:

// Транспонирование битмапа через ScanLine-кэш
procedure TransposeBitmap(Src: TBitmap; out Dst: TBitmap);
var
  X, Y, W, H: Integer;
  SrcCache, DstCache: TScanLineCache;
begin
  W := Src.Width;
  H := Src.Height;

  Dst := TBitmap.Create;
  Dst.PixelFormat := pf32bit;
  Dst.SetSize(H, W);  // ширина и высота меняются местами

  SrcCache := BuildScanLineCache(Src);
  DstCache := BuildScanLineCache(Dst);

  for Y := 0 to H - 1 do
    for X := 0 to W - 1 do
      DstCache[X][Y] := SrcCache[Y][X];
end;

// Транспонирование битмапа через ScanLine-кэш

procedure TransposeBitmap(Src: TBitmap; out Dst: TBitmap);

var

X, Y, W, H: Integer;

SrcCache, DstCache: TScanLineCache;

begin

W := Src.Width;

H := Src.Height;

Dst := TBitmap.Create;

Dst.PixelFormat := pf32bit;

Dst.SetSize(H, W); // ширина и высота меняются местами

SrcCache := BuildScanLineCache(Src);

DstCache := BuildScanLineCache(Dst);

for Y := 0 to H - 1 do

for X := 0 to W - 1 do

DstCache[X][Y] := SrcCache[Y][X];

end;

Горизонтальный проход

Это сердце оптимизации. Цикл разбит на три зоны:

// Горизонтальный проход с разделением краёв и середины
procedure HorzPassFast(SrcCache, DstCache: TScanLineCache;
  W, H: Integer; const K: TFastKernel);
var
  X, Y, I, R: Integer;
  SumR, SumG, SumB: Integer;
  Weight: Integer;
  SrcRow: PPixel32Array;
  DstRow: PPixel32Array;
begin
  R := K.Radius;

  for Y := 0 to H - 1 do
  begin
    SrcRow := PPixel32Array(SrcCache[Y]);
    DstRow := PPixel32Array(DstCache[Y]);

    //  Левый край: X = 0..R-1 (нужен Clamp)
    for X := 0 to R - 1 do
    begin
      // Центральный вес
      Weight := K.Weights[0];
      SumR := SrcRow[X].R * Weight;
      SumG := SrcRow[X].G * Weight;
      SumB := SrcRow[X].B * Weight;

      for I := 1 to R do
      begin
        Weight := K.Weights[I];

        // Левый сосед (может выйти за границу)
        if X - I >= 0 then
        begin
          SumR := SumR + SrcRow[X - I].R * Weight;
          SumG := SumG + SrcRow[X - I].G * Weight;
          SumB := SumB + SrcRow[X - I].B * Weight;
        end
        else
        begin
          SumR := SumR + SrcRow[0].R * Weight;
          SumG := SumG + SrcRow[0].G * Weight;
          SumB := SumB + SrcRow[0].B * Weight;
        end;

        // Правый сосед (всегда в диапазоне при X < R и R < W)
        if X + I < W then
        begin
          SumR := SumR + SrcRow[X + I].R * Weight;
          SumG := SumG + SrcRow[X + I].G * Weight;
          SumB := SumB + SrcRow[X + I].B * Weight;
        end
        else
        begin
          SumR := SumR + SrcRow[W - 1].R * Weight;
          SumG := SumG + SrcRow[W - 1].G * Weight;
          SumB := SumB + SrcRow[W - 1].B * Weight;
        end;
      end;

      DstRow[X].B := Clamp(SumB shr FIXED_SHIFT, 0, 255);
      DstRow[X].G := Clamp(SumG shr FIXED_SHIFT, 0, 255);
      DstRow[X].R := Clamp(SumR shr FIXED_SHIFT, 0, 255);
      DstRow[X].A := SrcRow[X].A;
    end;

    //  Середина: X = R..W-R-1 (без проверок границ!)
    for X := R to W - R - 1 do
    begin
      Weight := K.Weights[0];
      SumR := SrcRow[X].R * Weight;
      SumG := SrcRow[X].G * Weight;
      SumB := SrcRow[X].B * Weight;

      for I := 1 to R do
      begin
        Weight := K.Weights[I];
        // Гарантированно в диапазоне - никаких проверок!
        SumR := SumR + (SrcRow[X - I].R + SrcRow[X + I].R) * Weight;
        SumG := SumG + (SrcRow[X - I].G + SrcRow[X + I].G) * Weight;
        SumB := SumB + (SrcRow[X - I].B + SrcRow[X + I].B) * Weight;
      end;

      DstRow[X].B := SumB shr FIXED_SHIFT;
      DstRow[X].G := SumG shr FIXED_SHIFT;
      DstRow[X].R := SumR shr FIXED_SHIFT;
      DstRow[X].A := SrcRow[X].A;
    end;

    //  Правый край: X = W-R..W-1 (нужен Clamp)
    for X := Max(R, W - R) to W - 1 do
    begin
      Weight := K.Weights[0];
      SumR := SrcRow[X].R * Weight;
      SumG := SrcRow[X].G * Weight;
      SumB := SrcRow[X].B * Weight;

      for I := 1 to R do
      begin
        Weight := K.Weights[I];

        if X - I >= 0 then
        begin
          SumR := SumR + SrcRow[X - I].R * Weight;
          SumG := SumG + SrcRow[X - I].G * Weight;
          SumB := SumB + SrcRow[X - I].B * Weight;
        end
        else
        begin
          SumR := SumR + SrcRow[0].R * Weight;
          SumG := SumG + SrcRow[0].G * Weight;
          SumB := SumB + SrcRow[0].B * Weight;
        end;

        if X + I < W then
        begin
          SumR := SumR + SrcRow[X + I].R * Weight;
          SumG := SumG + SrcRow[X + I].G * Weight;
          SumB := SumB + SrcRow[X + I].B * Weight;
        end
        else
        begin
          SumR := SumR + SrcRow[W - 1].R * Weight;
          SumG := SumG + SrcRow[W - 1].G * Weight;
          SumB := SumB + SrcRow[W - 1].B * Weight;
        end;
      end;

      DstRow[X].B := Clamp(SumB shr FIXED_SHIFT, 0, 255);
      DstRow[X].G := Clamp(SumG shr FIXED_SHIFT, 0, 255);
      DstRow[X].R := Clamp(SumR shr FIXED_SHIFT, 0, 255);
      DstRow[X].A := SrcRow[X].A;
    end;
  end;
end;

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

// Горизонтальный проход с разделением краёв и середины

procedure HorzPassFast(SrcCache, DstCache: TScanLineCache;

W, H: Integer; const K: TFastKernel);

var

X, Y, I, R: Integer;

SumR, SumG, SumB: Integer;

Weight: Integer;

SrcRow: PPixel32Array;

DstRow: PPixel32Array;

begin

R := K.Radius;

for Y := 0 to H - 1 do

begin

SrcRow := PPixel32Array(SrcCache[Y]);

DstRow := PPixel32Array(DstCache[Y]);

// Левый край: X = 0..R-1 (нужен Clamp)

for X := 0 to R - 1 do

begin

// Центральный вес

Weight := K.Weights[0];

SumR := SrcRow[X].R * Weight;

SumG := SrcRow[X].G * Weight;

SumB := SrcRow[X].B * Weight;

for I := 1 to R do

begin

Weight := K.Weights[I];

// Левый сосед (может выйти за границу)

if X - I >= 0 then

begin

SumR := SumR + SrcRow[X - I].R * Weight;

SumG := SumG + SrcRow[X - I].G * Weight;

SumB := SumB + SrcRow[X - I].B * Weight;

end

else

begin

SumR := SumR + SrcRow[0].R * Weight;

SumG := SumG + SrcRow[0].G * Weight;

SumB := SumB + SrcRow[0].B * Weight;

end;

// Правый сосед (всегда в диапазоне при X < R и R < W)

if X + I < W then

begin

SumR := SumR + SrcRow[X + I].R * Weight;

SumG := SumG + SrcRow[X + I].G * Weight;

SumB := SumB + SrcRow[X + I].B * Weight;

end

else

begin

SumR := SumR + SrcRow[W - 1].R * Weight;

SumG := SumG + SrcRow[W - 1].G * Weight;

SumB := SumB + SrcRow[W - 1].B * Weight;

end;

DstRow[X].B := Clamp(SumB shr FIXED_SHIFT, 0, 255);

DstRow[X].G := Clamp(SumG shr FIXED_SHIFT, 0, 255);

DstRow[X].R := Clamp(SumR shr FIXED_SHIFT, 0, 255);

DstRow[X].A := SrcRow[X].A;

end;

// Середина: X = R..W-R-1 (без проверок границ!)

for X := R to W - R - 1 do

begin

Weight := K.Weights[0];

SumR := SrcRow[X].R * Weight;

SumG := SrcRow[X].G * Weight;

SumB := SrcRow[X].B * Weight;

for I := 1 to R do

begin

Weight := K.Weights[I];

// Гарантированно в диапазоне - никаких проверок!

SumR := SumR + (SrcRow[X - I].R + SrcRow[X + I].R) * Weight;

SumG := SumG + (SrcRow[X - I].G + SrcRow[X + I].G) * Weight;

SumB := SumB + (SrcRow[X - I].B + SrcRow[X + I].B) * Weight;

end;

DstRow[X].B := SumB shr FIXED_SHIFT;

DstRow[X].G := SumG shr FIXED_SHIFT;

DstRow[X].R := SumR shr FIXED_SHIFT;

DstRow[X].A := SrcRow[X].A;

end;

// Правый край: X = W-R..W-1 (нужен Clamp)

for X := Max(R, W - R) to W - 1 do

begin

Weight := K.Weights[0];

SumR := SrcRow[X].R * Weight;

SumG := SrcRow[X].G * Weight;

SumB := SrcRow[X].B * Weight;

for I := 1 to R do

begin

Weight := K.Weights[I];

if X - I >= 0 then

begin

SumR := SumR + SrcRow[X - I].R * Weight;

SumG := SumG + SrcRow[X - I].G * Weight;

SumB := SumB + SrcRow[X - I].B * Weight;

end

else

begin

SumR := SumR + SrcRow[0].R * Weight;

SumG := SumG + SrcRow[0].G * Weight;

SumB := SumB + SrcRow[0].B * Weight;

end;

if X + I < W then

begin

SumR := SumR + SrcRow[X + I].R * Weight;

SumG := SumG + SrcRow[X + I].G * Weight;

SumB := SumB + SrcRow[X + I].B * Weight;

end

else

begin

SumR := SumR + SrcRow[W - 1].R * Weight;

SumG := SumG + SrcRow[W - 1].G * Weight;

SumB := SumB + SrcRow[W - 1].B * Weight;

end;

DstRow[X].B := Clamp(SumB shr FIXED_SHIFT, 0, 255);

DstRow[X].G := Clamp(SumG shr FIXED_SHIFT, 0, 255);

DstRow[X].R := Clamp(SumR shr FIXED_SHIFT, 0, 255);

DstRow[X].A := SrcRow[X].A;

end;

Обратите внимание на середину — самую горячую зону. Здесь нет ни одного if, ни одного Abs, ни одного Clamp. Симметричные соседи складываются до умножения на вес: (SrcRow[X-I].R + SrcRow[X+I].R) * Weight — одно умножение вместо двух. А сумма весов ровно 65536, поэтому shr FIXED_SHIFT гарантированно даёт 0…255 без проверки.

Левый и правый края обрабатываются отдельно, с проверками — но это лишь 2R пикселей из каждой строки. При типичном изображении 1920 px и радиусе 30 — это 3% от строки.

Основная процедура

И наконец, основная процедура, которая собирает всё вместе:

procedure BlurGaussianSeparableFast(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  K: TFastKernel;
  W, H: Integer;
  Tmp, Transposed, BlurredT: TBitmap;
  SrcCache, TmpCache, TransCache, BlurTCache: TScanLineCache;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  BuildFastKernel(K, Sigma);

  Bitmap.PixelFormat := pf32bit;
  W := Bitmap.Width;
  H := Bitmap.Height;

  // --- Шаг 1: горизонтальный проход (Src -> Tmp) ---
  Tmp := TBitmap.Create;
  try
    Tmp.PixelFormat := pf32bit;
    Tmp.SetSize(W, H);

    SrcCache := BuildScanLineCache(Bitmap);
    TmpCache := BuildScanLineCache(Tmp);

    HorzPassFast(SrcCache, TmpCache, W, H, K);

    // --- Шаг 2: транспонировать ---
    TransposeBitmap(Tmp, Transposed);
  finally
    Tmp.Free;
  end;

  // --- Шаг 3: горизонтальный проход по транспонированному (= вертикальный) ---
  try
    BlurredT := TBitmap.Create;
    try
      BlurredT.PixelFormat := pf32bit;
      BlurredT.SetSize(H, W);  // размеры транспонированы

      TransCache := BuildScanLineCache(Transposed);
      BlurTCache := BuildScanLineCache(BlurredT);

      HorzPassFast(TransCache, BlurTCache, H, W, K);

      // --- Шаг 4: транспонировать обратно ---
      TransposeBitmap(BlurredT, Result);
    finally
      BlurredT.Free;
    end;
  finally
    Transposed.Free;
  end;
end;

procedure BlurGaussianSeparableFast(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Sigma: Double;

K: TFastKernel;

W, H: Integer;

Tmp, Transposed, BlurredT: TBitmap;

SrcCache, TmpCache, TransCache, BlurTCache: TScanLineCache;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

BuildFastKernel(K, Sigma);

Bitmap.PixelFormat := pf32bit;

W := Bitmap.Width;

H := Bitmap.Height;

// --- Шаг 1: горизонтальный проход (Src -> Tmp) ---

Tmp := TBitmap.Create;

try

Tmp.PixelFormat := pf32bit;

Tmp.SetSize(W, H);

SrcCache := BuildScanLineCache(Bitmap);

TmpCache := BuildScanLineCache(Tmp);

HorzPassFast(SrcCache, TmpCache, W, H, K);

// --- Шаг 2: транспонировать ---

TransposeBitmap(Tmp, Transposed);

finally

Tmp.Free;

end;

// --- Шаг 3: горизонтальный проход по транспонированному (= вертикальный) ---

try

BlurredT := TBitmap.Create;

try

BlurredT.PixelFormat := pf32bit;

BlurredT.SetSize(H, W); // размеры транспонированы

TransCache := BuildScanLineCache(Transposed);

BlurTCache := BuildScanLineCache(BlurredT);

HorzPassFast(TransCache, BlurTCache, H, W, K);

// --- Шаг 4: транспонировать обратно ---

TransposeBitmap(BlurredT, Result);

finally

BlurredT.Free;

end;

finally

Transposed.Free;

end;

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

Скрин не привожу — он один-в один с предыдущими. Время покажет бенчмарк.

Системные средства: GDI+ и Direct2D

Мы написали четыре реализации Gaussian blur, каждый раз ускоряя предыдущую. Но прежде чем смотреть бенчмарк — стоит задать вопрос: а что предлагает сама Windows?

В системе есть два готовых механизма размытия. Оба умеют делать Gaussian blur «из коробки», оба не требуют ни строчки математики от разработчика. Но у каждого есть свой характер, свои ограничения и свои сюрпризы.

GDI+ Blur

GDI+ версии 1.1 (Vista+) предоставляет функцию GdipBitmapApplyEffect с эффектом BlurEffect. API принимает параметр radius (судя по документации в диапазоне 0–255), но его связь с σ нигде не документирована. Эмпирически подобранный коэффициент σ × 2.1 даёт визуально похожий результат — но только при σ ≥ 10. При малых значениях размытие заметно отличается от эталонного Gaussian: ядро ведёт себя непредсказуемо. Дополнительная ложка дёгтя — без флага ExpandEdge изображение обрезается по краям.

За эффектом размытия по Гауссу в GDI+ сходим в телегу, там много всего такого, чего нет на сайте.

const
  BlurEffectGuid: TGUID = '{633C80A4-1843-482b-9EF2-BE2834C5FDD4}';

type
  GpEffect = Pointer;
  PGpBitmap = ^GpBitmap;

  TBlurParams = record
  public
    radius: Single; // 0..255
    expandEdge: BOOL; // TRUE
  end;
  PBlurParams = ^TBlurParams;
  TAccessGPBitmap = class(TGPBitmap);

function GdipCreateEffect(Guid: TGUID;
  out Effect: GpEffect): TStatus; stdcall; external WINGDIPDLL;
function GdipDeleteEffect(Effect: GpEffect): TStatus;
  stdcall; external WINGDIPDLL;
function GdipSetEffectParameters(Effect: GpEffect; const Params: Pointer;
  const Size: UINT): TStatus; stdcall; external WINGDIPDLL;
function GdipBitmapApplyEffect(Bitmap: GpBitmap; Effect: GpEffect; 
  roi: PRect; UseAuxData: Bool; out AuxData: Pointer; 
  out AuxDataSize: Integer): TStatus; stdcall; external WINGDIPDLL;

function SetBlurEffect(var ABitmap: TGPBitmap; ABlur: Single): Boolean;
var
  Effect: GpEffect;
  Params: TBlurParams;
  AuxData: Pointer;
  AuxDataSize: Integer;
begin
  Result := GdipCreateEffect(BlurEffectGuid, Effect) = Ok;
  if not Result then exit;
  try
    Params.radius := ABlur;
    Params.expandEdge := True;
    Result := GdipSetEffectParameters(Effect, @Params, SizeOf(Params)) = Ok;
    if not Result then exit;
    Result := GdipBitmapApplyEffect(TAccessGPBitmap(ABitmap).nativeImage,
      Effect, nil, False, AuxData, AuxDataSize) = Ok;
  finally
    GdipDeleteEffect(Effect);
  end;
end;

procedure BlurGaussianGDIP(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  Tmp : TGPBitmap;
  h: HBitmap;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;
  Bitmap.PixelFormat := pf32bit;

  Result := TBitmap.Create;
  Result.Assign(Bitmap);
  Result.PixelFormat := pf32bit;

  Tmp := TGPBitmap.Create(Result.Handle, Result.Palette);
  try
    if not SetBlurEffect(Tmp, Sigma*2.1) then
      exit;
    Tmp.GetHBITMAP(MakeColor(0, 0, 0, 0), H);
    Result.Handle := H;
  finally
    FreeAndNil(Tmp);
  end;
end;

const

BlurEffectGuid: TGUID = '{633C80A4-1843-482b-9EF2-BE2834C5FDD4}';

type

GpEffect = Pointer;

PGpBitmap = ^GpBitmap;

TBlurParams = record

public

radius: Single; // 0..255

expandEdge: BOOL; // TRUE

end;

PBlurParams = ^TBlurParams;

TAccessGPBitmap = class(TGPBitmap);

function GdipCreateEffect(Guid: TGUID;

out Effect: GpEffect): TStatus; stdcall; external WINGDIPDLL;

function GdipDeleteEffect(Effect: GpEffect): TStatus;

stdcall; external WINGDIPDLL;

function GdipSetEffectParameters(Effect: GpEffect; const Params: Pointer;

const Size: UINT): TStatus; stdcall; external WINGDIPDLL;

function GdipBitmapApplyEffect(Bitmap: GpBitmap; Effect: GpEffect;

roi: PRect; UseAuxData: Bool; out AuxData: Pointer;

out AuxDataSize: Integer): TStatus; stdcall; external WINGDIPDLL;

function SetBlurEffect(var ABitmap: TGPBitmap; ABlur: Single): Boolean;

var

Effect: GpEffect;

Params: TBlurParams;

AuxData: Pointer;

AuxDataSize: Integer;

begin

Result := GdipCreateEffect(BlurEffectGuid, Effect) = Ok;

if not Result then exit;

try

Params.radius := ABlur;

Params.expandEdge := True;

Result := GdipSetEffectParameters(Effect, @Params, SizeOf(Params)) = Ok;

if not Result then exit;

Result := GdipBitmapApplyEffect(TAccessGPBitmap(ABitmap).nativeImage,

Effect, nil, False, AuxData, AuxDataSize) = Ok;

finally

GdipDeleteEffect(Effect);

end;

procedure BlurGaussianGDIP(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Sigma: Double;

Tmp : TGPBitmap;

h: HBitmap;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

Bitmap.PixelFormat := pf32bit;

Result := TBitmap.Create;

Result.Assign(Bitmap);

Result.PixelFormat := pf32bit;

Tmp := TGPBitmap.Create(Result.Handle, Result.Palette);

try

if not SetBlurEffect(Tmp, Sigma*2.1) then

exit;

Tmp.GetHBITMAP(MakeColor(0, 0, 0, 0), H);

Result.Handle := H;

finally

FreeAndNil(Tmp);

end;

Время выполнения безусловно очень быстрое. Но размытие до сигмы ≤ 10 какое-то непонятное, смазанное. И чем больше сигма, тем быстрее выполняется. Скорее всего внутри GDI+ переключается на другой алгоритм или ядро вырождается. Существующий лимит 255 на radius не позволит задать сигму больше σ = 121.

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

Direct2D Gaussian Blur

Direct2D предлагает встроенный эффект CLSID_D2D1GaussianBlur, который выполняется на GPU. Параметр — честная сигма σ (standard deviation), без магических коэффициентов. Результат визуально идентичен эталонному Gaussian при любых значениях σ.

За исходниками снова сходим в телегу. Там выложен эффект стекла, в котором используется D2D-блюр.

procedure BlurGaussianD2D(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Item: TBlurD2DValue;
  Sigma: Double;
  Canvas: TDirect2DCanvas;
  Context: ID2D1DeviceContext;
  D2D1Bitmap: ID2D1Bitmap;
  Effect: ID2D1Effect;
begin
  Item := Value.AsType<TBlurD2DValue>;
  Sigma := Item.FSigma;
  Canvas := Item.FCanvas;

  Result := TBitmap.Create;
  try
    // Эффект работает с 32-битным изображением
    Bitmap.PixelFormat := pf32bit;
    // Просто переносим исходник в результат
    Result.Assign(Bitmap);
    Result.PixelFormat := pf32bit;
    // Указываем для отрисовки холст битмапа
    (Canvas.RenderTarget as ID2D1DCRenderTarget).BindDC(
      Result.Canvas.Handle, Result.Canvas.ClipRect);
    // Получаем контекст для хитрых отрисовок
    Canvas.RenderTarget.QueryInterface(ID2D1DeviceContext, Context);
    // Получаем битмап
    D2D1Bitmap := Canvas.CreateBitmap(Bitmap);
    // Стартуем сцену
    Canvas.BeginDraw;
    try
      // Применяем эффект (зря писали что ли в своё время)
      if GetGaussianBlurEffect(Context,
           D2D1Bitmap as ID2D1Image, Sigma,
           D2D1_GAUSSIANBLUR_OPTIMIZATION_BALANCED,
           D2D1_BORDER_MODE_HARD,
           Effect)
      then
        // Всё удачно - рисуем эффект на контексте
        Context.DrawImage(Effect as ID2D1Image);
    finally
      // Завершаем сцену
      Canvas.EndDraw;
    end;
  except
    FreeAndNil(Result);
  end;
end;

procedure BlurGaussianD2D(Bitmap: TBitmap; const Value: TValue;

out Result: TBitmap);

var

Item: TBlurD2DValue;

Sigma: Double;

Canvas: TDirect2DCanvas;

Context: ID2D1DeviceContext;

D2D1Bitmap: ID2D1Bitmap;

Effect: ID2D1Effect;

begin

Item := Value.AsType<TBlurD2DValue>;

Sigma := Item.FSigma;

Canvas := Item.FCanvas;

Result := TBitmap.Create;

try

// Эффект работает с 32-битным изображением

Bitmap.PixelFormat := pf32bit;

// Просто переносим исходник в результат

Result.Assign(Bitmap);

Result.PixelFormat := pf32bit;

// Указываем для отрисовки холст битмапа

(Canvas.RenderTarget as ID2D1DCRenderTarget).BindDC(

Result.Canvas.Handle, Result.Canvas.ClipRect);

// Получаем контекст для хитрых отрисовок

Canvas.RenderTarget.QueryInterface(ID2D1DeviceContext, Context);

// Получаем битмап

D2D1Bitmap := Canvas.CreateBitmap(Bitmap);

// Стартуем сцену

Canvas.BeginDraw;

try

// Применяем эффект (зря писали что ли в своё время)

if GetGaussianBlurEffect(Context,

D2D1Bitmap as ID2D1Image, Sigma,

D2D1_GAUSSIANBLUR_OPTIMIZATION_BALANCED,

D2D1_BORDER_MODE_HARD,

Effect)

then

// Всё удачно - рисуем эффект на контексте

Context.DrawImage(Effect as ID2D1Image);

finally

// Завершаем сцену

Canvas.EndDraw;

end;

except

FreeAndNil(Result);

end;

Как видим, этому монстру глубоко наплевать на значение сигмы. Аппаратное ускорение, графический процессор, параллельные вычисления. Всё летает, но привязано к железу и капризам ОС.

Цена — инициализация: D3D11 — DXGI — D2D Device — DeviceContext. Это десятки строк кода и обязательная привязка к COM. Однако если инициализация выполнена заранее, само размытие занимает единицы миллисекунд — и время не зависит от σ.

Зачем включать в сравнение

Оба метода участвуют в бенчмарке не для того, чтобы «выиграть», а чтобы показать шкалу. Наши CPU-реализации дают полный контроль над алгоритмом, работают везде и предсказуемо. Системные средства показывают, что вообще возможно в принципе, и за какую цену.

Бенчмарк

Мы не будем измерять два первых метода — наивный и сепарабельный. Они заведомо очень медленные и наш бенчмарк просто умрёт в ожидании конца. Да и на фоне их огромных значений не будут различимы детали остальных. Доступных для бенчмарка реализаций набралось пять: от Fixed-point Gaussian до GPU-ускоренного Direct2D.

Методика

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

Множество итераций. Каждый алгоритм прогоняется не один раз, а столько, сколько позволяет отведённое время. Результат — среднее по замерам.

Одинаковые условия. Все алгоритмы получают одно и то же изображение, одну и ту же сигму. Для GBlur2 константа MaxKernelSize расширена до 300, чтобы при больших σ ядро не обрезалось и замер был честным.

Возможность прерывания. Полный прогон на больших σ может занять минуты. Бенчмарк можно остановить в любой момент — результаты уже завершённых замеров сохраняются и отображаются на диаграмме. Прерванные алгоритмы помечаются отдельно, так что по графику всегда видно, какие данные полные, а какие — частичные.

Диаграммы

30 замеров на метод, изображение 600×600, два режима: σ = 5 (лёгкое размытие) и σ = 50 (сильное размытие). Время в миллисекундах, меньше — лучше.

32-bit

64-bit

Разбор результатов

Диаграмма показывает три вещи. Во-первых, каждая наша оптимизация давала измеримый результат: от GBlur2 до Gauss Fast — ×2.5 в 32-bit. Во-вторых, все CPU-реализации остаются O(r): десятикратный рост сигмы σ даёт десятикратный рост времени. В-третьих, системные средства живут по другим правилам: GDI+ при σ=50 быстрее, чем при σ=5 (вероятно, переключается на другой алгоритм), а Direct2D вообще не зависит от радиуса и по времени выполнения живёт где-то в другой вселенной.

Интереснее всего — разница между 32-bit и 64-bit:

32-bit vs 64-bit: кому помогает разрядность?

Самая интересная история начинается при сравнении двух диаграмм. Сведём в таблицу:

Метод	σ	32-bit	64-bit	Ускорение
Gauss Fixed	5	125.8	110.2	×1.1
	50	1209.2	1052.8	×1.1
GBlur2	5	176.8	78.7	×2.2
	50	1592.5	773.8	×2.1
Gauss Fast	5	71.4	44.1	×1.6
	50	719.0	329.4	×2.2
Gauss GDIP	5	38.9	40.8	×1.0
	50	13.4	12.4	×1.1
Gauss D2D	5	7.3	10.7	×0.7
	50	7.3	16.8	×0.4

GBlur2 получает двукратное ускорение и это целиком заслуга перехода с x87 FPU на SSE2. В 32-bit Delphi вещественная арифметика использует стековую архитектуру x87: каждая операция с Double сопровождается неявными зависимостями через стек регистров и невозможностью для процессора переупорядочивать инструкции. В 64-bit компилятор генерирует SSE2-инструкции с прямой адресацией регистров XMM0–XMM15, конвейер процессора раскрывается полностью. Мы это случайно открыли, когда делали кроссплатформенный MulDiv. Для алгоритма, который делает тысячи умножений Double × Double на каждую строку, это колоссальная разница.

Gauss Fixed ускоряется всего на 10–13%. Целочисленная арифметика и в 32-bit работает хорошо, shl/shr/imul не страдают от стековой архитектуры. Прирост даёт удвоенное количество регистров общего назначения (16 вместо 8): меньше промежуточных переменных вытесняется на стек.

Gauss Fast — от ×1.6 до ×2.2. Казалось бы, fixed-point, значит целые числа, значит выигрыш должен быть как у Gauss Fixed (~10%). Но Gauss Fast активно использует множество локальных переменных в горячем цикле: SumR, SumG, SumB, указатели на строки, веса ядра. В 32-bit с его восемью регистрами часть из них неизбежно живёт на стеке. В 64-bit всё помещается в регистры — и это даёт ощутимый прирост, особенно при больших ядрах (σ = 50), где горячий цикл длиннее.

GDI+ не меняется — работу делает системная DLL, разрядность приложения роли не играет.

Direct2D в 64-bit медленнее. GPU-часть идентична — видеокарта не знает о разрядности процесса. Но overhead копирования текстур между CPU и GPU, COM-вызовы и синхронизация — всё это проходит через другой путь в 64-bit (другой calling convention, другие размеры указателей в COM-структурах). Это мои предположения, я не знаю, как там реально обстоят дела.

Что из этого следует

Наши CPU-реализации дали ×2.5 ускорение от GBlur2 до Gauss Fast в 32-bit и ×1.8 в 64-bit. Это серьёзный результат. Переход на 64-bit добавляет ещё ×1.6–2.2 сверху бесплатно. Но все они остаются O(r), время линейно растёт с радиусом, и при σ = 50 даже лучший вариант в 64-bit тратит 329 мс. Системные средства показывают, что возможно принципиально иное: GDI+ каким-то образом укладывается в 12 мс, D2D — в 7–17 мс.

Если хотим полный контроль над алгоритмом, предсказуемое поведение и кроссплатформенность, то наш код работает везде, где есть TBitmap и ScanLine. Хотим скорость любой ценой, выбираем Direct2D и GPU.

Но есть и третий путь. Можно ли на CPU добиться O(1) по радиусу, сохранив качество Gaussian? Можно ли обойтись без GPU и чёрных ящиков, но при этом не платить за каждый пиксель ядра?

Об этом — во второй части: Box blur, Stack blur, тройное размытие и трюк с масштабированием. А если получится продолжить серию дальше, то есть много тем для рассмотрения: многопоточность, корректная обработка альфа-канала и рекурсивные (IIR) фильтры, о которых чуть подробней в следующем разделе.

А что насчёт рекурсивных (IIR) фильтров?

Все рассмотренные в статье алгоритмы — FIR-фильтры (Finite Impulse Response, фильтры с конечной импульсной характеристикой). Это значит, что для вычисления каждого выходного пикселя мы берём конечное окно соседей и суммируем их с весами ядра. Размер этого окна растёт линейно с σ, и именно это делает свёртку при больших радиусах такой дорогой, даже в сепарабельном варианте.

Но существует принципиально иной подход — IIR-фильтры (Infinite Impulse Response, фильтры с бесконечной импульсной характеристикой). Их ключевая идея: вместо того чтобы заглядывать в окно из десятков соседей, каждый выходной пиксель вычисляется на основе нескольких предыдущих выходных значений и нескольких входных. Это рекурсия — отсюда и название «рекурсивные фильтры».

Формально для одномерного случая это выглядит примерно так:

y[n] = b₀·x[n] + b₁·x[n-1] + ... - a₁·y[n-1] - a₂·y[n-2] - ...

1	y[n] = b₀·x[n] + b₁·x[n-1] + ... - a₁·y[n-1] - a₂·y[n-2] - ...

где x — входной сигнал, y — выходной, а коэффициенты b и a подбираются так, чтобы импульсная характеристика фильтра аппроксимировала гауссиану. Количество коэффициентов фиксировано (обычно 3–4) и не зависит от σ. Это даёт сложность O(1) на пиксель при любом радиусе размытия — хоть σ = 5, хоть σ = 500.

И, должен признаться, я ни разу не занимался построением IIR-фильтра. Знаю общую теорию. Хочется пощупать тему на практике.

Основные алгоритмы

Наиболее известны три подхода:

Deriche (1993). Рашид Дериш предложил аппроксимацию гауссианы экспоненциальными функциями. Фильтр четвёртого порядка, реализуется двумя проходами по строке (вперёд и назад — каузальная и антикаузальная части). Даёт хорошую точность, хорошо изучен. Именно этот алгоритм используется, например, в GIMP.

Young — van Vliet (1995). Ян Янг и Люк ван Влит предложили альтернативную рекурсивную схему третьего порядка с иным способом подбора коэффициентов. Чуть проще в реализации, чуть менее точен на малых σ, но на практике разница едва заметна.

Vliet — Young — Verbeek (1998). Доработка предыдущего подхода с улучшенной численной устойчивостью.

Все три алгоритма сепарабельны: сначала обрабатываются строки, потом столбцы — точно так же, как мы делали с FIR-свёрткой.

Почему они не вошли в эту статью?

Причин несколько, и все они существенные:

Это другая парадигма. FIR-свёртка концептуально проста: вот ядро, вот пиксели, перемножили-сложили. Рекурсивные фильтры требуют понимания z-преобразования, передаточных функций, устойчивости, это совсем другой разговор, и он заслуживает полноценного изложения, а не абзаца в конце и без того длинной статьи.

Нетривиальная обработка краёв. В FIR-фильтре краевые эффекты — это неудобство, которое решается парой условий. В IIR-фильтре неправильная инициализация рекурсии на границе изображения приводит к видимым артефактам, расползающимся на десятки пикселей. Корректная инициализация — это отдельная тема с несколькими подходами, каждый со своими компромиссами.

Числовая точность. Рекурсивные фильтры работают с числами с плавающей точкой и накапливают ошибку от пикселя к пикселю. Трюк с fixed-point арифметикой, который дал нам ощутимый прирост скорости в FIR-варианте, здесь применим далеко не так просто, так как потеря точности в рекурсии может привести к нестабильности фильтра.

Сложность корректного сравнения. Сравнивать FIR и IIR на малых σ (1–5) не очень честно: FIR и так быстр, а IIR несёт накладные расходы на инициализацию. На больших σ (50+) IIR выигрывает с разгромным счётом, но это другой сценарий использования. Для объективного бенчмарка нужно продумать методологию.

Так что рекурсивный Gaussian — это тема для отдельной, полноценной статьи. Если эта серия до неё доберётся, там будет что изучить и о чём поговорить.

Скачать

Друзья, спасибо за внимание!

Исходник (zip) 510 Кб. Delphi XE 7, 13.0 Florence

Исполняемый файл (zip) 1.54 Мб.

5 1 голос

Рейтинг статьи

0 комментариев

Старые

Новые Популярные

Межтекстовые Отзывы

Посмотреть все комментарии