K-means Clustering

in Japanese

Introduction

In this page, I will describe a brief explanation on the theory of the K-means clustering and implement a simple image segmentation by means of a function cv::kmeans the OpenCV provides to us.

Theory

Suppose we have a data set consisting of N points each of which is defined in the D-dimensional Euclidean space as
.
Our goal is to partition the data set into some number K of clusters, where the value of K is given.
When we introduce a vector to represent a center of each cluster, an objective function we should minimize is defined as

where, is a variable that equals 1 if the point belongs to the cluster k, otherwise 0.
The procedures to minimize J are as follows:

1. Initialize the cluster centers in some way.
2. Minimize J with respect to the , keeping the fixed.
3. Minimize J with respect to the , keeping the fixed.
4. Repeat these optimizations until both of the and the converge.

It's easy to realize the 2nd procedure. It's only necessary to choose which has the minimum distance to the point we now consider and set to 1.
In order to consider the 3rd procedure, we set J's derivative with respect to to zero,

It results in the following equation,

The denominator in this expression is equal to the number of points assigned to the cluster k. Namely, indicates the mean of all of the data points assigned to the cluster k. On the basis of these discussion, the procedures to minimize J are rewritten as follows:

1. Decide the initial centers of the clusters in some way.
2. Assign points to the closest clusters.
3. Calculate the mean of the points assigned to the cluster, and replace the center by the mean.
4. Repeat the 2nd and the 3rd procedures until the centers of the clusters converge.

A direct implementation of the 2nd procedure is to compute the distance between every center and every point and to find the minimum distance. As it is relatively slow, we have to introduce a data structure such as a tree.

Practice in OpenCV

A function to execute the K-means clustering is cv::kmeans. In the following program, the 3 dimensional space (RGB) is considered. A pixel on an image corresponds to a point in 3D space. The K-means clustering yields the K clusters each of which has a set of points with similar color.

1. 3rd-11th lines : Display an input image. A variable image indicates a H x W matrix with 3 channels. H and W are the height and the width of the image, respectively.
2. 13th-15th lines : The function cv::kmeans requires the matrix in which each row corresponds to a coordinate of a point. In this sample program, the coordinate of the point is (B,G,R). To meet the requirement, the matrix image must be modified by using a function reshape. After the modification, image translates into a T x 3 matrix with 1 channel (reshaped_image), where T=W * H.
3. 18th-21th lines : As the function cv::kmeans also requires the matrix consisting of floating-point elements, reshaped_image must be converted to the floating-point matrix (reshaped_image32f) by applying a function convertTo.
4. 23th-27th lines : After preparing two matrices for outputs (labels and centers) and creating an instance of cv::TermCriteria which represents conditions to quit the algorithm, we execute the function cv::kmeans. In this case, I chose the maximum number of iteration (100) as the condition to quit the algorithm. The final argument of the constructor of cv::TermCriteria is not used. The 6th argument of the function cv::kmeans is set to cv::KMEANS_RANDOM_CENTERS. This means that the initial centers of the clusters are randomly decided.
5. 29th line : Display result. A function show_result is as follows:

1. 3rd-7th lines : A variable labels indicates a T x 1 matrix. The each row of it has the value of the label. If cluster_number is 3, the value is either 0, 1, or 2. A variable centers is a cluster_number x 3 matrix. The each row of it has a coordinate of a center of a cluster. In this case, the coordinate corresponds to the vector specified in the RGB space.
2. 9th line : Prepare a matrix rgb_image in which the final image is stored.
3. 14th-16th lines : The matrix centers has floating-point values. These values are converted into std::uint8_t-type ones with the range [0,255]. Moreover, the number of channels is also translated into 3.
4. 18th-23th lines：The label is replaced by the corresponding (B,G,R) to make the final image rgb_image.

The results are as follows:
cluster_number=2

cluster_number=3

cluster_number=5

cluster_number=10

original image

Development Environment

1. Mac OS X 10.8.2
2. Processor：3.06 GHz Intel Core 2 Duo
3. Memory：4GB
4. Xcode4.6 with Apple LLVM 4.2（C++ Language Dialect → GNU++11, C++ Standard Library → libc++）
5. boost-1.51.0 built by the Apple LLVM 4.1. See here
6. opencv-2.4.3 built by the Apple LLVM 4.2. See here.

Reference

1. Pattern Recognition and Machine Learning, Christopher M. Bishop, Springer, p424-430
2. OpenCV Document

K-means クラスタリング

in English

はじめに

K-meansクラスタリングの理論とOpenCVを使った実践をまとめます。

理論

D次元ユークリッド空間内の点 を考える。

いま、N個の点 が与えられたとき、これらをK個のクラスタに分割することを考える。
各クラスタの中心座標を としたとき、最小にすべき目的関数Jを次式で定義する。

ここで、 は、点 がクラスタkに属するなら1を、そうでなければ0をとる変数である。
Jを最小にする手順は以下の通りである。

1. に適当な初期値を与える。
2. を固定し、 についてJを最小にする。
3. を固定し、 についてJを最小にする。
4. と が収束するまで、2,3の手順を繰り返す。

手順2を実現することは容易である。ある点 を考えたとき、それに最も近い を検索すればよい。
手順3を実現するため、Jを で偏微分する。

これより次式を得る。

上式の分母はクラスタkに属する点の総数である。すなわち、 は、クラスタkに属する点の重心となる。これらを踏まえて、Jを最小にする手順を書き直すと以下のようになる。

1. クラスタの中心座標を適当に決める。
2. 各点を一番近いクラスタに割り振る。
3. 割り振られた点の重心を計算し、それを新たなクラスタの中心座標とする。
4. クラスタの中心座標が変化しなくなるまで、2と3の手順を繰り返す。

手順2の最も直接的な実現方法は、各点と全てのクラスタの中心座標との距離を計算し、その最小値を探すというものである。これは大変時間がかかる。木構造などを導入し高速化すべきである。

OpenCVによる実践

K-meansクラスタリングを実行する関数はcv::kmeansである。以下の実装では、3次元色空間内の点を考え、類似色が同一のクラスタに属するような領域分割を行っている。

1. 3行目から11行目まで：入力画像を読み込み表示する。変数imageはチャンネル数3のH x Wの行列である。Hは画像の高さ、Wは画像の幅を表す。
2. 13行目から15行目まで：関数kmeansは、各行に点の座標を持つ行列を要求する。従って、入力画像を保持する行列imageの形状を、reshapeを使って変更する。この変更により、各行にB,G,Rの値の並ぶT x 3の行列になる。ここでT=W * Hである。
3. 18行目から21行目まで：関数kmeansは浮動小数点の行列を要求するので、convertToを使って値の型を変更する。
4. 23行目から27行目まで：関数kmeansを実行する。出力結果を格納する2つの行列labels, centersを用意する。さらに、アルゴリズムを停止させる条件を表現したクラス cv::TermCriteria のインスタンスを作成する。ここでは、最大繰り返し数を100とした。最後の引数の1はここではダミーである。kmeansの第6引数にcv::KMEANS_RANDOM_CENTERSを指定した。これは、クラスタの中心座標値の初期値を乱数で決めることを意味する。
5. 29行目：結果を表示する。その中身は以下の通りである。

1. 3行目から7行目まで：labelsはT x 1の行列である。各行にはラベルの値が収められている。cluster_numberを3とすれば、整数値0,1,2のいずれかとなる。centersはcluster_number x 3の行列である。各行にクラスタの中心座標が収められる。いまの場合、(B,G,R)の値が並ぶことになる。
2. 9行目から12行目まで：最終結果を収める行列rgb_imageを作る。
3. 14行目から16行目まで：centersの要素の型は浮動小数点である。これを0から255までのstd::uint8_tの型に変更する。さらに、チャンネル数を3に変更する。
4. 18行目から23行目まで：ラベルの値を(B,G,R)に置き換え、最終画像rgb_imageを作る。

結果を示す。
cluster_number=2

cluster_number=3

cluster_number=5

cluster_number=10

元画像

開発環境

1. Mac OS X 10.8.2
2. プロセッサ：3.06 GHz Intel Core 2 Duo
3. メモリ：4GB
4. Xcode4.6 with Apple LLVM 4.2（C++ Language Dialect → GNU++11, C++ Standard Library → libc++）
5. boost-1.51.0（Apple LLVM 4.1でコンパイルしたもの。こちら）
6. opencv-2.4.3（Apple LLVM 4.2でコンパイルしたもの。こちら）

参考文献

1. Pattern Recognition and Machine Learning, Christopher M. Bishop, Springer, p424-430
2. OpenCVドキュメント

Gallery of Level Set Method in 2D

in Japanese

I applied the previous implementation of the Level Set Method to various 2D images. See below.

gallery

Level Set Method -- Application to Arbitrary 3D Shapes --

in Japanese

Introduction

I have recently extended the previous implementation of the Level Set Method so that it can be applied to arbitrary 3D shapes. In this page, I would like to provide you with a brief explanation on the new implementation.
In order to accomplish the purpose, I used Open Asset Import Library (Assimp). The library allows us to load various 3D format files and to extract 3D information through the common interfaces the library comes with. I have made it possible for my engine to accept the STL format, which is one of 3D model formats, as an input.

Installation of Assimp Library

As MacPorts has it, we can easily install it as follows:

What's STL?

The following are a brief description:

1. It's an abbreviation for Stereo Lithography.
2. It describes a static 3D object.
3. It describes the object as a set of triangles.
4. The vertices of the triangle are listed in a counterclockwise order.
5. A normal vector for the triangle may also be listed. The vector is an unit one.
6. The STL file begins with solid, and ends with endsolid.
7. It has the format of either ASCII or BINARY.
8. It has no color information.
9. It has no compression.
10. It has 1 image.

For example, the STL file looks like this:

Extraction of 3D information

We can extract three vertices and one normal vector for each triangle through the interfaces Assimp Library provides.
In the 3D case, my implementation of the Level Set Method expects a std::uint8_t array with size of width x height x depth as input. The pixels on the surface of the object should be set to 128, and other pixels to 255. The procedures to convert the STL format into my requirement are as follows:

1. The scale conversion and the parallel shift are applied to an object to put it inside a box with the size of width x height x depth.
2. We find the 3D point (int-type) that has the distance less than 1 from each triangle and exists inside the triangle.
3. The pixel value is set to 128.

The adjacent pixels of the surface of the object has 128 and other pixels are set to 255.

Results

There are some demo movies below. All inputs have the size of 200x200x200 pixels.

gallery_of_level_set_method_in_3D

I downloaded some STL files from this site.

Development Environment

1. Mac OS X 10.8.2
2. Processor：3.06 GHz Intel Core 2 Duo
3. Memory：4GB
4. Xcode4.6 with Apple LLVM 4.2（C++ Language Dialect → C++11, C++ Standard Library → libc++）
5. boost-1.51.0 built by the Apple LLVM 4.1. See here
6. opencv-2.4.3 built by the Apple LLVM 4.2. See here.
7. assimp @3.0.1270 (MacPorts's package. Though MacPorts uses the different compiler from the one Xcode has, fortunately my implementation works well.)

My Source Code

Here is my source code. I tested it on the above environment.
The C++11 new features, such as lambda expression and thread library, are used to implement my engine.

Usage

See here.

Execution Command

The execution command looks like this: As soon as you type the return key, the application enters the pause mode. The calculation begins after typing 'p'.

Level Set 法 任意の3D形状への適用

in English

はじめに

先に実装したエンジンを、任意の３次元形状に適用できるよう拡張したので紹介します。
作業するにあたって、Open Asset Import Library（Assimpライブラリ）を利用しました。このライブラリは、さまざまな3Dモデルフォーマットを読み込み、ライブラリの提供するインターフェースを通して3D情報を取り出すことができます。今回は、STLフォーマットを入力値として受け取れるようLevel Set Methodのエンジンを拡張しました。

Assimpのインストール

macportsにあるのでこうします。

STLとは

簡単にまとめておきます。

1. Stereo Lithographyの略である。
2. 固定された３D物体を記述するものである。
3. 物体の表面を三角形に分割して表現する。
4. 三角形の頂点が反時計回りに列挙される。
5. 各三角形の法線ベクトルも記入されることがある。これは単位ベクトルでなければならない。
6. solid ではじまり、endsolidで終わる。
7. binaryとasciiがある。
8. 色情報はない。
9. 非圧縮である。
10. 1 imageである。

例えばこんな感じです。

3D情報の抽出

Assimpライブラリを使うと、頂点の座標と法線ベクトルを取り出すことができます。このとき頂点の数と法線ベクトルの数は同じです。上で見たように、本来であれば3つの頂点につき1つの法線ベクトルが対応するはずです。Assimpでは法線ベクトルが重複して取り出されます。
私のLevel Set Methodエンジンの入力値は、3Dの場合、サイズwidth x height x depthのstd::uint8_t型の配列です。オブジェクトに属していない画素の値は255、属している画素の値は128としています。入力値への変換手順は以下の通りです。

1. 適当なスケール変換と平行移動を行い、オブジェクトをwidth x height x depthの範囲に収める。
2. 各三角形（float型の頂点を持つ）から距離1以下の範囲に存在し、かつ、三角形内部に存在する座標値（int型）を求める。
3. その画素の値を128にする。

オブジェクトの表面近傍にある画素だけが128となり、その中身は255のままです。

結果

こちらで公開します。入力値のサイズは全て200x200x200 pixelsです。

gallery_of_level_set_method_in_3D

STLデータはこちらからダウンロードしました。

開発環境

1. Mac OS X 10.8.2
2. プロセッサ：3.06 GHz Intel Core 2 Duo
3. メモリ：4GB
4. Xcode4.6 with Apple LLVM 4.2（C++ Language Dialect → C++11, C++ Standard Library → libc++）
5. boost-1.51.0（Apple LLVM 4.1でコンパイルしたもの。こちら）
6. opencv-2.4.3（Apple LLVM 4.2でコンパイルしたもの。こちら）
7. assimp @3.0.1270（macportsでインストールしたもの。Xcodeで使われるコンパイラと違いますが動作しました。）

ソース

こちらです。上記の開発環境で動作確認しました。
注意：lambda式やthreadなど、C++11の新しい機能を使用しています。

実行方法

先のページと同じです。

コマンド例

起動後、すぐに一時停止するようにしています。「p」を打つと計算を始めます。