Mean Shift Filtering

Saya pernah menulis tentang Mean filter, lalu dilanjutkan dengan Bilateral Filter [1]. Kedua filter ini merupakan batu loncatan untuk memahami Mean Shift Filter. Kalau di Mean filter yang dipertimbangkan adalah seluruh piksel yang berada dalam suatu area ketetanggaan maka di Bilateral filter yang dipertimbangkan diseleksi lagi menjadi nilai piksel yang ‘berdekatan’ di ruang fitur (feature space) dengan piksel yang menjadi acuan (biasanya di tengah). Kriteria untuk menentukan suatu nilai ‘berdekatan’ bisa jadi berupa nilai keabuan (grayscale), jarak antara warna menggunakan metric euclid(atau metric lain seperti manhattan block distance) pada satu ruang warna (RGB, YUV, YCbCr, dll).

Untuk memahami mean-shift filter, perlu menambahkan konsep bahwa nilai rata-rata yang dipakai bisa jadi tidak terbatas di area ketetanggaan itu saja melainkan bisa saja ada di tempat lain. Oleh sebab itu, nilai rata-ratanya perlu ‘jalan-jalan’ atau bergeser (shift) sampai menemukan nilai rata-rata yang tepat. Kalau dalam mean filter dan bilateral filter yang dihitung hanya rata-rata nilai piksel, maka pada mean-shift filter koordinat piksel yang masuk kriteria seleksi juga dihitung nilai rata-ratanya. Dengan demikian didapat nilai rata-rata koordinat piksel sebagai koordinat baru. Nilai koordinat baru ini kemudian dijadikan titik acuan untuk menghitung nilai rata-rata di ruang fitur (intensitas piksel, warna, atau gradien). Jika titik koordinat baru ternyata tidak berubah maka nilai rata-rata di ruang fitur piksel digunakan untuk mengganti nilai piksel di seluruh koordinat yang pernah menjadi koordinat acuan sampai di koordinat semula.

Berbasis kode mean-shift filter di plugin ImageJ, saya buat versi delphi untuk mean-shift filtering. Kode aslinya dimodifikasi supaya lebih efisien dengan mengasumsikan koordinat yang pernah dikunjungi diganti dengan nilai akhir maka koordinat yang pernah dikunjungi tidak perlu lagi diproses. Efisiensi dilakukan dengan menggunakan citra biner yang berfungsi sebagai mask atau flag, yaitu penanda suatu piksel sudah pernah diproses (akibat penggeseran atau belum).

procedure TForm1.MeanShift1Click(Sender: TObject);
var
  rad               : array of TPoint;
  i, j, k           : integer;
  FeatureSpaceRadius: real;
  ImageSpaceRadius  : Integer;
  fsr2, isr2        : real;
  yiqimage          : array of array of array of real;
  p                 : array of PArrRGB;//lihat tulisan sebelumnya tentang mengakses piksel
  b                 : TBitmap;
  delta             : real;
  iter, maxiter     : integer;
  cur, old          : TPoint;
  sumx, sumy, sumf1, sumf2, sumf3: real;
  invnum            : real;
  xx, yy            : integer;
  f1, f2, f3        : real;
  of1, of2, of3     : real;
  df1, df2, df3     : real;
  dx, dy            : integer;
  bs                : TBitmap;
  mask              : array of array of boolean;
  path              : array of TPoint;

  function is_ok(x, y: integer): boolean;
  begin
    result := (x >= 0) and (y >= 0) and (x < b.Width) and (y < b.Height);
  end;
begin
  //transform image to feature space (YIQ color space)
  b := Image1.Picture.Bitmap;
  setlength(yiqimage, b.Height);
  setlength(mask, b.Height);
  setlength(p, b.Height);
  for j := 0 to high(yiqimage) do begin
    p[j] := b.ScanLine[j];
    setlength(yiqimage[j], b.Width);
    setlength(mask[j], b.Width);
    for i := 0 to b.Width - 1 do with p[j][i] do begin
        mask[j][i] := True;
        setlength(yiqimage[j][i], 3);
        yiqimage[j][i][0] := 0.299 * r + 0.587 * g + 0.114 * b; ;
        yiqimage[j][i][1] := 0.5957 * r - 0.2744 * g - 0.3212 * b; ;
        yiqimage[j][i][2] := 0.2114 * r - 0.5226 * g + 0.3111 * b;
      end;
  end;

  //precompute sampling displacement
  ImageSpaceRadius := 20;
  isr2 := ImageSpaceRadius * ImageSpaceRadius;
  for j := -ImageSpaceRadius to ImageSpaceRadius do begin
    for i := -ImageSpaceRadius to ImageSpaceRadius do begin
      if i * i + j * j <= isr2 then begin
        SetLength(rad, length(rad) + 1);
        rad[high(rad)] := Point(i, j);
      end;
    end;
  end;

  FeatureSpaceRadius := 40;
  fsr2 := FeatureSpaceRadius * FeatureSpaceRadius;
  maxiter := 100;
  for j := 0 to b.Height - 1 do begin
    for i := 0 to b.Width - 1 do begin
      if not mask[j][i] then continue;
      cur.X := i;
      cur.Y := j;
      iter := 0;
      f1 := yiqimage[j][i][0];
      f2 := yiqimage[j][i][1];
      f3 := yiqimage[j][i][2];
      setlength(path, 0);
      repeat
        old := cur;
        of1 := f1;
        of2 := f2;
        of3 := f3;

        invnum := 0;
        sumx := 0;
        sumy := 0;
        sumf1 := 0;
        sumf2 := 0;
        sumf3 := 0;
        for k := 0 to high(rad) do begin
          xx := cur.X + rad[k].X;
          yy := cur.Y + rad[k].Y;
          if is_ok(xx, yy) then begin
            df1 := f1 - yiqimage[yy][xx][0];
            df2 := f2 - yiqimage[yy][xx][1];
            df3 := f3 - yiqimage[yy][xx][2];
            if (df1 * df1 + df2 * df2 + df3 * df3 <= FSR2) then begin
              invnum := invnum + 1;
              sumx := sumx + xx;
              sumy := sumy + yy;
              sumf1 := sumf1 + yiqimage[yy][xx][0];
              sumf2 := sumf2 + yiqimage[yy][xx][1];
              sumf3 := sumf3 + yiqimage[yy][xx][2];
            end;
          end;
        end;
        if invnum < 1 then begin
          break;
        end;
        invnum := 1.0 / invnum;
        f1 := sumf1 * invnum;
        f2 := sumf2 * invnum;
        f3 := sumf3 * invnum;
        cur.X := round(sumx * invnum);
        cur.Y := round(sumy * invnum);
        if not mask[cur.Y][cur.X] then begin
          break;
        end;
        setlength(path, length(path) + 1);
        path[high(path)] := cur;
        dx := cur.X - old.X;
        dy := cur.Y - old.Y;
        df1 := f1 - of1;
        df2 := f2 - of2;
        df3 := f3 - of3;
        delta := dx * dx + dy * dy + df1 * df1 + df2 * df2 + df3 * df3;
        inc(iter);
      until (delta < 3) or (iter >= maxiter);
      p[j][i].r := round(f1 + 0.9563 * f2 + 0.6210 * f3);
      p[j][i].g := round(f1 - 0.2721 * f2 - 0.6473 * f3);
      p[j][i].b := round(f1 - 1.1070 * f2 + 1.7046 * f3);
      mask[j][i] := false;
      for k := 0 to high(path) do begin
        mask[j][i] := false;
        p[path[k].Y][path[k].X] := p[j][i];
      end;
    end;
  end;
  Image1.Refresh;
end;

Kode di atas menggunakan ruang warna YIQ untuk menghitung kedekatan antar piksel. Fungsi jarak yang digunakan adalah jarak Euclid. Selain fungsi jarak dan ruang fitur, filter mean-shift memiliki 2 parameter umum yaitu Radius dalam domain spasial (ImageSpaceRadius) yang dalam filter-filter umumnya menunjuk ke area ketetanggaan dan Radius dalam domain fitur (FeatureSpaceRadius) yaitu batas kedekatan suatu piksel pada ruang fitur. Contoh pemrosesan dengan filter mean-shift ditampilkan di gambar berikut.

2 comments

  1. ria · Desember 1, 2011

    terima kasih🙂.
    artikel anda sangat membantu.

    share ebook tentang mean shift dong.🙂
    kalo bisa yang bahasa indonesia.

    thanks y.

  2. Hend LuvNena · April 1, 2012

    bro, bisa jelasin secara sederhana gk bagaimana algoritma koding nya buat mean shift ?😦

    aku muter2 dari tadi akhirnya nemu blog punya bro, tapi masih kebingungan karena aku pake bahasa c# bro😦

    mohon bantuannya yah bro🙂

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s