Blur в Delphi. Часть II: Box, Stack, Downscale

В предыдущей части мы рассмотрели эталонный Gaussian Blur и способы ускорения этого процесса. Но там мы достигли потолка, сложность O(r) преодолеть не получилось. В этой части мы рассмотрим методы, приближающие нас к Гауссу, зато вычислительная сложность которых равна O(1).

В какой-то мере, это описание математических и алгоритмических фокусов. Попытался минимизировать математику и сконцентрироваться на понимании происходящего.

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

Box Blur

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

Скользящая сумма: O(1) на пиксель

Идея

Пример

Реализация

Стоит ли оно того?

Stack Blur

Алгоритм на ладони

Шаг 1. Исходное состояние

Шаг 2. Вычитаем слева

Шаг 3. Прибавляем справа

Идея

Приложение: Сравнение качества

Заключение

Скачать

Box Blur

Gaussian Blur даёт качественный результат, но даже оптимизированная версия остаётся O(r) на пиксель, время растёт с радиусом. Можно ли размывать быстрее? Можно, если сменить форму ядра.

Gaussian взвешивает соседей по колоколообразной кривой, ближние пиксели важнее, дальние почти не влияют. Box Blur проще — все соседи в окне равнозначны, результат — обычное среднее арифметическое. Ядро — прямоугольник (отсюда и название «box»):

Gaussian (R=3):  1  6  15  20  15  6  1
Box (R=3):       1  1   1   1   1  1  1

1 2	Gaussian (R=3): 1 6 15 20 15 6 1 Box (R=3): 1 1 1 1 1 1 1

Качество хуже (сильно хуже) — вместо плавного перехода на краях объектов появляются заметные ступеньки. Но одинаковые веса откроют нам трюк, невозможный с Gaussian. До него мы доберёмся через шаг, сначала напишем наивную версию, чтобы было с чем сравнивать.

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

Структура алгоритма один в один повторяет сепарабельный Gaussian, только вместо взвешенной суммы используем простое среднее (все веса одинаковы). Обратите внимание: параметр здесь радиус, не сигма σ:

Листинг: Наивный Box Blur

procedure BlurBoxNaive(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Radius: Integer;
  W, H, X, Y, I: Integer;
  SX, SY: Integer;
  SumR, SumG, SumB: Integer;
  Count: Integer;
  SrcCache, TmpCache, DstCache: TScanLineCache;
  Tmp: TBitmap;
  SrcPx: TPixel32;
begin
  Radius := Round(Value.AsType<Double>);
  if Radius < 1 then
    Radius := 1;

  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;
        Count := 0;

        for I := -Radius to Radius do
        begin
          SX := X + I;
          if (SX >= 0) and (SX < W) then
          begin
            SrcPx := SrcCache[Y][SX];
            SumB := SumB + SrcPx.B;
            SumG := SumG + SrcPx.G;
            SumR := SumR + SrcPx.R;
            Inc(Count);
          end;
        end;

        with TmpCache[Y][X] do
        begin
          B := SumB div Count;
          G := SumG div Count;
          R := SumR div Count;
          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;
        Count := 0;

        for I := -Radius to Radius do
        begin
          SY := Y + I;
          if (SY >= 0) and (SY < H) then
          begin
            SrcPx := TmpCache[SY][X];
            SumB := SumB + SrcPx.B;
            SumG := SumG + SrcPx.G;
            SumR := SumR + SrcPx.R;
            Inc(Count);
          end;
        end;

        with DstCache[Y][X] do
        begin
          B := SumB div Count;
          G := SumG div Count;
          R := SumR div Count;
          A := TmpCache[Y][X].A;
        end;
      end;
    end;

  finally
    Tmp.Free;
  end;
end;

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

out Result: TBitmap);

var

Radius: Integer;

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

SX, SY: Integer;

SumR, SumG, SumB: Integer;

Count: Integer;

SrcCache, TmpCache, DstCache: TScanLineCache;

Tmp: TBitmap;

SrcPx: TPixel32;

begin

Radius := Round(Value.AsType<Double>);

if Radius < 1 then

Radius := 1;

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;

Count := 0;

for I := -Radius to Radius do

begin

SX := X + I;

if (SX >= 0) and (SX < W) then

begin

SrcPx := SrcCache[Y][SX];

SumB := SumB + SrcPx.B;

SumG := SumG + SrcPx.G;

SumR := SumR + SrcPx.R;

Inc(Count);

end;

with TmpCache[Y][X] do

begin

B := SumB div Count;

G := SumG div Count;

R := SumR div Count;

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;

Count := 0;

for I := -Radius to Radius do

begin

SY := Y + I;

if (SY >= 0) and (SY < H) then

begin

SrcPx := TmpCache[SY][X];

SumB := SumB + SrcPx.B;

SumG := SumG + SrcPx.G;

SumR := SumR + SrcPx.R;

Inc(Count);

end;

with DstCache[Y][X] do

begin

B := SumB div Count;

G := SumG div Count;

R := SumR div Count;

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

end;

finally

Tmp.Free;

end;

[свернуть]

Все методы из первой части принимали сигму — стандартное отклонение нормального распределения. Для одиночного Box Blur это не самый естественный параметр: пользователь обычно задаёт радиус напрямую. Но у нас есть демо-приложение и бенчмарк, где все методы должны работать в одной системе координат. Если для одних мы указываем сигму, а для других — радиус, диаграмма теряет смысл: что значит «при σ=5 Box работает быстрее Gaussian», если Box на самом деле получил не σ=5, а произвольный радиус?

Поэтому напишем процедуру-обёртку, которая транслирует сигму в эквивалентный радиус Box и вызывает наивную реализацию. Для этого нам нужна формула пересчёта.

Дисперсия дискретной случайной величины, равномерно распределенной на отрезке из w целых чисел, вычисляется по формуле:

σ² = (w² - 1) / 12

1	σ² = (w² - 1) / 12

Примечание: число 12 в знаменателе — это фундаментальная константа для дисперсии равномерного распределения (как непрерывного, так и дискретного).

Выразим w через σ:

w = Sqrt(12 · σ² + 1)

1	w = Sqrt(12 · σ² + 1)

В нашем случае, w — это ширина окна. В графических алгоритмах под радиусом R обычно понимается количество пикселей в одну сторону от центрального. Ширина окна связана с радиусом формулой:

r = (w - 1) / 2

1	r = (w - 1) / 2

Подставив в эту формулу вместо w выражение чуть выше, получим:

r = (Sqrt(12 * σ² + 1) - 1) * 0.5

1	r = (Sqrt(12 * σ² + 1) - 1) * 0.5

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

procedure BlurBoxNaiveSigma(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var Sigma, Radius: Double;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;
  Radius := Max(1, (Sqrt(12.0 * Sigma * Sigma + 1.0) - 1) * 0.5);
  BlurBoxNaive(Bitmap, TValue.From<Double>(Radius), Result);
end;

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

out Result: TBitmap);

var Sigma, Radius: Double;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

Radius := Max(1, (Sqrt(12.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

BlurBoxNaive(Bitmap, TValue.From<Double>(Radius), Result);

end;

Теперь можем спокойно вызывать эту функцию в интерфейсе. В котором и наблюдаем обещанные прямоугольные артефакты.

Всё знакомо: два прохода, внутренний цикл по ядру, проверка границ. Сложность O(r) на пиксель, как и у сепарабельного Gaussian. Казалось бы, что тут такого революционного? Но именно благодаря тому, что все веса одинаковы, Box Blur делает возможным один очень впечатляющий трюк, который невозможен с Gaussian, где каждый пиксель имеет свой вес — метод Скользящей суммы.

Скользящая сумма: O(1) на пиксель

При сдвиге окна на один пиксель вправо сумма меняется предсказуемо: уходит один левый элемент, приходит один правый. Пересчитывать всё окно не нужно:

Sum(X+1) = Sum(X) - Pixel[X - R] + Pixel[X + R + 1]

1	Sum(X+1) = Sum(X) - Pixel[X - R] + Pixel[X + R + 1]

Комментарий к иллюстрации:

Допустим, мы уже знаем сумму для текущего положения (индексы 2-6): Sum(X)=10+20+30+40+50=150

Чтобы сдвинуть окно вправо, нам не нужно снова складывать 20+30+40+50+60, вместо этого мы делаем инкрементальное обновление:

Действие	Значение	Результат (Сумма)
Было	Текущая сумма	150
Вычитаем	Левый край (Pixel[2])	150-10=140
Добавляем	Новый правый (Pixel[7])	140+60=200

Два сложения и затем одно деление, и не важно, равен радиус 3 или 300:

Листинг: Скользящая сумма

procedure BlurBoxRunningSum(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Radius: Integer;
  W, H, X, Y: Integer;
  SumR, SumG, SumB: Integer;
  Count: Integer;
  LeftX, RightX: Integer;
  TopY, BottomY: Integer;
  SrcCache, TmpCache, DstCache: TScanLineCache;
  Tmp: TBitmap;
  AddPx, SubPx: TPixel32;
begin
  Radius := Round(Value.AsType<Double>);
  if Radius < 1 then
    Radius := 1;

  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
      // Инициализация: набираем первую сумму для X=0
      SumR := 0; SumG := 0; SumB := 0;
      Count := 0;

      for X := 0 to Min(Radius, W - 1) do
      begin
        AddPx := SrcCache[Y][X];
        SumR := SumR + AddPx.R;
        SumG := SumG + AddPx.G;
        SumB := SumB + AddPx.B;
        Inc(Count);
      end;

      // Записываем первый пиксель
      with TmpCache[Y][0] do
      begin
        R := SumR div Count;
        G := SumG div Count;
        B := SumB div Count;
        A := SrcCache[Y][0].A;
      end;

      // Скользим вправо
      for X := 1 to W - 1 do
      begin
        // Добавляем правый край
        RightX := X + Radius;
        if RightX < W then
        begin
          AddPx := SrcCache[Y][RightX];
          SumR := SumR + AddPx.R;
          SumG := SumG + AddPx.G;
          SumB := SumB + AddPx.B;
          Inc(Count);
        end;

        // Убираем левый край
        LeftX := X - Radius - 1;
        if LeftX >= 0 then
        begin
          SubPx := SrcCache[Y][LeftX];
          SumR := SumR - SubPx.R;
          SumG := SumG - SubPx.G;
          SumB := SumB - SubPx.B;
          Dec(Count);
        end;

        with TmpCache[Y][X] do
        begin
          R := SumR div Count;
          G := SumG div Count;
          B := SumB div Count;
          A := SrcCache[Y][X].A;
        end;
      end;
    end;

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

      for Y := 0 to Min(Radius, H - 1) do
      begin
        AddPx := TmpCache[Y][X];
        SumR := SumR + AddPx.R;
        SumG := SumG + AddPx.G;
        SumB := SumB + AddPx.B;
        Inc(Count);
      end;

      with DstCache[0][X] do
      begin
        R := SumR div Count;
        G := SumG div Count;
        B := SumB div Count;
        A := TmpCache[0][X].A;
      end;

      for Y := 1 to H - 1 do
      begin
        BottomY := Y + Radius;
        if BottomY < H then
        begin
          AddPx := TmpCache[BottomY][X];
          SumR := SumR + AddPx.R;
          SumG := SumG + AddPx.G;
          SumB := SumB + AddPx.B;
          Inc(Count);
        end;

        TopY := Y - Radius - 1;
        if TopY >= 0 then
        begin
          SubPx := TmpCache[TopY][X];
          SumR := SumR - SubPx.R;
          SumG := SumG - SubPx.G;
          SumB := SumB - SubPx.B;
          Dec(Count);
        end;

        with DstCache[Y][X] do
        begin
          R := SumR div Count;
          G := SumG div Count;
          B := SumB div Count;
          A := TmpCache[Y][X].A;
        end;
      end;
    end;

  finally
    Tmp.Free;
  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

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

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

out Result: TBitmap);

var

Radius: Integer;

W, H, X, Y: Integer;

SumR, SumG, SumB: Integer;

Count: Integer;

LeftX, RightX: Integer;

TopY, BottomY: Integer;

SrcCache, TmpCache, DstCache: TScanLineCache;

Tmp: TBitmap;

AddPx, SubPx: TPixel32;

begin

Radius := Round(Value.AsType<Double>);

if Radius < 1 then

Radius := 1;

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

// Инициализация: набираем первую сумму для X=0

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

Count := 0;

for X := 0 to Min(Radius, W - 1) do

begin

AddPx := SrcCache[Y][X];

SumR := SumR + AddPx.R;

SumG := SumG + AddPx.G;

SumB := SumB + AddPx.B;

Inc(Count);

end;

// Записываем первый пиксель

with TmpCache[Y][0] do

begin

R := SumR div Count;

G := SumG div Count;

B := SumB div Count;

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

end;

// Скользим вправо

for X := 1 to W - 1 do

begin

// Добавляем правый край

RightX := X + Radius;

if RightX < W then

begin

AddPx := SrcCache[Y][RightX];

SumR := SumR + AddPx.R;

SumG := SumG + AddPx.G;

SumB := SumB + AddPx.B;

Inc(Count);

end;

// Убираем левый край

LeftX := X - Radius - 1;

if LeftX >= 0 then

begin

SubPx := SrcCache[Y][LeftX];

SumR := SumR - SubPx.R;

SumG := SumG - SubPx.G;

SumB := SumB - SubPx.B;

Dec(Count);

end;

with TmpCache[Y][X] do

begin

R := SumR div Count;

G := SumG div Count;

B := SumB div Count;

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

end;

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

for X := 0 to W - 1 do

begin

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

Count := 0;

for Y := 0 to Min(Radius, H - 1) do

begin

AddPx := TmpCache[Y][X];

SumR := SumR + AddPx.R;

SumG := SumG + AddPx.G;

SumB := SumB + AddPx.B;

Inc(Count);

end;

with DstCache[0][X] do

begin

R := SumR div Count;

G := SumG div Count;

B := SumB div Count;

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

end;

for Y := 1 to H - 1 do

begin

BottomY := Y + Radius;

if BottomY < H then

begin

AddPx := TmpCache[BottomY][X];

SumR := SumR + AddPx.R;

SumG := SumG + AddPx.G;

SumB := SumB + AddPx.B;

Inc(Count);

end;

TopY := Y - Radius - 1;

if TopY >= 0 then

begin

SubPx := TmpCache[TopY][X];

SumR := SumR - SubPx.R;

SumG := SumG - SubPx.G;

SumB := SumB - SubPx.B;

Dec(Count);

end;

with DstCache[Y][X] do

begin

R := SumR div Count;

G := SumG div Count;

B := SumB div Count;

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

end;

finally

Tmp.Free;

end;

[свернуть]

Горизонтальный проход устроен просто: набираем начальную сумму для X=0, затем скользим вправо — на каждом шаге добавляем правый пиксель, убираем левый. Count отслеживает реальный размер окна у краёв изображения.

Вертикальный проход устроен также, но по столбцам. Здесь кэш-локальность хуже (соседние строки — разные адреса), однако сам факт O(1) на пиксель перевешивает: при радиусе 50 наивный Box делает 101 чтение на пиксель, скользящая сумма — всегда два.

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

procedure BlurBoxRunningSumSigma(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var Sigma, Radius: Double;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;
  Radius := Max(1, (Sqrt(12.0 * Sigma * Sigma + 1.0) - 1) * 0.5);
  BlurBoxRunningSum(Bitmap, TValue.From<Double>(Radius), Result);
end;

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

out Result: TBitmap);

var Sigma, Radius: Double;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

Radius := Max(1, (Sqrt(12.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

BlurBoxRunningSum(Bitmap, TValue.From<Double>(Radius), Result);

end;

В любом случае, качество одного прохода Box Blur оставляет желать лучшего — характерные «квадратные» артефакты видны невооружённым глазом. Но если повторить Box Blur три раза подряд, результат удивительно близок к Gaussian. Об этом — в следующем разделе.

Тройной Box Blur

Один проход Box Blur даёт характерные «квадратные» артефакты. Но если повторить его три раза — результат визуально почти неотличим от Gaussian.

Как это работает

Чтобы понять, как прямоугольники превращаются в колокол, нужно разобрать операцию свёртки (convolution). Именно она происходит, когда мы накладываем один фильтр на другой.

Математическое определение фильтра

Box Blur — это функция-прямоугольник. В одномерном виде она выглядит так:

Это «коробка» шириной 1 и высотой 1.

Первая свёртка: Прямоугольник * Прямоугольник

Когда мы применяем Box Blur второй раз, мы совершаем свертку функции саму с собой:

Визуально это можно представить так: один прямоугольник неподвижен, а второй «проезжает» сквозь него.

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

Источник: Свертка n квадратных импульсов в гауссовы

Результат — треугольная функция (Tent function):

Linear B-spline

Источник: Д. Лукашев Об оптимизации распределения узлов для параметризации диэлектрической функции с помощью B-сплайнов

[свернуть]

Вторая свёртка: Треугольник * Прямоугольник

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

На выходе получается B-сплайн второго порядка:

Это кусочно-параболическая кривая.
Она уже имеет гладкие «плечи» и закругленную вершину

Quadratic B-spline

[свернуть]

Это прямое следствие Центральной предельной теоремы: повторная свёртка любого распределения с самим собой стремится к гауссову, независимо от исходной формы.

Итак, тройной Box визуально похож на Gaussian. Осталось понять, какой радиус Box соответствует заданной σ.

Пересчёт радиуса.

Box и Gaussian параметризуются по-разному, поэтому для честного сравнения нужно пересчитать сигму Gaussian в радиус Box. Формула выводится из равенства дисперсий:

BoxRadius = (Sqrt(4 · σ² + 1) - 1) * 0.5

1	BoxRadius = (Sqrt(4 · σ² + 1) - 1) * 0.5

Вывод формулы Box-радиуса по Gaussian-сигме

Дисперсия одного Box Blur. Рассмотрим один проход фильтра Box Blur. В цифровой обработке сигналов ядро такого фильтра шириной w (где w = 2r + 1, r — радиус) представляет собой дискретное равномерное распределение.

Формула дисперсии дискретной случайной величины нам уже знакома:

Суммирование при трех проходах. Согласно свойствам дисперсии, при последовательной свёртке нескольких независимых фильтров их дисперсии складываются. Если мы применяем Box Blur трижды, общая дисперсия итогового распределения составит:

Подстановка и решение уравнения. Теперь подставим формулу дисперсии одного окна в уравнение суммы:

Наша цель — выразить ширину окна w через искомую дисперсию σ²_total (которую мы берем из желаемого Гауссова размытия):

Извлекаем корень, чтобы найти ширину окна w:

Переход от ширины (w) к радиусу (R). Знакомая нам формула по наивной реализации Box Blur:

Подставляем наше w в формулу для R:

Или:

R = Round( Sqrt(12 * σ² / 3 + 1) / 2 - 0.5 ) =>
BoxRadius = (Sqrt(4.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

1 2	R = Round( Sqrt(12 * σ² / 3 + 1) / 2 - 0.5 ) => BoxRadius = (Sqrt(4.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

[свернуть]

Формула гарантирует, что тройной Box Blur с этим радиусом даёт размытие, эквивалентное Gaussian с заданной σ.

procedure BlurBoxTriplePass(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma, BoxRadius: Double;
  Pass: Integer;
  Tmp, Src: TBitmap;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  // Пересчёт Sigma Gaussian -> radius Box для 3 проходов
  // Формула: r = Round((Sqrt((12 * Sigma^2 / 3) + 1) - 1) / 2)
  // Упрощённо: r примерно равен Round(Sigma * 1.15)
  // Это даёт 3*Box blur, максимально близкий к Gaussian с данной Sigma
  BoxRadius := Max(1, (Sqrt(4.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

  // Первый проход: из Bitmap -> Result
  BlurBoxRunningSum(Bitmap, TValue.From<Double>(BoxRadius), Result);

  // Второй и третий проходы: Result -> Tmp -> Result
  for Pass := 2 to 3 do
  begin
    Src := Result;
    try
      BlurBoxRunningSum(Src, TValue.From<Double>(BoxRadius), Tmp);
      Result := Tmp;
    finally
      Src.Free;
    end;
  end;
end;

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

out Result: TBitmap);

var

Sigma, BoxRadius: Double;

Pass: Integer;

Tmp, Src: TBitmap;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

// Пересчёт Sigma Gaussian -> radius Box для 3 проходов

// Формула: r = Round((Sqrt((12 * Sigma^2 / 3) + 1) - 1) / 2)

// Упрощённо: r примерно равен Round(Sigma * 1.15)

// Это даёт 3*Box blur, максимально близкий к Gaussian с данной Sigma

BoxRadius := Max(1, (Sqrt(4.0 * Sigma * Sigma + 1.0) - 1) * 0.5);

// Первый проход: из Bitmap -> Result

BlurBoxRunningSum(Bitmap, TValue.From<Double>(BoxRadius), Result);

// Второй и третий проходы: Result -> Tmp -> Result

for Pass := 2 to 3 do

begin

Src := Result;

try

BlurBoxRunningSum(Src, TValue.From<Double>(BoxRadius), Tmp);

Result := Tmp;

finally

Src.Free;

end;

Три вызова BlurBoxRunningSum, каждый — O(1) на пиксель. Итого — шесть линейных проходов (по два на каждый Box: горизонтальный + вертикальный). Скорость остаётся O(1) на пиксель, а качество — почти Gaussian.

Тройной Box Blur даёт хорошее приближение к Gaussian, но у него есть ограничение: радиус — целое число. Формула: R = Round(√(12σ²/n + 1) — 1) / 2

где n — количество проходов (итераций) фильтра, округляет результат, и фактическая σ размытия отличается от запрошенной. При малых σ ошибка особенно заметна: запросив σ=3, мы можем получить фактическую σ=2.7 или σ=3.4 — разница видна глазом. При больших σ округление некритично — разница в долях процента.

Можно ли точнее? Можно. Для этого существует метод Extended Box Blur — способ закрыть этот зазор, не увеличивая вычислительную нагрузку.

Extended Box Blur

Иван Кутскир, создатель онлайн-редактора Photopea, в своей статье предлагает элегантное решение: вместо одного округлённого радиуса использовать два соседних (в статье Кутскира ищем функцию boxesForGauss).

Идея

Часть проходов делается с радиусом r, часть — с r+1. Подбирая пропорцию, можно попасть в нужную σ с гораздо большей точностью.

Три прохода Box Blur с радиусами r₁, r₂, r₃ дают суммарную дисперсию:

σ² = (w₁² + w₂² + w₃² - 3) / 12

1	σ² = (w₁² + w₂² + w₃² - 3) / 12

где w_i = 2r_i + 1 — ширина ядра. Если мы ограничимся двумя значениями радиуса r и r+1, задача сводится к простому вопросу: сколько из трёх проходов делать с меньшим радиусом?

Кутскир даёт готовую формулу. Для n проходов:

Ideal width:  wIdeal = √(12σ²/n + 1)
Lower width:  wl = Floor(wIdeal)  (ближайшее нечётное снизу)
Upper width:  wu = wl + 2
Lower radius: r  = wl div 2
Upper radius: r + 1

Ideal width: wIdeal = √(12σ²/n + 1)

Lower width: wl = Floor(wIdeal) (ближайшее нечётное снизу)

Upper width: wu = wl + 2

Lower radius: r = wl div 2

Upper radius: r + 1

Количество проходов с меньшим радиусом (один в один как у автора метода):

m = Round((12σ² - n·wl² - 4n·wl - 3n) / (-4wl - 4))

1	m = Round((12σ² - n·wl² - 4n·wl - 3n) / (-4wl - 4))

В коде используется эквивалентная форма, выраженная через (wu² — wl²) — так она короче и нагляднее. Если в ней подставить wu = wl + 2 и сделать преобразования, то получится формула, приведённая выше. Здесь же формула дана в том виде, в котором она фигурирует в оригинальной статье Кутскира.

Оставшиеся n — m проходов используют r + 1.

Пример

Запрошено σ = 3, три прохода:

wIdeal = √(12·9/3 + 1) = √37 ≈ 6.08
wl = 5,  wu = 7
r  = 2,  r+1 = 3
m  = Round((108 - 3*25 - 60 - 9) / (-24)) = Round(-36 / -24) = Round(1.5) = 2

wIdeal = √(12·9/3 + 1) = √37 ≈ 6.08

wl = 5, wu = 7

r = 2, r+1 = 3

m = Round((108 - 3*25 - 60 - 9) / (-24)) = Round(-36 / -24) = Round(1.5) = 2

Результат: два прохода с r=2, один проход с r=3.

Фактическая σ:

σ² = (5² + 5² + 7² - 3) / 12 = (25 + 25 + 49 - 3) / 12 = 96/12 = 8.0
σ = √8 ≈ 2.83

1 2	σ² = (5² + 5² + 7² - 3) / 12 = (25 + 25 + 49 - 3) / 12 = 96/12 = 8.0 σ = √8 ≈ 2.83

Обычный тройной Box с округлённым r=2 даёт σ² = (75-3)/12 = 6.0, σ ≈ 2.45. С r=3: σ² = (147-3)/12 = 12.0, σ ≈ 3.46.

Extended попал в 2.83 — ближе к запрошенным 3.0, чем оба варианта с фиксированным радиусом.

Реализация

procedure BlurBoxExtended(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma, R: Double;
  N: Integer;
  wIdeal: Double;
  wl, wu, m, i: Integer;
  Radii: array of Integer;
  Tmp, Src: TBitmap;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  N := 3; // количество проходов

  // Шаг 1: идеальная ширина окна (вещественная)
  wIdeal := Sqrt(12.0 * Sigma * Sigma / N + 1.0);

  // Шаг 2: нижняя нечётная ширина и верхняя нечётная ширина
  wl := Trunc(wIdeal);
  if wl mod 2 = 0 then
    Dec(wl);
  if wl < 1 then
    wl := 1;
  wu := wl + 2;

  // Шаг 3: сколько проходов делать с меньшей шириной (m),
  //         остальные (N - m) — с большей
  //
  // Из равенства дисперсий:
  //   m * (wl² - 1)/12 + (N - m) * (wu² - 1)/12 = Sigma²
  //
  // Решаем относительно m:
  //   m = (N * wu² - N - 12 * Sigma²) / (wu² - wl²)
  // Если подставим wu = wl+2, то получим (формула в статье):
  //   m = Round((12σ² - n·wl² - 4n·wl - 3n) / (-4wl - 4))

  // Округляем и зажимаем в [0, N]
  if wu * wu <> wl * wl then // Деление на 0
    m := Round((N * wu * wu - N - 12.0 * Sigma * Sigma) / (wu * wu - wl * wl))
  else
    m := N;

  if m < 0 then m := 0;
  if m > N then m := N;

  // Шаг 4: собираем массив радиусов для каждого прохода
  SetLength(Radii, N);
  for i := 0 to N - 1 do
  begin
    if i < m then
      Radii[i] := Max(1, (wl - 1) div 2)
    else
      Radii[i] := Max(1, (wu - 1) div 2);
  end;

  // Шаг 5: последовательные проходы Box Blur (скользящая сумма)
  R := Radii[0];
  BlurBoxRunningSum(Bitmap, TValue.From<Double>(R), Result);

  for i := 1 to N - 1 do
  begin
    Src := Result;
    try
      R := Radii[i];
      BlurBoxRunningSum(Src, TValue.From<Double>(R), Tmp);
      Result := Tmp;
    finally
      Src.Free;
    end;
  end;
end;

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

out Result: TBitmap);

var

Sigma, R: Double;

N: Integer;

wIdeal: Double;

wl, wu, m, i: Integer;

Radii: array of Integer;

Tmp, Src: TBitmap;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

N := 3; // количество проходов

// Шаг 1: идеальная ширина окна (вещественная)

wIdeal := Sqrt(12.0 * Sigma * Sigma / N + 1.0);

// Шаг 2: нижняя нечётная ширина и верхняя нечётная ширина

wl := Trunc(wIdeal);

if wl mod 2 = 0 then

Dec(wl);

if wl < 1 then

wl := 1;

wu := wl + 2;

// Шаг 3: сколько проходов делать с меньшей шириной (m),

// остальные (N - m) — с большей

// Из равенства дисперсий:

// m * (wl² - 1)/12 + (N - m) * (wu² - 1)/12 = Sigma²

// Решаем относительно m:

// m = (N * wu² - N - 12 * Sigma²) / (wu² - wl²)

// Если подставим wu = wl+2, то получим (формула в статье):

// m = Round((12σ² - n·wl² - 4n·wl - 3n) / (-4wl - 4))

// Округляем и зажимаем в [0, N]

if wu * wu <> wl * wl then // Деление на 0

m := Round((N * wu * wu - N - 12.0 * Sigma * Sigma) / (wu * wu - wl * wl))

else

m := N;

if m < 0 then m := 0;

if m > N then m := N;

// Шаг 4: собираем массив радиусов для каждого прохода

SetLength(Radii, N);

for i := 0 to N - 1 do

begin

if i < m then

Radii[i] := Max(1, (wl - 1) div 2)

else

Radii[i] := Max(1, (wu - 1) div 2);

end;

// Шаг 5: последовательные проходы Box Blur (скользящая сумма)

R := Radii[0];

BlurBoxRunningSum(Bitmap, TValue.From<Double>(R), Result);

for i := 1 to N - 1 do

begin

Src := Result;

try

R := Radii[i];

BlurBoxRunningSum(Src, TValue.From<Double>(R), Tmp);

Result := Tmp;

finally

Src.Free;

end;

Код структурно идентичен тройному Box Blur. Единственное отличие — вместо одного радиуса для всех проходов мы вычисляем массив Radii, где часть элементов равна r, часть — r+1. Остальное без изменений.

На скрине сравниваем тройной проход Box Blur с эталонным Гауссом, который получили через D2D-эффект. Визуально — никакой разницы.

Стоит ли оно того?

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

Выгода — в точности. При σ=3 обычный тройной Box промахивается на ~18% (фактическая σ=2.45 или 3.46). Extended попадает в 2.83 — ошибка 5.7%. При больших σ разница между методами стирается: округление одного пикселя на фоне ядра шириной 100 — ничто.

Extended Box имеет смысл, если вам важна предсказуемость результата при малых σ, и не стоит ничего в плане производительности. Бесплатное улучшение — лучший вид улучшения.

Тройной Box Blur — это один из лучших компромиссов между скоростью и качеством в обработке изображений. Но Марио Клингеманн был им всегда недоволен — либо скорость низкая, либо качество невысокое и придумал свой блюр — Stack Blur, который быстрее тройного и значительно лучше одинарного.

Stack Blur

Stack Blur придумал Марио Клингеманн (Mario Klingemann) в 2004 году. Задача: получить качество лучше Box Blur, но сохранить скорость O(1) на пиксель. Решение — треугольное ядро, поддерживаемое через хитрую структуру данных.

Вспомним, как выглядят ядра трёх методов при радиусе 3:

Box:      1  1   1  1  1  1  1     — прямоугольник
Stack:    1  2   3  4  3  2  1     — треугольник
Gaussian: 1  6  15 20 15  6  1     — колокол

Box: 1 1 1 1 1 1 1 — прямоугольник

Stack: 1 2 3 4 3 2 1 — треугольник

Gaussian: 1 6 15 20 15 6 1 — колокол

Треугольник — это, по сути, два box-а, свёрнутых друг с другом (один проход Box Blur по результату другого). Качество заметно лучше одиночного Box, артефакты на краях мягче. А Клингеманн нашёл способ поддерживать такое взвешенное среднее инкрементально — без пересчёта всей суммы при каждом сдвиге.

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

Пусть вас не отпугивает выражение «хитрая структура данных». Всё на самом деле очень просто, надо только ухватить суть этого решения и весь алгоритм будет как на ладони.

Алгоритм на ладони

Вот хорошо было в Box Blur — все пиксели имеют один вес. Это позволило реализовать скользящую сумму — при каждом шаге вычитаем пиксель, оставшийся за окном слева и добавляем вошедший справа. Чтобы найти очередной цвет, мы просто делим на количество пикселей в скользящей сумме. А что делать, когда веса у пикселей разные? Оказывается, треугольные веса можно обновлять так же инкрементально — через два вспомогательных счётчика. Разберём по шагам.

Шаг 1. Исходное состояние

Давайте рассмотрим случай для R=2 строки из пикселей: p1, p2, p3, p4, p5, p6. Размер окна равен 5 (2R+1).

Для вычисления цвета в центре окна, нам необходима взвешенная сумма всех пикселей, входящих в окно. Вес пикселя можно определить по пирамидке над окном. Таким образом, взвешенная сумма для начального случая равна: Sum = p1 + 2*p2 + 3*p3 + 2*p4 + p5.

Шаг 2. Вычитаем слева

Теперь мы двигаем окно вправо, чтобы центром стал пиксель p4. Введём понятие суммы слева. Это сумма пикселей левых пикселей окна, включая центральный. В нашем случае, это: SumOut = p1 + p2 + p3. Давайте вычтем её из Sum:

Вычитая из взвешенной суммы пиксели из «левой» суммы, мы уменьшаем вес оставшихся пикселей на единицу. Sum′ = Sum — SumOut = p1 + 2*p2 + 3*p3 + 2*p4 + p5 — (p1 + p2 + p3) = p2 + 2*p3 + 2*p4 + p5.

Шаг 3. Прибавляем справа

Введём понятие суммы справа. Это сумма пикселей справа от центра в текущем окне. В нашем случае, это: SumIn = p4 + p5. Давайте прибавим её к Sum и добавим новый пиксель, заходящий в новое окно — p6:

Ого! А это действие добавляет единичку в пикселям во взвешенной сумме и пристыковывает новый пиксель справа с весом =1. Sum″ = Sum′ + SumIn + p_right = p2 + 2*p3 + 2*p4 + p5 + (p4 + p5) + p6 = p2 + 2*p3 + 3*p4 + 2*p5+p6. Полученное выражение — это ни что иное, как взвешенная сумма для пикселя p4 в окне (p2-p6).

Шаг 4. Результат

Таким образом, пирамидка весов переместилась вправо:

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

Как работает «стек»

Название метода — от структуры данных в его основе. Представьте кольцевой буфер размером 2R + 1, хранящий пиксели текущего окна. Буфер разделён на две зоны:

Out-зона (левая половина + центр) — пиксели, которые при следующем сдвиге потеряют один уровень веса или выпадут из окна.
In-зона (правая половина) — пиксели, которые при следующем сдвиге наберут один уровень веса.

Алгоритм поддерживает три суммы одновременно:

Сумма	Что хранит
Sum	Взвешенная сумма (для вычисления результата)
InSum	Сумма значений пикселей в In-зоне
OutSum	Сумма значений пикселей в Out-зоне

При сдвиге окна на один пиксель происходит восемь шагов — каждый за O(1):

Записать результат: Sum div Divisor, где Divisor = (R+1)².
Вычесть OutSum из Sum — все «уходящие» пиксели теряют по одному уровню веса.
Убрать самый старый элемент из OutSum — он полностью покидает окно.
Поместить новый пиксель в стек на освободившееся место.
Добавить новый пиксель в InSum — он входит в окно с весом 1.
Прибавить InSum к Sum — все «входящие» пиксели набирают по уровню веса.
Средний элемент переходит из In в Out — он достиг пика треугольника и начинает терять вес.
Сдвинуть указатель кольцевого буфера.

Почему Divisor = (R+1)²? Сумма весов треугольного ядра 1 + 2 + … + (R+1) + … + 2 + 1 равна (R+1)². Это удобное целое число — деление можно заменить на умножение с обратным значением или оставить как есть. Здесь мы оставим как есть, и возможно в будущем, в третьей части (вряд ли скоро), будем оптимизировать скорость.

Реализация

Треугольное ядро сепарабельно — двумерный треугольник раскладывается на два одномерных, поэтому горизонтальный и вертикальный проходы дают корректный результат.

Горизонтальный проход. Чётко придерживаемся плана, изложенного в разделе выше.

Листинг: Горизонтальный проход

procedure StackBlurHorz(SrcCache, DstCache: TScanLineCache; 
  W, H, Radius: Integer);
var
  X, Y, I, SrcIdx, StackIdx: Integer;
  Divisor, StackSize: Integer;
  SumR, SumG, SumB: Integer;
  InR, InG, InB: Integer;
  OutR, OutG, OutB: Integer;
  Px: TPixel32;
  Stack: array of TPixel32;
  Ri1: Integer; // Radius + 1
begin
  StackSize := Radius * 2 + 1;
  SetLength(Stack, StackSize);
  Ri1 := Radius + 1;
  Divisor := Ri1 * Ri1;

  for Y := 0 to H - 1 do
  begin
    SumR := 0; SumG := 0; SumB := 0;
    OutR := 0; OutG := 0; OutB := 0;
    InR := 0;  InG := 0;  InB := 0;

    // Инициализация стека для X=0
    // Левая часть стека (позиции 0..Radius): clamp к левому краю
    for I := 0 to Radius do
    begin
      SrcIdx := Clamp(I - Radius, 0, W - 1); // при R=5: -5,-4,...,0
      Px := SrcCache[Y][SrcIdx];
      Stack[I] := Px;
      // Вес = I + 1 (для позиции 0 вес=1, для Radius вес=Ri1)
      SumR := SumR + Px.R * (I + 1);
      SumG := SumG + Px.G * (I + 1);
      SumB := SumB + Px.B * (I + 1);
      OutR := OutR + Px.R;
      OutG := OutG + Px.G;
      OutB := OutB + Px.B;
    end;

    // Правая часть стека (позиции Radius+1..StackSize-1)
    for I := 1 to Radius do
    begin
      SrcIdx := Clamp(I, 0, W - 1);
      Px := SrcCache[Y][SrcIdx];
      Stack[Radius + I] := Px;
      // Вес = Ri1 - I
      SumR := SumR + Px.R * (Ri1 - I);
      SumG := SumG + Px.G * (Ri1 - I);
      SumB := SumB + Px.B * (Ri1 - I);
      InR := InR + Px.R;
      InG := InG + Px.G;
      InB := InB + Px.B;
    end;

    StackIdx := 0; // указатель на самый старый элемент

    //  Скольжение по строке
    for X := 0 to W - 1 do
    begin
      // 1. Записываем результат
      with DstCache[Y][X] do
      begin
        R := Clamp(SumR div Divisor, 0, 255);
        G := Clamp(SumG div Divisor, 0, 255);
        B := Clamp(SumB div Divisor, 0, 255);
        A := SrcCache[Y][X].A;
      end;

      // 2. Вычитаем OutSum из общей суммы
      SumR := SumR - OutR;
      SumG := SumG - OutG;
      SumB := SumB - OutB;

      // 3. Убираем самый старый элемент из OutSum
      Px := Stack[StackIdx];
      OutR := OutR - Px.R;
      OutG := OutG - Px.G;
      OutB := OutB - Px.B;

      // 4. Новый пиксель занимает место старого в стеке
      SrcIdx := Clamp(X + Ri1, 0, W - 1);
      Px := SrcCache[Y][SrcIdx];
      Stack[StackIdx] := Px;

      // 5. Добавляем новый пиксель в InSum
      InR := InR + Px.R;
      InG := InG + Px.G;
      InB := InB + Px.B;

      // 6. Прибавляем InSum к общей сумме
      SumR := SumR + InR;
      SumG := SumG + InG;
      SumB := SumB + InB;

      // 7. Средний элемент переходит из In -> Out
      //    Это элемент на позиции (StackIdx + Ri1) mod StackSize
      I := StackIdx + Ri1;
      if I >= StackSize then
        I := I - StackSize;

      Px := Stack[I];
      InR := InR - Px.R;
      InG := InG - Px.G;
      InB := InB - Px.B;
      OutR := OutR + Px.R;
      OutG := OutG + Px.G;
      OutB := OutB + Px.B;

      // 8. Сдвигаем указатель стека
      Inc(StackIdx);
      if StackIdx >= StackSize then
        StackIdx := 0;
    end;
  end;
end;

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

procedure StackBlurHorz(SrcCache, DstCache: TScanLineCache;

W, H, Radius: Integer);

var

X, Y, I, SrcIdx, StackIdx: Integer;

Divisor, StackSize: Integer;

SumR, SumG, SumB: Integer;

InR, InG, InB: Integer;

OutR, OutG, OutB: Integer;

Px: TPixel32;

Stack: array of TPixel32;

Ri1: Integer; // Radius + 1

begin

StackSize := Radius * 2 + 1;

SetLength(Stack, StackSize);

Ri1 := Radius + 1;

Divisor := Ri1 * Ri1;

for Y := 0 to H - 1 do

begin

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

OutR := 0; OutG := 0; OutB := 0;

InR := 0; InG := 0; InB := 0;

// Инициализация стека для X=0

// Левая часть стека (позиции 0..Radius): clamp к левому краю

for I := 0 to Radius do

begin

SrcIdx := Clamp(I - Radius, 0, W - 1); // при R=5: -5,-4,...,0

Px := SrcCache[Y][SrcIdx];

Stack[I] := Px;

// Вес = I + 1 (для позиции 0 вес=1, для Radius вес=Ri1)

SumR := SumR + Px.R * (I + 1);

SumG := SumG + Px.G * (I + 1);

SumB := SumB + Px.B * (I + 1);

OutR := OutR + Px.R;

OutG := OutG + Px.G;

OutB := OutB + Px.B;

end;

// Правая часть стека (позиции Radius+1..StackSize-1)

for I := 1 to Radius do

begin

SrcIdx := Clamp(I, 0, W - 1);

Px := SrcCache[Y][SrcIdx];

Stack[Radius + I] := Px;

// Вес = Ri1 - I

SumR := SumR + Px.R * (Ri1 - I);

SumG := SumG + Px.G * (Ri1 - I);

SumB := SumB + Px.B * (Ri1 - I);

InR := InR + Px.R;

InG := InG + Px.G;

InB := InB + Px.B;

end;

StackIdx := 0; // указатель на самый старый элемент

// Скольжение по строке

for X := 0 to W - 1 do

begin

// 1. Записываем результат

with DstCache[Y][X] do

begin

R := Clamp(SumR div Divisor, 0, 255);

G := Clamp(SumG div Divisor, 0, 255);

B := Clamp(SumB div Divisor, 0, 255);

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

end;

// 2. Вычитаем OutSum из общей суммы

SumR := SumR - OutR;

SumG := SumG - OutG;

SumB := SumB - OutB;

// 3. Убираем самый старый элемент из OutSum

Px := Stack[StackIdx];

OutR := OutR - Px.R;

OutG := OutG - Px.G;

OutB := OutB - Px.B;

// 4. Новый пиксель занимает место старого в стеке

SrcIdx := Clamp(X + Ri1, 0, W - 1);

Px := SrcCache[Y][SrcIdx];

Stack[StackIdx] := Px;

// 5. Добавляем новый пиксель в InSum

InR := InR + Px.R;

InG := InG + Px.G;

InB := InB + Px.B;

// 6. Прибавляем InSum к общей сумме

SumR := SumR + InR;

SumG := SumG + InG;

SumB := SumB + InB;

// 7. Средний элемент переходит из In -> Out

// Это элемент на позиции (StackIdx + Ri1) mod StackSize

I := StackIdx + Ri1;

if I >= StackSize then

I := I - StackSize;

Px := Stack[I];

InR := InR - Px.R;

InG := InG - Px.G;

InB := InB - Px.B;

OutR := OutR + Px.R;

OutG := OutG + Px.G;

OutB := OutB + Px.B;

// 8. Сдвигаем указатель стека

Inc(StackIdx);

if StackIdx >= StackSize then

StackIdx := 0;

end;

[свернуть]

Вертикальный проход — зеркальная копия горизонтального, где X и Y меняются ролями. Обходит память по столбцам, что менее эффективно для кэша. Оптимизация через транспонирование применима и здесь.

Листинг: Вертикальный проход

procedure StackBlurVert(SrcCache, DstCache: TScanLineCache; 
  W, H, Radius: Integer);
var
  X, Y, I, SrcIdx, StackIdx: Integer;
  Divisor, StackSize: Integer;
  SumR, SumG, SumB: Integer;
  InR, InG, InB: Integer;
  OutR, OutG, OutB: Integer;
  Px: TPixel32;
  Stack: array of TPixel32;
  Ri1: Integer;
begin
  StackSize := Radius * 2 + 1;
  SetLength(Stack, StackSize);
  Ri1 := Radius + 1;
  Divisor := Ri1 * Ri1;

  for X := 0 to W - 1 do
  begin
    SumR := 0; SumG := 0; SumB := 0;
    OutR := 0; OutG := 0; OutB := 0;
    InR := 0;  InG := 0;  InB := 0;

    for I := 0 to Radius do
    begin
      SrcIdx := Clamp(I - Radius, 0, H - 1);
      Px := SrcCache[SrcIdx][X];
      Stack[I] := Px;
      SumR := SumR + Px.R * (I + 1);
      SumG := SumG + Px.G * (I + 1);
      SumB := SumB + Px.B * (I + 1);
      OutR := OutR + Px.R;
      OutG := OutG + Px.G;
      OutB := OutB + Px.B;
    end;

    for I := 1 to Radius do
    begin
      SrcIdx := Clamp(I, 0, H - 1);
      Px := SrcCache[SrcIdx][X];
      Stack[Radius + I] := Px;
      SumR := SumR + Px.R * (Ri1 - I);
      SumG := SumG + Px.G * (Ri1 - I);
      SumB := SumB + Px.B * (Ri1 - I);
      InR := InR + Px.R;
      InG := InG + Px.G;
      InB := InB + Px.B;
    end;

    StackIdx := 0;

    for Y := 0 to H - 1 do
    begin
      with DstCache[Y][X] do
      begin
        R := Clamp(SumR div Divisor, 0, 255);
        G := Clamp(SumG div Divisor, 0, 255);
        B := Clamp(SumB div Divisor, 0, 255);
        A := SrcCache[Y][X].A;
      end;

      SumR := SumR - OutR;
      SumG := SumG - OutG;
      SumB := SumB - OutB;

      Px := Stack[StackIdx];
      OutR := OutR - Px.R;
      OutG := OutG - Px.G;
      OutB := OutB - Px.B;

      SrcIdx := Clamp(Y + Ri1, 0, H - 1);
      Px := SrcCache[SrcIdx][X];
      Stack[StackIdx] := Px;

      InR := InR + Px.R;
      InG := InG + Px.G;
      InB := InB + Px.B;

      SumR := SumR + InR;
      SumG := SumG + InG;
      SumB := SumB + InB;

      I := StackIdx + Ri1;
      if I >= StackSize then
        I := I - StackSize;

      Px := Stack[I];
      InR := InR - Px.R;
      InG := InG - Px.G;
      InB := InB - Px.B;
      OutR := OutR + Px.R;
      OutG := OutG + Px.G;
      OutB := OutB + Px.B;

      Inc(StackIdx);
      if StackIdx >= StackSize then
        StackIdx := 0;
    end;
  end;
end;

100

procedure StackBlurVert(SrcCache, DstCache: TScanLineCache;

W, H, Radius: Integer);

var

X, Y, I, SrcIdx, StackIdx: Integer;

Divisor, StackSize: Integer;

SumR, SumG, SumB: Integer;

InR, InG, InB: Integer;

OutR, OutG, OutB: Integer;

Px: TPixel32;

Stack: array of TPixel32;

Ri1: Integer;

begin

StackSize := Radius * 2 + 1;

SetLength(Stack, StackSize);

Ri1 := Radius + 1;

Divisor := Ri1 * Ri1;

for X := 0 to W - 1 do

begin

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

OutR := 0; OutG := 0; OutB := 0;

InR := 0; InG := 0; InB := 0;

for I := 0 to Radius do

begin

SrcIdx := Clamp(I - Radius, 0, H - 1);

Px := SrcCache[SrcIdx][X];

Stack[I] := Px;

SumR := SumR + Px.R * (I + 1);

SumG := SumG + Px.G * (I + 1);

SumB := SumB + Px.B * (I + 1);

OutR := OutR + Px.R;

OutG := OutG + Px.G;

OutB := OutB + Px.B;

end;

for I := 1 to Radius do

begin

SrcIdx := Clamp(I, 0, H - 1);

Px := SrcCache[SrcIdx][X];

Stack[Radius + I] := Px;

SumR := SumR + Px.R * (Ri1 - I);

SumG := SumG + Px.G * (Ri1 - I);

SumB := SumB + Px.B * (Ri1 - I);

InR := InR + Px.R;

InG := InG + Px.G;

InB := InB + Px.B;

end;

StackIdx := 0;

for Y := 0 to H - 1 do

begin

with DstCache[Y][X] do

begin

R := Clamp(SumR div Divisor, 0, 255);

G := Clamp(SumG div Divisor, 0, 255);

B := Clamp(SumB div Divisor, 0, 255);

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

end;

SumR := SumR - OutR;

SumG := SumG - OutG;

SumB := SumB - OutB;

Px := Stack[StackIdx];

OutR := OutR - Px.R;

OutG := OutG - Px.G;

OutB := OutB - Px.B;

SrcIdx := Clamp(Y + Ri1, 0, H - 1);

Px := SrcCache[SrcIdx][X];

Stack[StackIdx] := Px;

InR := InR + Px.R;

InG := InG + Px.G;

InB := InB + Px.B;

SumR := SumR + InR;

SumG := SumG + InG;

SumB := SumB + InB;

I := StackIdx + Ri1;

if I >= StackSize then

I := I - StackSize;

Px := Stack[I];

InR := InR - Px.R;

InG := InG - Px.G;

InB := InB - Px.B;

OutR := OutR + Px.R;

OutG := OutG + Px.G;

OutB := OutB + Px.B;

Inc(StackIdx);

if StackIdx >= StackSize then

StackIdx := 0;

end;

[свернуть]

Основная процедура связывает оба прохода. Как и в Box Blur, промежуточный буфер необходим — вертикальный проход не должен читать данные, которые сам же изменяет.

procedure BlurStackBlur(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Radius, W, H: Integer;
  Tmp: TBitmap;
  SrcCache, TmpCache, DstCache: TScanLineCache;
begin
  Radius := Round(Value.AsType<Double>);
  if Radius < 1 then
    Radius := 1;

  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);

    StackBlurHorz(SrcCache, TmpCache, W, H, Radius);
    StackBlurVert(TmpCache, DstCache, W, H, Radius);
  finally
    Tmp.Free;
  end;
end;

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

out Result: TBitmap);

var

Radius, W, H: Integer;

Tmp: TBitmap;

SrcCache, TmpCache, DstCache: TScanLineCache;

begin

Radius := Round(Value.AsType<Double>);

if Radius < 1 then

Radius := 1;

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);

StackBlurHorz(SrcCache, TmpCache, W, H, Radius);

StackBlurVert(TmpCache, DstCache, W, H, Radius);

finally

Tmp.Free;

end;

Два прохода, каждый — O(1) на пиксель. При этом, в отличие от тройного Box Blur, здесь всего два линейных прохода вместо шести.

По аналогии с предыдущими реализациями, где параметром выступает радиус, напишем аналог, принимающий сигму. Тут формула для радиуса может вызвать некоторое недоумение. Дело в том, что дисперсия треугольного ядра равна: σ² = ((R+1)² — 1)/ 6. Треугольный фильтр с радиусом R — это результат последовательного применения двух прямоугольных фильтров (Box Blur) с радиусом R/2 (или, точнее, со стороной w = R+1).

procedure BlurStackBlurSigma(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma, Radius: Double;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;
  // Пересчёт Sigma -> Radius для треугольного ядра
  // Дисперсия треугольного ядра: σ² = ((R+1)² - 1)/ 6
  // Решаем: R = Sqrt(6·σ² + 1) - 1
  Radius := Max(1, Sqrt(6.0 * Sigma * Sigma + 1) - 1.0);
  BlurStackBlur(Bitmap, Radius, Result);
end;

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

out Result: TBitmap);

var

Sigma, Radius: Double;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

// Пересчёт Sigma -> Radius для треугольного ядра

// Дисперсия треугольного ядра: σ² = ((R+1)² - 1)/ 6

// Решаем: R = Sqrt(6·σ² + 1) - 1

Radius := Max(1, Sqrt(6.0 * Sigma * Sigma + 1) - 1.0);

BlurStackBlur(Bitmap, Radius, Result);

end;

Появление этой процедуры обусловлено не только бенчмарком, но и тем, что в дальнейшем она будет использована в «хитром» методе размытия (несколькими разделами ниже).

По качеству Stack Blur (треугольное ядро) находится между одиночным Box и тройным Box (почти Gaussian). Тройной Box даёт колоколообразную кривую, более близкую к Gaussian, но требует шести линейных проходов. Stack Blur — компромисс: два прохода, качество лучше одиночного Box, но хуже тройного.

Stack Blur — один из самых популярных алгоритмов размытия в мобильных и веб-движках именно благодаря этому балансу.

Двойной Stack Blur

Раз тройной Box Blur приближает Гаусс благодаря Центральной предельной теореме, напрашивается вопрос: а что если применить Stack Blur дважды? Ядро Stack Blur — треугольник, а треугольник уже ближе к колоколу, чем прямоугольник. Значит, сходимость к Гауссу должна быть быстрее.

И это действительно так. Один Stack Blur — это, по сути, свёртка двух box-ядер, встроенная в один проход. Каждое повторное применение — ещё одна свёртка:

Метод	Эквивалент в box-свёртках	Форма ядра
1× Box	1	Прямоугольник
1× Stack	2	Треугольник
3× Box	3	Почти колокол
2× Stack	4	Гладкий колокол

Двойной Stack Blur эквивалентен четырём box-свёрткам и даёт B-сплайн третьего порядка — кубическую кривую, которая ближе к Гауссу, чем результат тройного Box. Звучит отлично. Но давайте посчитаем, во что это обходится.

Тройной Box Blur — это шесть линейных проходов, а двойной Stack — четыре. Проходов меньше, качество лучше — казалось бы, победа. Однако стоимость одного прохода у них принципиально разная:

Box Blur, один шаг скользящей суммы:
  Sum := Sum - Left + Right;          // 2 операции на канал

Stack Blur, один шаг:
  Sum := Sum - OutSum;                // обновить взвешенную сумму
  OutSum := OutSum - oldest;          // убрать старый из Out
  InSum := InSum + new;               // добавить новый в In
  Sum := Sum + InSum;                 // обновить взвешенную сумму
  // + переход In→Out, сдвиг указателя стека

Box Blur, один шаг скользящей суммы:

Sum := Sum - Left + Right; // 2 операции на канал

Stack Blur, один шаг:

Sum := Sum - OutSum; // обновить взвешенную сумму

OutSum := OutSum - oldest; // убрать старый из Out

InSum := InSum + new; // добавить новый в In

Sum := Sum + InSum; // обновить взвешенную сумму

// + переход In→Out, сдвиг указателя стека

Stack Blur выполняет примерно в 3–4 раза больше операций на пиксель. Умножаем на число проходов:

Метод	Проходов	Операций на проход	Итого	Качество
3× Box	6	~6	~36	≈ Гаусс
2× Stack	4	~20	~80	Чуть лучше

Двойной Stack Blur оказывается примерно вдвое дороже тройного Box при едва заметном выигрыше в качестве. А разницу между «почти Гаусс» и «ещё ближе к Гауссу» на реальных изображениях увидеть практически невозможно.

Поэтому, кажется, что двойной Stack Blur — приём теоретически корректный, но практически бессмысленный. Но давайте попробуем его реализовать и посмотрим в тестах, так ли это.

procedure BlurDoubleStackBlur(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma: Double;
  Radius: Integer;
  Tmp: TBitmap;
  RadiusValue: TValue;
begin
  Sigma := Value.AsType<Double>;
  Radius := Max(1, Round(Sqrt(3 * Sigma * Sigma + 1) - 1));
  RadiusValue := TValue.From<Double>(Radius);

  BlurStackBlur(Bitmap, RadiusValue, Tmp);
  try
    BlurStackBlur(Tmp, RadiusValue, Result);
  finally
    Tmp.Free;
  end;
end;

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

out Result: TBitmap);

var

Sigma: Double;

Radius: Integer;

Tmp: TBitmap;

RadiusValue: TValue;

begin

Sigma := Value.AsType<Double>;

Radius := Max(1, Round(Sqrt(3 * Sigma * Sigma + 1) - 1));

RadiusValue := TValue.From<Double>(Radius);

BlurStackBlur(Bitmap, RadiusValue, Tmp);

try

BlurStackBlur(Tmp, RadiusValue, Result);

finally

Tmp.Free;

end;

И вот, глядя на скрин, оценивая качество и время выполнения операции, кажется, что мы не зря потратили время на эту казалось бы бессмысленную затею.

До сих пор были математические фокусы, давайте теперь рассмотрим хитрость.

Хитрость вместо силы

Все предыдущие методы шли в лоб — обрабатывали каждый пиксель оригинала. Но если задуматься: зачем размывать изображение 4000×3000, если результат всё равно будет «кашей»? Мелкие детали, ради которых мы храним все 12 миллионов пикселей, при сильном размытии исчезнут. Значит, можно сначала уменьшить изображение, размыть маленькую копию (дёшево!), а затем растянуть обратно.

Идея

Четыре шага:

Downscale — уменьшаем изображение в Scale раз. Масштаб подбирается так, чтобы радиус размытия на уменьшенной копии составлял всего 3–5 пикселей.
Blur — размываем маленькое изображение любым методом (здесь — Stack Blur). На картинке 400×300 даже наивный алгоритм отработает мгновенно.
Upscale — растягиваем результат обратно до исходного размера. Билинейная интерполяция при увеличении сама по себе даёт дополнительное сглаживание.
Blur — очень лёгкий, сигма 1..2, сгладить оставшиеся артефакты.

При σ = 40 и изображении 4000×3000:

Scale = σ / 4 = 10 — уменьшаем в 10 раз.
Маленькая копия: 400×300 = 120 тысяч пикселей вместо 12 миллионов.
Радиус на маленькой копии: 40 / 10 = 4 пикселя — мгновенная свёртка.

Общая работа: два StretchBlt (выполняются аппаратно или высокооптимизированным кодом GDI) плюс размытие крошечной картинки. При больших σ этот подход оставляет далеко позади любую попиксельную обработку.

Реализация

procedure BlurDownscaleUpscale(Bitmap: TBitmap; const Value: TValue;
  out Result: TBitmap);
var
  Sigma, SmallSigma, Scale, FinalRadius: Double;
  SmallW, SmallH: Integer;
  SmallBmp, BlurredSmall, Tmp: TBitmap;
begin
  Sigma := Value.AsType<Double>;
  if Sigma < 0.5 then
    Sigma := 0.5;

  Scale := Max(1.0, Sigma / 4.0);
  SmallW := Max(2, Round(Bitmap.Width / Scale));
  SmallH := Max(2, Round(Bitmap.Height / Scale));
  SmallSigma := Sigma / Scale;
  FinalRadius := Max(1, Scale * 0.5);

  Bitmap.PixelFormat := pf32bit;
  // 1. Уменьшаем
  SmallBmp := TBitmap.Create;
  try
    SmallBmp.PixelFormat := pf32bit;
    SmallBmp.SetSize(SmallW, SmallH);

    SetStretchBltMode(SmallBmp.Canvas.Handle, HALFTONE);
    SetBrushOrgEx(SmallBmp.Canvas.Handle, 0, 0, nil);
    StretchBlt(
      SmallBmp.Canvas.Handle, 0, 0, SmallW, SmallH,
      Bitmap.Canvas.Handle, 0, 0, Bitmap.Width, Bitmap.Height,
      SRCCOPY
    );

    // 2. Размываем маленькое (Stack Blur)
    BlurStackBlurSigma(SmallBmp, TValue.From<Double>(SmallSigma), BlurredSmall);

    // 3. Увеличиваем обратно
    try
      Tmp := TBitmap.Create;
      try
        Tmp.PixelFormat := pf32bit;
        Tmp.SetSize(Bitmap.Width, Bitmap.Height);

        SetStretchBltMode(Tmp.Canvas.Handle, HALFTONE);
        SetBrushOrgEx(Tmp.Canvas.Handle, 0, 0, nil);
        StretchBlt(
          Tmp.Canvas.Handle, 0, 0, Bitmap.Width, Bitmap.Height,
          BlurredSmall.Canvas.Handle, 0, 0, SmallW, SmallH,
          SRCCOPY
        );

        // 4. Лёгкий проход по полноразмерному результату
        BlurBoxRunningSum(Tmp, TValue.From<Double>(FinalRadius), Result);
      finally
        Tmp.Free;
      end;
    finally
      BlurredSmall.Free;
    end;
  finally
    SmallBmp.Free;
  end;
end;

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

out Result: TBitmap);

var

Sigma, SmallSigma, Scale, FinalRadius: Double;

SmallW, SmallH: Integer;

SmallBmp, BlurredSmall, Tmp: TBitmap;

begin

Sigma := Value.AsType<Double>;

if Sigma < 0.5 then

Sigma := 0.5;

Scale := Max(1.0, Sigma / 4.0);

SmallW := Max(2, Round(Bitmap.Width / Scale));

SmallH := Max(2, Round(Bitmap.Height / Scale));

SmallSigma := Sigma / Scale;

FinalRadius := Max(1, Scale * 0.5);

Bitmap.PixelFormat := pf32bit;

// 1. Уменьшаем

SmallBmp := TBitmap.Create;

try

SmallBmp.PixelFormat := pf32bit;

SmallBmp.SetSize(SmallW, SmallH);

SetStretchBltMode(SmallBmp.Canvas.Handle, HALFTONE);

SetBrushOrgEx(SmallBmp.Canvas.Handle, 0, 0, nil);

StretchBlt(

SmallBmp.Canvas.Handle, 0, 0, SmallW, SmallH,

Bitmap.Canvas.Handle, 0, 0, Bitmap.Width, Bitmap.Height,

SRCCOPY

);

// 2. Размываем маленькое (Stack Blur)

BlurStackBlurSigma(SmallBmp, TValue.From<Double>(SmallSigma), BlurredSmall);

// 3. Увеличиваем обратно

try

Tmp := TBitmap.Create;

try

Tmp.PixelFormat := pf32bit;

Tmp.SetSize(Bitmap.Width, Bitmap.Height);

SetStretchBltMode(Tmp.Canvas.Handle, HALFTONE);

SetBrushOrgEx(Tmp.Canvas.Handle, 0, 0, nil);

StretchBlt(

Tmp.Canvas.Handle, 0, 0, Bitmap.Width, Bitmap.Height,

BlurredSmall.Canvas.Handle, 0, 0, SmallW, SmallH,

SRCCOPY

);

// 4. Лёгкий проход по полноразмерному результату

BlurBoxRunningSum(Tmp, TValue.From<Double>(FinalRadius), Result);

finally

Tmp.Free;

end;

finally

BlurredSmall.Free;

end;

finally

SmallBmp.Free;

end;

StretchBlt с режимом HALFTONE — качественное масштабирование силами GDI. При уменьшении он корректно усредняет пиксели (антиалиасинг), при увеличении — билинейная интерполяция. Никакого ручного ресэмплинга.

Ограничения

Метод не универсален. При малых σ (1–3) масштабирование бессмысленно — Scale получается 1, и мы просто делаем обычный Stack Blur. Выигрыш проявляется при σ ≥ 8–10, а по-настоящему впечатляет при σ ≥ 20, когда изображение сжимается в разы и размытие работает с крошечной копией.

Качество тоже не идеально — артефакты масштабирования могут быть заметны на резких границах. Но для большинства практических задач (фон под UI, эффект глубины резкости, превью) результат визуально неотличим от честного Gaussian.

Идея «размытие через масштабирование» доведена до предела в Dual Kawase Blur, где многоуровневая пирамида масштабирований сама по себе является размытием. Этот метод популярен в GPU-движках, но на CPU теряет смысл.

Бенчмарк

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

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

Диаграмма (64-bit)

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

Первое, что бросается в глаза, все наши методы выполнили обещание O(1) по радиусу. При десятикратном росте σ (с 5 до 50) время практически не меняется. Сравните с первой частью, где CPU-реализации Gaussian давали линейный рост с ростом сигмы.

Box Blur ×3 и Box Blur Extended — практически идентичны по скорости: 36-37 мс независимо от σ. Это ожидаемо: Extended отличается от обычного тройного прохода только тем, что часть проходов использует радиус r, а часть — r+1. Объём вычислений тот же, шесть проходов скользящей суммы (три горизонтальных + три вертикальных). Небольшая разница при σ=50 (37.5 мс против 36.8 мс) — статистический шум, а не реальная разница в производительности. Extended Box даёт более точное попадание в запрошенную σ бесплатно.

Stack Blur — самый быстрый из наших CPU-методов: 15-16 мс. Один проход с треугольным ядром вместо трёх проходов с прямоугольным. Вдвое меньше проходов по памяти — вдвое быстрее. Качество при этом сопоставимо с тройным Box (треугольник — хорошее приближение к Gaussian).

Stack Blur ×2 — двойной Stack Blur: 30-31 мс. Два прохода дают ядро в форме B-сплайна третьего порядка — самое гладкое приближение к Gaussian среди наших методов. Платим за это удвоением времени относительно одиночного Stack Blur, но всё ещё быстрее тройного Box.

Downscale Upscale — здесь проявляется природа метода. При σ=5 масштабирование минимально (Scale ≈ 1.25), и основное время тратится на Box Blur маленькой копии плюс финальный проход по полноразмерному изображению — 34.7 мс. При σ=50 изображение уменьшается в ~12 раз, размытие маленькой копии мгновенно, и всё время уходит на два StretchBlt — 14.9 мс. Парадокс: чем сильнее размытие, тем быстрее. Это единственный метод, который при больших σ обгоняет даже одиночный Stack Blur.

Gauss GDIP и Gauss D2D — системные реализации из первой части, добавлены для ориентира. GDI+ при σ=50 показывает 12.1 мс, что сопоставимо с нашим Downscale Upscale. Direct2D стабильно в районе 7 мс и GPU остаётся недосягаемым для CPU.

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

Если нужно…	Выбирайте
Максимальная скорость на CPU	Stack Blur (одиночный)
Лучшее качество приближения к Gaussian	Stack Blur ×2 или Box Blur Extended
Огромные σ (50+)	Downscale Upscale
Скорость любой ценой	Direct2D

Все методы из этой части подтвердили главное: время работы практически не зависит от величины σ. Десятикратный рост σ (с 5 до 50) не приводит к заметному росту времени — сравните синие и оранжевые столбцы. Все работают с одним TBitmap и ScanLine, без зависимостей. Разница между ними — в количестве проходов по памяти и в форме аппроксимирующего ядра.

В первой части мы дошли от 1593 мс (GBlur2, σ=50, 32-bit) до 329 мс (Gauss Fast, σ=50, 64-bit). Во второй — до 15 мс (Stack Blur) и 14.9 мс (Downscale Upscale). Это ×100 от наивной реализации, на чистом CPU, без GPU, без ассемблера, без внешних библиотек.

Приложение: Сравнение качества

Здесь представлены результаты размытия при сигме=5, левая картинка — эталон, полученный применением Direct2D-эффекта.

Сравнительная таблица

Картинки кликабельны. Чтобы легче было увидеть различия, можно кликнуть по одной из них, потом листать колесом мыши.

[свернуть]

Заключение

В первой части мы шли от наивного Gaussian к сепарабельному, выжимая скорость из прямой свёртки. Потолок был понятен: время линейно растёт с радиусом, и при σ=50 даже оптимизированная реализация тратила 329 мс.

Во второй части мы сломали эту линейную зависимость. Скользящая сумма, треугольное ядро Stack Blur, масштабирование — три разных способа уйти от O(r) к O(1). Результат: от 329 мс до 15 мс на том же изображении, том же CPU, без единой строки ассемблера.

Мы последовательно и подробно разобрали эти методы, каждый из которых опирается на предыдущие:

Box Blur — простейшая скользящая сумма. Один кирпичик, из которого складывается всё остальное.
Тройной Box — три кирпичика дают форму, близкую к колоколу. Центральная предельная теорема в действии.
Extended Box — бесплатная коррекция точности. Та же скорость, но σ попадает в цель.
Stack Blur — треугольное ядро за один проход. Два прохода по памяти вместо шести.
Downscale/Upscale — честный трюк: зачем размывать миллион пикселей, если можно размыть тысячу?

Каждый следующий метод не отменяет предыдущий — он добавляет инструмент. В реальном проекте выбор зависит от контекста: Stack Blur для превью в реальном времени, Extended Box когда важна точность, Downscale/Upscale для фоновых подложек с огромным σ.

Здесь намеренно не затронуты GPU, SIMD, многопоточность. Всё, что описано в этих двух частях, работает в одном потоке, с одним Bitmap и только средствами Delphi, без внешних библиотек. Это — базовый уровень. Потолок, от которого можно оттолкнуться дальше: оптимизировать или распараллелить проходы, корректно обработать альфу, применить другие стратегии для краевых пикселей.

Но это уже совсем другая история…

Скачать

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

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

Исполняемый файл x32 (zip) 1.55 Мб. Built in Delphi XE 7

Исполняемый файл x64 (zip) 2.08 Мб. Built in Delphi 13.0

5 1 голос

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

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

Старые

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

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

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