Particle Tracking: Least square fitting berbasis konvolusi

Ini adalah interpretasi bebas saya atas tutorial MDS di halaman granular material labs-nya. Idenya sederhana: menggunakan least square fitting berbasis konvolusi untuk menge-track pergerakan suatu partikel/image.

Least Square Fitting

Least square fitting adalah teknik kesesuaian atau kecocokan (godness fit) untuk mencari nilai terbaik yang paling mirip antara data observasi dengan data estimasi atau function fit. Jika data observasi dimisalkan dengan $y_i$ dan fungsi estimasi disimbolkan dengan $f(x_i)$ maka "least square fit" adalah jumlahan selisih nilai observasi dengan fitting function,

$$\sum\limits_{i=1}^n \left( {{y_i}-f ({x_i})}^2 \right)$$

Jika jika memiliki data dengan rentang error bar atau variansi, maka chi-square didefinisikan sebagai least square diatas dibagi dengan "error bar"nya, yakni, $\sigma$.

$$\chi^2=\sum\limits_{i=1}^n \left( \dfrac{{y_i}-f({x_i})}{\sigma} \right)^2$$

Nilai chi-square (χ2) berkisar antara 0 sampai dengan tak hingga (~), semakin kecil nilai chi-squared maka semakin mirip nilai observasi dengan nilai ekspektasi, dalam hal particle tracking, maka gambar bulatan particle akan semakin tipis sehingga semakin mudah dibedakan antar satu dengan yang lainnya. Pada tutorial kali ini, least square fitting yang digunakan adalah berbasis konvolusi yang akan dijabarkan pada tulisan dibawah.

Fungsi Partikel Ideal

Misalkan partikel yang ingin kita track posisinya memiliki fungsi ideal sebagai berikut, $$I_c(\vec{x})=\sum_{n=1}^{N} I_p(\vec{x}-\vec{x}_n(t);D,...),\;\;\;\;\;\;[1]$$ dengan $N$ adalah jumlah partikel dan $$I_p(\vec{x};D,...)$$ fungsi tersebut menggambarkan bentuk partikel ideal yang berada di tengah. Partikel ideal bergantung pada variabel diameter dan variabel lain yang digunakan dalam teknik pengolahan citra. Untuk demo ini kita menggunakan fungsi partikel ideal sbb,

$$I_p(\vec{x};D,w)= \dfrac{\bigl[1-tanh(\frac{|\vec{x}|-D/2}{w})\bigr]}{2}$$

Persamaan diatas diimplementasikan dengan matlab dalam script ipf.m, nilai bervariasi dari 1 di tengah partikel, 1/2 pada batas partikel dan 0 diluar partikel. Parameter W digunakan untuk menentukan seberapa tajam fungsi berubah (76% berubah pada rang +/- w). Pada prinsipnya nilai w akan berpengaruh seberapa focus sistem pencitraan yang dihasil melalui titik menyebar ke sekelilingnya.

Script untuk Partikel Ideal (Matlab)

x=-10:.1:10;
D=12;          % Diameter
w=1.3;         % Width
h=figure(2); set(h,'Position',[100 100 400 300],'Color',[1 1 1]);
plot(x,ipf(x,D,w),D/2*[1 1],[1/(1+exp(2)) 1/(1+exp(-2))],D/2*[1 1]-w,[0 1],D/2*[1 1]+w,[0 1],(w*[-1 1])+D/2,1/(1+exp(2))*[1 1],(w*[-1 1])+D/2,1/(1+exp(-2))*[1 1])
text(6,.9,'{\it 2w}','HorizontalAlignment','center');
text(6.5,.5,'76%');
xlabel('Position (Figure 1)');
ylabel('Intensity');

Hasilnya adalah sebagai berikut,

Plot image partikel ideal

D=12;          % Diameter
w=1.3;         % Width
ss=2*fix(D/2+4*w/2)-1;         % size of ideal particle image
os=(ss-1)/2;                   % (size-1)/2 of ideal particle image
[xx yy]=ndgrid(-os:os,-os:os);  % ideal particle image grid
r=abs(xx+i*yy);    % radial coordinate
h=figure(2); set(h,'Position',[100 100 400 400],'Color',[1 1 1]);
simage(ipf(r,D,w));
xlabel('Figure 2');

Hasilnya,

Fungsi Least-Squares Fit: Chi-squared

Posisi sebuah partikel dapat ditentukan dengan fungsi least-square-fitting. Kita definisikan chi-squared sebagai berikut, $$\chi^2(\vec{x}_0;D,w)=\int W(\vec{x}-\vec{x}_0)[I(\vec{x}) - I_p(\vec{x}-\vec{x}_0;D,w)]^2\; d\vec{x},\;\;\;\;\;[2]$$ dimana $W(x-x_0)$ adalah fungsi pembobot. Domain dari integrasi adalah pada luasan area dari citra eksperimental $I$. dimana domain $x_0$ adalah lebih besar. Secara umum, jika ukuran citra adalah Lx dan Ly dan ukuran $I_p$ adalah sx dan sy maka batas integrasi adalah [0 lx] dan [0 ly] namun range dari $x_0$ adalah [-sx Lx+sx] dan [-sy Ly+sy]. Ketika $x_0$ adalah posisi partikel, maka nilai chi-squarednya akan minimum. Karena itu, jika kita meminimalkan nilai chi-squared pada $x_0$ maka kita akan menemukan posisi partikel tersebut. Faktanya akan ada minimum chi-squared pada setiap partikel. Jika citra partikel adalah sama dengan persamaan [1] maka nilai chi-squared akan sama dengan nol pada saat $x_0=x_n$. Jadi, proses menemukan partikel sama dengan mencari nilai minimal dari chi-squared. Dengan konvolusi, proses ini akan menjadi mudah. Kembangkan persamaan [2] menjadi berikut, $$\chi^2(\vec{x}_0;D,w)=\int W(\vec{x}-\vec{x}_0)[I(\vec{x})^2 - 2I(\vec{x})I_p(\vec{x}- \vec{x}_0;D,w)+I_p(\vec{x}-\vec{x}_0;D,w)^2]\; d\vec{x},$$ $$\chi^2(\vec{x}_0;D,w)=I^2 \otimes W - 2I \otimes (W I_p) + \left\langle W I_p^2\right\rangle,\;\;\;\;\;[3]$$ dimana, $$f \otimes g = [f\otimes g](\vec{x}_0)=\int f(\vec{x})g(\vec{x}-\vec{x}_0)\;d\vec{x},\;\;\;\;\;[4]$$ $$\left\langle f \right\rangle=1 \otimes{f},$$ yang sayangnya bukan konstanta sederhana, namun fungsi $x_0$ karena domain integradi hanya pada area citra. Jika $W$ tidak simetris maka kita harus dengan hati-hati mengevaluasi persamaan [4] (berupa cross-correlation) menggunakan kovolusi, didefinisikan sebagai,
$$f*g = [f*g](\vec{x}_0)=\int f(\vec{x})g(\vec{x}_0-\vec{x})\;d\vec{x}=f(\vec{x}) \otimes g(-\vec{x}).$$

Cross-Correlation langsung

Beberapa pilihan yang mungkin untuk fungsi pembobot jika, $$W=1;\;\;\;\;\; \chi^2(\vec{x}_0;D,w)= \int I^2\;d\vec{x} - 2I \otimes I_p + \left\langle I_p^2 \right\rangle,$$ Dalam kasus ini, terminologi pertama ($W$) tidak bergantung pada $x_0$ dan terminologi kedua hanya begantung pada nilai $x_0$ didekat sudut citra, sehingga $$I \otimes I_p$$ akan menjadi maksimum untuk $x_0$ pada setiap posisi. Fakta bahwa cross-correlation adalah maksimum didekat pusat partikel menjadi titik awal untuk beberapa teknik particle tracking. Dari penjelasan sebelumnya, kita dapat melihat bahwa basis dari pendekatan ini adalah meminimalkan nilai chi-squared. Teknik ini berjalan baik pada sinyal citra dengan nilai signal-to-noise ratio yang baik dan terpisah secara cukup baik pula antar partikel satu dengan yang lain. Hal ini dikarenakan pilihan fungsi pembobot yang sangat lebar yang mengakibatkan fitting menjadi sensistif pada fakta bahwa ideal partikel memilik nilai nilai nol disekitarnya agar kontras pada partikel lain didekatnya. Pemilihan $W$ juga menghasilkan puncak yang lebar maksimal sehingga sensitif terhadap noise.

Meminimalkan chi-square secara penuh

Didapatkan bahwa pemilihan fungsi pembobot $W$ yang kecil menghasilkan hasil yang lebih baik. Biasanya digunakan, $$W=I_p;\;\;\;\;\; \chi^2(\vec{x}_0;D,w)=I^2 \otimes I_p - 2I \otimes I_p^2 + \left\langle I_p^3 \right\rangle.$$ dimana term chi-squared pertama menunjukkan area pada ukuran $I_p$ disekitar $x_0$ yang penting. Term kedua menunjukkan nilai konstan kecuali disekitar sudut citra. Disekitar citra dan semua term menjadi kecil karena fakta bahwa ada sedikit dan sedikit overlap diantara citra eksperimental dan citra partikel ideal. Jadi, untuk menormalisasi efek ini disekitar batas citra kita membagi dengan terminologi terakhir dan mendefiniskan chi-squared baru: $$\chi^2(\vec{x}_0;D,w)=\frac{I^2 \otimes I_p - 2I \otimes I_p^2}{ I_p^3 }+1.\;\;\;\;\;[5]$$ Persamaan [5] diimplementasikan pada file chiimg.m, yang juga menghitung nilai normalisasi dari persamaan [3]. $$\chi^2(\vec{x}_0;D,w)=\frac{I^2 \otimes W - 2I \otimes WI_p}{ WI_p^2 }+1.$$ Nilai chi-squared yang dinormalisasi masih diminimalkan pada posisi partikel, dan juga bisa digunakan untuk menemukan partikel yang memiliki pusat diluar citra.

Contoh:
Citra partikel pada kondisi padat dan Chi-squarednya

aw=imread('test.bmp');  % load image
[Nx Ny]=size(raw);       % image size

hi=250;  % hi and lo values come the image histogram
lo=10;   % hi/lo=typical pixel value outside/inside
ri=(hi-double(raw))/(hi-lo);  % normalize image

D=12;          % Diameter
w=1.3;         % Width
ss=2*fix(D/2+4*w/2)-1;         % size of ideal particle image
os=(ss-1)/2;                   % (size-1)/2 of ideal particle image
[xx yy]=ndgrid(-os:os,-os:os);  % ideal particle image grid
r=abs(xx+i*yy);    % radial coordinate

Cutoff=5;      % minimum peak intensity
MinSep=5;      % minimum separation between peaks
[Np px py]=findpeaks(1./chiimg(ri,ipf(r,D,w)),1,Cutoff,MinSep);  % find maxima

h=figure(2); set(h,'Position',[100 100 600 600],'Color',[1 1 1]);
simage([100*zerofill(ri,2*os,2*os) 100*zerofill(ri,2*os,2*os); 1./chiimg(ri,ipf(r,D,w)) 8*1./chiimg(ri,ipf(r,D,w))]); caxis([0 100])
hold on;
plot(py+2*os+Ny,px,'w.');
hold off;
xlabel('Figure 3.');

Hasil

Mencari pusat partikel

Sebuah citra monodispersi (berpencar dari titin tunggal) pada kira-kira 0.8 bagian area ditunjukkan pada gambar 3 kiri atas. Citra tersebut relatif berresolusi rendah dengan signal-to-noise ratio yang baik. Focus citra tersebut tidak begitu tajam, dan citra tersebut merepresentasikan contoh citra yang akan ditelusuri. Citra raw diambil dengan cahaya latar hitam sehingga partikel tampak gelap. Citra yang ditampilkan telah dinormalisasi dengan puncak terang dan gelap dari histogram citra sehingga memiliki nilai 1 didalam dan 0 diluar. Dibawahnya (kiri-bawah) ditunjukkan chi-squared dari persamaan [5] dengan $I_p$ dari gambar 2. Pada visualisasi invers dari chi-squared ditunjukkan bahwa pusat partikel tampak terang. Pada gambar bawah-kanan, citra denga 1/chi-squared ditunjukkan dengan delapan kali lebih terang untuk menunjukkan puncak yang lebih kecil. Fungsi Matlab findpeaks.m digunakan untuk mengekstrak semua peak (maxima) dalam $\dfrac{1}{\chi^2}$. Hal ini memungkinkan untuk menemukan pixel yang lebih besar dari 8 tetangga terdekat mereka, yang memiliki intensitas lebih lebih besar cari cut-off dan terpisah dari puncak (peak) lainnyca minimal sebesar "MinSep" pixel yang didefiniskan. Hasilnya adalah 122 puncak ditemukan. Titik-titik putih di gambar kanan atas menunjukkan partikel yang ditemukan. Dengan teknik ini, kita bisa menemukan semua partikel dan bahkan sebagian partikel pada kotak tersebut dengan resolusi yang relatif kecil dari 12 pixel/partikel dan fokus yang kurang. Hal ini memberikan pusat partikel kira-kira satu pixel akurasi.

Akurasi Sub-pixel dengan least square fitting

Kita juga dapat menggunakan least square fitting untuk mencari pusat partikel ke akurasi sub-pixel. Untuk melakukan hal tersebut, kita memodifikasi least-square sbb, $$\chi^2(\vec{x}_n;D,w)=\int [I(\vec{x}) - I_c(\vec{x},\vec{x}_n)]^2\; d\vec{x},\;\;\;\;\;[6]$$ dimana, $$I_c(\vec{x},\vec{x}_n)=\sum_n W_n(\vec{x})I_p(\vec{x}-\vec{x}_n;D,w),$$ dimana $W_n$ adalah fungsi dimana 1 didalam volume voronoi dari partikel $n$ dan 0 diluar. $W_n$ memasukkan partikel yang overlapp dalam perhitungannya. Kita mengimplementasikan $I_c$ menggunakan pgrid.m dan ipf.m . Berikut contoh script Matlab untuk perhitungan $I_c$ dan chi-squared.

[cxy over]=pgrid(px-os,py-os,Nx,Ny,[1 Nx 1 Ny],Np,2*os+3,0); % create local grid centered on each particle and overlap matrix
ci=ipf(cxy,D,w);  % create calculated image

di=ci-ri;                             % Calculate difference image
h=figure(2); set(h,'Position',[100 100 600 200],'Color',[1 1 1]);
simage([ri ci 4*di.^2]); caxis([0 1])   % Show the real, calculated, and chi-squared image
title('Real, Calculated, and \chi^2 images')
xlabel('figure 4.')

chi2=sum(di(:).^2);                   % Calculate Chi-Squared
fprintf('Chi-Squared=%f\n',chi2);
chi2o=chi2;di0=di;  % save for later

Hasilnya,

Meminimalkan Chi-squared

Untuk mencari pusat particle kita akan meminimalkan rata-rata pada persegi ketiga gambar diatas (figure 4), yang juga ekivalen dengan menyelesaikan persamaan berikut, $$\frac{\partial\chi^2(\vec{x}^*_n;D,w)}{\partial\vec{x}^*_n}=0,\;\;\;\;\;$$ untuk $x_{n}^{*}$. Karena kita memiliki akurasi pixel yang baik, maka kita dapat menyelesaikannya dengan menggunakan metode Newton. Satu langkah metode Newton diimplementasikan pada script cidp2.m. Fungsi ini mengembalikan nilai dalam $x_n$ yang dibutuhkan agar persamaan tersebut mendekati solusinya. Dengan memanggil fungsi tersebut beberapa kali maka solusinya akan dapat diketahui. Berikut adalah contohnya.

maxnr=5;       % maximum number of calls to Newton solver.
mindelchi2=1;  % minimum change in chi2 before stopping

nr=0;
delchi2=1e99;
while((abs(delchi2)>mindelchi2) && (nr < maxnr))
  [dpx dpy]=cidp2(cxy,over,di,Np,D,w);  %
  px=px+dpx;
  py=py+dpy;
  [cxy over]=pgrid(px-os,py-os,Nx,Ny,[1 Nx 1 Ny],Np,2*os+3,0); % create local grid centered on each particle and overlap matrix
  ci=ipf(cxy,D,w);  % create calculated image
  di=(ci-ri);
  delchi2=chi2-sum(di(:).^2); % Calculate change in Chi-Squared
  chi2=chi2-delchi2;          % Calculate Chi-Squared
  fprintf('.');
  nr=nr+1;
end
fprintf('\n');
chi2=sum(di(:).^2);                     % Calculate Chi-Squared
fprintf('Chi-Squared=%f\n',chi2);

h=figure(2); set(h,'Position',[100 100 600 400],'Color',[1 1 1]);
simage([di.^2 di0.^2]); caxis([0 .25]); % Show the chi-squared images
Title(sprintf('New \\chi^2=%6.2f           Original \\chi^2=%6.2f',chi2,chi2o),'fontsize',15)
xlabel('figure 5.')

Hasil,

Menemukan D dan W

Gambar 5 diatas menunjukkan peningkatan nilai chi-squared. Input dari user adalah diameter dan lebar dari citra partikel ideal. Dari perkiraan yang baik kita dapat menghitung nilai yang lebih baik dengan meminimalkan nilai chi-squared pada $D$ dan $w$, yakni menemukan $D*$ dan $w*$ dari persamaan berikut, $$\frac{\partial\chi^2(\vec{x}_0;D^*,w^*)}{\partial D^*}=0;\;\;\;\;\; \frac{\partial\chi^2(\vec{x}_0;D^*,w^*)}{\partial w^*}=0$$ Dengan menggunakan metode Newton lagi kita dapat menyelesaikan persamaan tersebut. Satu langkah metode Newton diimplementasikan pada script ciDw.m. Fungsi tersebut akan menghasilkan nilai $D$ dan $w$ agar lebih dekat dengan solusinya. Dengan memanggil fungsi tersebut beberapa kali maka kita akan mendapatkan solusi-nya.

Attachement:
Script Matlab di Github (TeknikOptik/code/imageProc)

Bersambung disini, Particle Tracking: Elongated particle and biological application.