k-means clustering using custom distance measuring method

ကြိုက်တဲ့ အကွာအဝေး တိုင်းတဲ့ နည်းသုံးပြီး k-means clustering လုပ်ပေးတဲ့ MATLAB ဖန်ရှင် လေး တစ်ခု ရေးကြည့် ထားပါတယ်။ ဥပမာ histograms တွေကို ခွဲခြားဖို့ဆိုရင် chi-square distance ကို သုံးတာက ပိုကောင်း ချင် ကောင်းနိုင်ပါတယ်။ အောက်က နမူနာမှာ chi-square distance သုံးထားတာကို ပြထားပေမယ့် အဲဒီ distance ဖန်ရှင်နေရာမှာ ကိုယ်ကြိုက်တဲ့ နည်း ကို အစားထိုးပြီး သုံးလို့ရပါတယ်။

---

%k-means test program
X = [randn(100,2)+ones(100,2);...
     randn(100,2)-ones(100,2)];

[idx,ctrs] = KMeansCustom(X,2);
%[idx,ctrs] = kmeans(X,2);

plot(X(idx==1,1),X(idx==1,2),'r.','MarkerSize',12)
hold on
plot(X(idx==2,1),X(idx==2,2),'b.','MarkerSize',12)
plot(ctrs(:,1),ctrs(:,2),'kx',...
     'MarkerSize',12,'LineWidth',2)
legend('Cluster 1','Cluster 2','Centroids',...
       'Location','NW')

---

function [Idx,C]=KMeansCustom(X,k)
%KMeansCustom partitions the points in the n-by-d data matrix X into k clusters.
%[Idx,C]= KMeansCustom(X,k) returns 
%n-by-1 vector IDX containing the cluster indices of each point and 
%k-by-d matrix C containing the k cluster centroid locations.
%For n sample points with d dimensions in each point, X has n rows and d columns.
%File name: KMeanCustom.m
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Define maximum number of iterations
MaxIter=500;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,d]=size(X);
k=round(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step1 :arbitrarily choose k samples as the initial cluster centers
p=randperm(n);
Mu=X(p(1:k),:);
D=zeros(k,d);
for t=1:MaxIter
 %step2:distribute the samples X  to the clusters 
 for j=1:k
        for i=1:n
            D(j,i)=ChiDist(X(i,:),Mu(j,:));%Use custom distance
        end
 end
 [ValMin,IndexMin]=min(D);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step 3: update the cluster centers
    OldMu=Mu;
 for i=1:k
        Mu(i,:)=mean(X(IndexMin==i,:));
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step4 :check convergence
 if sum(sum(abs(OldMu-Mu))) == 0 %< 1e-9
        break
 end
end
Idx=IndexMin';
C=Mu;

---

function d=ChiDist(v1,v2)
    dv=(v1-v2).^2;
    sv=abs(v1)+abs(v2);
    %------------------------------------------------------
    %eliminate zero denominator
    sv(sv==0)=1e-9;
    %------------------------------------------------------
    d=sum(dv./sv)./2;    
end

Amplitude spectrum of a signal using Fourier transform

signal တစ်ခုရဲ့ amplitude spectrum ကို Fourier transform သုံးပြီး ရှာတဲ့ MATLAB ဖန်ရှင်လေး တစ်ခု ရေးကြည့်ထားပါတယ်။ အဲဒီ ဖန်ရှင်ကို နမူနာ သုံးတဲ့ ပုံစံ ရယ်၊ ဖန်ရှင် ရေးထားတဲ့ ကုဒ်တွေ ရယ်ကို အောက်မှာ ဖော်ပြထားပါတယ်။ MATLAB မှာ ပါတဲ့ ဖန်ရှင် fft ကို သုံးတဲ့ ဟာ ကို လည်း နောက်ထပ် နမူနာ တစ်ခု ပြထားပါတယ်။ ဒီဟာလေးတွေက amplitude spectrum တို့၊ power spectrum တို့ ရှာချင်တဲ့ သူတွေ အတွက် အသုံးဝင်မယ်လို့ မျှော်လင့်ပါတယ်။

---

%Example MATLAB code to test function FX
Fs=100;
t=(0:1/Fs:1-1/Fs)';
x=2*cos(2*pi*5*t)+sin(2*pi*10*t);
[f a]=FX(x,Fs);
plot(f,a);

---

function [f a]=FX(x,fs)
%Calculates amplitude spectrum of a signal x using Fourier transform
%[f a]=FX(x,fs)
%Input: x=signal, fs = sampling rate
%Output: f = frequency axis, a = amplitude spectrum
%File name: FX.m
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/

dT = 1/fs; 
N=length(x);
dF=1/(N*dT);
NF=fs/2;%Nyquist Freq
%You can always limit freq range for faster performance
%NF=20;
t=(0:dT:(N-1)*dT)';
f=(0:dF:NF)';
a=f;%initialize a
a(1)=mean(x);
for i=2 : length(a)
    b=(2*mean(x.*cos(2*pi*dF*(i-1)*t)));
    c=(2*mean(x.*sin(2*pi*dF*(i-1)*t)));
    a(i)=sqrt(b^2+c^2);
end

---

%MATLAB program to calculates amplitude spectrum of a signal x using fft
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/

Fs=100;
t=(0:1/Fs:1-1/Fs)';
x=2*cos(2*pi*5*t)+sin(2*pi*10*t);
L=length(x);

A=2*abs(fft(x)/L);
A=A(1:Fs/2+1);
f = Fs/2*linspace(0,1,Fs/2+1);
plot(f,A);

---

Changing sampling rate using quadratic regression

signal တစ်ခုရဲ့ sampling rate ကို 333Hz ကနေ 5000Hz ကို ပြောင်းဖို့ အတွက် လုပ်ကြည့်ထားပါတယ်။ ပိုမို ချောမွေ့တဲ့ ရလာဒ် ရချင်တာမို့ ပုံမှန် zero order hold (ZOH) ကို သုံးမယ့်အစား second order hold ကို သုံးကြည့်ထားတာပါ။ ဒီမှာ sampling rate ကို တင်ဖို့လုပ်ထားပေမယ့်၊ ဒီနည်းကို down sampling လုပ်ဖို့လည်း သုံးလို့ရပါတယ်။ အခြား ဖြစ်နိုင်ချေတွေကတော့ noise ကို phase delay မရှိပဲ filter လုပ်တဲ့ဟာမျိုးတွေ အတွက်လည်း သုံးလို့ ရနိုင်ပါတယ်။ real-time application တွေအတွက်ပေါ့။
Approach 1: Quadratic regression
Quadratic ဖန်ရှင်တစ်ခုကို အောက်ပါအတိုင်း သတ်မှတ်နိုင်ပါတယ်။

$$f=w_0+w_1 x + w_2 x^2$$

အနည်းဆုံး ဖြစ်အောင် လုပ်မယ့် cost function ကို အောက်ကလို သတ်မှတ်မှာပါ။

$$J(w)=\frac{1}{2}\sum_{i=1}^{n}(y_i-f_i)^2$$
$$J(w)=\frac{1}{2}\sum_{i=1}^{n}(y_i-w_0+w_1 x_i + w_2 x_i^2)^2$$

ဒါဆို Optimal weights တွေကို differentiate လုပ်ပြီးသုည ညီလိုက်ပြီး ရှာနိုင်ပါတယ်။

$$\frac{\partial J(w)}{\partial w_0}=0$$
$$-\sum_{i=1}^{n}(y_i-w_0+w_1 x_i + w_2 x_i^2)=0$$
$$w_0\sum_{i=1}^{n}1+w_1\sum_{i=1}^{n}x_i+w_2\sum_{i=1}^{n}x_i^2=\sum_{i=1}^{n}y_i$$

w₁ နဲ့ w₂ ကိုလည်း အဲဒီလိုပဲ ရှာလိုက်မယ် ဆိုရင် အောက်အတိုင်း ထပ်ရလာပါမယ်။

$$w_0\sum_{i=1}^{n}x_i+w_1\sum_{i=1}^{n}x_i^2+w_2\sum_{i=1}^{n}x_i^3=\sum_{i=1}^{n}y_i x_i$$
$$w_0\sum_{i=1}^{n}x_i^2+w_1\sum_{i=1}^{n}x_i^3+w_2\sum_{i=1}^{n}x_i^4=\sum_{i=1}^{n}y_i x_i^2$$

အဲဒီ equation သုံးခုကို မေ့ထရစ် ပုံစံ နဲ့ အောက်ပါအတိုင်း ရေးလို့ရပါတယ်။

$$\begin{bmatrix} \sum_{i=1}^{n}1 & \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 \\ \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 & \sum_{i=1}^{n}x_i^3 \\ \sum_{i=1}^{n}x_i^2 & \sum_{i=1}^{n}x_i^3 & \sum_{i=1}^{n}x_i^4 \end{bmatrix} \begin{bmatrix} w_0 \\ w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{n}y_i \\ \sum_{i=1}^{n}y_i x_i \\ \sum_{i=1}^{n}y_i x_i^2 \end{bmatrix}$$
$$\mathbf{A}\mathbf{W}=\mathbf{B}$$
$$\mathbf{W}=\mathbf{A}^{-1}\mathbf{B}$$

Approach 2: Three equations
တကယ်တော့ ဒီကိစ္စမှာ နောက်ဆုံး အမှတ် ၃ ခုကိုပဲ သုံးလိုက်ရင်ရတဲ့ အတွက်၊ ညီမျှခြင်း သုံးကြောင်း ကို အောက်က အတိုင်း တစ်ခါတည်း ချရေးလို့လည်း ရပါတယ်။

$$y_1=w_0+w_1 x_1 + w_2 x_1^2$$
$$y_2=w_0+w_1 x_2 + w_2 x_2^2$$
$$y_3=w_0+w_1 x_3 + w_2 x_3^2$$
$$\begin{bmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ 1 & x_3 & x_3^2 \end{bmatrix} \begin{bmatrix} w_0 \\ w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}$$
$$\mathbf{A}\mathbf{W}=\mathbf{B}$$
$$\mathbf{W}=\mathbf{A}^{-1}\mathbf{B}$$

Approach 3: Two equations အကယ်၍ x₁=-1, x₂=0, နဲ့ x₃=1 လို့ သတ်မှတ်လို့ ရတယ်ဆိုရင် y₂ က w₀ နဲ့ ညီသွားပါမယ်။ နောက်ကျန်တဲ့ အမှတ် ၂ ခုကနေ ညီမျှခြင်း နှစ်ကြောင်းပဲ ရှင်းဖို့ လိုပါတော့တယ်။

$$y_1=y_2-w_1+ w_2$$
$$y_3=y_2+w_1+ w_2$$
$$\begin{bmatrix} -1 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} y_1-y_2 \\ y_3-y2 \end{bmatrix}$$
$$\begin{bmatrix} w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} -0.5 & 0.5 \\ 0.5 & 0.5 \end{bmatrix} \begin{bmatrix} y_1-y_2 \\ y_3-y2 \end{bmatrix}$$

နောက် အောက်ကအတိုင်းရလာပါမယ်။

$$w_0=y_2$$
$$w_1=-0.5(y_1-y_2)+0.5(y_3-y2)$$
$$w_2=0.5(y_1-y_2)+0.5(y_3-y2)$$

အထက်ကနည်းတွေကို MatLab နဲ့စမ်းကြည့်ထားတာကို အောက်က ကုဒ်တွေမှာ တွေ့နိုင်ပါတယ်။

%-------------------------------------------------------------------------
clc;
close all;
clear all;
%-------------------------------------------------------------------------
% y= w0 + w1*x + w2* x^2;
%-------------------------------------------------------------------------
%Got x and y
x=[-1 0 1]';
Wo=[4 3 2]';
y=Wo(1)+Wo(2)*x+Wo(3).*x.*x;

%-------------------------------------------------------------------------
%Approach 1
%Polynomial regression of order 2
%For n=3
S1=3;
Sx=sum(x);
Sx2=sum(x.*x);
Sx3=sum(x.*x.*x);
Sx4=sum(x.*x.*x.*x);
Sy=sum(y);
Syx=sum(y.*x);
Syx2=sum(y.*x.*x);

P=[S1 Sx Sx2; Sx Sx2 Sx3; Sx2 Sx3 Sx4];
B=[Sy Syx Syx2]';
%P1=P^(-1);
W1=P\B

%-------------------------------------------------------------------------
%Approach 2
%Linear equations
A=[1 x(1) x(1)*x(1);1 x(2) x(2)*x(2); 1 x(3) x(3)*x(3)];
W2=A\y
%-------------------------------------------------------------------------
%Approach 3
%Only 2 linear equations
w0=y(2);
w1=-0.5*( y(1)- y(2))+0.5*( y(3)- y(2));
w2=0.5*( y(1)- y(2))+0.5*(y(3)- y(2));
W3=[w0 w1 w2]'
%-------------------------------------------------------------------------

အောက်ကပုံမှာ ဒီနည်းကို သုံးလို့ရတဲ့ ရလဒ် (အပြာ) နဲ့ ပုံမှန် zero order hold သုံးရင် ရမယ့် ရလဒ် (အနက်) တို့ကို ယှဉ်ပြထားပါတယ်။ ဒီနည်းက ပိုပြီး ညက်ညော ချောမွေ့တဲ့ ရလဒ်ကို ရပေမယ့် one sample delay ရှိသွားတာကိုတော့ သတိထားဖို့ လိုပါတယ်။

ဒီနည်းလမ်းကို LabVIEW မှာ C code သုံးပြီး လုပ်ထားတာကို အောက်ကပုံမှာ တွေ့နိုင်ပါတယ်။

ပထမ နှစ်ခုမှာ 3x3 matrix ကို inverse ရှာဖို့ လိုတဲ့အတွက် inverse ရှာတဲ့ C program လေးတစ်ခုကို ရေးပြီး စမ်းကြည့်ထားပါတယ်။

#include
#include
main()
{
    float M[3][3]={{3,0,2},{0,2,0},{2,0,2}}; //initialize a 3x3 matrix
    float N[3][3]={{0,0,0},{0,0,0},{0,0,0}}; //allocate for inverse
    int i,j;
    float d;
    //-------------------------------------------------------------------------
    N[0][0]=(M[1][1]*M[2][2]-M[2][1]*M[1][2]);
    N[1][0]=-(M[1][0]*M[2][2]-M[2][0]*M[1][2]);
    N[2][0]=(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
    d=M[0][0]*N[0][0]+M[0][1]*N[1][0]+M[0][2]*N[2][0];
    N[0][0]/=d;
    N[1][0]/=d;
    N[2][0]/=d;
    N[0][1]=-(M[0][1]*M[2][2]-M[0][2]*M[2][1])/d;
    N[1][1]=(M[0][0]*M[2][2]-M[0][2]*M[2][0])/d;
    N[2][1]=-(M[0][0]*M[2][1]-M[0][1]*M[2][0])/d;
    N[0][2]=(M[0][1]*M[1][2]-M[0][2]*M[1][1])/d;
    N[1][2]=-(M[0][0]*M[1][2]-M[0][2]*M[1][0])/d;
    N[2][2]=(M[0][0]*M[1][1]-M[0][1]*M[1][0])/d;
    //-------------------------------------------------------------------------
    //print 3x3 matrix
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++) printf("%3.4f ",N[i][j]);
        printf("\n");
    }
    getch();
    return 0;
}

Rotations in 3D space using Euler angles

ကိရိယာ တစ်ခုရဲ့ body reference frame နဲ့ world reference frame ရဲ့ ဆက်သွယ်မှု ကို ရှာဖို့အတွက် Euler angle sequence တစ်ခုကို သုံးပြီး သရီးဒီ စပေ့ ထဲမှာ rotation လုပ်တဲ့ အကြောင်း ပြောချင်ပါတယ်။ အွိုင်လာ တင်ပြခဲ့တဲ့ သီအိုရမ်တစ်ခု မှာ ပြောထားတာက အောက်ပါအတိုင်း ဖြစ်ပါတယ်။
Any two independent orthonormal coordinate frames can be related by a sequence of rotations (not more than three) about coordinate axes, where no two successive rotations may be about the same axis.
{M} နဲ့ ကိုယ်စားပြု ခေါ်မယ့် world reference frame နဲ့ {B} နဲ့ ကိုယ်စားပြုမယ့် body reference frame တို့ကို အောက်မှာ ပုံနဲ့ နမူနာ ပြထားပါတယ်။

Pan angle rotation ကို α လို့ ခေါ်မှာ ဖြစ်ပြီး သူ့ကို computer vision ကနေ ရနိုင်ပါတယ်။ နောက် tilt angle rotation - β နဲ့ roll angle rotation - γ တို့ကိုတော့ ကိရိယာ ထဲမှာ တပ်ထားမယ့် accelerometers တွေကနေ မြေဆွဲအား (gravity- g ) ကို အာရုံခံ ပြီး သိရှိနိုင်ပါတယ်။ အရင်ဆုံး pan and tilt မရှိဘူး လို့ ယာယီ ယူဆပြီး roll angle ကို စဉ်းစားကြည့်ပါမယ်။ γ က သုည ဖြစ်နေချိန်မှာ Y ဝင်ရိုးက အထက်ကို ညွှန်ပြီး၊ X ဝင်ရိုးက right handed rule အတိုင်း ရေပြင်ညီ ရှိမယ်လို့ သတ်မှတ် ပါမယ်။ X ဝင်ရိုး accelerometer နဲ့ Y ဝင်ရိုး accelerometer တို့က မြေဆွဲအားကို အာရုံခံ တိုင်းလို့ရတဲ့ တန်ဖိုးတွေကို gx နဲ့ gy လို့ခေါ်ပါမယ်။ roll angle ကို သတ်မှတ် ထားပုံကို အောက် က ပုံမှာ ဖော်ပြထားပါတယ်။

Roll angle - γ ကို အောက်က ညီမျှခြင်း ကိုသုံးပြီး တွက်ယူလို့ ရပါတယ်။

$$\gamma = \tan^{-1}\frac{gx}{gy}$$

နောက် second နဲ့ third quadrants တွေ မှာ ရှိမရှိကို ပြန်စစ်ပါမယ်။

if(gy<0)
  if(gx>=0)
    γ=π+γ;
  else
    γ=-π+γ;
  end
end

ဒါဆို roll angle rotation matrix - R_γ ကို အောက်က အတိုင်းရနိုင်ပါတယ်။

$$\mathbf{R}_{\gamma}= \begin{bmatrix} \cos(\gamma) & -\sin(\gamma) & 0 \\ \sin(\gamma) & \cos(\gamma) & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

Roll angle ကို ပြင်ပြီးတဲ့နောက် tilt angle ကိုဆက်စဉ်းစားပါမယ်။ Tilt angle ရဲ့ သတ်မှတ်ပုံ ကို အောက်မှာ ပြထားပါတယ်။

Tilt angle - β ကို Y နဲ့ Z မြေဆွဲအား တိုင်းတာမှုတွေ ကနေ တွက်လို့ရပါတယ်။ သတိထား ဖို့ လိုတာကတော့ အဲဒီတန်ဖိုးတွေက roll angle ပြင်ထားတဲ့ တန်ဖိုးတွေ ဖြစ်ရပါမယ်။ ပြောင်းသွားတဲ့ Y နဲ့ Z မြေဆွဲအား တိုင်းတာမှုတွေ ကို အောက်ပါအတိုင်း တွက်ယူနိုင်ပါတယ်။

$$\begin{bmatrix} gx2 \\ gy2 \\ gz2 \end{bmatrix} = \mathbf{R}_{\gamma} \begin{bmatrix} gx \\ gy \\ gz \end{bmatrix}$$

ဒါဆို β ကို အောက်ပါ ညီမျှခြင်း ကနေရှာလို့ ရပါပြီ။

$$\beta = \tan^{-1}\frac{gy2}{gz2}$$ .

နောက် second and third quadrants တွေ အတွက် စစ်ပါမယ်။

if(gz2>=0)
  if(gy2>=0)
    β=-π+β;
  else
    β=π+β;
  end
end

အဲဒီအခါ Tilt angle rotation matrix -R_β က အောက်ပါအတိုင်းရလာပါမယ်။

$$\mathbf{R}_{\beta}= \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos(\gamma) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \cos(\gamma) \end{bmatrix}$$

အဲဒီလိုပါပဲ၊ pan angle rotation matrix - R_α နဲ့ pan angle ရဲ့ သတ်မှတ်ပုံတွေ ကိုလည်း အောက်ပါအတိုင်း တွေနိုင်ပါတယ်။

$$\mathbf{R}_{\alpha}= \begin{bmatrix} \cos(\alpha) & -\sin(\alpha) & 0 \\ \sin(\alpha) & \cos(\alpha) & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

အဲဒီ rotation matrices တွေရဲ့ မြှောက်လဒ် က rotation matrix တစ်ခု ပါပဲ။ $$\mathbf{R}=\mathbf{R}_{\alpha} \mathbf{R}_{\beta} \mathbf{R}_{\gamma}$$
အဲဒီအခါ အောက်ပါအတိုင်း ရရှိလာမှာပါ။
M=RB

R က rotation matrix တစ်ခု ဖြစ်တာမို့ သူ့ရဲ့ inverse ကို ရှာဖို့ လိုခဲ့ရင် အဲဒါက R transpose နဲ့ အတူတူပါပဲ။ ဒီ 3D rotation ကို LabVIEW မှာ MabLab code သုံးပြီး လုပ်ထားတာကို အောက်ကပုံမှာ မြင်နိုင်ပါတယ်။

အခုလို 3D transformation တွေအတွက် Quaternion algebra ကလည်း လွယ်ကူပြီး ခေတ်စားနဲ့ နည်းတစ်နည်း ပါပဲ။

Reference:
Kuipers, Jack B., Quaternions and rotation sequences : a primer with applications to orbits, aerospace, and virtual reality, Princeton University Press, 1999, ISBN: 0691058725.