Myanmar Astrological Calendar Days

A typical Myanmar calendar usually describes astrological days such as auspicious days and ill fated days. In this post, we will present about the following days and javascript code for them.

Algorithm, Program and Calculation of Myanmar Calendar

မြန်မာဘာသာ ဖြင့်ဖတ်ရန်

An easier and faster modern method to calculate a past, present or future date in the Myanmar calendar (the Burmese calendar) is presented. The constants, formulas and steps are clearly defined to derive all the essential elements of the Myanmar calendar such as Myanmar year, Myanmar month, waxing or waning moon and day. Unlike the existing methods, this method is easy even for a person who is not familiar with the Myanmar calendar and its terms.

k-means clustering using custom distance measuring method

I have developed a MATLAB function to perform k-means clustering which enables custom distance measuring method. For example, to sort out histograms, chi-square distance may be more suitable. The following example function uses chi-square distance and you can replace it with any distance measurement method.

---

%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

Rotations in 3D space using Euler angles

Rotations in 3D space to relate body reference frame of a device to world reference frame using an Euler angle sequence are discussed. Euler stated that

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.
The illustration of the world reference frame, {M}, and the body reference frame, {B}, is shown in the following figure.

Pan angle rotation, α, is acquired using computer vision. And, tilt angle rotation, β, and roll angle rotation, γ, are acquired by sensing the gravity, g, using onboard accelerometers. Let us consider the roll angle first, assuming there are no pan and tilt angles. When the roll angle, γ, is zero Y-axis is upward vertical and X-axis is horizontal pointing according to right handed rule. The gravity sensed by X-axis accelerometer and Y-axis accelerometer are denoted by gx and gy respectively. The definition of the roll angle is illustrated in the following figure.

The roll angle, γ, is calculated as

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

Then, it is checked for second and third quadrants as

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

The roll angle rotation matrix, R_γ, is obtained as follows.

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

After the roll angle is corrected, let us consider for tilt angle. The definition of the tilt angle is illustrated in the following figure.

The tilt angle, β, can be calculated from Y and Z components of the gravity. It is important to note that these gravity component values should be in the roll angle rotated frame. The new values for the gravity components are obtained as

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

The tilt angle, β, is calculated as
$$\beta = \tan^{-1}\frac{gy2}{gz2} $$.
Then, it is checked for second and third quadrants as

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

The tilt angle rotation matrix, R_β, is obtained as follows.

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

Similarly, calculation of the pan angle rotation matrix, R_α, and description of pan angle definition are as follows.

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

The product of rotation matrices is itself a rotation matrix. $$ \mathbf{R}=\mathbf{R}_{\alpha} \mathbf{R}_{\beta} \mathbf{R}_{\gamma}$$
And, M=RB
Since R is a rotation matrix, inverse of R is obtained by transposing R. The implementation of this 3D rotation in LabVIEW using MabLab code is shown in the following figure.

Quaternion algebra is also easy and popular approach for such 3D transformation.

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

Skeleton of Image Region

The structural shape of an image region can be represented by a graph. It is achieved by using a thinning algorithm. The resulting graph is called the skeleton of the region. It can be a useful preprocessing step for some applications such as optical character recognition. The following figures shows an original structural shape and its skeleton.

The skeleton of a region is obtained using various methods. Two popular methods among them are Medial Axis Transformation (MAT) , and Two Step Thinning. In Medial Axis Transformation, each point in the region finds its closest neighbour on the boundary of the region. If it finds more than one such neighbour, it belongs to the skeleton of the region. Although MAT is easy to understand, it needs a lot of calculation. We find the distance transform of the region first. The distance transform replaces each pixel value in the region with the distance of the pixel from the nearest neighbour on the boundary of the region. The nearer pixels from the boundary have the lower distance values (lower intensity) and the farther pixels from the boundary have the larger distance values (higher intensity). From the distance values, the local maximums along row or column are searched and defined as pixels in the skeleton. The example MATLAB code for MAT can be seen here (MATeg.m). The original image, the distance transformed image, and the skeleton using MAT are shown in the following figures.

Two Step Thinning algorithm is more efficient and faster than MAT. Two step thinning algorithm iterately delete edge points of a region subject to the constraints that deletion of these points does not remove end points, does not break connectivity, and does not cause excessive erosion of the region. The region points are assumed to have value 1 and background points to have value 0. The method consists of successive passes of two basic steps applied to the given region.

With reference to the 8 neighborhood notation shown in the figure, step 1 flags a contour point p1 for deletion if the following conditions are satisfied: (a) 2 ≤ N(p1) ≤ 6 (b) T(p1) = 1 (c) p2.p4.p6 = 0 (d) p4.p6.p8 = 0 where N(p1) is the number of nonzero neighbors of p1, and T(p1) is the number of 0 to 1 transitions in the ordered sequence p2, p3, ... , p8, p9, p2. In step 2, conditions (a) and (b) remain the same, but conditions (c) and (d) are changed to (c') p2.p4.p8= 0 (d') p2.p6.p8= 0 Step 1 is applied to every border pixel in the binary region under consideration. If one or more of conditions (a)-(d) are violated, the value of the point in question is not changed. If all conditions are satisfied the point is flagged for deletion. However, the point is not deleted until all border points have been processed. This delay prevents changing the structure of the data during execution of the algorithm. After step 1 has been applied to all border points, those that were flagged are deleted (changed to 0). Then step 2 is applied to the resulting data in exactly the same manner as step 1. Thus one iteration of the thinning algorithm consists of (1) applying step 1 to flag border points for deletion; (2) deleting the flagged points; (3) applying step 2 to flag the remaining border points for deletion; and (4) deleting the flagged points. This basic procedure is applied iteratively until no further points are deleted, at which time the algorithm terminates, yielding the skeleton of the region. Condition (c) p2.p4.p6 = 0 and (d) p4.p6.p8 = 0 are satisfied simultaneously by the minimum set of values; (p4 = 0 or p6= 0) or (p2=0 and p8=0). Similarly, conditions (c') and (d') are satisfied simultaneously by the following minimum set of values: (p2=0 or p8=0) or (p4=0 and p6=0). The example MATLAB code for Two Step Thinning algorithm can be seen here (TSTeg.m).

The above codes may be useful to understand the methods and to implement them in other programming languages. In MATLAB, the skeleton of all regions in a binary image is generated via function bwmorph. sk=bwmorph(bwImg,'skel',Inf); The result of the above operation is shown below.

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

Monte Carlo integration

As I am reading about Particle Filter, just for fun, I wrote a Matlab program that performs Monte Carlo integration using uniform distribution. The equation and the program are as follows.

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

CRC Calculation in VB and C

Just to share a few software modules that were written in Visual Basic 2005 and C for the calculation of CRC -Cyclic Redundancy Check.

CRC Calculation in VB2005

The followings are the source code for various CRC calculations in VB2005. To make the calculation faster, they use CRC tables.


CRC Calculation - GitHub

An example usage for calculation of CRC16 CCITT is shown below.

Dim StrIn as String= "String to calculate CRC"
Dim CRCVal16 As UInt16 = 0
Dim crc As String
CRCVal16 = CRC16_CCITT.Calculate(StrIn)
crc = CRC16_CCITT.ToString(CRCVal16)

Initial value for CRC16 CCITT is 0xFFFF. The following example calculate CRC for Str1 and use that CRC value as initial value to calculate Str2.

CRCVal16 = CRC16_CCITT.Calculate(Str1)
CRCVal16 = CRC16_CCITT.Calculate(Str2, CRCVal16)
crc = CRC16_CCITT.ToString(CRCVal16)

CRC Calculation in C

The followings are the source code for various CRC calculations in C. To save storage, they do not use CRC tables .


CRC Calculation - GitHub

An example usage for calculation of CRC16 CCITT is shown below.

#define STRLEN 4
char str[STRLEN]={0x01,0x01,0x00,0x0B};
unsigned char c[2];
unsigned int crc;
//Calculate CRC16 CCITT
crc=CRC16CCITT_InitialValue();
crc=CRC16CCITT_Calculate(str,STRLEN,crc);
CRC16CCITT_ToString(crc,c);
printf("CRC16 CCITT = %02X %02X \n",c[0],c[1]);

Online checksum calculator such as the following one may be useful to debug the code.

Online Checksum Calculator