Integration of Accelerometers

3-axis accelerometer တွေကို ကိုယ်ထည်အမာ တစ်ခုမှာ လိုအပ်ချက်နဲ့ ကိုက်ညီသလို တပ်ထားပါတယ်။ အဲဒီ ကိုယ်ထည်အမာ ကလည်း ရွေ့ချင်တဲ့ အဖြောင့် (linear) နဲ့ အလှည့် (angular) လှုပ်ရှားမှုတွေနဲ့ ရွေ့နေပါတယ်။ တပ်ထားတဲ့ accelerometer တွေက linear acceleration တွေကိုပဲ အာရုံခံနိုင်တဲ့ အတွက် angular ရွေ့လျားမှု တွေကို မသိဘူး ဖြစ်နေပါတယ်။

သိချင်တာ တစ်ခုက accelerometer တွေ အများကြီးက တိုင်းတာလို့ ရတဲ့ တန်ဖိုးတွေကို ပေါင်းပြီး ထပ်တူညီတဲ့ စုပေါင်း integrated equivalent accelerometer တစ်ခုအနေနဲ့ ရနိုင်မလား ဆိုတာပါပဲ။
စဉ်းစားကြည့်လိုက်တော့ ရနိုင်တယ် လို့ ထင်ပါတယ်။

Low Pass and High Pass Filter

Low pass နဲ့ high pass filter လေးတွေက အမြဲသုံးနေ၊ တွေ့နေ ရလို့ သူတို့ နဲ့ ပတ်သက်တဲ့ note အချို့ ရေးချ ထားချင်ပါတယ်။ ရိုးရှင်းတဲ့ first order RC circuit လေးတွေ အတွက်ပါပဲ။

Figure. Low pass and high pass filter.

Fitting a curve to a function

ပေးထားတဲ့ curve တခုနဲ့ အကိုက်ညီဆံုး piecewise linear function တခုကို ရှာတဲ့ နမူနာ MatLab code တခုပါ။

%Fitting a curve to a piecewise linear function using least square method
%with the fixed start and end points

%clear command windows
clc;

%clear workspace
clear all;
%--------------------------------------------------------------------------
%Read curves
load Curve.dat;

%find lower left and upper right corners 
xD=Curve(1,:);
yD=Curve(2,:);
x0=min(xD); x1=max(xD); 
y0=min(yD);  y1=max(yD);
n=5;%order n i.e, the number of segments
FS=x1-x0;
global xs;
global yB;
global yE;
xs=(0:FS/n:FS);
yB=y0;
yE=y1;
%initialize coefficients (y values) to find 
%excluding the first and the last one
ys_initial=xs(2:n);
% options = optimset('MaxFunEvals',1000,'MaxIter',1000);
% LowerBoundC=y0;
% UpperBoundC=y1;
% ys =lsqcurvefit(@LineSeg,ys_initial,xD,yD,LowerBoundC,UpperBoundC,options);
ys =lsqcurvefit(@LineSeg,ys_initial,xD,yD);
%fixed start and end points
ys=[y0 ys y1];
%--------------------------------------------------------------------------
%Plot 
hFig1 = figure(1);
set(hFig1, 'Position', [600 100 500 300])
plot(xD,yD,':g','LineWidth',3,...
                'MarkerEdgeColor','b',...
                'MarkerFaceColor','b',...
                'MarkerSize',2)
hold on;             
plot(xs,ys,'-rs','LineWidth',1,...
     'MarkerEdgeColor','r',...
     'MarkerFaceColor','r',...
     'MarkerSize',4)
hold off;        
grid on;
%--------------------------------------------------------------------------

LineSeg.m

function [y] = LineSeg(ys,x)

global xs
global yB
global yE
L=length(x);
y=zeros(1,L);
ys=[yB ys yE];
LS=length(ys);
for i=1:L
    xt=x(i);
    %--------------------------
    %find segment
    for j=2:LS
        if(xs(j)>=xt) 
            break;
        end
    end    
    %interpolate
    x2=xs(j);
    y2=ys(j);
    y1=ys(j-1);
    x1=xs(j-1);
    y(i)=y1+(y2-y1)/(x2-x1)*(xt-x1);
    %--------------------------
end 

Image Retrieval

ဓာတ်ပုံ တပုံ (test image) ပေးလိုက်ရင် သိမ်းထားတဲ့ ပုံတွေထဲက အတူဆုံး ပုံတွေ ကို ပြန်ရွေးထုတ်ပေးတဲ့ program ကို MatLab နဲ့ ရေးကြည့်ထားပါတယ်။

Selective Color Modification

အောက်မှာ ပြထားတဲ့ ပုံထဲက လမ်းအပြာလေးကိုပဲ ရွေးပြီး အစိမ်းရောင် ပြောင်းကြည့် ချင်တာနဲ့ MatLab program လေး ရေးကြည့်ပြီး Octave သုံးပြီး run ကြည့်လိုက်တယ်။ အပြာရောင် ကို အစိမ်းပြောင်းချင်တာဆိုတော့ rgb color model မှာပဲ blue နဲ့ green ကို နေရာချင်း လဲလိုက်ရင် ရပေမယ့် သိပ် ယေဘူယျ မဆန်ဘူးထင်တာနဲ့ hsv color model သုံးလိုက်တယ်။ အဲဒါဆို hue ကို ကြိုက်တဲ့ အရောင်ကနေ လိုချင်တဲ့ အရောင်ကို အလွယ်တကူ ပြောင်းလို့ရသွားတာပေါ့။

GNU Octave as an Alternative to MatLab

GNU Octave၊ FreeMat နဲ့ Scilab တို့က MATLAB လိုမျိုး တွက်ချက်မှုတွေ လုပ်ပေးနိုင်တဲ့ အလကားရတဲ့ free open source software တွေဖြစ်ပါတယ်။ လိုက်ရှာ ဖတ်ကြည့်တော့ GNU Octave က အကောင်းဆုံးလို့ ကြားတာနဲ့ သူ့ကို စက်ထဲမှာ ထည့်ပြီး စမ်းကြည့်တော့ တော်တော် သဘောကျ သွားပါတယ်။ GNU Octave ကို ဝင်းဒိုး ပေါ်မှာတင်ရတာ တော့ တော်တော်လွယ်ပါတယ်။ သူ့ရဲ့ Download page က ဖိုင်တွေကို ယူပြီး စက်ထဲမှာ folder တစ်ခု အနေနဲ့ ထုတ်ယူလိုက်ပါမယ်။ ပြီးရင် shortcut ဖိုင်ကို ကြိုက်တဲ့နေရာမှာ ထားပြီး သူ့ရဲ့ properties မှာ path နဲ့ icon ကို ပြင်ပေး လိုက်ရုံပါပဲ။ သူရဲ့ shortcut ကနေ Octave ကို ဖွင့်လိုက်ရင် အောက်က ပုံမှာ ပြထားတဲ့ အတိုင်း ကွန်မန်း ပေးလို့ရတဲ့ ဝင်းဒိုး ပွင့်လာပါမယ်။ အဲဒီမှာ MATLAB ကွန်မန်း တွေ တိုက်ရိုက် ရိုက်ထည့်ရင် လည်းရပါတယ်။

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;
}

Skeleton of Image Region

ပုံရိပ် တစ်ခုရဲ့ တည်ဆောက်ပုံ ပုံသဏ္ဌာန် ကို မျဉ်းကြောင်းများနဲ့ တုတ်ချောင်းပုံ လုပ်ပြီး ကိုယ်စားပြု ဖော်ပြနိုင်ပါတယ်။ အဲဒီလို ဖော်ပြဖို့အတွက် ပုံသဏ္ဌာန် တစ်ခုကို တဖြည်းဖြည်း ပိန်ပါးသွားအောင် လုပ်နည်း တစ်ခုခု ကို သုံးလို့ရပါတယ်။ ရလာတဲ့ မျဉ်းကြောင်းများနဲ့ တုတ်ချောင်းပုံ ကို အဲဒီ ပုံရဲ့ အရိုး (skeleton) လို့ ခေါ်ပါမယ်။ ဒီနည်းဟာ စာလုံးများကို စက်က သိရှိနားလည် မှတ်မိ အောင် လုပ်တဲ့ နေရာ (optical character recognition) တွေမှာ ပဏာမ အဆင့် အနေနဲ့ အသုံးကျကောင်း ကျနိုင်ပါတယ်။ အောက်ကပုံတွေမှာ နမူနာ မူရင်းပုံနဲ့ သူ့ရဲ့ အရိုးကို ဖော်ပြထားပါတယ်.

ပုံရိပ်တစ်ခု ရဲ့ အရိုးကို နည်းအမျိုးမျိုးနဲ့ ရယူနိုင်ပြီး အဲဒီ နည်းတွေထဲက လူသိများတဲ့ နည်းနှစ်နည်းက အလယ်ဝင်ရိုးပြောင်းနည်း Medial Axis Transformation (MAT) နဲ့ နှစ်ဆင့်ပါးနည်း Two Step Thinning တို့ပါ။ MAT မှာ ပုံသဏ္ဌာန် ထဲမှာ ရှိတဲ့ အစက် တစ်စက် ချင်းစီ အတွက် သူတို့ ရဲ့ အနီးဆုံး အနားသတ် ကို လိုက်ရှာပါတယ်။ အဲဒီလို အကွာအဝေးတူ အနီးဆုံး အနားသတ် တစ်ခုမက ရှိတဲ့ အစက်ကို အရိုး လို့သတ်မှတ်ပါတယ်။ MAT ကနားလည်ဖို့ လွယ်ပေမယ့် သူ့ကိုရဖို့ တွက်ချက်မှုတွေ အများကြီး လိုပါတယ်။ သူ့ကို တွက်ဖို့ အကွာအဝေးပုံစံ ပြောင်းခြင်း (distance transform) ကို အရင် လုပ်ပါတယ်။ အကွာအဝေးပုံစံ ပြောင်းတယ်ဆိုတာ ပုံသဏ္ဌာန် ထဲမှာ ရှိတဲ့ အစက် တစ်စက် ချင်းစီ ကို အနီးဆုံး အနားသတ်အထိ အကွာအဝေးနဲ့ အစားထိုးတာပါ။ ရလာတဲ့ ပုံမှာ အနားသတ်နဲ့ နီးတဲ့ အစက်က အကွာအဝေး တန်ဖိုးနည်းတော့ အလင်းအားပျော့ ပြီး၊ အနားသတ်နဲ့ ဝေးတဲ့ အစက်က အကွာအဝေးတန်ဖိုးကြီးတော့ ပိုလင်းပါမယ်။ နောက် အစက်တစ်စက် ခြင်းစီကို ဘေးတိုက် ဒါမှမဟုတ် ဒေါင်လိုက် ကပ်ရပ် အစက်တွေနဲ့ ယှဉ်ကြည့်ပြီး အလယ်မှာ အမြင့်ဆုံး ဖြစ်နေတဲ့ အစက်တွေကို အရိုးလို့ သတ်မှတ်ပါမယ်။ နမူနာ MATLAB ကုဒ် ကို ဒီမှာ (MATeg.m) ကြည့်နိုင်ပါတယ်။ မူရင်း ပုံသဏ္ဌာန်ရယ်၊ အကွာအဝေးပုံစံ ရယ်၊ ရလာတဲ့ အရိုးရယ် ကို အောက်က ပုံတွေမှာ ကြည့်နိုင်ပါတယ်။

နှစ်ဆင့်ပါးနည်း ကတော့ MAT နဲ့ယှဉ်ရင် ပိုမြန်၊ ပိုထိရောက်ပါတယ်။ နှစ်ဆင့်ပါးနည်း က ပုံသဏ္ဌာန် တစ်ခုရဲ့ အနားတွေကို ထပ်တလဲလဲ ဖျက်သွားတာပါ။ အဲဒီလို ဖျက်တဲ့ အခါမှာ အဆုံးမှတ်တွေကို မဖျက်ပါဘူး။ တနေရာနဲ့ တနေရာ ဆက်စပ်နေတဲ့ အမှတ်တွေကို မဖျက်ပါဘူး။ ပုံသဏ္ဌာန် ထဲမှာ ရှိတဲ့ အစက် တွေကို တန်ဖိုး တစ်ရှိတယ်လို့ မှတ်ယူ မှာ ဖြစ်ပြီး နောက်ခံ အစက်တွေကို တော့ တန်ဖိုး သုည ရှိတယ် လို့ မှတ်ယူပါမယ်။ ဒီနည်းမှာ အဆင့် နှစ်ဆင့် ကို ပုံပေါ်မှာ အထပ်ထပ် အဖန်ဖန် နောက်ဆုံး မပြောင်းလဲတဲ့ အရိုး ရတဲ့အထိ သုံးသွားမှာပါ။

အပေါ်ပုံမှာ အလယ် အစက်ရဲ့ ကပ်ရပ် အစက်တွေကို အခေါ်အဝေါ်သတ်မှတ်ပုံကို ပြထားပါတယ်။ ပထမ အဆင့်မှာ အစက် p1 ကို အောက်က အချက်လေးချက်နဲ့ ကိုက်ညီရင် ဖျက်ပစ် ပါမယ်။ (a) 2 ≤ N(p1) ≤ 6 (b) T(p1) = 1 (c) p2.p4.p6 = 0 (d) p4.p6.p8 = 0 အဲဒီမှာ N(p1) ဆိုတာက p1 နဲ့ ကပ်ရပ် အစက်တွေထဲက သုညမဟုတ်တဲ့ အစက် အရေအတွက်ဖြစ်ပြီး၊ T(p1) ဆိုတာက p2, p3, ... , p8, p9, p2 ဆိုတဲ့ အစီအစဉ်အတိုင်း ကြည့်သွားရင် ၀ ကနေ ၁ ကို ပြောင်းသွားတဲ့ အရေ အတွက်ပါ။ ဒုတိယ အဆင့်မှာ (a) နဲ့ (b) ကအတူတူပါပဲ။ ဒါပေမယ့် (c) နဲ့ (d) ကတော့ အောက်မှာ ပြထားတဲ့ အတိုင်း ပြောင်းသွားပါတယ်။ (c') p2.p4.p8= 0 (d') p2.p6.p8= 0 ပထမ အဆင့် ကို ပုံသဏ္ဌာန်တစ်ခုရဲ့ အနားသတ်ပေါ်မှာရှိတဲ့ အစက်တိုင်းအတွက် စဉ်းစားပါမယ်။ အချက် (a) ကနေ (d) ထဲက တစ်ချက်ချက် မကိုက်တာနဲ့ အဲဒီ အစက်ကို မဖျက်တော့ ပါဘူး။ အချက်အားလုံး ကိုက်ညီရင် တော့ ဖျက်ဖို့ အတွက် မှတ်ထားလိုက် ပါမယ်။ ဒါပေမယ့် အဲဒီ အစက် ကို အနားသတ် ပေါ်မှာရှိတဲ့ အစက် အားလုံးကို စဉ်းစားလို့ မပြီးမချင်း တကယ် မဖျက်သေးပါဘူး။ အဲဒီလို နောက်ကျမှ ဖျက်တဲ့ အတွက် တစ်ခြား အစက် တွေကို စဉ်းစား နေတဲ့ အချိန်မှာ ပုံသဏ္ဌာန် ပြောင်းသွားတာကို ကာကွယ်နိုင်ပါတယ်။ ပထမ အဆင့် ကို အနားသတ်ပေါ်မှာရှိတဲ့ အစက် အားလုံးကို စဉ်းစာပြီးပြီဆိုတာနဲ့ ဖျက်ဖို့မှတ်ထားတဲ့ အစက် အားလုံးကို ဖျက်လိုက်ပါမယ်။ (အစက်ရဲ့ တန်ဖိုးကို သုည ပြောင်းလိုက်တာပါ။) ရလာတဲ့ ရလဒ် ပေါ်မှာ ဒုတိယ အဆင့် ကိုလဲ ပထမအဆင့် အတိုင်း ထပ်သုံးပါမယ်။ ဒီ အဆင့်နှစ်ဆင့်ကို အကြိမ်ကြိမ် အဖန်ဖန် ပြန်သုံးမှာ ဖြစ်ပြီး အကြိမ်တစ်ကြိမ်တိုင်းမှာ အောက်က လုပ်ငန်းလေးခုကို လုပ်ပါတယ်။ (၁) ပထမ အဆင့်ကို သုံးပြီး ဖျက်ရမယ့် အစက်တွေကို မှတ်ပါတယ်။ (၂) မှတ်ထားတဲ့ အစက်တွေကို ဖျက်ပါတယ်။ (၃) ဒုတိယ အဆင့်ကို သုံးပြီး ဖျက်ရမယ့် အစက်တွေကို မှတ်ပါတယ်။ (၄) မှတ်ထားတဲ့ အစက်တွေကို ဖျက်ပါတယ်။ ဒီအဆင့်တွေကို သုံးပြီး ပါးပါးသွားရင်း နောက်ဆုံး ဖျက်စရာ အစက် မကျန်တော့ဘူးဆိုရင် ရပ်လိုက်ပါမယ်။ ရလာတဲ့ ပုံက အရိုးဖြစ်ပါတယ်။ အချက် (c) p2.p4.p6 = 0 နဲ့ (d) p4.p6.p8 = 0 ကို တစ်ခါထဲ ပေါင်းပြီး (p4 = 0 or p6= 0) or (p2=0 and p8=0) လို့ တန်းစစ်လို့လဲရပါတယ်။ အဲဒီလိုပဲ (c') နဲ့ (d') ကိုလဲ (p2=0 or p8=0) or (p4=0 and p6=0) လို့ ပေါင်းစစ်လို့ ရပါတယ်။ Two Step Thinning အတွက်နမူနာ MATLAB ကုဒ်ကို ဒီမှာ (TSTeg.m) တွေ့နိုင်ပါတယ်။

ရှေ့မှာ ဖော်ပြခဲ့တဲ့ ကုဒ်တွေက နည်းတွေကိုနားလည်ဖို့ ရယ်၊ တစ်ခြား ပရိုဂရမ် ဘာသာတွေနဲ့ ပြန်ရေးဖို့ရယ် အတွက်တော့ အထောက် အကူ ဖြစ်မယ်ထင်ပါတယ်။ တကယ်တော့ MATLAB မှာ bwmorph ဆိုတဲ့ ဖန်ရှင်ကို သုံးလိုက်တာနဲ့ အရိုးကို တန်းရနိုင်ပါတယ်။ sk=bwmorph(bwImg,'skel',Inf); အဲဒီလိုခိုင်းလိုက်ရင် ထွက်လာတဲ့ ရလဒ် ကို အောက်မှာ ပြထားပါတယ်။

Reference: Rafael C. Gonzalez, Richard E. Woods, Steven L. Eddins, "Digital Image Processing Using MATLAB", Second Edition, Mc Graw Hill (Asia), 2011.

Testing fsolve and anonymous function

တစ်ခါက coefficient တွေ ပြောင်းလဲ နေတဲ့ equation တွေကို ဖြေရှင်းဖို့ လို လာပါတယ်။ ဒါနဲ့ MatLab ဖန်ရှင်တစ်ခု ဖြစ်တဲ့ fsolve ရယ်၊ function handle နဲ့ anonymous function တွေကို သုံးပြီး စမ်းကြည့်ထားတာပါ။

s=[3 7 11];
d=[1 3 5];
%ans= (2,1), (5,2), (8,3)
for i=1:3
x0 = [-5; -5]; % Make a starting guess at the solution
myfun =@(x) [x(1) + x(2) - s(i);
x(1) - x(2) - d(i)];
[x,fval] = fsolve(myfun,x0)  % Call 
end

Curve Fitting in Matlab

မိတ်ဆွေ တစ်ယောက်က data sample တွေက နေ curve equation ထုတ်ဖို့ လိုနေ ပါတယ်။ Excel မှာ polynomial ကို order 6 အထိပဲ ရလို့ higher order ရအောင် Matlab program ကူရေးပေး ဖို့ ဆိုလာ ပါတယ်။ ဂရပ် ပုံကို ပါ ဆွဲပေး ဖို့ program ရယ်၊ နမူနာ ဒေတာ ဖိုင်ရယ် ကို အောက်မှာ ရေးကြည့်ထား ပါတယ်။ ကျွန်တော် ကိုယ်တိုင်လဲ မကြာခဏ data fit လုပ်ဖို့ လိုတာ ကြောင့် ကျွန်တော့် အတွက်လဲ ဒါဟာ အသုံးဝင်ပါတယ်။ See d.txt and cf.m.

%--------------------------------------
%Curve fitting
%2011-Jun-03
%Programmer: Yan Naing Aye
%--------------------------------------
clc;
close all;
clear all;
%--------------------------------------
%get data from file that can be edited
load('d.txt');
X=d(:,1);
Y=d(:,2);
%--------------------------------------
%define order
N=2;
%--------------------------------------
%curve fitting
p=polyfit(X,Y,N);
%generate curve fitted data
[nr nc]=size(X);
YA=zeros(nr,1);
for i=0:N
    YA=YA+p(i+1)*X.^(N-i);
end
%--------------------------------------
%Plot Y vs X
figure;
plot(X,Y,' rd','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',6)
hold on;
plot(X,YA,'b','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',6)
hold off;            
grid on;
Title('Y vs X')
xlabel('X');
ylabel('Y');
legend('Sampled data','Fitted curve',...
       'Location','NW')
%--------------------------------------

Monte Carlo integration

Particle Filter အကြောင်း ဖတ်ကြည့်ရင်း uniform distribution နဲ့ Monte Carlo integration လုပ်တဲ့ Matlab program လေး တစ်ပုဒ် ကို အပျော် ရေးကြည့် ထားပါတယ်။ အဲဒီ equation ရယ်၊ program ရယ် ကို အောက်မှာ ကြည့်နိုင်ပါတယ်။

clc;
close all;
clear all;
%--------------------------------------
%Monte Carlo Integration
%--------------------------------------
%--------------------------------------
%limits
LowerLimit=1;
UpperLimit=3;
%--------------------------------------
%generate uniformly distributed random 
%numbers within the range
n=10000;
x=LowerLimit+(UpperLimit-LowerLimit)*rand(n,1);
%--------------------------------------
%Probability distribution function for 
%uniformly distributed random numbers
q=1/(UpperLimit-LowerLimit);
%--------------------------------------
%function to integrate
f=3.*x.*x;
i=mean(f./q)
%--------------------------------------