Размытие изображения

← →
anonimos (2009-05-02 21:55) [0]

Возникла необходимость в процедуре быстрого размытия изображения (найденные мною в инете были медленными или имели другие проблемы). Если не затруднит, подскажите как можно еще ускорить мною написанную процедуру (или где найти нормальную ^_^).

вот код, который на компе 2,2ГГц и при радиусе размытия 255 точек, размывает изображение 1000х700 за 22 сек (долговато)

... procedure Smooth(Dest, Source: TBitmap; S: Byte; X: Integer = 0; Y: Integer = 0; MinCasualOpacity: Byte = 0; MaxCasualOpacity: Byte = 0); implementation {----------------------------------------------------------------} function GetBounds(VLeft,VTop,VWidth,VHeight:Integer):TBounds; begin with Result do begin Left:=VLeft; Top:=VTop; Width:=VWidth; Height:=VHeight; end; end; {----------------------------------------------------------------} function BoundsAnd(B1,B2:TBounds):TBounds; var LStart,LLength:Integer; {------------------------------------------------------} procedure AndFuncAtOneDirection(Start1,Length1:Integer; Start2,Length2:Integer); begin LStart:=0; LLength:=0; if (Start1<=Start2)AND(Start1+Length1<=Start2+Length2) then begin LStart:=Start2; LLength:=Start1+Length1-Start2; end; if (Start1<=Start2)AND(Start1+Length1>Start2+Length2) then begin LStart:=Start2; LLength:=Length2; end; if (Start1>=Start2)AND(Start1+Length1>=Start2+Length2) then begin LStart:=Start1; LLength:=Start2+Length2-Start1; end; if (Start1>=Start2)AND(Start1+Length1<Start2+Length2) then begin LStart:=Start1; LLength:=Length1; end; if LLength<0 then LLength:=0; end; {------------------------------------------------------} begin AndFuncAtOneDirection(B1.Left,B1.Width,B2.Left,B2.Width); Result.Left:=LStart; Result.Width:=LLength; AndFuncAtOneDirection(B1.Top,B1.Height,B2.Top,B2.Height); Result.Top:=LStart; Result.Height:=LLength; end; {----------------------------------------------------------------} procedure Smooth(Dest, Source: TBitmap; S: Byte; X: Integer = 0; Y: Integer = 0; MinCasualOpacity: Byte = 0; MaxCasualOpacity: Byte = 0); var DestP, SourceP, Temp: PAC; DestRow, SourceRow: Integer; i1, i2, j1, j2, n1, n2, m1, m2: Integer; MinWidth, MinHeight: Integer; Bounds1, Bounds2: TBounds; {------------------------------------------------------} k: Integer; CAY: array [0..3] of ACD; M: ACD; MIndex: Integer; CAX: array [0..3] of Cardinal; StartX, FinX, StartY, FinY: Integer; DivisorY, DivisorX: Cardinal; {------------------------------------------------------} function Borders(Value, MinV, MaxV: Integer): Integer; begin Result := Value; if Result < MinV then Result := MinV; if Result > MaxV then Result := MaxV; end; {------------------------------------------------------} begin {------------------------------------------------------} Bounds1 := BoundsAnd(GetBounds(0, 0, Dest.Width, Dest.Height), GetBounds(X, Y, Source.Width, Source.Height)); Bounds2 := GetBounds(Bounds1.Left - X - S, Bounds1.Top - Y - S, Bounds1.Width + 2 * S, Bounds1.Height + 2 * S); Bounds2 := BoundsAnd(Bounds2, GetBounds(0, 0, Source.Width, Source.Height)); if (Bounds1.Width <= 0) OR (Bounds1.Height <= 0) then Exit; Dest.PixelFormat:=pf32bit; DestP:=Dest.ScanLine[Bounds1.Top]; DestRow:=Dest.Width*4; Source.PixelFormat:=pf32bit; SourceP:=Source.ScanLine[Bounds1.Top - Y]; SourceRow:=Source.Width*4; {------------------------------------------------------} for k := 0 to 3 do begin SetLength(CAY[k], Bounds2.Width); end; SetLength(M, (S * 2) + 1); i1 := (S * 2); for MIndex := 0 to i1 do begin if MIndex > S then M[MIndex] := M[2 * S - MIndex] else M[MIndex] := (255 * (MIndex + 1)) DIV (S + 1); end; {------------------------------------------------------} i1 := Bounds1.Top; while i1 < (Bounds1.Top + Bounds1.Height) do begin i2 := i1 - Y; {------------------------------------------------------} StartY := Borders(i2 - S, Bounds2.Top, Bounds2.Top + Bounds2.Height - 1); FinY := Borders(i2 + S, Bounds2.Top, Bounds2.Top + Bounds2.Height - 1); DivisorY := 0; m2 := 0; while m2 < (Bounds2.Width) do begin {------------------------------------------------------} j2 := m2 + Bounds2.Left; for k := 0 to 3 do begin CAY[k][m2] := 0; end; n2 := StartY; DivisorY := 0; while n2 <= FinY do begin Integer(Temp) := Integer(SourceP) - SourceRow * (n2 - i2); {------------------------------------------------------} MIndex := n2 - i2 + S; DivisorY := DivisorY + M[MIndex]; for k := 0 to 3 do begin CAY[k][m2] := CAY[k][m2] + RDWord(Temp^[j2]).B[k] * M[MIndex]; end; {------------------------------------------------------} n2 := n2 + 1; end; {------------------------------------------------------} m2 := m2 + 1; end; {------------------------------------------------------} j1 := Bounds1.Left; while j1 < (Bounds1.Left + Bounds1.Width) do begin j2 := j1 - X; {------------------------------------------------------} StartX := Borders(j2 - S, Bounds2.Left, Bounds2.Left + Bounds2.Width - 1) - Bounds2.Left; FinX := Borders(j2 + S, Bounds2.Left, Bounds2.Left + Bounds2.Width - 1) - Bounds2.Left; for k := 0 to 3 do begin CAX[k] := 0; end; DivisorX := 0; m2 := StartX; while m2 <= FinX do begin {------------------------------------------------------} MIndex := m2 - j2 + Bounds2.Left + S; DivisorX := DivisorX + M[MIndex]; for k := 0 to 3 do begin CAX[k] := CAX[k] + (CAY[k][m2] * M[MIndex]) DIV (DivisorY); end; {------------------------------------------------------} m2 := m2 + 1; end; if MaxCasualOpacity = MinCasualOpacity then begin for k := 0 to 3 do begin RDWord(DestP^[j1]).B[k] := (CAX[k]) DIV (DivisorX); end; end else begin m1:= RandomRange(MinCasualOpacity, MaxCasualOpacity); for k := 0 to 3 do begin RDWord(DestP^[j1]).B[k] := (RDWord(SourceP^[j2]).B[k] * Cardinal(m1) + (CAX[k]) DIV (DivisorX) * (255 - Cardinal(m1))) DIV (255); end; end; {------------------------------------------------------} j1 := j1 + 1; end; {------------------------------------------------------} i1 := i1 + 1; Integer(SourceP) := Integer(SourceP) - SourceRow; Integer(DestP) := Integer(DestP) - DestRow; end; {------------------------------------------------------} for k := 0 to 3 do begin SetLength(CAY[k], 0); end; SetLength(M, 0); end; end.

← →
Б (2009-05-03 09:49) [1]

См. DRKB "Об ускорении работы с графикой" и др.

← →
Sapersky (2009-05-03 11:02) [2]

А какое размытие требуется?
Насколько я понял, в M генерируется некий kernel, но какой конкретно... что-то вроде gauss, когда центральные пиксели имеют большее влияние, чем крайние? И так ли вообще нужен этот gauss, может, подойдёт mean (коэфф. для всех пикселей одинаков, просто сложить и разделить на кол-во). Вроде как нужно сильное размытие (судя по используемому радиусу), а mean как раз мажет сильнее, gauss даёт более тонкий эффект.
Mean с большим радиусом можно оптимизировать за счёт хранения общей суммы пикселей, для каждого нового пикселя из неё вычитается сумма крайнего левого столбца и добавляется сумма крайнего правого (для цикла по горизонтали), вместо того, чтобы каждый раз считать сумму всего квадрата 255*255.
Вообще при таком сильном размытии можно попробовать сначала уменьшить картинку раза в 2, размазать, потом увеличить обратно - может, разница в качестве и не будет заметна.

Конкретно по коду - не рекомендую использовать div и мелкие внутренние циклы for k := 0 to 3 do.

← →
anonimos (2009-05-03 12:58) [3]

> центральные пиксели имеют большее влияние, чем крайние?

да, М - матрица весовых коэффициентов, задается в начале и может быть любая.
Если использовать Mean то, например, черная линия, на белом фоне, станет широкой серой полосой без плавного перехода (хотя это действительно быстрее).

> не рекомендую использовать div и мелкие внутренние циклы
> for k := 0 to 3 do

спасибо за совет.

Кстати, можно ли подобными вычислениями "запрячь" ММХ ?

← →
Sapersky (2009-05-03 14:05) [4]

Кстати, можно ли подобными вычислениями "запрячь" ММХ ?

Можно, во всяком случае если ограничиться 32-битными битмапами.
Все цветовые компоненты будут обрабатываться параллельно, так что пресловутый цикл for k := 0 to 3 do выкинется естественным образом. Хотя потребуются доп. команды для преобразования данных - распаковка,упаковка... так что особых чудес от MMX тоже ждать не стоит - как показывает практика, он даёт ускорение раза в 2 максимум по сравнению с качественным (!) не-MMX кодом.

Ещё можно посмотреть готовые библиотеки, у Intel что-то такое было, кажется IPL называется. Наверняка там MMX уже задействован.
Хотя последние версии всех интеловских библиотек платные, бесплатно можно найти только что-то совсем старое (2000 г. или около того) и не лучшим образом оптимизированное под современные CPU. Хотя вот IJL1.5 (загрузка jpeg), например, до сих пор "живее всех живых".

← →
DVM © (2009-05-03 17:21) [5]

> Sapersky (03.05.09 14:05) [4]

> Хотя вот IJL1.5 (загрузка jpeg), например, до сих пор "живее
> всех живых".

Новым версиям она уже проигрывает, особенно на больших изображениях.

← →
Pavia © (2009-05-03 23:11) [6]

Размытие для произвольного радиуса лучше делать через FIR или рекурсивный фильтр. Твое ядро хорошо ложится.

Да и фильтр можно применять сначало по строчкам потом по столбцам. Что уменьшает умножение с N^2 до 2*N.

> Б (03.05.09 09:49) [1]
> См. DRKB "Об ускорении работы с графикой" и др.

Насколько я помню там был ужасный код.

← →
MBo © (2009-05-04 05:24) [7]

Быстрая свертка (convolution)

← →
Б (2009-05-04 07:36) [8]

> Насколько я помню там был ужасный код.

Но довольно быстрый.

← →
anonimos (2009-06-20 22:41) [9]

я немного чайник и не знаю что такое

> Быстрая свертка (convolution)

Попытка присобачить MMX для выполнения более простой задачи дала отрицательный результат (увеличение времени вычислений).

> См. DRKB "Об ускорении работы с графикой" и др.

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

Intel IPL (или IPP) - бесплатного не нашел (если честно почти ничего норм. не нашел)

Можно ли каким нибудь образом ускорить вычисления за счет применения видеокарты (прошу прощения если вопрос абсурден, никогда не работал с видеоадаптером)?

← →
Б (2009-06-21 08:03) [10]

> Можно ли каким нибудь образом ускорить вычисления за счет
> применения видеокарты (прошу прощения если вопрос абсурден,
> никогда не работал с видеоадаптером)?

Только D3D, DirectDraw и OpenGL.

← →
DVM © (2009-06-21 11:33) [11]

> Б (21.06.09 08:03) [10]

> Только D3D, DirectDraw и OpenGL.

Не только.

Есть еще CUDA, есть OpenGV (OpenCV c поддержкой аппаратного ускорения).

← →
Б (2009-06-21 15:12) [12]

Будем знать. ;)

← →
Sapersky (2009-06-22 13:34) [13]

Попытка присобачить MMX для выполнения более простой задачи дала отрицательный результат (увеличение времени вычислений).

Я сейчас тоже пробую прикрутить MMX к Gauss, результат пока отчасти положительный - на P3, как и положено, в 1.5-2 раза быстрее, на P4 - немного медленнее. Надо будет похимичить с выравниванием kernel.
Ещё выяснился такой момент. Фильтрация делается в два прохода, сначала по горизонтали, потом по вертикали, это обычная реализация подобных фильтров, в [6] ещё писали. На больших радиусах вертикальный фильтр начинает сильно тормозить, независимо от использования MMX - видимо, из-за того, что потребные для обработки скан-линии данные перестают влезать в кэш. Так что с определенного радиуса (зависит от ширины картинки и кэша) выгоднее вместо вертикального фильтра заранее развернуть на 90 град, отфильтровать горизонтальным фильтром, потом обратно.
С применением этой оптимизации gauss у меня работает значительно быстрее варианта [0], даже без MMX:
1024*768 R=255, Cel2.8 (24 bpp, noMMX) -> 2.4 c, P3-933 (MMX) -> 4.7 c
И ещё, глядя на результаты gauss, возникают сомнения по поводу необходимости радиуса 255. Gauss c радиусом 50 уже превращает картинку в набор невнятных цветных пятен, а 255 - вообще одно пятно доминирующего цвета. Может, в [0] просто фильтр неэффективный?

Другой вариант - покопать в сторону Frequency Filter:

One difference is that the computational cost of the spatial filter increases with the standard deviation (i.e. with the size of the filter mask), whereas the costs for a frequency filter are independent of the filter function. Hence, the spatial Gaussian filter is more appropriate for narrow lowpass filters, while the Butterworth filter is a better implementation for wide lowpass filters.
(http://homepages.inf.ed.ac.uk/rbf/HIPR2/freqfilt.htm)

Я так понял - скорость обработки в frequency domain не зависит от размера фильтра, хотя сами по себе FFT/обратное FFT (= БПФ, быстрое преобразование Фурье) довольно тяжёлые операции.

Только D3D, DirectDraw и OpenGL.

Точнее, шейдеры (D3D, OGL) или GPGPU-средства (OpenGV, CUDA, у ATI наверное тоже есть свой API).
C шейдерами получить большой радиус фильтра за один проход вряд ли возможно, у них ограниченное число чтений из текстуры, но можно вместо этого последовательно применить к картинке несколько фильтров с меньшим радиусом (можно попробовать этот метод и на софтвере, кстати). Математически это не вполне корректно, но визуально может быть и приемлемо.
Для получения представления о шейдерах можно посмотреть DX9 SDK (clootie.ru), пример PostProcess, там демонстрируются применение различных фильтров, в т.ч. размытия. Применяются они к 3D-сцене, но вместо 3D может быть и просто картинка.

Я всетаки за IIR там требуется всего 12умножений на пиксель, а то и 8. Не зависит от радиуса и быстрее чем FFT.
Пока непонял где косяк. Но преблизительно время работы 1024x1024x24 0.4 с на 2 ГГц

> есть OpenGV (OpenCV c поддержкой аппаратного ускорения).

Я хотел сказать GpuCV:

https://picoforge.int-evry.fr/cgi-bin/twiki/view/Gpucv/Web/WebHome

← →
anonimos (2009-06-22 22:11) [16]

Если я правильно понял что написал Sapersky, то возникают вопросы:
для БПФ вроде надо количество точек кратное 2^n. Тогда, если точек 1026 - придется взять матрицу 65536 у которой всё остальное равно 0. А это дополнительное увеличение кол-ва расчетов, плюс "реактивный" расход памяти, т.к. нужно временное хранилище.

Пока я только слегка ускорил алгоритм откинув обработку компоненты непрозрачности (opacity) - визуально некрасиво получалось.

Спасибо за вектор поиска (копать придется долго)

> для БПФ вроде надо количество точек кратное 2^n. Тогда,
> если точек 1026 - придется взять матрицу 65536 у которой
> всё остальное равно 0.

Достаточно 2048. А вообще есть быстры алгоритмы для призвольного n. В любом случии тут надо уже брать готовую библиотеку где код уже проаптимизирован.

Размытие изображения Найти похожие ветки