Implementasi Gaussian Blur

Sebagai pengguna perangkat lunak manipulasi citra/foto teknik gaussian blur merupakan salah satu yang sering digunakan. Belum lama, ada yang menanyakan tentang dekomposisi wavelet (LL,LH,HL,HH) dan saya pikir perlu ada tulisan ini sebagai pengantar tulisan selanjutnya.

Gaussian blur sebetulnya merupakan konvolusi citra dengan fungsi gaussian. Mengapa perlu nama khusus? hal ini karena gaussian memiliki beberapa sifat yang cukup unik. Pertama karena hasil transformasi fourier fungsi gaussian ternyata menghasilkan gaussian juga. Kedua yang menarik dari fungsi gaussian adalah untuk implementasi fungsi gaussian 2D. Fungsi gaussian 2D secara matematis memiliki sifat dapat dipisahkan (separable, silakan cari sendiri bagi yang minat utak-atik simbol matematis :D ). Sifat inilah yang sangat membantu dalam mengimplementasi proses pengaburan (blurring).

Sebagaimana telah ditulis sebelumnya, dalam operasi perata-rataan tetangga (konvolusi) digunakan jendela atau area yang menyatakan hubungan ketetanggaan yang dilibatkan dalam operasi. Semakin besar ukuran jendela tersebut maka gambar akan semakin kabur dan juga waktu yang diperlukan untuk melakukan operasi tersebut. Dengan adanya operasi yang bersifat dapat dipisahkan, proses konvolusi bisa dilakukan dua kali namun hanya menggunakan jendela berdimensi satu. proses konvolusi pertama kali dilakukan untuk tiap baris dan dilanjutkan dengan konvolusi untuk tiap kolom sebagai data berdimensi satu juga. Akibat partisi ini operasi konvolusi bisa mengeksploitasi paralelisme baik menggunakan thread ataupun prosesor grafik.

Kode dibawah sudah didekomposisi konvolusi untuk elemen baris dan kolom sebagai upafungsi tersendiri sehingga dapat memudahkan untuk memproses sebagian saja (seperti dekomposisi LL,LH,HL,HH di atas).
view source
print?
01 function gaussian( b: TBitmap; sigma: real ): TBitmap;
02 var
03 mid, ksize : integer;
04 ker : array of real;
05 p : array of PArrRGB;
06 b2 : TBitmap;
07
08 procedure gauss( sigma: real );
09 var
10 i : integer;
11 sum, w : real;
12 norm : real;
13 ediv : real;
14 begin
15 norm := 1 / ( sigma * sqrt( 2 * pi ) );
16 ediv := 0.5 / ( sigma * sigma );
17 ksize := 1 + 2 * round( 3 * sigma );
18 if ksize < 3 then ksize := 3;
19 if ( ( ksize and 1 ) <> 1 ) then ksize := ksize + 1;
20 setlength( ker, ksize );
21 mid := ksize div 2;
22 sum := 0;
23 for i := 0 to ksize - 1 do begin
24 w := norm * exp( -( i - mid ) * ( i - mid ) * ediv );
25 ker[i] := w;
26 sum := sum + w;
27 end;
28 for i := 0 to high( ker ) do
29 ker[i] := ker[i] / sum;
30 end;
31
32 function Hconv( b: TBitmap ): TBitmap;
33 var
34 j, i, k : integer;
35 corr : real;
36 sum : real;
37 pp : ParrRGB;
38 begin
39 result := citra_clone( b );
40 for j := 0 to high( p ) do
41 p[j] := b.ScanLine[j];
42 for j := 0 to b.Height - 1 do begin
43 pp := result.ScanLine[j];
44 for i := 0 to b.Width - 1 do begin
45 sum := 0;
46 corr := 0;
47 for k := -mid to mid do begin
48 if ( i + k < 0 ) or ( i + k >= b.Width ) then begin
49 corr := corr + ker[mid + k];
50 continue;
51 end;
52 sum := sum + ker[mid + k] * p[j][i + k].r;
53 end;
54
55 if corr > 0 then sum := sum * ( 1 / ( 1 - corr ) ); // memo1.lines.add(format('corr %.8f', [corr]));
56 pp[i] := warna_create( clamp( round( sum ) ), clamp( round( sum ) ), clamp( round( sum ) ) );
57 end;
58 end;
59 end;
60
61 function Vconv( b: TBitmap ): TBitmap;
62 var
63 j, i, k : integer;
64 corr : real;
65 sum : real;
66 pp : ParrRGB;
67 begin
68 result := citra_clone( b );
69 for j := 0 to high( p ) do
70 p[j] := b.ScanLine[j];
71 for j := 0 to b.Width - 1 do begin
72 for i := 0 to b.Height - 1 do begin
73 pp := result.ScanLine[i];
74 sum := 0;
75 corr := 0;
76 for k := -mid to mid do begin
77 if ( i + k < 0 ) or ( i + k >= b.Height ) then begin
78 corr := corr + ker[mid + k];
79 continue;
80 end;
81 sum := sum + ker[mid + k] * p[i + k][j].r;
82 end;
83 if corr > 0 then sum := sum * ( 1 / ( 1 - corr ) );
84 pp[j] := warna_create( clamp( round( sum ) ), clamp( round( sum ) ), clamp( round( sum ) ) );
85 end;
86 end;
87 end;
88 begin
89 setlength( p, b.Height );
90 gauss( sigma );
91 b2 := hconv( b );
92 result := vconv( b2 );
93 b2.Free;
94 end;

coming up next.. edge sharpening/edge-preserving smoothing berbasis dekomposisi LH (lowpass-highpass filtering) dan canny edge detection ;)