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

function gaussian( b: TBitmap; sigma: real ): TBitmap;
var
  mid, ksize        : integer;
  ker               : array of real;
  p                 : array of PArrRGB;
  b2                : TBitmap;

  procedure gauss( sigma: real );
  var
    i               : integer;
    sum, w          : real;
    norm            : real;
    ediv            : real;
  begin
    norm := 1 / ( sigma * sqrt( 2 * pi ) );
    ediv := 0.5 / ( sigma * sigma );
    ksize := 1 + 2 * round( 3 * sigma );
    if ksize < 3 then ksize := 3;
    if ( ( ksize and 1 ) <> 1 ) then ksize := ksize + 1;
    setlength( ker, ksize );
    mid := ksize div 2;
    sum := 0;
    for i := 0 to ksize - 1 do begin
      w := norm * exp( -( i - mid ) * ( i - mid ) * ediv );
      ker[i] := w;
      sum := sum + w;
    end;
    for i := 0 to high( ker ) do
      ker[i] := ker[i] / sum;
  end;

  function Hconv( b: TBitmap ): TBitmap;
  var
    j, i, k         : integer;
    corr            : real;
    sum             : real;
    pp              : ParrRGB;
  begin
    result := citra_clone( b );
    for j := 0 to high( p ) do
      p[j] := b.ScanLine[j];
    for j := 0 to b.Height - 1 do begin
      pp := result.ScanLine[j];
      for i := 0 to b.Width - 1 do begin
        sum := 0;
        corr := 0;
        for k := -mid to mid do begin
          if ( i + k < 0 ) or ( i + k >= b.Width ) then begin
            corr := corr + ker[mid + k];
            continue;
          end;
          sum := sum + ker[mid + k] * p[j][i + k].r;
        end;

        if corr > 0 then sum := sum * ( 1 / ( 1 - corr ) ); //        	memo1.lines.add(format('corr %.8f', [corr]));
        pp[i] := warna_create( clamp( round( sum ) ), clamp( round( sum ) ), clamp( round( sum ) ) );
      end;
    end;
  end;

  function Vconv( b: TBitmap ): TBitmap;
  var
    j, i, k         : integer;
    corr            : real;
    sum             : real;
    pp              : ParrRGB;
  begin
    result := citra_clone( b );
    for j := 0 to high( p ) do
      p[j] := b.ScanLine[j];
    for j := 0 to b.Width - 1 do begin
      for i := 0 to b.Height - 1 do begin
        pp := result.ScanLine[i];
        sum := 0;
        corr := 0;
        for k := -mid to mid do begin
          if ( i + k < 0 ) or ( i + k >= b.Height ) then begin
            corr := corr + ker[mid + k];
            continue;
          end;
          sum := sum + ker[mid + k] * p[i + k][j].r;
        end;
        if corr > 0 then sum := sum * ( 1 / ( 1 - corr ) );
        pp[j] := warna_create( clamp( round( sum ) ), clamp( round( sum ) ), clamp( round( sum ) ) );
      end;
    end;
  end;
begin
  setlength( p, b.Height );
  gauss( sigma );
  b2 := hconv( b );
  result := vconv( b2 );
  b2.Free;
end;

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